Commit ec40bf59 authored by Marco Selig's avatar Marco Selig

'explicit' update.

parent e241cb20
......@@ -486,7 +486,7 @@ class _about(object): ## nifty support class for global settings
"""
## version
self._version = "0.6.2"
self._version = "0.6.5"
## switches and notifications
self._errors = notification(default=True,ccode=notification._code)
......@@ -3279,12 +3279,17 @@ class rg_space(space):
elif(naxes==2):
if(np.iscomplexobj(x)):
about.infos.cprint("INFO: absolute values, real and imaginary part plotted.")
about.infos.cprint("INFO: absolute values and phases are plotted.")
if(title):
title += " "
self.get_plot(np.absolute(x),title=title+"(absolute)",vmin=vmin,vmax=vmax,power=False,unit=unit,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
self.get_plot(np.real(x),title=title+"(real part)",vmin=vmin,vmax=vmax,power=False,unit=unit,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
self.get_plot(np.imag(x),title=title+"(imaginary part)",vmin=vmin,vmax=vmax,power=False,unit=unit,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
# self.get_plot(np.real(x),title=title+"(real part)",vmin=vmin,vmax=vmax,power=False,unit=unit,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
# self.get_plot(np.imag(x),title=title+"(imaginary part)",vmin=vmin,vmax=vmax,power=False,unit=unit,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
if(unit):
unit = "rad"
if(cmap is None):
cmap = pl.cm.hsv_r
self.get_plot(np.angle(x,deg=False),title=title+"(phase)",vmin=0,vmax=6.28319,power=False,unit=unit,norm=None,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs) ## vmax == 2 pi
else:
if(vmin is None):
vmin = np.min(x,axis=None,out=None)
......@@ -4122,9 +4127,9 @@ class lm_space(space):
ax0.legend()
ax0.set_xlim(xaxes[1],xaxes[-1])
ax0.set_xlabel(r"$l$")
ax0.set_xlabel(r"$\ell$")
ax0.set_ylim(vmin,vmax)
ax0.set_ylabel(r"$l(2l+1) C_l$")
ax0.set_ylabel(r"$\ell(2\ell+1) C_\ell$")
ax0.set_title(title)
else:
......@@ -4133,8 +4138,11 @@ class lm_space(space):
if(title):
title += " "
self.get_plot(np.absolute(x),title=title+"(absolute)",vmin=vmin,vmax=vmax,power=False,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
self.get_plot(np.real(x),title=title+"(real part)",vmin=vmin,vmax=vmax,power=False,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
self.get_plot(np.imag(x),title=title+"(imaginary part)",vmin=vmin,vmax=vmax,power=False,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
# self.get_plot(np.real(x),title=title+"(real part)",vmin=vmin,vmax=vmax,power=False,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
# self.get_plot(np.imag(x),title=title+"(imaginary part)",vmin=vmin,vmax=vmax,power=False,norm=norm,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs)
if(cmap is None):
cmap = pl.cm.hsv_r
self.get_plot(np.angle(x,deg=False),title=title+"(phase)",vmin=0,vmax=6.28319,power=False,norm=None,cmap=cmap,cbar=cbar,other=None,legend=False,**kwargs) ## vmax == 2 pi
else:
if(vmin is None):
vmin = np.min(x,axis=None,out=None)
......@@ -4165,7 +4173,7 @@ class lm_space(space):
sub = ax0.pcolormesh(xaxes,yaxes,np.ma.masked_where(np.isnan(xmesh),xmesh),cmap=cmap,norm=n_,vmin=vmin,vmax=vmax,clim=(vmin,vmax))
ax0.set_xlim(xaxes[0],xaxes[-1])
ax0.set_xticks([0],minor=False)
ax0.set_xlabel(r"$l$")
ax0.set_xlabel(r"$\ell$")
ax0.set_ylim(yaxes[0],yaxes[-1])
ax0.set_yticks([0],minor=False)
ax0.set_ylabel(r"$m$")
......@@ -7698,17 +7706,17 @@ class operator(object):
self.sym = bool(sym)
self.uni = bool(uni)
if(self.domain.discrete):
self.imp = True
else:
self.imp = bool(imp)
if(target is None)or(self.sym)or(self.uni):
target = self.domain
elif(not isinstance(target,space)):
raise TypeError(about._errors.cstring("ERROR: invalid input."))
self.target = target
if(self.domain.discrete)and(self.target.discrete):
self.imp = True
else:
self.imp = bool(imp)
if(para is not None):
self.para = para
......@@ -9023,6 +9031,39 @@ class diagonal_operator(operator):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def det(self):
"""
Computes the determinant of the matrix.
Returns
-------
det : float
The determinant
"""
if(self.domain.dim(split=False)<self.domain.dof()): ## hidden degrees of freedom
return np.exp(self.domain.calc_dot(np.ones(self.domain.dim(split=True),dtype=self.domain.datatype,order='C'),np.log(self.val)))
else:
return np.prod(self.val,axis=None,dtype=None,out=None)
def inverse_det(self):
"""
Computes the determinant of the inverse operator.
Returns
-------
det : float
The determinant
"""
det = self.det()
if(det<>0):
return 1/det
else:
raise ValueError(about._errors.cstring("ERROR: singular operator."))
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def get_random_field(self,domain=None,target=None,**kwargs):
"""
Generates a Gaussian random field with variance equal to the
......@@ -10275,21 +10316,13 @@ class response_operator(operator):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def _multiply(self,x,**kwargs): ## > applies the operator to a given field
# ## weight ## TODO: remove in future version
# if(self.domain.discrete)and(self.den):
# x_ = self.domain.calc_weight(x.val,power=1)
# elif(not self.domain.discrete)and(not self.den):
# x_ = self.domain.calc_weight(x.val,power=-1)
# else:
# x_ = x.val
# ## smooth
# x_ = self.domain.calc_smooth(x_.val,sigma=self.sigma)
## smooth
x_ = self.domain.calc_smooth(x.val,sigma=self.sigma)
## mask
x_ = self.mask*x_
## assign
return x_[self.assign.tolist()]
#return x_[self.assign.tolist()]
return field(self.target,val=x_[self.assign.tolist()],target=kwargs.get("target",None))
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')
......@@ -10299,12 +10332,8 @@ class response_operator(operator):
x_ = self.mask*x_
## smooth
x_ = self.domain.calc_smooth(x_,sigma=self.sigma)
# ## weight ## TODO: remove in future version
# if(self.domain.discrete)and(self.den):
# x_ = self.domain.calc_weight(x_,power=1)
# elif(not self.domain.discrete)and(not self.den):
# x_ = self.domain.calc_weight(x_,power=-1)
return x_
#return x_
return field(self.domain,val=x_,target=kwargs.get("target",None))
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
......@@ -10352,8 +10381,6 @@ class response_operator(operator):
##-----------------------------------------------------------------------------
## IDEA: explicit_operator
......@@ -12064,15 +12091,14 @@ class trace_probing(probing):
else:
domain = op.domain
else:
if(function in [op.inverse_times,op.adjoint_times]):
if(not op.target.check_codomain(domain)): ## restrictive
raise ValueError(about._errors.cstring("ERROR: incompatible domains."))
if(target is None)or(not isinstance(target,space)):
if(not op.domain.check_codomain(domain)): ## restrictive
raise ValueError(about._errors.cstring("ERROR: incompatible domains."))
if(not op.target.check_codomain(domain)): ## restrictive
raise ValueError(about._errors.cstring("ERROR: incompatible domains."))
if(target is None)or(not isinstance(target,space)):
if(function in [op.inverse_times,op.adjoint_times]):
target = op.target
else:
if(not op.domain.check_codomain(domain)): ## restrictive
raise ValueError(about._errors.cstring("ERROR: incompatible domains."))
if(target is None)or(not isinstance(target,space)):
else:
target = op.domain
self.domain = domain
......@@ -12124,6 +12150,8 @@ class trace_probing(probing):
if(f is None):
return None
else:
if(f.domain!=self.domain):
f.transform(target=self.domain,overwrite=True)
return self.domain.calc_dot(probe.val,f.val) ## discrete inner product
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
......@@ -12184,7 +12212,7 @@ class trace_probing(probing):
os.kill()
else:
if(result is not None):
if(np.iscomplexobj(result)):
if(isinstance(_share.sum,tuple)):
_share.sum[0].acquire(block=True,timeout=None)
_share.sum[0].value += np.real(result)
_share.sum[0].release()
......@@ -12454,15 +12482,14 @@ class diagonal_probing(probing):
else:
domain = op.domain
else:
if(function in [op.inverse_times,op.adjoint_times]):
if(not op.target.check_codomain(domain)): ## restrictive
raise ValueError(about._errors.cstring("ERROR: incompatible domains."))
if(target is None)or(not isinstance(target,space)):
if(not op.domain.check_codomain(domain)): ## restrictive
raise ValueError(about._errors.cstring("ERROR: incompatible domains."))
if(not op.target.check_codomain(domain)): ## restrictive
raise ValueError(about._errors.cstring("ERROR: incompatible domains."))
if(target is None)or(not isinstance(target,space)):
if(function in [op.inverse_times,op.adjoint_times]):
target = op.target
else:
if(not op.domain.check_codomain(domain)): ## restrictive
raise ValueError(about._errors.cstring("ERROR: incompatible domains."))
if(target is None)or(not isinstance(target,space)):
else:
target = op.domain
self.domain = domain
......@@ -12594,6 +12621,8 @@ class diagonal_probing(probing):
if(f is None):
return None
else:
if(f.domain!=self.domain):
f.transform(target=self.domain,overwrite=True)
result = np.conjugate(probe.val)*f.val
if(self.save is not None):
np.save(self.save+"%08u"%idnum,result)
......@@ -12609,7 +12638,7 @@ class diagonal_probing(probing):
os.kill()
else:
if(result is not None):
if(np.iscomplexobj(result)):
if(isinstance(_share.sum,tuple)):
_share.sum[0].acquire(block=True,timeout=None)
_share.sum[0][:] += np.real(result)
_share.sum[0].release()
......
......@@ -20,13 +20,13 @@
## along with this program. If not, see <http://www.gnu.org/licenses/>.
"""
.. __ ____ __
.. /__/ / _/ / /_
.. __ ___ __ / /_ / _/ __ __
.. / _ | / / / _/ / / / / / /
.. / / / / / / / / / /_ / /_/ /
.. /__/ /__/ /__/ /__/ \___/ \___ / explicit
.. /______/
.. __ ____ __
.. /__/ / _/ / /_
.. __ ___ __ / /_ / _/ __ __
.. / _ | / / / _/ / / / / / /
.. / / / / / / / / / /_ / /_/ /
.. /__/ /__/ /__/ /__/ \___/ \___ / explicit
.. /______/
TODO: documentation
......@@ -38,14 +38,24 @@ from nifty_core import *
##-----------------------------------------------------------------------------
class matrix_operator(operator):
class explicit_operator(operator):
"""
..
..
.. __ __ __ __
.. / / /__/ /__/ / /_
.. _______ __ __ ______ / / __ _______ __ / _/
.. / __ / \ \/ / / _ | / / / / / ____/ / / / /
.. / /____/ / / / /_/ / / /_ / / / /____ / / / /_
.. \______/ /__/\__\ / ____/ \___/ /__/ \______/ /__/ \___/ operator
.. /__/
TODO: documentation
"""
epsilon = 1E-12 ## absolute precision for comparisons to identity
def __init__(self,domain,matrix,bare=True,sym=None,uni=None,target=None,para=None):
def __init__(self,domain,matrix,bare=True,sym=None,uni=None,target=None):
"""
TODO: documentation
......@@ -58,9 +68,9 @@ class matrix_operator(operator):
self.domain = domain
## check shape
val = np.array(matrix).reshape(-1,self.domain.dim(split=False))
val = np.array(matrix).reshape((-1,self.domain.dim(split=False)))
if(val.size>1048576):
about.warnings.cprint("WARNING: matrix size > 2 ** 20.")
about.infos.cprint("INFO: matrix size > 2 ** 20.")
## check target
if(target is not None):
......@@ -77,21 +87,98 @@ class matrix_operator(operator):
raise TypeError(about._errors.cstring("ERROR: insufficient input."))
self.target = target
## check datatype
if(np.any(np.iscomplex(val))):
datatype,purelyreal = max(min(val.dtype,self.domain.datatype),min(val.dtype,self.target.datatype)),False
else:
datatype,purelyreal = max(min(val.dtype,self.domain.vol.dtype),min(val.dtype,self.target.vol.dtype)),True
## weight if ... (given `domain` and `target`)
if(isinstance(bare,tuple)):
if(len(bare)!=2):
raise ValueError(about._errors.cstring("ERROR: invalid input."))
else:
val = self._calc_weight_rows(val,-bool(bare[0]))
val = self._calc_weight_cols(val,-bool(bare[1]))
elif(not bare):
val = self._calc_weight_rows(val,-1)
if(purelyreal):
self.val = np.real(val).astype(datatype)
else:
self.val = val.astype(datatype)
# if(np.iscomplexobj(val)):
# if(np.all(np.imag(val)==0)):
# val = np.real(val).astype(min(val.dtype,self.domain.vol.dtype,self.target.vol.dtype))
# else:
# val = val.astype(min(val.dtype,self.domain.datatype,self.target.datatype))
# else:
# val = val.astype(min(val.dtype,self.domain.vol.dtype,self.target.vol.dtype))
# ## weight if ... (given `domain` and `target`)
# if(isinstance(bare,tuple)):
# if(len(bare)!=2):
# raise ValueError(about._errors.cstring("ERROR: invalid input."))
# else:
# val = self._calc_weight_rows(val,-bool(bare[0]))
# val = self._calc_weight_cols(val,-bool(bare[1]))
# elif(not bare):
# val = self._calc_weight_rows(val,-1)
# self.val = val
## check hidden degrees of freedom
self._hidden = np.array([self.domain.dim(split=False)<self.domain.dof(),self.target.dim(split=False)<self.target.dof()],dtype=np.bool)
# if(np.any(self._hidden)):
# about.infos.cprint("INFO: inappropriate space.")
## check flags
if(val.shape[0]==val.shape[1]):
self.sym,self.uni = self._check_flags(sym=sym,uni=uni)
if(self.domain.discrete)and(self.target.discrete):
self.imp = True
else:
self.imp = False ## bare matrix is stored for efficiency
self._inv = None ## defined when needed
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def _check_flags(self,sym=None,uni=None): ## > determine `sym` and `uni`
if(self.val.shape[0]==self.val.shape[1]):
if(sym is None):
adj = np.conjugate(val.T)
sym = np.all(np.absolute(val-adj)<self.epsilon)
adj = np.conjugate(self.val.T)
sym = np.all(np.absolute(self.val-adj)<self.epsilon)
if(uni is None):
uni = np.all(np.absolute(np.dot(val,adj)-np.diag(np.ones(len(val),dtype=np.int,order='C'),k=0))<self.epsilon)
uni = np.all(np.absolute(self._calc_mul(adj,0)-np.diag(1/self.target.get_meta_volume(total=False),k=0))<self.epsilon)
elif(uni is None):
adj = np.conjugate(val.T)
uni = np.all(np.absolute(np.dot(val,adj)-np.diag(np.ones(len(val),dtype=np.int,order='C'),k=0))<self.epsilon)
adj = np.conjugate(self.val.T)
uni = np.all(np.absolute(self._calc_mul(adj,0)-np.diag(1/self.target.get_meta_volume(total=False),k=0))<self.epsilon)
return bool(sym),bool(uni)
else:
sym = False
uni = False
self.sym = bool(sym)
self.uni = bool(uni)
return False,False
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def _set_inverse(self): ## > define inverse matrix
if(self._inv is None):
if(np.any(self._hidden)):
about.warnings.cprint("WARNING: inappropriate inversion.")
self._inv = np.linalg.inv(self.weight(rowpower=1,colpower=1,overwrite=False))
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def set_matrix(self,newmatrix,bare=True,sym=None,uni=None):
"""
TODO: documentation
"""
## check shape
val = np.array(newmatrix).reshape((-1,self.domain.dim(split=False)))
if(val.size>1048576):
about.warnings.cprint("WARNING: matrix size > 2 ** 20.")
## check datatype
if(np.iscomplexobj(val)):
......@@ -101,37 +188,38 @@ class matrix_operator(operator):
val = val.astype(min(val.dtype,self.domain.datatype,self.target.datatype))
else:
val = val.astype(min(val.dtype,self.domain.vol.dtype,self.target.vol.dtype))
## weight if ... (given `domain`, `target` and `val`)
## weight if ... (given `domain` and `target`)
if(isinstance(bare,tuple)):
if(len(bare)!=2):
raise ValueError(about._errors.cstring("ERROR: invalid input."))
else:
val = self._calc_weight_rows(val,-bool(bare[0]))
val = self._calc_weight_cols(val,-bool(bare[1]))
val = self._calc_weight_rows(val,power=-bool(bare[0]))
val = self._calc_weight_cols(val,power=-bool(bare[1]))
elif(not bare):
val = self._calc_weight_rows(val,-1)
self.val = val
if(self.domain.discrete)and(self.target.discrete): ## FIXME: operator ???
self.imp = True
else:
self.imp = False
## set parameters
if(para is not None):
self.para = para
self.inv = None ## defined when needed
## check flags
self.sym,self.uni = self._check_flags(sym=sym,uni=uni)
self._inv = None ## reset
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def set_matrix(self,newmatrix,bare=True,sym=None,uni=None,):
def get_matrix(self,bare=True):
"""
TODO: documentation
"""
self.val = newmatrix
self.inv = None
if(bare==True)or(self.imp):
return self.val
## weight if ...
elif(isinstance(bare,tuple)):
if(len(bare)!=2):
raise ValueError(about._errors.cstring("ERROR: invalid input."))
else:
return self.weight(rowpower=bool(bare[0]),colpower=bool(bare[1]),overwrite=False)
elif(not bare):
return self.weight(rowpower=bool(bare),colpower=0,overwrite=False)
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
......@@ -142,10 +230,10 @@ class matrix_operator(operator):
return self.target.calc_weight(x,power=power).flatten(order='C')
def _calc_weight_rows(self,X,power=1): ## > weight all rows
return np.apply_along_axis(self._calc_weight_row,0,X,power)
return np.apply_along_axis(self._calc_weight_row,1,X,power)
def _calc_weight_cols(self,X,power=1): ## > weight all columns
return np.apply_along_axis(self._calc_weight_col,1,X,power)
return np.apply_along_axis(self._calc_weight_col,0,X,power)
def weight(self,rowpower=0,colpower=0,overwrite=False):
"""
......@@ -153,34 +241,46 @@ class matrix_operator(operator):
"""
if(overwrite):
if(rowpower): ## rowpower <> 0
if(not self.domain.discrete)and(rowpower): ## rowpower <> 0
self.val = self._calc_weight_rows(self.val,rowpower)
if(colpower): ## colpower <> 0
if(not self.target.discrete)and(colpower): ## colpower <> 0
self.val = self._calc_weight_cols(self.val,colpower)
else:
X = np.copy(self.val)
if(rowpower): ## rowpower <> 0
if(not self.domain.discrete)and(rowpower): ## rowpower <> 0
X = self._calc_weight_rows(X,rowpower)
if(colpower): ## colpower <> 0
if(not self.target.discrete)and(colpower): ## colpower <> 0
X = self._calc_weight_cols(X,colpower)
return X
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def _multiply(self,x,**kwargs): ## > applies the matirx to a given field
x_ = field(self.target,val=np.dot(self.val,x.val.flatten(order='C'),out=None),target=x.target)
if(self._hidden[0]): ## hidden degrees of freedom
x_ = np.apply_along_axis(self.domain.calc_dot,1,self.val,np.conjugate(x.val))
else:
x_ = np.dot(self.val,x.val.flatten(order='C'),out=None)
return x_
def _adjoint_multiply(self,x,**kwargs): ## > applies the adjoint operator to a given field
x_ = field(self.domain,val=np.dot(np.conjugate(self.val.T),x.val.flatten(order='C'),out=None),target=x.target)
if(self._hidden[1]): ## hidden degrees of freedom
x_ = np.apply_along_axis(self.target.calc_dot,0,np.conjugate(self.val),np.conjugate(x.val))
else:
x_ = np.dot(np.conjugate(self.val.T),x.val.flatten(order='C'),out=None)
return x_
def _inverse_multiply(self,x,**kwargs): ## > applies the inverse operator to a given field
x_ = field(self.domain,val=np.dot(self.inv,x.val.flatten(order='C'),out=None),target=x.target)
if(self._hidden[1]): ## hidden degrees of freedom
x_ = np.apply_along_axis(self.target.calc_dot,1,self._inv,np.conjugate(x.val))
else:
x_ = np.dot(self._inv,x.val.flatten(order='C'),out=None)
return x_
def _inverse_adjoint_multiply(self,x,**kwargs): ## > applies the adjoint inverse operator to a given field
x_ = field(self.target,val=np.dot(np.conjugate(self.inv.T),x.val.flatten(order='C'),out=None),target=x.target)
if(self._hidden[0]): ## hidden degrees of freedom
x_ = np.apply_along_axis(self.domain.calc_dot,0,np.conjugate(self.val),np.conjugate(x.val))
else:
x_ = np.dot(np.conjugate(self._inv.T),x.val.flatten(order='C'),out=None)
return x_
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
......@@ -212,10 +312,9 @@ class matrix_operator(operator):
## check whether square matrix
elif(self.nrow()!=self.ncol()):
raise ValueError(about._errors.cstring("ERROR: inverse ill-defined for "+str(op.nrow())+" x "+str(op.ncol())+" matrices."))
## set inverse if ...
elif(self.inv is None):
self.inv = np.linalg.inv(self.val)
raise ValueError(about._errors.cstring("ERROR: inverse ill-defined for "+str(self.nrow())+" x "+str(self.ncol())+" matrices."))
self._set_inverse()
## prepare
x_ = self._briefing(x,self.target,True)
......@@ -250,7 +349,7 @@ class matrix_operator(operator):
If it is no square matrix.
"""
return self.inverse_adjoint_times(x,**kwargs)
return self._inverse_adjoint_times(x,**kwargs)
def inverse_adjoint_times(self,x,**kwargs):
"""
......@@ -280,18 +379,16 @@ class matrix_operator(operator):
"""
## check whether square matrix
if(self.nrow()!=self.ncol()):
raise ValueError(about._errors.cstring("ERROR: inverse ill-defined for "+str(op.nrow())+" x "+str(op.ncol())+" matrices."))
raise ValueError(about._errors.cstring("ERROR: inverse ill-defined for "+str(self.nrow())+" x "+str(self.ncol())+" matrices."))
## check whether self-adjoint
if(self.sym):
return self.inverse_times(x,**kwargs)
return self._inverse_times(x,**kwargs)
## check whether unitary
if(self.uni):
return self.times(x,**kwargs)
## set inverse if ...
elif(self.inv is None):
self.inv = np.linalg.inv(self.val)
self._set_inverse()
## prepare
x_ = self._briefing(x,self.domain,True)
......@@ -302,15 +399,960 @@ class matrix_operator(operator):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def tr(self,domain=None,**kwargs):
"""
Computes the trace of the operator.
Parameters
----------
domain : space, *optional*
space wherein the probes live (default: self.domain)
target : space, *optional*
space wherein the transform of the probes live
(default: None, applies target of the domain)
random : string, *optional*
Specifies the pseudo random number generator. Valid
options are "pm1" for a random vector of +/-1, or "gau"
for a random vector with entries drawn from a Gaussian
distribution with zero mean and unit variance.
(default: "pm1")
ncpu : int, *optional*
number of used CPUs to use. (default: 2)
nrun : int, *optional*
total number of probes (default: 8)
nper : int, *optional*
number of tasks performed by one process (default: 1)
var : bool, *optional*
Indicates whether to additionally return the probing variance
or not (default: False).
loop : bool, *optional*
Indicates whether or not to perform a loop i.e., to
parallelise (default: False)
##-----------------------------------------------------------------------------
Returns
-------
tr : float
Trace of the operator
delta : float, *optional*
Probing variance of the trace. Returned if `var` is True in
of probing case.
"""
if(self.domain!=self.target):
raise ValueError(about._errors.cstring("ERROR: trace ill-defined."))
if(domain is None)or(domain==self.domain):
diag = self.diag(bare=False,domain=self.domain)
if(self._hidden[0]): ## hidden degrees of freedom
return self.domain.calc_dot(np.ones(self.domain.dim(split=True),dtype=self.domain.datatype,order='C'),diag) ## discrete inner product
else:
return np.sum(diag,axis=None,dtype=None,out=None)
else:
return super(explicit_operator,self).tr(domain=domain,**kwargs) ## probing
def inverse_tr(self,domain=None,**kwargs):
"""
Computes the trace of the inverse operator.
Parameters
----------
domain : space, *optional*
space wherein the probes live (default: self.domain)
target : space, *optional*
space wherein the transform of the probes live
(default: None, applies target of the domain)
random : string, *optional*
Specifies the pseudo random number generator. Valid
options are "pm1" for a random vector of +/-1, or "gau"
for a random vector with entries drawn from a Gaussian
distribution with zero mean and unit variance.
(default: "pm1")
ncpu : int, *optional*
number of used CPUs to use. (default: 2)
nrun : int, *optional*
total number of probes (default: 8)
nper : int, *optional*
number of tasks performed by one process (default: 1)