Commit f1e913f5 authored by Marco Selig's avatar Marco Selig
Browse files

response double measurement adjoint bug fixed; alpha retrieval implemented.

parent 1a308962
...@@ -163,7 +163,7 @@ import powerspectrum as gp ...@@ -163,7 +163,7 @@ import powerspectrum as gp
pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679 pi = 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
__version__ = "0.9.5" __version__ = "0.9.6"
##----------------------------------------------------------------------------- ##-----------------------------------------------------------------------------
...@@ -10429,11 +10429,11 @@ class response_operator(operator): ...@@ -10429,11 +10429,11 @@ class response_operator(operator):
for ii in xrange(np.size(assign,axis=1)): for ii in xrange(np.size(assign,axis=1)):
if(np.any(assign[:,ii]>=self.domain.dim(split=True)[ii]))or(np.any(assign[:,ii]<-self.domain.dim(split=True)[ii])): if(np.any(assign[:,ii]>=self.domain.dim(split=True)[ii]))or(np.any(assign[:,ii]<-self.domain.dim(split=True)[ii])):
raise IndexError(about._errors.cstring("ERROR: invalid bounds.")) raise IndexError(about._errors.cstring("ERROR: invalid bounds."))
self.assign = assign.T ## transpose self.assign = assign ## transpose
if(target is None): if(target is None):
## set target ## set target
target = point_space(np.size(assign,axis=0),datatype=self.domain.datatype) target = point_space(np.size(self.assign,axis=0),datatype=self.domain.datatype)
else: else:
## check target ## check target
if(not isinstance(target,space)): if(not isinstance(target,space)):
...@@ -10442,8 +10442,8 @@ class response_operator(operator): ...@@ -10442,8 +10442,8 @@ class response_operator(operator):
raise ValueError(about._errors.cstring("ERROR: continuous codomain.")) ## discrete(!) raise ValueError(about._errors.cstring("ERROR: continuous codomain.")) ## discrete(!)
elif(np.size(target.dim(split=True))!=1): elif(np.size(target.dim(split=True))!=1):
raise ValueError(about._errors.cstring("ERROR: structured codomain.")) ## unstructured(!) raise ValueError(about._errors.cstring("ERROR: structured codomain.")) ## unstructured(!)
elif(np.size(assign,axis=0)!=target.dim(split=False)): elif(np.size(self.assign,axis=0)!=target.dim(split=False)):
raise ValueError(about._errors.cstring("ERROR: dimension mismatch ( "+str(np.size(assign,axis=0))+" <> "+str(target.dim(split=False))+" ).")) raise ValueError(about._errors.cstring("ERROR: dimension mismatch ( "+str(np.size(self.assign,axis=0))+" <> "+str(target.dim(split=False))+" )."))
self.target = target self.target = target
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
...@@ -10493,13 +10493,14 @@ class response_operator(operator): ...@@ -10493,13 +10493,14 @@ class response_operator(operator):
## mask ## mask
x_ *= self.mask x_ *= self.mask
## assign ## assign
#return x_[self.assign.tolist()] #return x_[self.assign.T.tolist()]
return field(self.target,val=x_[self.assign.tolist()],target=kwargs.get("target",None)) return field(self.target,val=x_[self.assign.T.tolist()],target=kwargs.get("target",None))
def _adjoint_multiply(self,x,**kwargs): ## > applies the adjoint operator to a given field def _adjoint_multiply(self,x,**kwargs): ## > applies the adjoint operator to a given field
x_ = np.zeros(self.domain.dim(split=True),dtype=self.domain.datatype,order='C') x_ = np.zeros(self.domain.dim(split=True),dtype=self.domain.datatype,order='C')
## assign (transposed) ## assign (transposed)
x_[self.assign.tolist()] += x.val.flatten(order='C') for ii in xrange(np.size(self.assign,axis=0)):
x_[np.array([self.assign[ii]]).T.tolist()] += x[ii]
## mask ## mask
x_ *= self.mask x_ *= self.mask
## smooth ## smooth
......
...@@ -1006,6 +1006,8 @@ class steepest_descent(object): ...@@ -1006,6 +1006,8 @@ class steepest_descent(object):
self.c = c ## 0 < c1 < c2 < 1 self.c = c ## 0 < c1 < c2 < 1
self.note = notification(default=bool(note)) self.note = notification(default=bool(note))
self._alpha = None ## last alpha
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def __call__(self,x0,alpha=1,tol=1E-4,clevel=8,limii=100000): def __call__(self,x0,alpha=1,tol=1E-4,clevel=8,limii=100000):
...@@ -1041,6 +1043,10 @@ class steepest_descent(object): ...@@ -1041,6 +1043,10 @@ class steepest_descent(object):
raise TypeError(about._errors.cstring("ERROR: invalid input.")) raise TypeError(about._errors.cstring("ERROR: invalid input."))
self.x = x0 self.x = x0
## check for exsisting alpha
if(alpha is None)and(self._alpha is not None):
alpha = self._alpha
clevel = max(1,int(clevel)) clevel = max(1,int(clevel))
limii = int(limii) limii = int(limii)
...@@ -1091,6 +1097,10 @@ class steepest_descent(object): ...@@ -1091,6 +1097,10 @@ class steepest_descent(object):
if(self.spam is not None): if(self.spam is not None):
self.spam(self.x,ii) self.spam(self.x,ii)
## memorise last alpha
if(alpha is not None):
self._alpha = alpha/a ## undo update
return self.x,convergence return self.x,convergence
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment