diff --git a/.gitignore b/.gitignore
index 368117a8f077d1bd83c1623cc95c65a9ecc0d308..15bb48beb44c8d1e41def32f99ca7c13ff1215eb 100644
--- a/.gitignore
+++ b/.gitignore
@@ -13,6 +13,13 @@
.document
build
+operators/*
+!operators/__init__py
+!operators/nifty_explicit.py
+!operators/nifty_operators.py
+!operators/nifty_probing.py
+
+
rg/*
!rg/__init__.py
!rg/fft_rg.py
diff --git a/operators/nifty_explicit.py b/operators/nifty_explicit.py
new file mode 100644
index 0000000000000000000000000000000000000000..099875959750c5f8acb2bdea33c1995985a03361
--- /dev/null
+++ b/operators/nifty_explicit.py
@@ -0,0 +1,2379 @@
+## NIFTY (Numerical Information Field Theory) has been developed at the
+## Max-Planck-Institute for Astrophysics.
+##
+## Copyright (C) 2013 Max-Planck-Society
+##
+## Author: Marco Selig
+## Project homepage: <http://www.mpa-garching.mpg.de/ift/nifty/>
+##
+## This program is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+## See the GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+"""
+ .. __ ____ __
+ .. /__/ / _/ / /_
+ .. __ ___ __ / /_ / _/ __ __
+ .. / _ | / / / _/ / / / / / /
+ .. / / / / / / / / / /_ / /_/ /
+ .. /__/ /__/ /__/ /__/ \___/ \___ / explicit
+ .. /______/
+
+ This module extends NIFTY's versatility to the usage of explicit matrix
+ representations of linear operator by the :py:class:`explicit_operator`.
+ In order to access explicit operators, this module provides the
+ :py:class:`explicit_probing` class and the :py:func:`explicify` function.
+ Those objects are supposed to support the user in solving information field
+ theoretical problems in low (or moderate) dimensions, or in debugging
+ algorithms by studying operators in detail.
+
+"""
+from __future__ import division
+from sys import stdout as so
+import numpy as np
+import pylab as pl
+from matplotlib.colors import LogNorm as ln
+from matplotlib.ticker import LogFormatter as lf
+from multiprocessing import Pool as mp
+from multiprocessing import Value as mv
+from multiprocessing import Array as ma
+
+from nifty.nifty_about import about
+from nifty.nifty_core import space, \
+ field
+from nifty_operators import operator,\
+ diagonal_operator,\
+ identity,\
+ vecvec_operator
+from nifty_probing import probing
+
+
+##-----------------------------------------------------------------------------
+
+class explicit_operator(operator):
+ """
+ ..
+ ..
+ .. __ __ __ __
+ .. / / /__/ /__/ / /_
+ .. _______ __ __ ______ / / __ _______ __ / _/
+ .. / __ / \ \/ / / _ | / / / / / ____/ / / / /
+ .. / /____/ / / / /_/ / / /_ / / / /____ / / / /_
+ .. \______/ /__/\__\ / ____/ \___/ /__/ \______/ /__/ \___/ operator class
+ .. /__/
+
+ NIFTY subclass for explicit linear operators.
+
+ This class essentially supports linear operators with explicit matrix
+ representation in the NIFTY framework. Note that this class is not
+ suited for handling huge dimensionalities.
+
+ Parameters
+ ----------
+ domain : space
+ The space wherein valid arguments live.
+ matrix : {list, array}
+ The matrix representation of the operator, ideally shaped in 2D
+ according to dimensionality of target and domain (default: None).
+ bare : {bool, 2-tuple}, *optional*
+ Whether the matrix entries are `bare` or not
+ (mandatory for the correct incorporation of volume weights)
+ (default: True)
+ sym : bool, *optional*
+ Indicates whether the operator is self-adjoint or not
+ (default: False).
+ uni : bool, *optional*
+ Indicates whether the operator is unitary or not (default: False).
+ target : space, *optional*
+ The space wherein the operator output lives (default: domain).
+
+ See Also
+ --------
+ explicify : A function that returns an explicit oparator given an
+ implicit one.
+
+ Notes
+ -----
+ The ambiguity of `bare` or non-bare diagonal entries is based
+ on the choice of a matrix representation of the operator in
+ question. The naive choice of absorbing the volume weights
+ into the matrix leads to a matrix-vector calculus with the
+ non-bare entries which seems intuitive, though. The choice of
+ keeping matrix entries and volume weights separate deals with the
+ bare entries that allow for correct interpretation of the matrix
+ entries; e.g., as variance in case of an covariance operator.
+
+ Examples
+ --------
+ >>> x_space = rg_space(2) # default `dist` == 0.5
+ >>> A = explicit_operator(x_space, matrix=[[2, 0], [1, 1]], bare=False)
+ >>> A.get_matrix(bare=False)
+ array([[2, 0],
+ [1, 1]])
+ >>> A.get_matrix(bare=True)
+ array([[4, 0],
+ [2, 2]])
+ >>> c = field(x_space, val=[3, 5])
+ >>> A(c).val
+ array([ 6., 8.])
+ >>> A.inverse()
+ <nifty_explicit.explicit_operator>
+ >>> (A * A.inverse()).get_matrix(bare=False) # == identity
+ array([[ 1., 0.],
+ [ 0., 1.]])
+ >>> B = A + diagonal_operator(x_space, diag=2, bare=False)
+ >>> B.get_matrix(bare=False)
+ array([[ 4., 0.],
+ [ 1., 3.]])
+ >>> B(c).val
+ array([ 12., 18.])
+ >>> B.tr()
+ 7.0
+ >>> B.det()
+ 12.000000000000005
+
+ Attributes
+ ----------
+ domain : space
+ The space wherein valid arguments live.
+ val : array
+ An array containing the `bare` matrix entries.
+ sym : bool
+ Indicates whether the operator is self-adjoint or not
+ uni : bool
+ Indicates whether the operator is unitary or not
+ imp : bool
+ Indicates whether the incorporation of volume weights in
+ multiplications is already implemented in the `multiply`
+ instance methods or not.
+ target : space
+ The space wherein the operator output lives.
+ _hidden : array
+ An array specifying hidden degrees of freedom in domain and target.
+ _inv : array
+ An array containing the inverse matrix; set when needed.
+
+ """
+ epsilon = 1E-12 ## absolute precision for comparisons to identity
+
+ def __init__(self,domain,matrix=None,bare=True,sym=None,uni=None,target=None):
+ """
+ Initializes the explicit operator and sets the standard operator
+ properties as well as `values`.
+
+ Parameters
+ ----------
+ domain : space
+ The space wherein valid arguments live.
+ matrix : {list, array}
+ The matrix representation of the operator, ideally shaped in 2D
+ according to dimensionality of target and domain (default: None).
+ bare : {bool, 2-tuple}, *optional*
+ Whether the matrix entries are `bare` or not
+ (mandatory for the correct incorporation of volume weights)
+ (default: True).
+ sym : bool, *optional*
+ Indicates whether the operator is self-adjoint or not
+ (default: False).
+ uni : bool, *optional*
+ Indicates whether the operator is unitary or not (default: False).
+ target : space, *optional*
+ The space wherein the operator output lives (default: domain).
+
+ Notes
+ -----
+ The ambiguity of `bare` or non-bare diagonal entries is based
+ on the choice of a matrix representation of the operator in
+ question. The naive choice of absorbing the volume weights
+ into the matrix leads to a matrix-vector calculus with the
+ non-bare entries which seems intuitive, though. The choice of
+ keeping matrix entries and volume weights separate deals with the
+ bare entries that allow for correct interpretation of the matrix
+ entries; e.g., as variance in case of an covariance operator.
+
+ Raises
+ ------
+ TypeError
+ If invalid input is given.
+ ValueError
+ If dimensions of `domain`, `target`, and `matrix` mismatch;
+ or if `bare` is invalid.
+
+ """
+ ## check domain
+ if(not isinstance(domain,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ self.domain = domain
+
+ ## check matrix and target
+ if(matrix is None):
+ if(target is None):
+ val = np.zeros((self.domain.dim(split=False),self.domain.dim(split=False)),dtype=np.int,order='C')
+ target = self.domain
+ else:
+ if(not isinstance(target,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ elif(target!=self.domain):
+ sym = False
+ uni = False
+ val = np.zeros((target.dim(split=False),self.domain.dim(split=False)),dtype=np.int,order='C')
+ elif(np.size(matrix,axis=None)%self.domain.dim(split=False)==0):
+ val = np.array(matrix).reshape((-1,self.domain.dim(split=False)))
+ if(target is not None):
+ if(not isinstance(target,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ elif(val.shape[0]!=target.dim(split=False)):
+ raise ValueError(about._errors.cstring("ERROR: dimension mismatch ( "+str(val.shape[0])+" <> "+str(target.dim(split=False))+" )."))
+ elif(target!=self.domain):
+ sym = False
+ uni = False
+ elif(val.shape[0]==val.shape[1]):
+ target = self.domain
+ else:
+ raise TypeError(about._errors.cstring("ERROR: insufficient input."))
+ else:
+ raise ValueError(about._errors.cstring("ERROR: dimension mismatch ( "+str(np.size(matrix,axis=None))+" <> "+str(self.domain.dim(split=False))+" )."))
+# if(val.size>1048576):
+# about.infos.cprint("INFO: matrix size > 2 ** 20.")
+ 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,power=-int(not bare[0]))
+ val = self._calc_weight_cols(val,power=-int(not 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)
+
+ ## 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
+ 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(self.val.T)
+ sym = np.all(np.absolute(self.val-adj)<self.epsilon)
+ if(uni is None):
+ uni = np.all(np.absolute(self._calc_mul(adj,0)-np.diag(1/self.target.get_meta_volume(total=False).flatten(order='C'),k=0))<self.epsilon)
+ elif(uni is None):
+ 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).flatten(order='C'),k=0))<self.epsilon)
+ return bool(sym),bool(uni)
+ else:
+ 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 cast_domain(self,newdomain):
+ """
+ Casts the domain of the operator.
+
+ Parameters
+ ----------
+ newdomain : space
+ New space wherein valid argument lives.
+
+ Returns
+ -------
+ None
+
+ Raises
+ ------
+ TypeError
+ If `newdomain` is no instance of the nifty_core.space class
+ ValueError
+ If `newdomain` does not match the (unsplit) dimensionality of
+ the current domain.
+
+ """
+ if(not isinstance(newdomain,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ elif(newdomain.dim(split=False)!=self.domain.dim(split=False)):
+ raise ValueError(about._errors.cstring("ERROR: dimension mismatch ( "+str(newdomain.dim(split=False))+" <> "+str(self.domain.dim(split=False))+" )."))
+ self.domain = newdomain
+
+ def cast_target(self,newtarget):
+ """
+ Casts the target of the operator.
+
+ Parameters
+ ----------
+ newtarget : space
+ New space wherein the operator output lives.
+
+ Returns
+ -------
+ None
+
+ Raises
+ ------
+ TypeError
+ If `newtarget` is no instance of the nifty_core.space class
+ ValueError
+ If `newtarget` does not match the (unsplit) dimensionality of
+ the current target.
+
+ """
+ if(not isinstance(newtarget,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ elif(newtarget.dim(split=False)!=self.target.dim(split=False)):
+ raise ValueError(about._errors.cstring("ERROR: dimension mismatch ( "+str(newtarget.dim(split=False))+" <> "+str(self.target.dim(split=False))+" )."))
+ self.target = newtarget
+
+ def cast_spaces(self,newdomain=None,newtarget=None):
+ """
+ Casts the domain and/or the target of the operator.
+
+ Parameters
+ ----------
+ newdomain : space, *optional*
+ New space wherein valid argument lives (default: None).
+ newtarget : space, *optional*
+ New space wherein the operator output lives (default: None).
+
+ Returns
+ -------
+ None
+
+ """
+ if(newdomain is not None):
+ self.cast_domain(newdomain)
+ if(newtarget is not None):
+ self.cast_target(newtarget)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def set_matrix(self,newmatrix,bare=True,sym=None,uni=None):
+ """
+ Resets the entire matrix.
+
+ Parameters
+ ----------
+ matrix : {list, array}
+ New matrix representation of the operator, ideally shaped in 2D
+ according to dimensionality of target and domain (default: None).
+ bare : {bool, 2-tuple}, *optional*
+ Whether the new matrix entries are `bare` or not
+ (mandatory for the correct incorporation of volume weights)
+ (default: True).
+ sym : bool, *optional*
+ Indicates whether the new operator is self-adjoint or not
+ (default: False).
+ uni : bool, *optional*
+ Indicates whether the new operator is unitary or not
+ (default: False).
+
+ Notes
+ -----
+ The ambiguity of `bare` or non-bare diagonal entries is based
+ on the choice of a matrix representation of the operator in
+ question. The naive choice of absorbing the volume weights
+ into the matrix leads to a matrix-vector calculus with the
+ non-bare entries which seems intuitive, though. The choice of
+ keeping matrix entries and volume weights separate deals with the
+ bare entries that allow for correct interpretation of the matrix
+ entries; e.g., as variance in case of an covariance operator.
+
+ Returns
+ -------
+ None
+
+ Raises
+ ------
+ ValueError
+ If matrix' dimensions mismatch;
+ or if `bare` is invalid.
+
+ """
+ ## check matrix
+ if(np.size(newmatrix,axis=None)==self.domain.dim(split=False)*self.target.dim(split=False)):
+ val = np.array(newmatrix).reshape((self.target.dim(split=False),self.domain.dim(split=False)))
+ if(self.target!=self.domain):
+ sym = False
+ uni = False
+# if(val.size>1048576):
+# about.infos.cprint("INFO: matrix size > 2 ** 20.")
+ else:
+ raise ValueError(about._errors.cstring("ERROR: dimension mismatch ( "+str(np.size(newmatrix,axis=None))+" <> "+str(self.nrow())+" x "+str(self.ncol())+" )."))
+
+ ## 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,power=-int(not bare[0]))
+ val = self._calc_weight_cols(val,power=-int(not 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)
+
+ ## check flags
+ self.sym,self.uni = self._check_flags(sym=sym,uni=uni)
+ self._inv = None ## reset
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def get_matrix(self,bare=True):
+ """
+ Returns the entire matrix.
+
+ Parameters
+ ----------
+ bare : {bool, 2-tuple}, *optional*
+ Whether the returned matrix entries are `bare` or not
+ (default: True).
+
+ Notes
+ -----
+ The ambiguity of `bare` or non-bare diagonal entries is based
+ on the choice of a matrix representation of the operator in
+ question. The naive choice of absorbing the volume weights
+ into the matrix leads to a matrix-vector calculus with the
+ non-bare entries which seems intuitive, though. The choice of
+ keeping matrix entries and volume weights separate deals with the
+ bare entries that allow for correct interpretation of the matrix
+ entries; e.g., as variance in case of an covariance operator.
+
+ Returns
+ -------
+ matrix : numpy.array
+ The matrix representation of the operator, shaped in 2D
+ according to dimensionality of target and domain.
+
+ Raises
+ ------
+ ValueError
+ If `bare` is invalid.
+
+ """
+ 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=int(not bare[0]),colpower=int(not bare[1]),overwrite=False)
+ elif(not bare):
+ return self.weight(rowpower=int(not bare),colpower=0,overwrite=False)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def _calc_weight_row(self,x,power): ## > weight row and flatten
+ return self.domain.calc_weight(x,power=power).flatten(order='C')
+
+ def _calc_weight_col(self,x,power): ## > weight column and flatten
+ return self.target.calc_weight(x,power=power).flatten(order='C')
+
+ def _calc_weight_rows(self,X,power=1): ## > weight all rows
+ if(np.any(np.iscomplex(X)))and(not issubclass(self.domain.datatype,np.complexfloating)):
+ return (np.apply_along_axis(self._calc_weight_row,1,np.real(X),power)+np.apply_along_axis(self._calc_weight_row,1,np.imag(X),power)*1j)
+ else:
+ return np.apply_along_axis(self._calc_weight_row,1,X,power)
+
+ def _calc_weight_cols(self,X,power=1): ## > weight all columns
+ if(np.any(np.iscomplex(X)))and(not issubclass(self.target.datatype,np.complexfloating)):
+ return (np.apply_along_axis(self._calc_weight_col,0,np.real(X),power)+np.apply_along_axis(self._calc_weight_col,0,np.imag(X),power)*1j)
+ else:
+ return np.apply_along_axis(self._calc_weight_col,0,X,power)
+
+ def weight(self,rowpower=0,colpower=0,overwrite=False):
+ """
+ Returns the entire matrix, weighted with the volume factors to a
+ given power. The matrix entries will optionally be overwritten.
+
+ Parameters
+ ----------
+ rowpower : scalar, *optional*
+ Specifies the power of the volume factors applied to the rows
+ of the matrix (default: 0).
+ rowpower : scalar, *optional*
+ Specifies the power of the volume factors applied to the columns
+ of the matrix (default: 0).
+ overwrite : bool, *optional*
+ Whether to overwrite the matrix or not (default: False).
+
+ Returns
+ -------
+ field : field, *optional*
+ If overwrite is ``False``, the weighted matrix is returned;
+ otherwise, nothing is returned.
+
+ """
+ if(overwrite):
+ if(not self.domain.discrete)and(rowpower): ## rowpower <> 0
+ self.val = self._calc_weight_rows(self.val,rowpower)
+ 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(not self.domain.discrete)and(rowpower): ## rowpower <> 0
+ X = self._calc_weight_rows(X,rowpower)
+ 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
+ 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
+ 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
+ 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
+ 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_
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def _briefing(self,x,domain,inverse): ## > prepares x for `multiply`
+ ## inspect x
+ if(not isinstance(x,field)):
+ x_ = field(domain,val=x,target=None)
+ else:
+ ## check x.domain
+ if(domain==x.domain):
+ x_ = x
+ ## transform
+ else:
+ x_ = x.transform(target=domain,overwrite=False)
+ ## weight if ...
+ if(not self.imp)and(not domain.discrete):
+ x_ = x_.weight(power=1,overwrite=False)
+ return x_
+
+ def _debriefing(self,x,x_,target,inverse): ## > evaluates x and x_ after `multiply`
+ if(x_ is None):
+ return None
+ else:
+ ## inspect x_
+ if(not isinstance(x_,field)):
+ x_ = field(target,val=x_,target=None)
+ elif(x_.domain!=target):
+ raise ValueError(about._errors.cstring("ERROR: invalid output domain."))
+ ## inspect x
+ if(isinstance(x,field)):
+ ## repair ...
+ if(self.domain==self.target!=x.domain):
+ x_ = x_.transform(target=x.domain,overwrite=False) ## ... domain
+ if(x_.domain==x.domain)and(x_.target!=x.target):
+ x_.set_target(newtarget=x.target) ## ... codomain
+ return x_
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def inverse_times(self,x,**kwargs):
+ """
+ Applies the inverse operator to a given object.
+
+ Parameters
+ ----------
+ x : {scalar, ndarray, field}
+ Scalars are interpreted as constant arrays, and an array will
+ be interpreted as a field on the domain space of the operator.
+
+ Returns
+ -------
+ OIx : field
+ Mapped field on the target space of the operator.
+
+ Raises
+ ------
+ ValueError
+ If it is no square matrix.
+
+ """
+ ## check whether self-inverse
+ if(self.sym)and(self.uni):
+ return self.times(x,**kwargs)
+
+ ## check whether square matrix
+ elif(self.nrow()!=self.ncol()):
+ 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)
+ ## apply operator
+ x_ = self._inverse_multiply(x_,**kwargs)
+ ## evaluate
+ return self._debriefing(x,x_,self.domain,True)
+
+ def adjoint_inverse_times(self,x,**kwargs):
+ """
+ Applies the inverse adjoint operator to a given object.
+
+ Parameters
+ ----------
+ x : {scalar, ndarray, field}
+ Scalars are interpreted as constant arrays, and an array will
+ be interpreted as a field on the target space of the operator.
+
+ Returns
+ -------
+ OAIx : field
+ Mapped field on the domain of the operator.
+
+ Notes
+ -----
+ For linear operators represented by square matrices, inversion and
+ adjungation and inversion commute.
+
+ Raises
+ ------
+ ValueError
+ If it is no square matrix.
+
+ """
+ return self.inverse_adjoint_times(x,**kwargs)
+
+ def inverse_adjoint_times(self,x,**kwargs):
+ """
+ Applies the adjoint inverse operator to a given object.
+
+ Parameters
+ ----------
+ x : {scalar, ndarray, field}
+ Scalars are interpreted as constant arrays, and an array will
+ be interpreted as a field on the target space of the operator.
+
+ Returns
+ -------
+ OIAx : field
+ Mapped field on the domain of the operator.
+
+ Notes
+ -----
+ For linear operators represented by square matrices, inversion and
+ adjungation and inversion commute.
+
+ Raises
+ ------
+ ValueError
+ If it is no square matrix.
+
+ """
+ ## check whether square matrix
+ if(self.nrow()!=self.ncol()):
+ 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)
+ ## check whether unitary
+ if(self.uni):
+ return self.times(x,**kwargs)
+
+ self._set_inverse()
+
+ ## prepare
+ x_ = self._briefing(x,self.domain,True)
+ ## apply operator
+ x_ = self._inverse_adjoint_multiply(x_,**kwargs)
+ ## evaluate
+ return self._debriefing(x,x_,self.target,True)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ 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).
+ random : string, *optional*
+ Specifies the pseudo random number generator (default: "pm1");
+ supported distributions are:
+
+ - "pm1" (uniform distribution over {+1,-1} or {+1,+i,-1,-i}
+ - "gau" (normal distribution with zero-mean and a given standard deviation or variance)
+
+ ncpu : integer, *optional*
+ Number of used CPUs to use (default: 2).
+ nrun : integer, *optional*
+ Total number of probes (default: 8).
+ nper : integer, *optional*
+ Number of tasks performed by one worker (default: 1).
+ var : bool, *optional*
+ Whether to additionally return the probing variance or not
+ (default: False).
+ loop : bool, *optional*
+ Whether to perform a loop or 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
+ case of probing.
+
+ See Also
+ --------
+ trace_probing
+
+ Raises
+ ------
+ ValueError
+ If `domain` and `target` are unequal.
+
+ """
+ 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).
+ random : string, *optional*
+ Specifies the pseudo random number generator (default: "pm1");
+ supported distributions are:
+
+ - "pm1" (uniform distribution over {+1,-1} or {+1,+i,-1,-i}
+ - "gau" (normal distribution with zero-mean and a given standard deviation or variance)
+
+ ncpu : integer, *optional*
+ Number of used CPUs to use (default: 2).
+ nrun : integer, *optional*
+ Total number of probes (default: 8).
+ nper : integer, *optional*
+ Number of tasks performed by one worker (default: 1).
+ var : bool, *optional*
+ Whether to additionally return the probing variance or not
+ (default: False).
+ loop : bool, *optional*
+ Whether to perform a loop or 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
+ case of probing.
+
+ See Also
+ --------
+ trace_probing
+
+ Raises
+ ------
+ ValueError
+ If `domain` and `target` are unequal.
+
+ """
+ if(self.domain!=self.target):
+ raise ValueError(about._errors.cstring("ERROR: trace ill-defined."))
+
+ if(domain is None)or(domain==self.domain):
+ diag = self.inverse_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 diag(self,bare=False,domain=None,**kwargs):
+ """
+ Computes the diagonal of the operator.
+
+ Parameters
+ ----------
+ bare : bool, *optional*
+ Whether the returned diagonal entries are `bare` or not
+ (default: True).
+ domain : space, *optional*
+ Space wherein the probes live (default: self.domain).
+ target : space, *optional*
+ Space wherein the transform of the probes live (default: None).
+ random : string, *optional*
+ Specifies the pseudo random number generator (default: "pm1");
+ supported distributions are:
+
+ - "pm1" (uniform distribution over {+1,-1} or {+1,+i,-1,-i}
+ - "gau" (normal distribution with zero-mean and a given standard deviation or variance)
+
+ ncpu : integer, *optional*
+ Number of used CPUs to use (default: 2).
+ nrun : integer, *optional*
+ Total number of probes (default: 8).
+ nper : integer, *optional*
+ Number of tasks performed by one worker (default: 1).
+ var : bool, *optional*
+ Whether to additionally return the probing variance or not
+ (default: False).
+ save : bool, *optional*
+ Whether all individual probing results are saved or not
+ (default: False).
+ path : string, *optional*
+ Path wherein the results are saved (default: "tmp").
+ prefix : string, *optional*
+ Prefix for all saved files (default: "")
+ loop : bool, *optional*
+ Whether to perform a loop or to parallelise (default: False).
+
+ Returns
+ -------
+ diag : ndarray
+ Diagonal of the operator.
+ delta : float, *optional*
+ Probing variance of the trace; returned if `var` is ``True`` in
+ case of probing.
+
+ Notes
+ -----
+ The ambiguity of `bare` or non-bare diagonal entries is based
+ on the choice of a matrix representation of the operator in
+ question. The naive choice of absorbing the volume weights
+ into the matrix leads to a matrix-vector calculus with the
+ non-bare entries which seems intuitive, though. The choice of
+ keeping matrix entries and volume weights separate deals with the
+ bare entries that allow for correct interpretation of the matrix
+ entries; e.g., as variance in case of an covariance operator.
+
+ See Also
+ --------
+ diagonal_probing
+
+ Raises
+ ------
+ ValueError
+ If it is no square matrix.
+
+ """
+ if(self.val.shape[0]!=self.val.shape[1]):
+ raise ValueError(about._errors.cstring("ERROR: diagonal ill-defined for "+str(self.val.shape[0])+" x "+str(self.val.shape[1])+" matrices."))
+ if(self.domain!=self.target)and(not bare):
+ about.warnings.cprint("WARNING: ambiguous non-bare diagonal.")
+
+ if(domain is None)or(domain==self.domain):
+ diag = np.diagonal(self.val,offset=0,axis1=0,axis2=1)
+ ## weight if ...
+ if(not self.domain.discrete)and(not bare):
+ diag = self.domain.calc_weight(diag,power=1)
+ ## check complexity
+ if(np.all(np.imag(diag)==0)):
+ diag = np.real(diag)
+ return diag
+ elif(domain==self.target):
+ diag = np.diagonal(np.conjugate(self.val.T),offset=0,axis1=0,axis2=1)
+ ## weight if ...
+ if(not self.target.discrete)and(not bare):
+ diag = self.target.calc_weight(diag,power=1)
+ ## check complexity
+ if(np.all(np.imag(diag)==0)):
+ diag = np.real(diag)
+ return diag
+ else:
+ return super(explicit_operator,self).diag(bare=bare,domain=domain,**kwargs) ## probing
+
+ def inverse_diag(self,bare=False,domain=None,**kwargs):
+ """
+ Computes the diagonal of the inverse operator.
+
+ Parameters
+ ----------
+ bare : bool, *optional*
+ Whether the returned diagonal entries are `bare` or not
+ (default: True).
+ domain : space, *optional*
+ Space wherein the probes live (default: self.domain).
+ target : space, *optional*
+ Space wherein the transform of the probes live (default: None).
+ random : string, *optional*
+ Specifies the pseudo random number generator (default: "pm1");
+ supported distributions are:
+
+ - "pm1" (uniform distribution over {+1,-1} or {+1,+i,-1,-i}
+ - "gau" (normal distribution with zero-mean and a given standard deviation or variance)
+
+ ncpu : integer, *optional*
+ Number of used CPUs to use (default: 2).
+ nrun : integer, *optional*
+ Total number of probes (default: 8).
+ nper : integer, *optional*
+ Number of tasks performed by one worker (default: 1).
+ var : bool, *optional*
+ Whether to additionally return the probing variance or not
+ (default: False).
+ save : bool, *optional*
+ Whether all individual probing results are saved or not
+ (default: False).
+ path : string, *optional*
+ Path wherein the results are saved (default: "tmp").
+ prefix : string, *optional*
+ Prefix for all saved files (default: "")
+ loop : bool, *optional*
+ Whether to perform a loop or to parallelise (default: False).
+
+ Returns
+ -------
+ diag : ndarray
+ Diagonal of the inverse operator.
+ delta : float, *optional*
+ Probing variance of the trace; returned if `var` is ``True`` in
+ case of probing.
+
+ Notes
+ -----
+ The ambiguity of `bare` or non-bare diagonal entries is based
+ on the choice of a matrix representation of the operator in
+ question. The naive choice of absorbing the volume weights
+ into the matrix leads to a matrix-vector calculus with the
+ non-bare entries which seems intuitive, though. The choice of
+ keeping matrix entries and volume weights separate deals with the
+ bare entries that allow for correct interpretation of the matrix
+ entries; e.g., as variance in case of an covariance operator.
+
+ See Also
+ --------
+ diagonal_probing
+
+ Raises
+ ------
+ ValueError
+ If it is no square matrix.
+
+ """
+ if(self.val.shape[0]!=self.val.shape[1]):
+ raise ValueError(about._errors.cstring("ERROR: inverse ill-defined for "+str(self.val.shape[0])+" x "+str(self.val.shape[1])+" matrices."))
+ if(self.domain!=self.target)and(not bare):
+ about.warnings.cprint("WARNING: ambiguous non-bare diagonal.")
+
+ if(domain is None)or(domain==self.target):
+ self._set_inverse()
+ diag = np.diagonal(self._inv,offset=0,axis1=0,axis2=1)
+ ## weight if ...
+ if(not self.target.discrete)and(not bare):
+ diag = self.target.calc_weight(diag,power=1)
+ ## check complexity
+ if(np.all(np.imag(diag)==0)):
+ diag = np.real(diag)
+ return diag
+ elif(domain==self.domain):
+ self._set_inverse()
+ diag = np.diagonal(np.conjugate(self._inv.T),offset=0,axis1=0,axis2=1)
+ ## weight if ...
+ if(not self.domain.discrete)and(not bare):
+ diag = self.domain.calc_weight(diag,power=1)
+ ## check complexity
+ if(np.all(np.imag(diag)==0)):
+ diag = np.real(diag)
+ return diag
+ else:
+ return super(explicit_operator,self).inverse_diag(bare=bare,domain=domain,**kwargs) ## probing
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def det(self):
+ """
+ Computes the determinant of the operator.
+
+ Parameters
+ ----------
+ None
+
+ Returns
+ -------
+ det : float
+ Determinant of the operator.
+
+ Raises
+ ------
+ ValueError
+ If `domain` and `target` are unequal.
+
+ """
+ if(self.domain!=self.target):
+ raise ValueError(about._errors.cstring("ERROR: determinant ill-defined."))
+
+ if(np.any(self._hidden)):
+ about.warnings.cprint("WARNING: inappropriate determinant calculation.")
+ det = np.linalg.det(self.weight(rowpower=0.5,colpower=0.5,overwrite=False))
+ if(np.isreal(det)):
+ return np.asscalar(np.real(det))
+ else:
+ return det
+
+ def inverse_det(self):
+ """
+ Computes the determinant of the inverse operator.
+
+ Parameters
+ ----------
+ None
+
+ Returns
+ -------
+ det : float
+ Determinant of the inverse operator.
+
+ Raises
+ ------
+ ValueError
+ If it is a singular matrix.
+
+ """
+ det = self.det()
+ if(det<>0):
+ return 1/det
+ else:
+ raise ValueError(about._errors.cstring("ERROR: singular matrix."))
+
+ def log_det(self):
+ """
+ Computes the logarithm of the determinant of the operator
+ (if applicable).
+
+ Returns
+ -------
+ logdet : float
+ The logarithm of the determinant
+
+ See Also
+ --------
+ numpy.linalg.slogdet
+
+ Raises
+ ------
+ ValueError
+ If `domain` and `target` are unequal or it is non-positive
+ definite matrix.
+
+ """
+ if(self.domain!=self.target):
+ raise ValueError(about._errors.cstring("ERROR: determinant ill-defined."))
+
+ if(np.any(self._hidden)):
+ about.warnings.cprint("WARNING: inappropriate determinant calculation.")
+ sign,logdet = np.linalg.slogdet(self.weight(rowpower=0.5,colpower=0.5,overwrite=False))
+ if(abs(sign)<0.1): ## abs(sign) << 1
+ raise ValueError(about._errors.cstring("ERROR: singular matrix."))
+ if(sign==-1):
+ raise ValueError(about._errors.cstring("ERROR: non-positive definite matrix."))
+ else:
+ logdet += np.log(sign)
+ if(np.isreal(logdet)):
+ return np.asscalar(np.real(logdet))
+ else:
+ return logdet
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def __len__(self):
+ return int(self.nrow()[0])
+
+ def __getitem__(self,key):
+ return self.val[key]
+
+ def __setitem__(self,key,value):
+ self.val[key] = self.val.dtype(value)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def __pos__(self):
+ return explicit_operator(self.domain,self.val,bare=True,sym=self.sym,uni=self.uni,target=self.target)
+
+ def __neg__(self):
+ return explicit_operator(self.domain,-self.val,bare=True,sym=self.sym,uni=self.uni,target=self.target)
+
+ def __abs__(self):
+ return explicit_operator(self.domain,np.absolute(self.val),bare=True,sym=self.sym,uni=self.uni,target=self.target)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def transpose(self):
+ """
+ Computes the transposed matrix.
+
+ Returns
+ -------
+ T : explicit_operator
+ The transposed matrix.
+
+ """
+ return explicit_operator(self.target,self.val.T,bare=True,sym=self.sym,uni=self.uni,target=self.domain)
+
+ def conjugate(self):
+ """
+ Computes the complex conjugated matrix.
+
+ Returns
+ -------
+ CC : explicit_operator
+ The complex conjugated matrix.
+
+ """
+ return explicit_operator(self.domain,np.conjugate(self.val),bare=True,sym=self.sym,uni=self.uni,target=self.target)
+
+ def adjoint(self):
+ """
+ Computes the adjoint matrix.
+
+ Returns
+ -------
+ A : explicit_operator
+ The adjoint matrix.
+
+ """
+ return explicit_operator(self.target,np.conjugate(self.val.T),bare=True,sym=self.sym,uni=self.uni,target=self.domain)
+
+ def inverse(self):
+ """
+ Computes the inverted matrix.
+
+ Returns
+ -------
+ I : explicit_operator
+ The inverted matrix.
+
+ Raises
+ ------
+ ValueError
+ If it is no square matrix.
+
+ """
+ ## check whether square matrix
+ if(self.nrow()!=self.ncol()):
+ raise ValueError(about._errors.cstring("ERROR: inverse ill-defined for "+str(self.nrow())+" x "+str(self.ncol())+" matrices."))
+ self._set_inverse()
+ return explicit_operator(self.target,self._inv,bare=True,sym=self.sym,uni=self.uni,target=self.domain)
+
+ def adjoint_inverse(self):
+ """
+ Computes the adjoint inverted matrix.
+
+ Returns
+ -------
+ AI : explicit_operator
+ The adjoint inverted matrix.
+
+ Notes
+ -----
+ For linear operators represented by square matrices, inversion and
+ adjungation and inversion commute.
+
+ Raises
+ ------
+ ValueError
+ If it is no square matrix.
+
+ """
+ return self.inverse_adjoint()
+
+ def inverse_adjoint(self):
+ """
+ Computes the inverted adjoint matrix.
+
+ Returns
+ -------
+ IA : explicit_operator
+ The inverted adjoint matrix.
+
+ Notes
+ -----
+ For linear operators represented by square matrices, inversion and
+ adjungation and inversion commute.
+
+ Raises
+ ------
+ ValueError
+ If it is no square matrix.
+
+ """
+ ## check whether square matrix
+ if(self.nrow()!=self.ncol()):
+ raise ValueError(about._errors.cstring("ERROR: inverse ill-defined for "+str(self.nrow())+" x "+str(self.ncol())+" matrices."))
+ self._set_inverse()
+ return explicit_operator(self.target,np.conjugate(self._inv.T),bare=True,sym=self.sym,uni=self.uni,target=self.domain)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def __add__(self,X): ## __add__ : self + X
+ if(isinstance(X,operator)):
+ if(self.domain!=X.domain):
+ raise ValueError(about._errors.cstring("ERROR: inequal domains."))
+ sym = (self.sym and X.sym)
+ uni = None
+ if(isinstance(X,explicit_operator)):
+ if(self.target!=X.target):
+ raise ValueError(about._errors.cstring("ERROR: inequal codomains."))
+ matrix = self.val+X.val
+ elif(isinstance(X,diagonal_operator)):
+ if(self.target.dim(split=False)!=X.target.dim(split=False))or(not self.target.check_codomain(X.target)):
+ raise ValueError(about._errors.cstring("ERROR: incompatible codomains."))
+ matrix = self.val+np.diag(X.diag(bare=True,domain=None),k=0) ## domain == X.domain
+ elif(isinstance(X,vecvec_operator)):
+ if(self.target!=X.target):
+ raise ValueError(about._errors.cstring("ERROR: inequal codomains."))
+ matrix = self.val+np.tensordot(X.val,X.val,axes=0)
+ else:
+ raise TypeError(about._errors.cstring("ERROR: unsupported or incompatible operator."))
+ elif(np.size(X)==1):
+ if(self.nrow()!=self.ncol()):
+ raise ValueError(about._errors.cstring("ERROR: identity ill-defined for "+str(self.nrow())+" x "+str(self.ncol())+" matrices."))
+ sym = self.sym
+ uni = None
+ X = X*np.ones(self.domain.dim(split=False),dtype=np.int,order='C')
+ X = self.domain.calc_weight(X,power=-1).astype(max(min(X.dtype,self.domain.datatype),min(X.dtype,self.target.datatype)))
+ matrix = self.val+np.diag(X.flatten(order='C'),k=0)
+ elif(np.size(X)==np.size(self.val)):
+ sym = None
+ uni = None
+ X = np.array(X).reshape(self.val.shape)
+ if(np.all(np.isreal(X))):
+ X = np.real(X)
+ X = X.astype(max(min(X.dtype,self.domain.datatype),min(X.dtype,self.target.datatype)))
+ matrix = self.val+X
+ else:
+ raise ValueError(about._errors.cstring("ERROR: dimension mismatch ( "+str(np.size(X))+" <> "+str(self.nrow())+" x "+str(self.ncol())+" )."))
+ return explicit_operator(self.domain,matrix,bare=True,sym=sym,uni=uni,target=self.target)
+
+ __radd__ = __add__ ## __add__ : X + self
+
+ def __iadd__(self,X): ## __iadd__ : self += X
+ if(isinstance(X,operator)):
+ if(self.domain!=X.domain):
+ raise ValueError(about._errors.cstring("ERROR: inequal domains."))
+ self.sym = (self.sym and X.sym)
+ self.uni = None
+ if(isinstance(X,explicit_operator)):
+ if(self.target!=X.target):
+ raise ValueError(about._errors.cstring("ERROR: inequal codomains."))
+ X = X.val
+ elif(isinstance(X,diagonal_operator)):
+ if(self.target.dim(split=False)!=X.target.dim(split=False))or(not self.target.check_codomain(X.target)):
+ raise ValueError(about._errors.cstring("ERROR: incompatible codomains."))
+ X = np.diag(X.diag(bare=True,domain=None),k=0) ## domain == X.domain
+ elif(isinstance(X,vecvec_operator)):
+ if(self.target!=X.target):
+ raise ValueError(about._errors.cstring("ERROR: inequal codomains."))
+ X = np.tensordot(X.val,X.val,axes=0)
+ else:
+ raise TypeError(about._errors.cstring("ERROR: unsupported or incompatible operator."))
+ elif(np.size(X)==1):
+ if(self.nrow()!=self.ncol()):
+ raise ValueError(about._errors.cstring("ERROR: identity ill-defined for "+str(self.nrow())+" x "+str(self.ncol())+" matrices."))
+ self.uni = None
+ X = X*np.ones(self.domain.dim(split=False),dtype=np.int,order='C')
+ X = np.diag(self.domain.calc_weight(X,power=-1).astype(max(min(X.dtype,self.domain.datatype),min(X.dtype,self.target.datatype))).flatten(order='C'),k=0)
+ elif(np.size(X)==np.size(self.val)):
+ self.sym = None
+ self.uni = None
+ X = np.array(X).reshape(self.val.shape)
+ if(np.all(np.isreal(X))):
+ X = np.real(X)
+ X = X.astype(max(min(X.dtype,self.domain.datatype),min(X.dtype,self.target.datatype)))
+ else:
+ raise ValueError(about._errors.cstring("ERROR: dimension mismatch ( "+str(np.size(X))+" <> "+str(self.nrow())+" x "+str(self.ncol())+" )."))
+
+ ## add
+ if(X.dtype>self.val.dtype):
+ about.warnings.cprint("WARNING: datatype reset.")
+ self.val += X
+
+ ## check flags
+ self.sym,self.uni = self._check_flags(sym=self.sym,uni=self.uni)
+
+ self._inv = None ## reset
+
+ return self
+
+ def __sub__(self,X): ## __sub__ : self - X
+ if(isinstance(X,operator)):
+ if(self.domain!=X.domain):
+ raise ValueError(about._errors.cstring("ERROR: inequal domains."))
+ sym = (self.sym and X.sym)
+ uni = None
+ if(isinstance(X,explicit_operator)):
+ if(self.target!=X.target):
+ raise ValueError(about._errors.cstring("ERROR: inequal codomains."))
+ matrix = self.val-X.val
+ elif(isinstance(X,diagonal_operator)):
+ if(self.target.dim(split=False)!=X.target.dim(split=False))or(not self.target.check_codomain(X.target)):
+ raise ValueError(about._errors.cstring("ERROR: incompatible codomains."))
+ matrix = self.val-np.diag(X.diag(bare=True,domain=None),k=0) ## domain == X.domain
+ elif(isinstance(X,vecvec_operator)):
+ if(self.target!=X.target):
+ raise ValueError(about._errors.cstring("ERROR: inequal codomains."))
+ matrix = self.val-np.tensordot(X.val,X.val,axes=0)
+ else:
+ raise TypeError(about._errors.cstring("ERROR: unsupported or incompatible operator."))
+ elif(np.size(X)==1):
+ if(self.nrow()!=self.ncol()):
+ raise ValueError(about._errors.cstring("ERROR: identity ill-defined for "+str(self.nrow())+" x "+str(self.ncol())+" matrices."))
+ sym = self.sym
+ uni = None
+ X = X*np.ones(self.domain.dim(split=False),dtype=np.int,order='C')
+ X = self.domain.calc_weight(X,power=-1).astype(max(min(X.dtype,self.domain.datatype),min(X.dtype,self.target.datatype)))
+ matrix = self.val-np.diag(X.flatten(order='C'),k=0)
+ elif(np.size(X)==np.size(self.val)):
+ sym = None
+ uni = None
+ X = np.array(X).reshape(self.val.shape)
+ if(np.all(np.isreal(X))):
+ X = np.real(X)
+ X = X.astype(max(min(X.dtype,self.domain.datatype),min(X.dtype,self.target.datatype)))
+ matrix = self.val-X
+ else:
+ raise ValueError(about._errors.cstring("ERROR: dimension mismatch ( "+str(np.size(X))+" <> "+str(self.nrow())+" x "+str(self.ncol())+" )."))
+ return explicit_operator(self.domain,matrix,bare=True,sym=sym,uni=uni,target=self.target)
+
+ def __rsub__(self,X): ## __rsub__ : X - self
+ if(isinstance(X,operator)):
+ if(self.domain!=X.domain):
+ raise ValueError(about._errors.cstring("ERROR: inequal domains."))
+ sym = (self.sym and X.sym)
+ uni = None
+ if(isinstance(X,explicit_operator)):
+ if(self.target!=X.target):
+ raise ValueError(about._errors.cstring("ERROR: inequal codomains."))
+ matrix = X.val-self.val
+ elif(isinstance(X,diagonal_operator)):
+ if(self.target.dim(split=False)!=X.target.dim(split=False))or(not self.target.check_codomain(X.target)):
+ raise ValueError(about._errors.cstring("ERROR: incompatible codomains."))
+ matrix = np.diag(X.diag(bare=True,domain=None),k=0)-self.val ## domain == X.domain
+ elif(isinstance(X,vecvec_operator)):
+ if(self.target!=X.target):
+ raise ValueError(about._errors.cstring("ERROR: inequal codomains."))
+ matrix = np.tensordot(X.val,X.val,axes=0)-self.val
+ else:
+ raise TypeError(about._errors.cstring("ERROR: unsupported or incompatible operator."))
+ elif(np.size(X)==1):
+ if(self.nrow()!=self.ncol()):
+ raise ValueError(about._errors.cstring("ERROR: identity ill-defined for "+str(self.nrow())+" x "+str(self.ncol())+" matrices."))
+ sym = self.sym
+ uni = None
+ X = X*np.ones(self.domain.dim(split=False),dtype=np.int,order='C')
+ X = self.domain.calc_weight(X,power=-1).astype(max(min(X.dtype,self.domain.datatype),min(X.dtype,self.target.datatype)))
+ matrix = np.diag(X.flatten(order='C'),k=0)-self.val
+ elif(np.size(X)==np.size(self.val)):
+ sym = None
+ uni = None
+ X = np.array(X).reshape(self.val.shape)
+ if(np.all(np.isreal(X))):
+ X = np.real(X)
+ X = X.astype(max(min(X.dtype,self.domain.datatype),min(X.dtype,self.target.datatype)))
+ matrix = X-self.val
+ else:
+ raise ValueError(about._errors.cstring("ERROR: dimension mismatch ( "+str(np.size(X))+" <> "+str(self.nrow())+" x "+str(self.ncol())+" )."))
+ return explicit_operator(self.domain,matrix,bare=True,sym=sym,uni=uni,target=self.target)
+
+ def __isub__(self,X): ## __isub__ : self -= X
+ if(isinstance(X,operator)):
+ if(self.domain!=X.domain):
+ raise ValueError(about._errors.cstring("ERROR: inequal domains."))
+ self.sym = (self.sym and X.sym)
+ self.uni = None
+ if(isinstance(X,explicit_operator)):
+ if(self.target!=X.target):
+ raise ValueError(about._errors.cstring("ERROR: inequal codomains."))
+ X = X.val
+ elif(isinstance(X,diagonal_operator)):
+ if(self.target.dim(split=False)!=X.target.dim(split=False))or(not self.target.check_codomain(X.target)):
+ raise ValueError(about._errors.cstring("ERROR: incompatible codomains."))
+ X = np.diag(X.diag(bare=True,domain=None),k=0) ## domain == X.domain
+ elif(isinstance(X,vecvec_operator)):
+ if(self.target!=X.target):
+ raise ValueError(about._errors.cstring("ERROR: inequal codomains."))
+ X = np.tensordot(X.val,X.val,axes=0)
+ else:
+ raise TypeError(about._errors.cstring("ERROR: unsupported or incompatible operator."))
+ elif(np.size(X)==1):
+ if(self.nrow()!=self.ncol()):
+ raise ValueError(about._errors.cstring("ERROR: identity ill-defined for "+str(self.nrow())+" x "+str(self.ncol())+" matrices."))
+ self.uni = None
+ X = X*np.ones(self.domain.dim(split=False),dtype=np.int,order='C')
+ X = np.diag(self.domain.calc_weight(X,power=-1).astype(max(min(X.dtype,self.domain.datatype),min(X.dtype,self.target.datatype))).flatten(order='C'),k=0)
+ elif(np.size(X)==np.size(self.val)):
+ self.sym = None
+ self.uni = None
+ X = np.array(X).reshape(self.val.shape)
+ if(np.all(np.isreal(X))):
+ X = np.real(X)
+ X = X.astype(max(min(X.dtype,self.domain.datatype),min(X.dtype,self.target.datatype)))
+ else:
+ raise ValueError(about._errors.cstring("ERROR: dimension mismatch ( "+str(np.size(X))+" <> "+str(self.nrow())+" x "+str(self.ncol())+" )."))
+
+ ## subtract
+ if(X.dtype>self.val.dtype):
+ about.warnings.cprint("WARNING: datatype reset.")
+ self.val -= X
+
+ ## check flags
+ self.sym,self.uni = self._check_flags(sym=self.sym,uni=self.uni)
+
+ self._inv = None ## reset
+
+ return self
+
+ def _calc_mul(self,X,side): ## > multiplies self with X ...
+ if(side==0): ## ... from right
+ if(self._hidden[0]): ## hidden degrees of freedom
+ return np.array([np.apply_along_axis(self.domain.calc_dot,0,X,np.conjugate(rr)) for rr in self.weight(rowpower=1,colpower=0,overwrite=False)])
+ else:
+ return np.dot(self.weight(rowpower=1,colpower=0,overwrite=False),X,out=None)
+ elif(side==1): ## ... from left
+ if(self._hidden[1]): ## hidden degrees of freedom
+ return np.array([np.apply_along_axis(self.target.calc_dot,0,self.weight(rowpower=0,colpower=1,overwrite=False),np.conjugate(rr)) for rr in X])
+ else:
+ return np.dot(X,self.weight(rowpower=0,colpower=1,overwrite=False),out=None)
+ else:
+ raise ValueError(about._errors.cstring("ERROR: invalid input."))
+
+ def __mul__(self,X): ## __mul__ : self * X
+ if(isinstance(X,operator)):
+ if(self.domain!=X.target):
+ raise ValueError(about._errors.cstring("ERROR: incompatible spaces."))
+ newdomain = X.domain
+ sym = None
+ uni = (self.uni and X.uni)
+ if(isinstance(X,explicit_operator)):
+ X = X.val
+ elif(isinstance(X,diagonal_operator)):
+ X = np.diag(X.diag(bare=True,domain=None),k=0) ## domain == X.domain
+ elif(isinstance(X,vecvec_operator)):
+ X = np.tensordot(X.val,X.val,axes=0)
+ else:
+ raise TypeError(about._errors.cstring("ERROR: unsupported or incompatible operator."))
+ elif(np.size(X)==1):
+ newdomain = self.domain
+ sym = None
+ uni = None
+ X = X*np.ones(self.domain.dim(split=False),dtype=np.int,order='C')
+ X = np.diag(self.domain.calc_weight(X,power=-1).astype(max(min(X.dtype,self.domain.datatype),min(X.dtype,self.target.datatype))).flatten(order='C'),k=0)
+ elif(np.size(X)==self.val.shape[1]**2):
+ newdomain = self.domain
+ sym = None
+ uni = None
+ X = np.array(X).reshape((self.val.shape[1],self.val.shape[1]))
+ if(np.all(np.isreal(X))):
+ X = np.real(X)
+ X = X.astype(max(min(X.dtype,self.domain.datatype),min(X.dtype,self.target.datatype)))
+ else:
+ raise ValueError(about._errors.cstring("ERROR: dimension mismatch ( "+str(np.size(X))+" <> "+str(self.nrow())+" x "+str(self.nrow())+" )."))
+ return explicit_operator(newdomain,self._calc_mul(X,0),bare=True,sym=sym,uni=uni,target=self.target)
+
+ def __rmul__(self,X): ## __mul__ : X * self
+ if(isinstance(X,operator)):
+ if(X.domain!=self.target):
+ raise ValueError(about._errors.cstring("ERROR: incompatible spaces."))
+ newtarget = X.target
+ sym = None
+ uni = (self.uni and X.uni)
+ if(isinstance(X,explicit_operator)):
+ X = X.val
+ elif(isinstance(X,diagonal_operator)):
+ X = np.diag(X.diag(bare=True,domain=None),k=0) ## domain == X.domain
+ elif(isinstance(X,vecvec_operator)):
+ X = np.tensordot(X.val,X.val,axes=0)
+ else:
+ raise TypeError(about._errors.cstring("ERROR: unsupported or incompatible operator."))
+ elif(np.size(X)==1):
+ newtarget = self.target
+ sym = None
+ uni = None
+ X = X*np.ones(self.domain.dim(split=False),dtype=np.int,order='C')
+ X = np.diag(self.domain.calc_weight(X,power=-1).astype(max(min(X.dtype,self.domain.datatype),min(X.dtype,self.target.datatype))).flatten(order='C'),k=0)
+ elif(np.size(X)==self.val.shape[0]**2):
+ newtarget = self.target
+ sym = None
+ uni = None
+ X = np.array(X).reshape((self.val.shape[0],self.val.shape[0]))
+ if(np.all(np.isreal(X))):
+ X = np.real(X)
+ X = X.astype(max(min(X.dtype,self.domain.datatype),min(X.dtype,self.target.datatype)))
+ else:
+ raise ValueError(about._errors.cstring("ERROR: dimension mismatch ( "+str(np.size(X))+" <> "+str(self.ncol())+" x "+str(self.ncol())+" )."))
+ return explicit_operator(self.domain,self._calc_mul(X,1),bare=True,sym=sym,uni=uni,target=newtarget)
+
+ def __imul__(self,X): ## __imul__ : self *= X
+ if(isinstance(X,operator)):
+ if(self.domain!=X.target):
+ raise ValueError(about._errors.cstring("ERROR: incompatible spaces."))
+ if(self.domain!=X.domain):
+ raise ValueError(about._errors.cstring("ERROR: incompatible operator."))
+ self.sym = None
+ self.uni = (self.uni and X.uni)
+ if(isinstance(X,explicit_operator)):
+ X = X.val
+ elif(isinstance(X,diagonal_operator)):
+ X = np.diag(X.diag(bare=True,domain=None),k=0) ## domain == X.domain
+ elif(isinstance(X,vecvec_operator)):
+ X = np.tensordot(X.val,X.val,axes=0)
+ else:
+ raise TypeError(about._errors.cstring("ERROR: unsupported or incompatible operator."))
+ elif(np.size(X)==1):
+ self.sym = None
+ self.uni = None
+ X = X*np.ones(self.domain.dim(split=False),dtype=np.int,order='C')
+ X = np.diag(self.domain.calc_weight(X,power=-1).astype(max(min(X.dtype,self.domain.datatype),min(X.dtype,self.target.datatype))).flatten(order='C'),k=0)
+ elif(np.size(X)==self.val.shape[1]**2):
+ self.sym = None
+ self.uni = None
+ X = np.array(X).reshape((self.val.shape[1],self.val.shape[1]))
+ if(np.all(np.isreal(X))):
+ X = np.real(X)
+ X = X.astype(max(min(X.dtype,self.domain.datatype),min(X.dtype,self.target.datatype)))
+ else:
+ raise ValueError(about._errors.cstring("ERROR: dimension mismatch ( "+str(np.size(X))+" <> "+str(self.nrow())+" x "+str(self.nrow())+" )."))
+
+ ## multiply
+ if(X.dtype>self.val.dtype):
+ about.warnings.cprint("WARNING: datatype reset.")
+ self.val = self._calc_mul(X,0)
+
+ ## check flags
+ self.sym,self.uni = self._check_flags(sym=self.sym,uni=self.uni)
+
+ self._inv = None ## reset
+
+ return self
+
+ def __div__(self,X):
+ raise Exception(about._errors.cstring("ERROR: matrix division ill-defined."))
+
+ __rdiv__ = __div__
+ __idiv__ = __div__
+ __truediv__ = __div__
+ __rtruediv__ = __rdiv__
+ __itruediv__ = __idiv__
+
+ def __pow__(self,x): ## __pow__(): self ** x
+ if(not isinstance(x,(int,long))):
+ raise TypeError(about._errors.cstring("ERROR: non-integer exponent."))
+ elif(self.domain<>self.target):
+ raise ValueError(about._errors.cstring("ERROR: incompatible spaces."))
+ elif(x==0):
+ return identity(self.domain)
+ elif(x<0):
+ return self.inverse().__pow__(-x)
+ elif(x==1):
+ return self
+ else:
+ matrix = self._calc_mul(self.val,0)
+ for ii in xrange(x-1):
+ matrix = self._calc_mul(matrix,0)
+ return explicit_operator(self.domain,matrix,bare=True,sym=None,uni=self.uni,target=self.target)
+
+ def __rpow__(self,X): ## __pow__(): X ** self
+ raise Exception(about._errors.cstring("ERROR: matrix exponential ill-defined."))
+
+ def __ipow__(self,x): ## __pow__(): self **= x
+ if(not isinstance(x,(int,long))):
+ raise TypeError(about._errors.cstring("ERROR: non-integer exponent."))
+ elif(self.domain<>self.target):
+ raise ValueError(about._errors.cstring("ERROR: incompatible spaces."))
+ elif(x==0):
+ self.val = np.diag(self.domain.calc_weight(np.ones(self.domain.dim(split=False),dtype=np.int,order='C'),power=-1).astype(self.val.dtype),k=0)
+ elif(x<0):
+ self.val = (self.inverse().__pow__(-x)).val
+ elif(x==1):
+ pass
+ else:
+ matrix = self._calc_mul(self.val,0)
+ for ii in xrange(x-1):
+ matrix = self._calc_mul(matrix,0)
+ self.val = matrix
+
+ ## check flags
+ self.sym,self.uni = self._check_flags(sym=None,uni=self.uni)
+
+ self._inv = None ## reset
+
+ return self
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ @staticmethod
+ def get_plot(X,title="",vmin=None,vmax=None,unit="",norm=None,cmap=None,cbar=True,**kwargs):
+ """
+ Creates a plot of the matrix according to the given specifications.
+
+ Parameters
+ ----------
+ X : numpy.ndarray
+ Array containing the matrix.
+
+ Returns
+ -------
+ None
+
+ Other parameters
+ ----------------
+ title : string, *optional*
+ Title of the plot (default: "").
+ vmin : float, *optional*
+ Minimum value to be displayed (default: ``min(x)``).
+ vmax : float, *optional*
+ Maximum value to be displayed (default: ``max(x)``).
+ unit : string, *optional*
+ Unit of the field values (default: "").
+ norm : string, *optional*
+ Scaling of the field values before plotting (default: None).
+ cmap : matplotlib.colors.LinearSegmentedColormap, *optional*
+ Color map to be used for two-dimensional plots (default: None).
+ cbar : bool, *optional*
+ Whether to show the color bar or not (default: True).
+ save : string, *optional*
+ Valid file name where the figure is to be stored, by default
+ the figure is not saved (default: False).
+
+ Notes
+ -----
+ This method is static and thus also plot a valid numpy.ndarray.
+
+ See Also
+ --------
+ explicit_operator.plot
+
+ Raises
+ ------
+ ValueError
+ If the matrix `X` is not two-dimensional;
+ or if the logarithmic normalisation encounters negative values.
+
+ """
+ if(not pl.isinteractive())and(not bool(kwargs.get("save",False))):
+ about.warnings.cprint("WARNING: interactive mode off.")
+
+ if(len(X.shape)!=2):
+ raise ValueError(about._errors.cstring("ERROR: invalid matirx."))
+
+ if(vmin is None):
+ vmin = np.min(X,axis=None,out=None)
+ if(vmax is None):
+ vmax = np.max(X,axis=None,out=None)
+ if(norm=="log")and(vmin<=0):
+ raise ValueError(about._errors.cstring("ERROR: nonpositive value(s)."))
+
+ s_ = np.array([X.shape[1]/max(X.shape),X.shape[0]/max(X.shape)*(1.0+0.159*bool(cbar))])
+ fig = pl.figure(num=None,figsize=(6.4*s_[0],6.4*s_[1]),dpi=None,facecolor=None,edgecolor=None,frameon=False,FigureClass=pl.Figure)
+ ax0 = fig.add_axes([0.06/s_[0],0.06/s_[1],1.0-0.12/s_[0],1.0-0.12/s_[1]])
+ if(norm=="log"):
+ n_ = ln(vmin=vmin,vmax=vmax)
+ else:
+ n_ = None
+ sub = ax0.pcolormesh(X[::-1,:],cmap=cmap,norm=n_,vmin=vmin,vmax=vmax)
+ ax0.set_xlim(0,X.shape[1])
+ ax0.set_xticks([],minor=False)
+ ax0.set_ylim(0,X.shape[0])
+ ax0.set_yticks([],minor=False)
+ ax0.set_aspect("equal")
+ if(cbar):
+ if(norm=="log"):
+ f_ = lf(10,labelOnlyBase=False)
+ b_ = sub.norm.inverse(np.linspace(0,1,sub.cmap.N+1))
+ v_ = np.linspace(sub.norm.vmin,sub.norm.vmax,sub.cmap.N)
+ else:
+ f_ = None
+ b_ = None
+ v_ = None
+ cb0 = fig.colorbar(sub,ax=ax0,orientation="horizontal",fraction=0.1,pad=0.05,shrink=0.75,aspect=20,ticks=[vmin,vmax],format=f_,drawedges=False,boundaries=b_,values=v_)
+ cb0.ax.text(0.5,-1.0,unit,fontdict=None,withdash=False,transform=cb0.ax.transAxes,horizontalalignment="center",verticalalignment="center")
+ ax0.set_title(title)
+
+ if(bool(kwargs.get("save",False))):
+ fig.savefig(str(kwargs.get("save")),dpi=None,facecolor=None,edgecolor=None,orientation='portrait',papertype=None,format=None,transparent=False,bbox_inches=None,pad_inches=0.1)
+ pl.close(fig)
+ else:
+ fig.canvas.draw()
+
+ def plot(self,bare=True,**kwargs):
+ """
+ Creates a plot of the matrix according to the given specifications.
+
+ Parameters
+ ----------
+ bare : {bool, 2-tuple}, *optional*
+ Whether the returned matrix entries are `bare` or not
+ (default: True).
+
+ Returns
+ -------
+ None
+
+ Other parameters
+ ----------------
+ title : string, *optional*
+ Title of the plot (default: "").
+ vmin : float, *optional*
+ Minimum value to be displayed (default: ``min(x)``).
+ vmax : float, *optional*
+ Maximum value to be displayed (default: ``max(x)``).
+ unit : string, *optional*
+ Unit of the field values (default: "").
+ norm : string, *optional*
+ Scaling of the field values before plotting (default: None).
+ cmap : matplotlib.colors.LinearSegmentedColormap, *optional*
+ Color map to be used for two-dimensional plots (default: None).
+ cbar : bool, *optional*
+ Whether to show the color bar or not (default: True).
+ save : string, *optional*
+ Valid file name where the figure is to be stored, by default
+ the figure is not saved (default: False).
+
+ Notes
+ -----
+ The ambiguity of `bare` or non-bare diagonal entries is based
+ on the choice of a matrix representation of the operator in
+ question. The naive choice of absorbing the volume weights
+ into the matrix leads to a matrix-vector calculus with the
+ non-bare entries which seems intuitive, though. The choice of
+ keeping matrix entries and volume weights separate deals with the
+ bare entries that allow for correct interpretation of the matrix
+ entries; e.g., as variance in case of an covariance operator.
+
+ See Also
+ --------
+ explicit_operator.get_plot
+
+ """
+ interactive = pl.isinteractive()
+ pl.matplotlib.interactive(not bool(kwargs.get("save",False)))
+
+ X = self.get_matrix(bare=bare)
+
+ if(np.any(np.iscomplex(X))):
+ about.infos.cprint("INFO: absolute values and phases are plotted.")
+ if(kwargs.has_key("title")):
+ title = kwargs.get("title")+" "
+ kwargs.__delitem__("title")
+ else:
+ title = ""
+ self.get_plot(np.absolute(X),title=title+"(absolute)",**kwargs)
+ if(kwargs.has_key("vmin")):
+ kwargs.__delitem__("vmin")
+ if(kwargs.has_key("vmin")):
+ kwargs.__delitem__("vmax")
+ if(kwargs.has_key("unit")):
+ kwargs["unit"] = "rad"
+ if(kwargs.has_key("norm")):
+ kwargs["norm"] = None
+ if(not kwargs.has_key("cmap")):
+ kwargs["cmap"] = pl.cm.hsv_r
+ self.get_plot(np.angle(X,deg=False),title=title+"(phase)",vmin=-3.1416,vmax=3.1416,**kwargs) ## values in [-pi,pi]
+ else:
+ self.get_plot(np.real(X),**kwargs)
+
+ pl.matplotlib.interactive(interactive)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def __repr__(self):
+ return "<nifty_explicit.explicit_operator>"
+
+##-----------------------------------------------------------------------------
+
+
+
+
+
+##-----------------------------------------------------------------------------
+
+class _share(object):
+
+ __init__ = None
+
+ @staticmethod
+ def _init_share(_mat,_num):
+ _share.mat = _mat
+ _share.num = _num
+
+##-----------------------------------------------------------------------------
+
+##-----------------------------------------------------------------------------
+
+class explicit_probing(probing):
+ """
+ ..
+ ..
+ .. __ __ __ __
+ .. / / /__/ /__/ / /_
+ .. _______ __ __ ______ / / __ _______ __ / _/
+ .. / __ / \ \/ / / _ | / / / / / ____/ / / / /
+ .. / /____/ / / / /_/ / / /_ / / / /____ / / / /_
+ .. \______/ /__/\__\ / ____/ \___/ /__/ \______/ /__/ \___/ probing class
+ .. /__/
+
+
+ NIFTY subclass for explicit probing (using multiprocessing)
+
+ Called after initialization, an explicit matrix representation of a
+ linear operator or function is sampled by applying (weighted) canonical
+ vectors in a specified basis.
+
+ Parameters
+ ----------
+ op : operator, *optional*
+ Operator to be probed; if not given, the probing will resort to
+ `function` (default: None).
+ function : function, *optional*
+ Function applied to the probes; either `op` or `function` must
+ be given (default: `op.times`).
+ domain : space, *optional*
+ Space wherein the probes live, defines the domain of the
+ explicified operator (default: `op.domain`).
+ codomain : space, *optional*
+ Space wherein the output of the explicified operator lives
+ (default: `op.target`).
+ target : space, *optional*
+ Space wherein the transform of the probes live (default: None).
+ ncpu : integer, *optional*
+ Number of used CPUs to use (default: 2).
+ nper : integer, *optional*
+ Number of tasks performed by one worker (default: 1).
+
+ See Also
+ --------
+ probing
+
+ Examples
+ --------
+ >>> v = field(point_space(3), val=[1, 2, 3])
+ >>> W = vecvec_operator(val=v) # implicit operator
+ >>> W_ij = explicit_probing(Wim)(loop=True) # explicit operator
+ >>> W_ij.get_matrix()
+ array([[ 1., 2., 3.],
+ [ 2., 4., 6.],
+ [ 3., 6., 9.]])
+
+ Attributes
+ ----------
+ function : function
+ Function applied to the probes.
+ domain : space, *optional*
+ Space wherein the probes live, defines the domain of the
+ explicified operator.
+ codomain : space, *optional*
+ Space wherein the output of the explicified operator lives.
+ target : space
+ Space wherein the transform of the probes live.
+ ncpu : integer
+ Number of used CPUs to use.
+ nper : integer
+ Number of tasks performed by one worker.
+ quargs : dict
+ Keyword arguments passed to `function` in each call.
+
+ """
+ def __init__(self,op=None,function=None,domain=None,codomain=None,target=None,ncpu=2,nper=1,**quargs):
+ """
+ Initializes the explicit probing and sets the standard probing
+ attributes except for `random`, `nrun`, and `var`.
+
+ Parameters
+ ----------
+ op : operator, *optional*
+ Operator to be probed; if not given, the probing will resort to
+ `function` (default: None).
+ function : function, *optional*
+ Function applied to the probes; either `op` or `function` must
+ be given (default: `op.times`).
+ domain : space, *optional*
+ Space wherein the probes live, defines the domain of the
+ explicified operator (default: `op.domain`).
+ codomain : space, *optional*
+ Space wherein the output of the explicified operator lives
+ (default: `op.target`).
+ target : space, *optional*
+ Space wherein the transform of the probes live (default: None).
+ ncpu : integer, *optional*
+ Number of used CPUs to use (default: 2).
+ nper : integer, *optional*
+ Number of tasks performed by one worker (default: 1).
+
+ Other Parameters
+ ----------------
+ quargs : dict
+ Keyword arguments passed to `function` in each call.
+
+ Raises
+ ------
+ TypeError
+ If input is invalid or insufficient.
+ NameError
+ If `function` is not an attribute of `op`.
+ ValueError
+ If spaces are incompatible.
+
+ """
+ if(op is None):
+ ## check whether callable
+ if(function is None)or(not hasattr(function,"__call__")):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ ## check given domain
+ if(domain is None):
+ raise TypeError(about._errors.cstring("ERROR: insufficient input."))
+ elif(not isinstance(domain,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input domain."))
+ ## check given codomain
+ if(codomain is None):
+ raise TypeError(about._errors.cstring("ERROR: insufficient input."))
+ elif(not isinstance(codomain,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input codomain."))
+ else:
+ if(not isinstance(op,operator)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ ## check whether callable
+ if(function is None)or(not hasattr(function,"__call__")):
+ function = op.times
+ elif(op==function):
+ function = op.times
+ ## check whether correctly bound
+ if(op!=function.im_self):
+ raise NameError(about._errors.cstring("ERROR: invalid input function."))
+ ## check given domain
+ if(domain is None)or(not isinstance(domain,space)):
+ if(function in [op.inverse_times,op.adjoint_times]):
+ domain = op.target
+ else:
+ domain = op.domain
+ else:
+ if(function in [op.inverse_times,op.adjoint_times]):
+ op.target.check_codomain(domain) ## a bit pointless
+ else:
+ op.domain.check_codomain(domain) ## a bit pointless
+ ## check given codomain
+ if(codomain is None)or(not isinstance(codomain,space)):
+ if(function in [op.inverse_times,op.adjoint_times]):
+ codomain = op.domain
+ else:
+ codomain = op.target
+ else:
+ if(function in [op.inverse_times,op.adjoint_times]):
+ if(not op.domain.check_codomain(domain))and(op.domain.dim(split=False)!=codomain.dim(split=False)):
+ raise ValueError(about._errors.cstring("ERROR: incompatible input codomain."))
+ else:
+ if(not op.target.check_codomain(domain))and(op.target.dim(split=False)!=codomain.dim(split=False)):
+ raise ValueError(about._errors.cstring("ERROR: incompatible input codomain."))
+
+ self.function = function
+ self.domain = domain
+ self.codomain = codomain
+
+ ## check target
+ if(target is None):
+ target = self.domain.get_codomain()
+ else:
+ self.domain.check_codomain(target) ## a bit pointless
+ self.target = target
+
+ ## check shape
+ if(self.domain.dim(split=False)*self.codomain.dim(split=False)>1048576):
+ about.warnings.cprint("WARNING: matrix size > 2 ** 20.")
+
+ self.ncpu = int(max(1,ncpu))
+ self.nrun = self.domain.dim(split=False)
+ if(nper is None):
+ self.nper = None
+ else:
+ self.nper = int(max(1,min(self.nrun//self.ncpu,nper)))
+
+ self.quargs = quargs
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def configure(self,**kwargs):
+ """
+ Changes the performance attributes of the instance.
+
+ Parameters
+ ----------
+ ncpu : integer, *optional*
+ Number of used CPUs to use.
+ nper : integer, *optional*
+ Number of tasks performed by one worker.
+
+ """
+ if("ncpu" in kwargs):
+ self.ncpu = int(max(1,kwargs.get("ncpu")))
+ if("nper" in kwargs):
+ if(kwargs.get("nper") is None):
+ self.nper = None
+ else:
+ self.nper = int(max(1,min(self.nrun//self.ncpu,kwargs.get("nper"))))
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def gen_probe(self,index,value):
+ """
+ Generates a probe.
+
+ For explicit probing, each probe is a (weighted) canonical base.
+
+ Parameters
+ ----------
+ index : int
+ Index where to put the value.
+ value : scalar
+ Weighted 1.
+
+ Returns
+ -------
+ probe : field
+ Weighted canonical base.
+
+ """
+ probe = field(self.domain,val=None,target=self.target)
+ probe[np.unravel_index(index,self.domain.dim(split=True),order='C')] = value
+ return probe
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def probing(self,idnum,probe):
+ """
+ Computes a single probing result meaning a single column of the
+ linear operators matrix representation.
+
+ Parameters
+ ----------
+ probe : field
+ Weighted canonical base on which `function` will be applied.
+ idnum : int
+ Identification number; obsolete.
+
+ Returns
+ -------
+ result : ndarray
+ Result of function evaluation (equaling a column of the matrix).
+
+ """
+ f = self.function(probe,**self.quargs) ## field
+ if(f is None):
+ return None
+ elif(isinstance(f,field)):
+ if(f.domain!=self.codomain):
+ try:
+ f.transform(target=self.codomain,overwrite=True)
+ except(ValueError): ## unsupported transformation
+ pass ## checkless version of f.cast_domain(self.codomain,newtarget=None,force=True)
+ return f.val.flatten(order='C')
+ else:
+ return f.flatten(order='C')
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def evaluate(self,matrix,num):
+ """
+ Evaluates the probing results.
+
+ Parameters
+ ----------
+ matrix : numpy.array
+ Matrix representation of the probed linear operator.
+ num : int
+ Number of successful probings (not returning ``None``).
+
+ Returns
+ -------
+ result : explicit_operator
+ The probed linear operator as explicit operator instance.
+
+ """
+ if(num<self.nrun):
+ about.infos.cflush(" ( %u probe(s) failed, effectiveness == %.1f%% )\n"%(self.nrun-num,100*num/self.nrun))
+ if(num==0):
+ about.warnings.cprint("WARNING: probing failed.")
+ return None
+ else:
+ about.infos.cflush("\n")
+
+ return explicit_operator(self.domain,matrix,bare=True,sym=None,uni=None,target=self.codomain)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def _single_probing(self,zipped): ## > performs one probing operation
+ ## generate probe
+ probe = self.gen_probe(*zipped)
+ ## do the actual probing
+ return self.probing(zipped[0],probe)
+
+ def _serial_probing(self,zipped): ## > performs the probing operation serially
+ try:
+ result = self._single_probing(zipped)
+ except Exception as exception:
+ raise exception
+ except BaseException: ## capture system-exiting exception including KeyboardInterrupt
+ raise Exception(about._errors.cstring("ERROR: unknown."))
+ else:
+ if(result is not None):
+ result = np.array(result).flatten(order='C')
+ rindex = zipped[0]*self.codomain.dim(split=False)
+ if(isinstance(_share.mat,tuple)):
+ _share.mat[0].acquire(block=True,timeout=None)
+ _share.mat[0][rindex:rindex+self.codomain.dim(split=False)] = np.real(result)
+ _share.mat[0].release()
+ _share.mat[1].acquire(block=True,timeout=None)
+ _share.mat[1][rindex:rindex+self.codomain.dim(split=False)] = np.imag(result)
+ _share.mat[1].release()
+ else:
+ _share.mat.acquire(block=True,timeout=None)
+ _share.mat[rindex:rindex+self.codomain.dim(split=False)] = result
+ _share.mat.release()
+ _share.num.acquire(block=True,timeout=None)
+ _share.num.value += 1
+ _share.num.release()
+ self._progress(_share.num.value)
+
+ def _parallel_probing(self): ## > performs the probing operations in parallel
+ ## define weighted canonical base
+ base = self.domain.calc_weight(self.domain.enforce_values(1,extend=True),power=-1).flatten(order='C')
+ ## define shared objects
+ if(issubclass(self.codomain.datatype,np.complexfloating)):
+ _mat = (ma('d',np.empty(self.nrun*self.codomain.dim(split=False),dtype=np.float64,order='C'),lock=True),ma('d',np.empty(self.nrun*self.codomain.dim(split=False),dtype=np.float64,order='C'),lock=True)) ## tuple(real,imag)
+ else:
+ _mat = ma('d',np.empty(self.nrun*self.codomain.dim(split=False),dtype=np.float64,order='C'),lock=True)
+ _num = mv('I',0,lock=True)
+ ## build pool
+ if(about.infos.status):
+ so.write(about.infos.cstring("INFO: multiprocessing "+(' ')*10))
+ so.flush()
+ pool = mp(processes=self.ncpu,initializer=_share._init_share,initargs=(_mat,_num),maxtasksperchild=self.nper)
+ try:
+ ## pooling
+ pool.map_async(self._serial_probing,zip(np.arange(self.nrun,dtype=np.int),base),chunksize=None,callback=None).get(timeout=None)
+ ## close and join pool
+ about.infos.cflush(" done.")
+ pool.close()
+ pool.join()
+ except BaseException as exception:
+ ## terminate and join pool
+ about._errors.cprint("\nERROR: terminating pool.")
+ pool.terminate()
+ pool.join()
+ ## re-raise exception
+ raise exception ## traceback by looping
+ ## evaluate
+ if(issubclass(self.codomain.datatype,np.complexfloating)):
+ _mat = (np.array(_mat[0][:])+np.array(_mat[1][:])*1j) ## comlpex array
+ else:
+ _mat = np.array(_mat[:])
+ _mat = _mat.reshape((self.nrun,self.codomain.dim(split=False))).T
+ return self.evaluate(_mat,_num.value)
+
+ def _nonparallel_probing(self): ## > performs the probing operations one after another
+ ## define weighted canonical base
+ base = self.domain.calc_weight(self.domain.enforce_values(1,extend=True),power=-1).flatten(order='C')
+ ## looping
+ if(about.infos.status):
+ so.write(about.infos.cstring("INFO: looping "+(' ')*10))
+ so.flush()
+ _mat = np.empty((self.nrun,self.codomain.dim(split=False)),dtype=self.codomain.datatype,order='C')
+ _num = 0
+ for ii in xrange(self.nrun):
+ result = self._single_probing((ii,base[ii])) ## zipped tuple
+ if(result is None):
+ _mat[ii] = np.zeros(self.codomain.dim(split=False),dtype=self.codomain.datatype)
+ else:
+ _mat[ii] = np.array(result,dtype=self.codomain.datatype)
+ _num += 1
+ self._progress(_num)
+ about.infos.cflush(" done.")
+ ## evaluate
+ return self.evaluate(_mat.T,_num)
+
+ def __call__(self,loop=False,**kwargs):
+ """
+ Start the explicit probing procedure.
+
+ Parameters
+ ----------
+ loop : bool, *optional*
+ Whether to perform a loop or to parallelise (default: False).
+
+ Returns
+ -------
+ result : explicit_operator
+ The probed linear operator as explicit operator instance.
+
+ Other Parameters
+ ----------------
+ ncpu : integer, *optional*
+ Number of used CPUs to use.
+ nper : integer, *optional*
+ Number of tasks performed by one worker.
+
+
+
+ """
+ self.configure(**kwargs)
+ if(not about.multiprocessing.status)or(loop):
+ return self._nonparallel_probing()
+ else:
+ return self._parallel_probing()
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def __repr__(self):
+ return "<nifty_explicit.explicit_probing>"
+
+##-----------------------------------------------------------------------------
+
+##-----------------------------------------------------------------------------
+
+def explicify(op,newdomain=None,newtarget=None,ncpu=2,nper=1,loop=False,**quargs):
+ """
+ Explicifys an (implicit) operator.
+
+ This function wraps the :py:class:`explicit_probing` class with a more
+ user-friendly interface.
+
+ Parameters
+ ----------
+ op : operator
+ Operator to be explicified.
+ newdomain : space, *optional*
+ Space wherein the probes live, defines the domain of the
+ explicified operator (default: `op.domain`).
+ newtarget : space, *optional*
+ Space wherein the output of the explicified operator lives
+ (default: `op.target`).
+ ncpu : integer, *optional*
+ Number of used CPUs to use (default: 2).
+ nper : integer, *optional*
+ Number of tasks performed by one worker (default: 1).
+
+ Returns
+ -------
+ EO : explicit_operator
+ The explicified linear operator as explicit operator instance.
+
+ Other Parameters
+ ----------------
+ quargs : dict
+ Keyword arguments passed to `function` in each call.
+
+ See Also
+ --------
+ explicit_probing
+
+ Examples
+ --------
+ >>> x_space = rg_space(4)
+ >>> k_space = x_space.get_codomain()
+ >>> S = power_operator(k_space, spec=[9, 3, 1])
+ >>> S.diag() # implicit operator in k_space
+ array([ 1., 3., 9., 3.])
+ >>> S_kq = explicify(S)
+ >>> S_kq.get_matrix() # explicit operator in k_space
+ array([[ 1., 0., 0., 0.],
+ [ 0., 3., 0., 0.],
+ [ 0., 0., 9., 0.],
+ [ 0., 0., 0., 3.]])
+ >>> S_xy = explicify(S, newdomain=x_space)
+ >>> S_xy.get_matrix() # explicit operator in x_space
+ array([[ 16., 8., 4., 8.],
+ [ 8., 16., 8., 4.],
+ [ 4., 8., 16., 8.],
+ [ 8., 4., 8., 16.]])
+
+ Raises
+ ------
+ TypeError
+ If `op` is no operator instance.
+
+ """
+ if(not isinstance(op,operator)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ elif(newdomain is not None)and(newtarget is None)and(op.domain==op.target):
+ newtarget = newdomain
+ if(newdomain is None)or(newdomain==op.domain):
+ target_ = None
+ else:
+ target_ = op.domain
+ return explicit_probing(op=op,function=op.times,domain=newdomain,codomain=newtarget,target=target_,ncpu=ncpu,nper=nper,**quargs)(loop=loop)
+
+##-----------------------------------------------------------------------------
+
diff --git a/operators/nifty_operators.py b/operators/nifty_operators.py
new file mode 100644
index 0000000000000000000000000000000000000000..894466601259186e243a4278725a8f6cee4d64cf
--- /dev/null
+++ b/operators/nifty_operators.py
@@ -0,0 +1,3477 @@
+## NIFTY (Numerical Information Field Theory) has been developed at the
+## Max-Planck-Institute for Astrophysics.
+##
+## Copyright (C) 2015 Max-Planck-Society
+##
+## Author: Marco Selig
+## Project homepage: <http://www.mpa-garching.mpg.de/ift/nifty/>
+##
+## This program is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+## See the GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+from __future__ import division
+import numpy as np
+from nifty.nifty_about import about
+from nifty.nifty_core import space, \
+ point_space, \
+ nested_space, \
+ field
+from nifty.nifty_tools import conjugate_gradient
+from nifty_probing import trace_probing, \
+ diagonal_probing
+
+##=============================================================================
+
+class operator(object):
+ """
+ .. __
+ .. / /_
+ .. ______ ______ _______ _____ ____ __ / _/ ______ _____
+ .. / _ | / _ | / __ / / __/ / _ / / / / _ | / __/
+ .. / /_/ / / /_/ / / /____/ / / / /_/ / / /_ / /_/ / / /
+ .. \______/ / ____/ \______/ /__/ \______| \___/ \______/ /__/ class
+ .. /__/
+
+ NIFTY base class for (linear) operators
+
+ The base NIFTY operator class is an abstract class from which other
+ specific operator subclasses, including those preimplemented in NIFTY
+ (e.g. the diagonal operator class) must be derived.
+
+ Parameters
+ ----------
+ domain : space
+ The space wherein valid arguments live.
+ sym : bool, *optional*
+ Indicates whether the operator is self-adjoint or not
+ (default: False)
+ uni : bool, *optional*
+ Indicates whether the operator is unitary or not
+ (default: False)
+ imp : bool, *optional*
+ Indicates whether the incorporation of volume weights in
+ multiplications is already implemented in the `multiply`
+ instance methods or not (default: False)
+ target : space, *optional*
+ The space wherein the operator output lives (default: domain)
+ para : {single object, list of objects}, *optional*
+ This is a freeform list of parameters that derivatives of the
+ operator class can use. Not used in the base operators.
+ (default: None)
+
+ See Also
+ --------
+ diagonal_operator : An operator class for handling purely diagonal
+ operators.
+ power_operator : Similar to diagonal_operator but with handy features
+ for dealing with diagonal operators whose diagonal
+ consists of a power spectrum.
+ vecvec_operator : Operators constructed from the outer product of two
+ fields
+ response_operator : Implements a modeled instrument response which
+ translates a signal into data space.
+ projection_operator : An operator that projects out one or more
+ components in a basis, e.g. a spectral band
+ of Fourier components.
+
+ Attributes
+ ----------
+ domain : space
+ The space wherein valid arguments live.
+ sym : bool
+ Indicates whether the operator is self-adjoint or not
+ uni : bool
+ Indicates whether the operator is unitary or not
+ imp : bool
+ Indicates whether the incorporation of volume weights in
+ multiplications is already implemented in the `multiply`
+ instance methods or not
+ target : space
+ The space wherein the operator output lives
+ para : {single object, list of objects}
+ This is a freeform list of parameters that derivatives of the
+ operator class can use. Not used in the base operators.
+ """
+ def __init__(self,domain,sym=False,uni=False,imp=False,target=None,para=None):
+ """
+ Sets the attributes for an operator class instance.
+
+ Parameters
+ ----------
+ domain : space
+ The space wherein valid arguments live.
+ sym : bool, *optional*
+ Indicates whether the operator is self-adjoint or not
+ (default: False)
+ uni : bool, *optional*
+ Indicates whether the operator is unitary or not
+ (default: False)
+ imp : bool, *optional*
+ Indicates whether the incorporation of volume weights in
+ multiplications is already implemented in the `multiply`
+ instance methods or not (default: False)
+ target : space, *optional*
+ The space wherein the operator output lives (default: domain)
+ para : {object, list of objects}, *optional*
+ This is a freeform list of parameters that derivatives of the
+ operator class can use. Not used in the base operators.
+ (default: None)
+
+ Returns
+ -------
+ None
+ """
+ if(not isinstance(domain,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ self.domain = domain
+ self.sym = bool(sym)
+ self.uni = bool(uni)
+
+ 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
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def nrow(self):
+ """
+ Computes the number of rows.
+
+ Returns
+ -------
+ nrow : int
+ number of rows (equal to the dimension of the codomain)
+ """
+ return self.target.dim(split=False)
+
+ def ncol(self):
+ """
+ Computes the number of columns
+
+ Returns
+ -------
+ ncol : int
+ number of columns (equal to the dimension of the domain)
+ """
+ return self.domain.dim(split=False)
+
+ def dim(self,axis=None):
+ """
+ Computes the dimension of the space
+
+ Parameters
+ ----------
+ axis : int, *optional*
+ Axis along which the dimension is to be calculated.
+ (default: None)
+
+ Returns
+ -------
+ dim : {int, ndarray}
+ The dimension(s) of the operator.
+
+ """
+ if(axis is None):
+ return np.array([self.nrow(),self.ncol()])
+ elif(axis==0):
+ return self.nrow()
+ elif(axis==1):
+ return self.ncol()
+ else:
+ raise ValueError(about._errors.cstring("ERROR: invalid input axis."))
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def set_para(self,newpara):
+ """
+ Sets the parameters and creates the `para` property if it does
+ not exist
+
+ Parameters
+ ----------
+ newpara : {object, list of objects}
+ A single parameter or a list of parameters.
+
+ Returns
+ -------
+ None
+
+ """
+ self.para = newpara
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def _multiply(self,x,**kwargs): ## > applies the operator to a given field
+ raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'multiply'."))
+
+ def _adjoint_multiply(self,x,**kwargs): ## > applies the adjoint operator to a given field
+ raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'adjoint_multiply'."))
+
+ def _inverse_multiply(self,x,**kwargs): ## > applies the inverse operator to a given field
+ raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'inverse_multiply'."))
+
+ def _adjoint_inverse_multiply(self,x,**kwargs): ## > applies the inverse adjoint operator to a given field
+ raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'adjoint_inverse_multiply'."))
+
+ def _inverse_adjoint_multiply(self,x,**kwargs): ## > applies the adjoint inverse operator to a given field
+ raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'inverse_adjoint_multiply'."))
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def _briefing(self,x,domain,inverse): ## > prepares x for `multiply`
+ ## inspect x
+ if(not isinstance(x,field)):
+ x_ = field(domain,val=x,target=None)
+ else:
+ ## check x.domain
+ if(domain==x.domain):
+ x_ = x
+ ## transform
+ else:
+ x_ = x.transform(target=domain,overwrite=False)
+ ## weight if ...
+ if(not self.imp)and(not domain.discrete)and(not inverse):
+ x_ = x_.weight(power=1,overwrite=False)
+ return x_
+
+ def _debriefing(self,x,x_,target,inverse): ## > evaluates x and x_ after `multiply`
+ if(x_ is None):
+ return None
+ else:
+ ## inspect x_
+ if(not isinstance(x_,field)):
+ x_ = field(target,val=x_,target=None)
+ elif(x_.domain!=target):
+ raise ValueError(about._errors.cstring("ERROR: invalid output domain."))
+ ## weight if ...
+ if(not self.imp)and(not target.discrete)and(inverse):
+ x_ = x_.weight(power=-1,overwrite=False)
+ ## inspect x
+ if(isinstance(x,field)):
+ ## repair ...
+ if(self.domain==self.target!=x.domain):
+ x_ = x_.transform(target=x.domain,overwrite=False) ## ... domain
+ if(x_.domain==x.domain)and(x_.target!=x.target):
+ x_.set_target(newtarget=x.target) ## ... codomain
+ return x_
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def times(self,x,**kwargs):
+ """
+ Applies the operator to a given object
+
+ Parameters
+ ----------
+ x : {scalar, ndarray, field}
+ Scalars are interpreted as constant arrays, and an array will
+ be interpreted as a field on the domain of the operator.
+
+ Returns
+ -------
+ Ox : field
+ Mapped field on the target domain of the operator.
+ """
+ ## prepare
+ x_ = self._briefing(x,self.domain,False)
+ ## apply operator
+ x_ = self._multiply(x_,**kwargs)
+ ## evaluate
+ return self._debriefing(x,x_,self.target,False)
+
+ def __call__(self,x,**kwargs):
+ return self.times(x,**kwargs)
+
+ def adjoint_times(self,x,**kwargs):
+ """
+ Applies the adjoint operator to a given object.
+
+ Parameters
+ ----------
+ x : {scalar, ndarray, field}
+ Scalars are interpreted as constant arrays, and an array will
+ be interpreted as a field on the target space of the operator.
+
+ Returns
+ -------
+ OAx : field
+ Mapped field on the domain of the operator.
+
+ """
+ ## check whether self-adjoint
+ if(self.sym):
+ return self.times(x,**kwargs)
+ ## check whether unitary
+ if(self.uni):
+ return self.inverse_times(x,**kwargs)
+
+ ## prepare
+ x_ = self._briefing(x,self.target,False)
+ ## apply operator
+ x_ = self._adjoint_multiply(x_,**kwargs)
+ ## evaluate
+ return self._debriefing(x,x_,self.domain,False)
+
+ def inverse_times(self,x,**kwargs):
+ """
+ Applies the inverse operator to a given object.
+
+ Parameters
+ ----------
+ x : {scalar, ndarray, field}
+ Scalars are interpreted as constant arrays, and an array will
+ be interpreted as a field on the domain space of the operator.
+
+ Returns
+ -------
+ OIx : field
+ Mapped field on the target space of the operator.
+ """
+ ## check whether self-inverse
+ if(self.sym)and(self.uni):
+ return self.times(x,**kwargs)
+
+ ## prepare
+ x_ = self._briefing(x,self.target,True)
+ ## apply operator
+ x_ = self._inverse_multiply(x_,**kwargs)
+ ## evaluate
+ return self._debriefing(x,x_,self.domain,True)
+
+ def adjoint_inverse_times(self,x,**kwargs):
+ """
+ Applies the inverse adjoint operator to a given object.
+
+ Parameters
+ ----------
+ x : {scalar, ndarray, field}
+ Scalars are interpreted as constant arrays, and an array will
+ be interpreted as a field on the target space of the operator.
+
+ Returns
+ -------
+ OAIx : field
+ Mapped field on the domain of the operator.
+
+ """
+ ## check whether self-adjoint
+ if(self.sym):
+ return self.inverse_times(x,**kwargs)
+ ## check whether unitary
+ if(self.uni):
+ return self.times(x,**kwargs)
+
+ ## prepare
+ x_ = self._briefing(x,self.domain,True)
+ ## apply operator
+ x_ = self._adjoint_inverse_multiply(x_,**kwargs)
+ ## evaluate
+ return self._debriefing(x,x_,self.target,True)
+
+ def inverse_adjoint_times(self,x,**kwargs):
+ """
+ Applies the adjoint inverse operator to a given object.
+
+ Parameters
+ ----------
+ x : {scalar, ndarray, field}
+ Scalars are interpreted as constant arrays, and an array will
+ be interpreted as a field on the target space of the operator.
+
+ Returns
+ -------
+ OIAx : field
+ Mapped field on the domain of the operator.
+
+ """
+ ## check whether self-adjoint
+ if(self.sym):
+ return self.inverse_times(x,**kwargs)
+ ## check whether unitary
+ if(self.uni):
+ return self.times(x,**kwargs)
+
+ ## prepare
+ x_ = self._briefing(x,self.domain,True)
+ ## apply operator
+ x_ = self._inverse_adjoint_multiply(x_,**kwargs)
+ ## evaluate
+ return self._debriefing(x,x_,self.target,True)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def tr(self,domain=None,target=None,random="pm1",ncpu=2,nrun=8,nper=1,var=False,loop=False,**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.
+
+ See Also
+ --------
+ probing : The class used to perform probing operations
+ """
+ if(domain is None):
+ domain = self.domain
+ return trace_probing(self,function=self.times,domain=domain,target=target,random=random,ncpu=(ncpu,1)[bool(loop)],nrun=nrun,nper=nper,var=var,**kwargs)(loop=loop)
+
+ def inverse_tr(self,domain=None,target=None,random="pm1",ncpu=2,nrun=8,nper=1,var=False,loop=False,**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)
+ 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 inverse operator
+ delta : float, *optional*
+ Probing variance of the trace. Returned if `var` is True in
+ of probing case.
+
+ See Also
+ --------
+ probing : The class used to perform probing operations
+ """
+ if(domain is None):
+ domain = self.target
+ return trace_probing(self,function=self.inverse_times,domain=domain,target=target,random=random,ncpu=(ncpu,1)[bool(loop)],nrun=nrun,nper=nper,var=var,**kwargs)(loop=loop)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def diag(self,bare=False,domain=None,target=None,random="pm1",ncpu=2,nrun=8,nper=1,var=False,save=False,path="tmp",prefix="",loop=False,**kwargs):
+ """
+ Computes the diagonal of the operator via probing.
+
+ Parameters
+ ----------
+ bare : bool, *optional*
+ Indicatese whether the diagonal entries are `bare` or not
+ (mandatory for the correct incorporation of volume weights)
+ (default: False)
+ 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).
+ save : bool, *optional*
+ whether all individual probing results are saved or not
+ (default: False)
+ path : string, *optional*
+ path wherein the results are saved (default: "tmp")
+ prefix : string, *optional*
+ prefix for all saved files (default: "")
+ loop : bool, *optional*
+ Indicates whether or not to perform a loop i.e., to
+ parallelise (default: False)
+
+ Returns
+ -------
+ diag : ndarray
+ The matrix diagonal
+ delta : float, *optional*
+ Probing variance of the trace. Returned if `var` is True in
+ of probing case.
+
+ See Also
+ --------
+ probing : The class used to perform probing operations
+
+ Notes
+ -----
+ The ambiguity of `bare` or non-bare diagonal entries is based
+ on the choice of a matrix representation of the operator in
+ question. The naive choice of absorbing the volume weights
+ into the matrix leads to a matrix-vector calculus with the
+ non-bare entries which seems intuitive, though. The choice of
+ keeping matrix entries and volume weights separate deals with the
+ bare entries that allow for correct interpretation of the matrix
+ entries; e.g., as variance in case of an covariance operator.
+
+ """
+ if(domain is None):
+ domain = self.domain
+ diag = diagonal_probing(self,function=self.times,domain=domain,target=target,random=random,ncpu=(ncpu,1)[bool(loop)],nrun=nrun,nper=nper,var=var,save=save,path=path,prefix=prefix,**kwargs)(loop=loop)
+ if(diag is None):
+# about.warnings.cprint("WARNING: forwarding 'NoneType'.")
+ return None
+ ## weight if ...
+ elif(not domain.discrete)and(bare):
+ if(isinstance(diag,tuple)): ## diag == (diag,variance)
+ return domain.calc_weight(diag[0],power=-1),domain.calc_weight(diag[1],power=-1)
+ else:
+ return domain.calc_weight(diag,power=-1)
+ else:
+ return diag
+
+ def inverse_diag(self,bare=False,domain=None,target=None,random="pm1",ncpu=2,nrun=8,nper=1,var=False,save=False,path="tmp",prefix="",loop=False,**kwargs):
+ """
+ Computes the diagonal of the inverse operator via probing.
+
+ Parameters
+ ----------
+ bare : bool, *optional*
+ Indicatese whether the diagonal entries are `bare` or not
+ (mandatory for the correct incorporation of volume weights)
+ (default: False)
+ 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).
+ save : bool, *optional*
+ whether all individual probing results are saved or not
+ (default: False)
+ path : string, *optional*
+ path wherein the results are saved (default: "tmp")
+ prefix : string, *optional*
+ prefix for all saved files (default: "")
+ loop : bool, *optional*
+ Indicates whether or not to perform a loop i.e., to
+ parallelise (default: False)
+
+ Returns
+ -------
+ diag : ndarray
+ The diagonal of the inverse matrix
+ delta : float, *optional*
+ Probing variance of the trace. Returned if `var` is True in
+ of probing case.
+
+ See Also
+ --------
+ probing : The class used to perform probing operations
+
+ Notes
+ -----
+ The ambiguity of `bare` or non-bare diagonal entries is based
+ on the choice of a matrix representation of the operator in
+ question. The naive choice of absorbing the volume weights
+ into the matrix leads to a matrix-vector calculus with the
+ non-bare entries which seems intuitive, though. The choice of
+ keeping matrix entries and volume weights separate deals with the
+ bare entries that allow for correct interpretation of the matrix
+ entries; e.g., as variance in case of an covariance operator.
+
+ """
+ if(domain is None):
+ domain = self.target
+ diag = diagonal_probing(self,function=self.inverse_times,domain=domain,target=target,random=random,ncpu=(ncpu,1)[bool(loop)],nrun=nrun,nper=nper,var=var,save=save,path=path,prefix=prefix,**kwargs)(loop=loop)
+ if(diag is None):
+# about.warnings.cprint("WARNING: forwarding 'NoneType'.")
+ return None
+ ## weight if ...
+ elif(not domain.discrete)and(bare):
+ if(isinstance(diag,tuple)): ## diag == (diag,variance)
+ return domain.calc_weight(diag[0],power=-1),domain.calc_weight(diag[1],power=-1)
+ else:
+ return domain.calc_weight(diag,power=-1)
+ else:
+ return diag
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def det(self):
+ """
+ Computes the determinant of the operator.
+
+ Returns
+ -------
+ det : float
+ The determinant
+
+ """
+ raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'det'."))
+
+ def inverse_det(self):
+ """
+ Computes the determinant of the inverse operator.
+
+ Returns
+ -------
+ det : float
+ The determinant
+
+ """
+ raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'inverse_det'."))
+
+ def log_det(self):
+ """
+ Computes the logarithm of the determinant of the operator (if applicable).
+
+ Returns
+ -------
+ logdet : float
+ The logarithm of the determinant
+
+ """
+ raise NotImplementedError(about._errors.cstring("ERROR: no generic instance method 'log_det'."))
+
+ def tr_log(self):
+ """
+ Computes the trace of the logarithm of the operator (if applicable).
+
+ Returns
+ -------
+ logdet : float
+ The trace of the logarithm
+
+ """
+ return self.log_det()
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def hat(self,bare=False,domain=None,target=None,**kwargs):
+ """
+ Translates the operator's diagonal into a field
+
+ Parameters
+ ----------
+ bare : bool, *optional*
+ Indicatese whether the diagonal entries are `bare` or not
+ (mandatory for the correct incorporation of volume weights)
+ (default: False)
+ 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)
+ save : bool, *optional*
+ whether all individual probing results are saved or not
+ (default: False)
+ path : string, *optional*
+ path wherein the results are saved (default: "tmp")
+ prefix : string, *optional*
+ prefix for all saved files (default: "")
+ loop : bool, *optional*
+ Indicates whether or not to perform a loop i.e., to
+ parallelise (default: False)
+
+ Returns
+ -------
+ x : field
+ The matrix diagonal as a field living on the operator
+ domain space
+
+ See Also
+ --------
+ probing : The class used to perform probing operations
+
+ Notes
+ -----
+ The ambiguity of `bare` or non-bare diagonal entries is based
+ on the choice of a matrix representation of the operator in
+ question. The naive choice of absorbing the volume weights
+ into the matrix leads to a matrix-vector calculus with the
+ non-bare entries which seems intuitive, though. The choice of
+ keeping matrix entries and volume weights separate deals with the
+ bare entries that allow for correct interpretation of the matrix
+ entries; e.g., as variance in case of an covariance operator.
+
+ """
+ if(domain is None):
+ domain = self.domain
+ diag = self.diag(bare=bare,domain=domain,target=target,var=False,**kwargs)
+ if(diag is None):
+ about.warnings.cprint("WARNING: forwarding 'NoneType'.")
+ return None
+ return field(domain,val=diag,target=target)
+
+ def inverse_hat(self,bare=False,domain=None,target=None,**kwargs):
+ """
+ Translates the inverse operator's diagonal into a field
+
+ Parameters
+ ----------
+ bare : bool, *optional*
+ Indicatese whether the diagonal entries are `bare` or not
+ (mandatory for the correct incorporation of volume weights)
+ (default: False)
+ 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)
+ save : bool, *optional*
+ whether all individual probing results are saved or not
+ (default: False)
+ path : string, *optional*
+ path wherein the results are saved (default: "tmp")
+ prefix : string, *optional*
+ prefix for all saved files (default: "")
+ loop : bool, *optional*
+ Indicates whether or not to perform a loop i.e., to
+ parallelise (default: False)
+
+ Returns
+ -------
+ x : field
+ The matrix diagonal as a field living on the operator
+ domain space
+
+ See Also
+ --------
+ probing : The class used to perform probing operations
+
+ Notes
+ -----
+ The ambiguity of `bare` or non-bare diagonal entries is based
+ on the choice of a matrix representation of the operator in
+ question. The naive choice of absorbing the volume weights
+ into the matrix leads to a matrix-vector calculus with the
+ non-bare entries which seems intuitive, though. The choice of
+ keeping matrix entries and volume weights separate deals with the
+ bare entries that allow for correct interpretation of the matrix
+ entries; e.g., as variance in case of an covariance operator.
+
+ """
+ if(domain is None):
+ domain = self.target
+ diag = self.inverse_diag(bare=bare,domain=domain,target=target,var=False,**kwargs)
+ if(diag is None):
+ about.warnings.cprint("WARNING: forwarding 'NoneType'.")
+ return None
+ return field(domain,val=diag,target=target)
+
+ def hathat(self,domain=None,**kwargs):
+ """
+ Translates the operator's diagonal into a diagonal 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)
+ save : bool, *optional*
+ whether all individual probing results are saved or not
+ (default: False)
+ path : string, *optional*
+ path wherein the results are saved (default: "tmp")
+ prefix : string, *optional*
+ prefix for all saved files (default: "")
+ loop : bool, *optional*
+ Indicates whether or not to perform a loop i.e., to
+ parallelise (default: False)
+
+ Returns
+ -------
+ D : diagonal_operator
+ The matrix diagonal as an operator
+
+ See Also
+ --------
+ probing : The class used to perform probing operations
+
+ Notes
+ -----
+ The ambiguity of `bare` or non-bare diagonal entries is based
+ on the choice of a matrix representation of the operator in
+ question. The naive choice of absorbing the volume weights
+ into the matrix leads to a matrix-vector calculus with the
+ non-bare entries which seems intuitive, though. The choice of
+ keeping matrix entries and volume weights separate deals with the
+ bare entries that allow for correct interpretation of the matrix
+ entries; e.g., as variance in case of an covariance operator.
+
+ """
+ if(domain is None):
+ domain = self.domain
+ diag = self.diag(bare=False,domain=domain,var=False,**kwargs)
+ if(diag is None):
+ about.warnings.cprint("WARNING: forwarding 'NoneType'.")
+ return None
+ return diagonal_operator(domain=domain,diag=diag,bare=False)
+
+ def inverse_hathat(self,domain=None,**kwargs):
+ """
+ Translates the inverse operator's diagonal into a diagonal
+ 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)
+ save : bool, *optional*
+ whether all individual probing results are saved or not
+ (default: False)
+ path : string, *optional*
+ path wherein the results are saved (default: "tmp")
+ prefix : string, *optional*
+ prefix for all saved files (default: "")
+ loop : bool, *optional*
+ Indicates whether or not to perform a loop i.e., to
+ parallelise (default: False)
+
+ Returns
+ -------
+ D : diagonal_operator
+ The diagonal of the inverse matrix as an operator
+
+ See Also
+ --------
+ probing : The class used to perform probing operations
+
+ Notes
+ -----
+ The ambiguity of `bare` or non-bare diagonal entries is based
+ on the choice of a matrix representation of the operator in
+ question. The naive choice of absorbing the volume weights
+ into the matrix leads to a matrix-vector calculus with the
+ non-bare entries which seems intuitive, though. The choice of
+ keeping matrix entries and volume weights separate deals with the
+ bare entries that allow for correct interpretation of the matrix
+ entries; e.g., as variance in case of an covariance operator.
+
+ """
+ if(domain is None):
+ domain = self.target
+ diag = self.inverse_diag(bare=False,domain=domain,var=False,**kwargs)
+ if(diag is None):
+ about.warnings.cprint("WARNING: forwarding 'NoneType'.")
+ return None
+ return diagonal_operator(domain=domain,diag=diag,bare=False)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def __repr__(self):
+ return "<nifty_core.operator>"
+
+##=============================================================================
+
+
+
+##-----------------------------------------------------------------------------
+
+class diagonal_operator(operator):
+ """
+ .. __ __ __
+ .. / / /__/ / /
+ .. ____/ / __ ____ __ ____ __ ______ __ ___ ____ __ / /
+ .. / _ / / / / _ / / _ / / _ | / _ | / _ / / /
+ .. / /_/ / / / / /_/ / / /_/ / / /_/ / / / / / / /_/ / / /_
+ .. \______| /__/ \______| \___ / \______/ /__/ /__/ \______| \___/ operator class
+ .. /______/
+
+ NIFTY subclass for diagonal operators
+
+ Parameters
+ ----------
+ domain : space, *optional*
+ The space wherein valid arguments live. If no domain is given
+ then the diag parameter *must* be a field and the domain
+ of that field is assumed. (default: None)
+ diag : {scalar, ndarray, field}
+ The diagonal entries of the operator. For a scalar, a constant
+ diagonal is defined having the value provided. If no domain
+ is given, diag must be a field. (default: 1)
+ bare : bool, *optional*
+ whether the diagonal entries are `bare` or not
+ (mandatory for the correct incorporation of volume weights)
+ (default: False)
+
+ Notes
+ -----
+ The ambiguity of `bare` or non-bare diagonal entries is based
+ on the choice of a matrix representation of the operator in
+ question. The naive choice of absorbing the volume weights
+ into the matrix leads to a matrix-vector calculus with the
+ non-bare entries which seems intuitive, though. The choice of
+ keeping matrix entries and volume weights separate deals with the
+ bare entries that allow for correct interpretation of the matrix
+ entries; e.g., as variance in case of an covariance operator.
+
+ The inverse applications of the diagonal operator feature a ``pseudo``
+ flag indicating if zero divison shall be ignored and return zero
+ instead of causing an error.
+
+ Attributes
+ ----------
+ domain : space
+ The space wherein valid arguments live.
+ val : ndarray
+ A field containing the diagonal entries of the matrix.
+ sym : bool
+ Indicates whether the operator is self-adjoint or not
+ uni : bool
+ Indicates whether the operator is unitary or not
+ imp : bool
+ Indicates whether the incorporation of volume weights in
+ multiplications is already implemented in the `multiply`
+ instance methods or not
+ target : space
+ The space wherein the operator output lives
+ """
+ def __init__(self,domain=None,diag=1,bare=False):
+ """
+ Sets the standard operator properties and `values`.
+
+ Parameters
+ ----------
+ domain : space, *optional*
+ The space wherein valid arguments live. If no domain is given
+ then the diag parameter *must* be a field and the domain
+ of that field is assumed. (default: None)
+ diag : {scalar, ndarray, field}, *optional*
+ The diagonal entries of the operator. For a scalar, a constant
+ diagonal is defined having the value provided. If no domain
+ is given, diag must be a field. (default: 1)
+ bare : bool, *optional*
+ whether the diagonal entries are `bare` or not
+ (mandatory for the correct incorporation of volume weights)
+ (default: False)
+
+ Returns
+ -------
+ None
+
+ Notes
+ -----
+ The ambiguity of `bare` or non-bare diagonal entries is based
+ on the choice of a matrix representation of the operator in
+ question. The naive choice of absorbing the volume weights
+ into the matrix leads to a matrix-vector calculus with the
+ non-bare entries which seems intuitive, though. The choice of
+ keeping matrix entries and volume weights separate deals with the
+ bare entries that allow for correct interpretation of the matrix
+ entries; e.g., as variance in case of an covariance operator.
+
+ """
+ if(domain is None)and(isinstance(diag,field)):
+ domain = diag.domain
+ if(not isinstance(domain,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ self.domain = domain
+
+ diag = self.domain.enforce_values(diag,extend=True)
+ ## weight if ...
+ if(not self.domain.discrete)and(bare):
+ diag = self.domain.calc_weight(diag,power=1)
+ ## check complexity
+ if(np.all(np.imag(diag)==0)):
+ self.val = np.real(diag)
+ self.sym = True
+ else:
+ self.val = diag
+# about.infos.cprint("INFO: non-self-adjoint complex diagonal operator.")
+ self.sym = False
+
+ ## check whether identity
+ if(np.all(diag==1)):
+ self.uni = True
+ else:
+ self.uni = False
+
+ self.imp = True ## correctly implemented for efficiency
+ self.target = self.domain
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def set_diag(self,newdiag,bare=False):
+ """
+ Sets the diagonal of the diagonal operator
+
+ Parameters
+ ----------
+ newdiag : {scalar, ndarray, field}
+ The new diagonal entries of the operator. For a scalar, a
+ constant diagonal is defined having the value provided. If
+ no domain is given, diag must be a field.
+
+ bare : bool, *optional*
+ whether the diagonal entries are `bare` or not
+ (mandatory for the correct incorporation of volume weights)
+ (default: False)
+
+ Returns
+ -------
+ None
+ """
+ newdiag = self.domain.enforce_values(newdiag,extend=True)
+ ## weight if ...
+ if(not self.domain.discrete)and(bare):
+ newdiag = self.domain.calc_weight(newdiag,power=1)
+ ## check complexity
+ if(np.all(np.imag(newdiag)==0)):
+ self.val = np.real(newdiag)
+ self.sym = True
+ else:
+ self.val = newdiag
+# about.infos.cprint("INFO: non-self-adjoint complex diagonal operator.")
+ self.sym = False
+
+ ## check whether identity
+ if(np.all(newdiag==1)):
+ self.uni = True
+ else:
+ self.uni = False
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def _multiply(self,x,**kwargs): ## > applies the operator to a given field
+ x_ = field(self.target,val=None,target=x.target)
+ x_.val = x.val*self.val ## bypasses self.domain.enforce_values
+ return x_
+
+ def _adjoint_multiply(self,x,**kwargs): ## > applies the adjoint operator to a given field
+ x_ = field(self.domain,val=None,target=x.target)
+ x_.val = x.val*np.conjugate(self.val) ## bypasses self.domain.enforce_values
+ return x_
+
+ def _inverse_multiply(self,x,pseudo=False,**kwargs): ## > applies the inverse operator to a given field
+ if(np.any(self.val==0)):
+ if(pseudo):
+ x_ = field(self.domain,val=None,target=x.target)
+ x_.val = np.ma.filled(x.val/np.ma.masked_where(self.val==0,self.val,copy=False),fill_value=0) ## bypasses self.domain.enforce_values
+ return x_
+ else:
+ raise AttributeError(about._errors.cstring("ERROR: singular operator."))
+ else:
+ x_ = field(self.domain,val=None,target=x.target)
+ x_.val = x.val/self.val ## bypasses self.domain.enforce_values
+ return x_
+
+ def _adjoint_inverse_multiply(self,x,pseudo=False,**kwargs): ## > applies the inverse adjoint operator to a given field
+ if(np.any(self.val==0)):
+ if(pseudo):
+ x_ = field(self.domain,val=None,target=x.target)
+ x_.val = np.ma.filled(x.val/np.ma.masked_where(self.val==0,np.conjugate(self.val),copy=False),fill_value=0) ## bypasses self.domain.enforce_values
+ return x_
+ else:
+ raise AttributeError(about._errors.cstring("ERROR: singular operator."))
+ else:
+ x_ = field(self.target,val=None,target=x.target)
+ x_.val = x.val/np.conjugate(self.val) ## bypasses self.domain.enforce_values
+ return x_
+
+ def _inverse_adjoint_multiply(self,x,pseudo=False,**kwargs): ## > applies the adjoint inverse operator to a given field
+ if(np.any(self.val==0)):
+ if(pseudo):
+ x_ = field(self.domain,val=None,target=x.target)
+ x_.val = np.ma.filled(x.val/np.conjugate(np.ma.masked_where(self.val==0,self.val,copy=False)),fill_value=0) ## bypasses self.domain.enforce_values
+ return x_
+ else:
+ raise AttributeError(about._errors.cstring("ERROR: singular operator."))
+
+ else:
+ x_ = field(self.target,val=None,target=x.target)
+ x_.val = x.val*np.conjugate(1/self.val) ## bypasses self.domain.enforce_values
+ return x_
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ 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(domain is None)or(domain==self.domain):
+ if(self.uni): ## identity
+ return (self.domain.datatype(self.domain.dof())).real
+ elif(self.domain.dim(split=False)<self.domain.dof()): ## hidden degrees of freedom
+ return self.domain.calc_dot(np.ones(self.domain.dim(split=True),dtype=self.domain.datatype,order='C'),self.val) ## discrete inner product
+ else:
+ return np.sum(self.val,axis=None,dtype=None,out=None)
+ else:
+ if(self.uni): ## identity
+ if(not isinstance(domain,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ ## check degrees of freedom
+ if(self.domain.dof()>domain.dof()):
+ about.infos.cprint("INFO: variant numbers of degrees of freedom ( "+str(self.domain.dof())+" / "+str(domain.dof())+" ).")
+ return (domain.datatype(domain.dof())).real
+ else:
+ return super(diagonal_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)
+ 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 inverse operator
+ delta : float, *optional*
+ Probing variance of the trace. Returned if `var` is True in
+ of probing case.
+
+ """
+ if(np.any(self.val==0)):
+ raise AttributeError(about._errors.cstring("ERROR: singular operator."))
+
+ if(domain is None)or(domain==self.target):
+ if(self.uni): ## identity
+ return np.real(self.domain.datatype(self.domain.dof()))
+ elif(self.domain.dim(split=False)<self.domain.dof()): ## hidden degrees of freedom
+ return self.domain.calc_dot(np.ones(self.domain.dim(split=True),dtype=self.domain.datatype,order='C'),1/self.val) ## discrete inner product
+ else:
+ return np.sum(1/self.val,axis=None,dtype=None,out=None)
+ else:
+ if(self.uni): ## identity
+ if(not isinstance(domain,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ ## check degrees of freedom
+ if(self.domain.dof()>domain.dof()):
+ about.infos.cprint("INFO: variant numbers of degrees of freedom ( "+str(self.domain.dof())+" / "+str(domain.dof())+" ).")
+ return np.real(domain.datatype(domain.dof()))
+ else:
+ return super(diagonal_operator,self).inverse_tr(domain=domain,**kwargs) ## probing
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def diag(self,bare=False,domain=None,**kwargs):
+ """
+ Computes the diagonal of the operator.
+
+ Parameters
+ ----------
+ bare : bool, *optional*
+ Indicatese whether the diagonal entries are `bare` or not
+ (mandatory for the correct incorporation of volume weights)
+ (default: False)
+ 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).
+ save : bool, *optional*
+ whether all individual probing results are saved or not
+ (default: False)
+ path : string, *optional*
+ path wherein the results are saved (default: "tmp")
+ prefix : string, *optional*
+ prefix for all saved files (default: "")
+ loop : bool, *optional*
+ Indicates whether or not to perform a loop i.e., to
+ parallelise (default: False)
+
+ Returns
+ -------
+ diag : ndarray
+ The matrix diagonal
+ delta : float, *optional*
+ Probing variance of the trace. Returned if `var` is True in
+ of probing case.
+
+ Notes
+ -----
+ The ambiguity of `bare` or non-bare diagonal entries is based
+ on the choice of a matrix representation of the operator in
+ question. The naive choice of absorbing the volume weights
+ into the matrix leads to a matrix-vector calculus with the
+ non-bare entries which seems intuitive, though. The choice of
+ keeping matrix entries and volume weights separate deals with the
+ bare entries that allow for correct interpretation of the matrix
+ entries; e.g., as variance in case of an covariance operator.
+
+ """
+ if(domain is None)or(domain==self.domain):
+ ## weight if ...
+ if(not self.domain.discrete)and(bare):
+ diag = self.domain.calc_weight(self.val,power=-1)
+ ## check complexity
+ if(np.all(np.imag(diag)==0)):
+ diag = np.real(diag)
+ return diag
+ else:
+ return self.val
+ else:
+ if(self.uni): ## identity
+ if(not isinstance(domain,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ ## weight if ...
+ if(not domain.discrete)and(bare):
+ return np.real(domain.calc_weight(domain.enforce_values(1,extend=True),power=-1))
+ else:
+ return np.real(domain.enforce_values(1,extend=True))
+ else:
+ return super(diagonal_operator,self).diag(bare=bare,domain=domain,**kwargs) ## probing
+
+ def inverse_diag(self,bare=False,domain=None,**kwargs):
+ """
+ Computes the diagonal of the inverse operator.
+
+ Parameters
+ ----------
+ bare : bool, *optional*
+ Indicatese whether the diagonal entries are `bare` or not
+ (mandatory for the correct incorporation of volume weights)
+ (default: False)
+ 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).
+ save : bool, *optional*
+ whether all individual probing results are saved or not
+ (default: False)
+ path : string, *optional*
+ path wherein the results are saved (default: "tmp")
+ prefix : string, *optional*
+ prefix for all saved files (default: "")
+ loop : bool, *optional*
+ Indicates whether or not to perform a loop i.e., to
+ parallelise (default: False)
+
+ Returns
+ -------
+ diag : ndarray
+ The diagonal of the inverse matrix
+ delta : float, *optional*
+ Probing variance of the trace. Returned if `var` is True in
+ of probing case.
+
+ See Also
+ --------
+ probing : The class used to perform probing operations
+
+ Notes
+ -----
+ The ambiguity of `bare` or non-bare diagonal entries is based
+ on the choice of a matrix representation of the operator in
+ question. The naive choice of absorbing the volume weights
+ into the matrix leads to a matrix-vector calculus with the
+ non-bare entries which seems intuitive, though. The choice of
+ keeping matrix entries and volume weights separate deals with the
+ bare entries that allow for correct interpretation of the matrix
+ entries; e.g., as variance in case of an covariance operator.
+
+ """
+ if(domain is None)or(domain==self.target):
+ ## weight if ...
+ if(not self.domain.discrete)and(bare):
+ diag = self.domain.calc_weight(1/self.val,power=-1)
+ ## check complexity
+ if(np.all(np.imag(diag)==0)):
+ diag = np.real(diag)
+ return diag
+ else:
+ return 1/self.val
+ else:
+ if(self.uni): ## identity
+ if(not isinstance(domain,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ ## weight if ...
+ if(not domain.discrete)and(bare):
+ return np.real(domain.calc_weight(domain.enforce_values(1,extend=True),power=-1))
+ else:
+ return np.real(domain.enforce_values(1,extend=True))
+ else:
+ return super(diagonal_operator,self).inverse_diag(bare=bare,domain=domain,**kwargs) ## probing
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def det(self):
+ """
+ Computes the determinant of the matrix.
+
+ Returns
+ -------
+ det : float
+ The determinant
+
+ """
+ if(self.uni): ## identity
+ return 1
+ elif(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
+
+ """
+ if(self.uni): ## identity
+ return 1
+ det = self.det()
+ if(det<>0):
+ return 1/det
+ else:
+ raise ValueError(about._errors.cstring("ERROR: singular operator."))
+
+ def log_det(self):
+ """
+ Computes the logarithm of the determinant of the operator.
+
+ Returns
+ -------
+ logdet : float
+ The logarithm of the determinant
+
+ """
+ if(self.uni): ## identity
+ return 0
+ elif(self.domain.dim(split=False)<self.domain.dof()): ## hidden degrees of freedom
+ return self.domain.calc_dot(np.ones(self.domain.dim(split=True),dtype=self.domain.datatype,order='C'),np.log(self.val))
+ else:
+ return np.sum(np.log(self.val),axis=None,dtype=None,out=None)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def get_random_field(self,domain=None,target=None):
+ """
+ Generates a Gaussian random field with variance equal to the
+ diagonal.
+
+ Parameters
+ ----------
+ domain : space, *optional*
+ space wherein the field lives (default: None, indicates
+ to use self.domain)
+ target : space, *optional*
+ space wherein the transform of the field lives
+ (default: None, indicates to use target of domain)
+
+ Returns
+ -------
+ x : field
+ Random field.
+
+ """
+ ## weight if ...
+ if(not self.domain.discrete):
+ diag = self.domain.calc_weight(self.val,power=-1)
+ ## check complexity
+ if(np.all(np.imag(diag)==0)):
+ diag = np.real(diag)
+ else:
+ diag = self.val
+
+ if(domain is None)or(domain==self.domain):
+ return field(self.domain,val=None,target=target,random="gau",var=self.diag(bare=True,domain=self.domain))
+ else:
+ return field(self.domain,val=None,target=domain,random="gau",var=self.diag(bare=True,domain=self.domain)).transform(target=domain,overwrite=False)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def __repr__(self):
+ return "<nifty_core.diagonal_operator>"
+
+##-----------------------------------------------------------------------------
+
+##-----------------------------------------------------------------------------
+
+def identity(domain):
+ """
+ Returns an identity operator.
+
+ The identity operator is represented by a `diagonal_operator` instance,
+ which is applicable to a field-like object; i.e., a scalar, list,
+ array or field. (The identity operator is unrelated to PYTHON's
+ built-in function :py:func:`id`.)
+
+ Parameters
+ ----------
+ domain : space
+ The space wherein valid arguments live.
+
+ Returns
+ -------
+ id : diagonal_operator
+ The identity operator as a `diagonal_operator` instance.
+
+ See Also
+ --------
+ diagonal_operator
+
+ Examples
+ --------
+ >>> I = identity(rg_space(8,dist=0.2))
+ >>> I.diag()
+ array([ 1., 1., 1., 1., 1., 1., 1., 1.])
+ >>> I.diag(bare=True)
+ array([ 5., 5., 5., 5., 5., 5., 5., 5.])
+ >>> I.tr()
+ 8.0
+ >>> I(3)
+ <nifty.field>
+ >>> I(3).val
+ array([ 3., 3., 3., 3., 3., 3., 3., 3.])
+ >>> I(np.arange(8))[:]
+ array([ 0., 1., 2., 3., 4., 5., 6., 7.])
+ >>> f = I.get_random_field()
+ >>> print(I(f) - f)
+ nifty.field instance
+ - domain = <nifty.rg_space>
+ - val = [...]
+ - min.,max. = [0.0, 0.0]
+ - med.,mean = [0.0, 0.0]
+ - target = <nifty.rg_space>
+ >>> I.times(f) ## equal to I(f)
+ <nifty.field>
+ >>> I.inverse_times(f)
+ <nifty.field>
+
+ """
+ return diagonal_operator(domain=domain,diag=1,bare=False)
+
+##-----------------------------------------------------------------------------
+
+##-----------------------------------------------------------------------------
+
+class power_operator(diagonal_operator):
+ """
+ .. ______ ______ __ __ _______ _____
+ .. / _ | / _ | | |/\/ / / __ / / __/
+ .. / /_/ / / /_/ / | / / /____/ / /
+ .. / ____/ \______/ |__/\__/ \______/ /__/ operator class
+ .. /__/
+
+ NIFTY subclass for (signal-covariance-type) diagonal operators containing a power spectrum
+
+ Parameters
+ ----------
+ domain : space, *optional*
+ The space wherein valid arguments live. If no domain is given
+ then the diag parameter *must* be a field and the domain
+ of that field is assumed. (default: None)
+ spec : {scalar, list, array, field, function}
+ The power spectrum. For a scalar, a constant power
+ spectrum is defined having the value provided. If no domain
+ is given, diag must be a field. (default: 1)
+ bare : bool, *optional*
+ whether the entries are `bare` or not
+ (mandatory for the correct incorporation of volume weights)
+ (default: True)
+ pindex : ndarray, *optional*
+ indexing array, obtainable from domain.get_power_indices
+ (default: None)
+
+ Other Parameters
+ ----------------
+ log : bool, *optional*
+ Flag specifying if the spectral binning is performed on logarithmic
+ scale or not; if set, the number of used bins is set
+ automatically (if not given otherwise); by default no binning
+ is done (default: None).
+ nbin : integer, *optional*
+ Number of used spectral bins; if given `log` is set to ``False``;
+ integers below the minimum of 3 induce an automatic setting;
+ by default no binning is done (default: None).
+ binbounds : {list, array}, *optional*
+ User specific inner boundaries of the bins, which are preferred
+ over the above parameters; by default no binning is done
+ (default: None). vmin : {scalar, list, ndarray, field}, *optional*
+ Lower limit of the uniform distribution if ``random == "uni"``
+ (default: 0).
+
+ Notes
+ -----
+ The ambiguity of `bare` or non-bare diagonal entries is based
+ on the choice of a matrix representation of the operator in
+ question. The naive choice of absorbing the volume weights
+ into the matrix leads to a matrix-vector calculus with the
+ non-bare entries which seems intuitive, though. The choice of
+ keeping matrix entries and volume weights separate deals with the
+ bare entries that allow for correct interpretation of the matrix
+ entries; e.g., as variance in case of an covariance operator.
+
+ Attributes
+ ----------
+ domain : space
+ The space wherein valid arguments live.
+ val : ndarray
+ A field containing the diagonal entries of the matrix.
+ sym : bool
+ Indicates whether the operator is self-adjoint or not
+ uni : bool
+ Indicates whether the operator is unitary or not
+ imp : bool
+ Indicates whether the incorporation of volume weights in
+ multiplications is already implemented in the `multiply`
+ instance methods or not
+ target : space
+ The space wherein the operator output lives
+
+ """
+ def __init__(self,domain,spec=1,bare=True,pindex=None,**kwargs):
+ """
+ Sets the diagonal operator's standard properties
+
+ Parameters
+ ----------
+ domain : space, *optional*
+ The space wherein valid arguments live. If no domain is given
+ then the diag parameter *must* be a field and the domain
+ of that field is assumed. (default: None)
+ spec : {scalar, list, array, field, function}
+ The power spectrum. For a scalar, a constant power
+ spectrum is defined having the value provided. If no domain
+ is given, diag must be a field. (default: 1)
+ bare : bool, *optional*
+ whether the entries are `bare` or not
+ (mandatory for the correct incorporation of volume weights)
+ (default: True)
+ pindex : ndarray, *optional*
+ indexing array, obtainable from domain.get_power_indices
+ (default: None)
+
+ Returns
+ -------
+ None
+
+ Other Parameters
+ ----------------
+ log : bool, *optional*
+ Flag specifying if the spectral binning is performed on logarithmic
+ scale or not; if set, the number of used bins is set
+ automatically (if not given otherwise); by default no binning
+ is done (default: None).
+ nbin : integer, *optional*
+ Number of used spectral bins; if given `log` is set to ``False``;
+ integers below the minimum of 3 induce an automatic setting;
+ by default no binning is done (default: None).
+ binbounds : {list, array}, *optional*
+ User specific inner boundaries of the bins, which are preferred
+ over the above parameters; by default no binning is done
+ (default: None). vmin : {scalar, list, ndarray, field}, *optional*
+ Lower limit of the uniform distribution if ``random == "uni"``
+ (default: 0).
+
+ """
+ if(not isinstance(domain,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ self.domain = domain
+ ## check implicit pindex
+ if(pindex is None):
+ try:
+ self.domain.set_power_indices(**kwargs)
+ except:
+ raise ValueError(about._errors.cstring("ERROR: invalid input."))
+ else:
+ pindex = self.domain.power_indices.get("pindex")
+ ## check explicit pindex
+ else:
+ pindex = np.array(pindex,dtype=np.int)
+ if(not np.all(np.array(np.shape(pindex))==self.domain.dim(split=True))):
+ raise ValueError(about._errors.cstring("ERROR: shape mismatch ( "+str(np.array(np.shape(pindex)))+" <> "+str(self.domain.dim(split=True))+" )."))
+ ## set diagonal
+ try:
+ #diag = self.domain.enforce_power(spec,size=np.max(pindex,axis=None,out=None)+1)[pindex]
+ temp_spec = self.domain.enforce_power(
+ spec,size=np.max(pindex,axis=None,out=None)+1)
+ diag = pindex.apply_scalar_function(lambda x: temp_spec[x])
+ except(AttributeError):
+ raise ValueError(about._errors.cstring("ERROR: invalid input."))
+ ## weight if ...
+ if(not self.domain.discrete)and(bare):
+ self.val = np.real(self.domain.calc_weight(diag,power=1))
+ else:
+ self.val = diag
+
+ self.sym = True
+
+ ## check whether identity
+ if(np.all(spec==1)):
+ self.uni = True
+ else:
+ self.uni = False
+
+ self.imp = True
+ self.target = self.domain
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def set_power(self,newspec,bare=True,pindex=None,**kwargs):
+ """
+ Sets the power spectrum of the diagonal operator
+
+ Parameters
+ ----------
+ newspec : {scalar, list, array, field, function}
+ The entries of the operator. For a scalar, a constant
+ diagonal is defined having the value provided. If no domain
+ is given, diag must be a field. (default: 1)
+ bare : bool
+ whether the entries are `bare` or not
+ (mandatory for the correct incorporation of volume weights)
+ pindex : ndarray, *optional*
+ indexing array, obtainable from domain.get_power_indices
+ (default: None)
+
+ Returns
+ -------
+ None
+
+ Other Parameters
+ ----------------
+ log : bool, *optional*
+ Flag specifying if the spectral binning is performed on logarithmic
+ scale or not; if set, the number of used bins is set
+ automatically (if not given otherwise); by default no binning
+ is done (default: None).
+ nbin : integer, *optional*
+ Number of used spectral bins; if given `log` is set to ``False``;
+ integers below the minimum of 3 induce an automatic setting;
+ by default no binning is done (default: None).
+ binbounds : {list, array}, *optional*
+ User specific inner boundaries of the bins, which are preferred
+ over the above parameters; by default no binning is done
+ (default: None).
+
+ """
+# if(bare is None):
+# about.warnings.cprint("WARNING: bare keyword set to default.")
+# bare = True
+ ## check implicit pindex
+ if(pindex is None):
+ try:
+ self.domain.set_power_indices(**kwargs)
+ except:
+ raise ValueError(about._errors.cstring("ERROR: invalid domain."))
+ else:
+ pindex = self.domain.power_indices.get("pindex")
+ ## check explicit pindex
+ else:
+ pindex = np.array(pindex,dtype=np.int)
+ if(not np.all(np.array(np.shape(pindex))==self.domain.dim(split=True))):
+ raise ValueError(about._errors.cstring("ERROR: shape mismatch ( "+str(np.array(np.shape(pindex)))+" <> "+str(self.domain.dim(split=True))+" )."))
+ ## set diagonal
+ try:
+ diag = self.domain.enforce_power(newspec,size=np.max(pindex,axis=None,out=None)+1)[pindex]
+ except(AttributeError):
+ raise ValueError(about._errors.cstring("ERROR: invalid input."))
+ ## weight if ...
+ if(not self.domain.discrete)and(bare):
+ self.val = np.real(self.domain.calc_weight(diag,power=1))
+ else:
+ self.val = diag
+
+ ## check whether identity
+ if(np.all(newspec==1)):
+ self.uni = True
+ else:
+ self.uni = False
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def get_power(self,bare=True,pundex=None,pindex=None,**kwargs):
+ """
+ Computes the power spectrum
+
+ Parameters
+ ----------
+ bare : bool, *optional*
+ whether the entries are `bare` or not
+ (mandatory for the correct incorporation of volume weights)
+ (default: True)
+ pundex : ndarray, *optional*
+ unindexing array, obtainable from domain.get_power_indices
+ (default: None)
+ pindex : ndarray, *optional*
+ indexing array, obtainable from domain.get_power_indices
+ (default: None)
+
+ Returns
+ -------
+ spec : ndarray
+ The power spectrum
+
+ Other Parameters
+ ----------------
+ log : bool, *optional*
+ Flag specifying if the spectral binning is performed on logarithmic
+ scale or not; if set, the number of used bins is set
+ automatically (if not given otherwise); by default no binning
+ is done (default: None).
+ nbin : integer, *optional*
+ Number of used spectral bins; if given `log` is set to ``False``;
+ integers below the minimum of 3 induce an automatic setting;
+ by default no binning is done (default: None).
+ binbounds : {list, array}, *optional*
+ User specific inner boundaries of the bins, which are preferred
+ over the above parameters; by default no binning is done
+ (default: None). vmin : {scalar, list, ndarray, field}, *optional*
+ Lower limit of the uniform distribution if ``random == "uni"``
+ (default: 0).
+
+ """
+ ## weight if ...
+ if(not self.domain.discrete)and(bare):
+ diag = np.real(self.domain.calc_weight(self.val,power=-1))
+ else:
+ diag = self.val
+ ## check implicit pundex
+ if(pundex is None):
+ if(pindex is None):
+ try:
+ self.domain.set_power_indices(**kwargs)
+ except:
+ raise ValueError(about._errors.cstring("ERROR: invalid domain."))
+ else:
+ pundex = self.domain.power_indices.get("pundex")
+ else:
+ pindex = np.array(pindex,dtype=np.int)
+ if(not np.all(np.array(np.shape(pindex))==self.domain.dim(split=True))):
+ raise ValueError(about._errors.cstring("ERROR: shape mismatch ( "+str(np.array(np.shape(pindex)))+" <> "+str(self.domain.dim(split=True))+" )."))
+ ## quick pundex
+ pundex = np.unique(pindex,return_index=True,return_inverse=False)[1]
+ ## check explicit pundex
+ else:
+ pundex = np.array(pundex,dtype=np.int)
+
+ return diag.flatten(order='C')[pundex]
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def get_projection_operator(self,pindex=None,**kwargs):
+ """
+ Generates a spectral projection operator
+
+ Parameters
+ ----------
+ pindex : ndarray
+ indexing array obtainable from domain.get_power_indices
+ (default: None)
+
+ Returns
+ -------
+ P : projection_operator
+
+ Other Parameters
+ ----------------
+ log : bool, *optional*
+ Flag specifying if the spectral binning is performed on logarithmic
+ scale or not; if set, the number of used bins is set
+ automatically (if not given otherwise); by default no binning
+ is done (default: None).
+ nbin : integer, *optional*
+ Number of used spectral bins; if given `log` is set to ``False``;
+ integers below the minimum of 3 induce an automatic setting;
+ by default no binning is done (default: None).
+ binbounds : {list, array}, *optional*
+ User specific inner boundaries of the bins, which are preferred
+ over the above parameters; by default no binning is done
+ (default: None). vmin : {scalar, list, ndarray, field}, *optional*
+ Lower limit of the uniform distribution if ``random == "uni"``
+ (default: 0).
+
+ """
+ ## check implicit pindex
+ if(pindex is None):
+ try:
+ self.domain.set_power_indices(**kwargs)
+ except:
+ raise ValueError(about._errors.cstring("ERROR: invalid domain."))
+ else:
+ pindex = self.domain.power_indices.get("pindex")
+ ## check explicit pindex
+ else:
+ pindex = np.array(pindex,dtype=np.int)
+ if(not np.all(np.array(np.shape(pindex))==self.domain.dim(split=True))):
+ raise ValueError(about._errors.cstring("ERROR: shape mismatch ( "+str(np.array(np.shape(pindex)))+" <> "+str(self.domain.dim(split=True))+" )."))
+
+ return projection_operator(self.domain,assign=pindex)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def __repr__(self):
+ return "<nifty_core.power_operator>"
+
+##-----------------------------------------------------------------------------
+
+
+
+##-----------------------------------------------------------------------------
+
+class projection_operator(operator):
+ """
+ .. __ __ __
+ .. /__/ / /_ /__/
+ .. ______ _____ ______ __ _______ _______ / _/ __ ______ __ ___
+ .. / _ | / __/ / _ | / / / __ / / ____/ / / / / / _ | / _ |
+ .. / /_/ / / / / /_/ / / / / /____/ / /____ / /_ / / / /_/ / / / / /
+ .. / ____/ /__/ \______/ / / \______/ \______/ \___/ /__/ \______/ /__/ /__/ operator class
+ .. /__/ /___/
+
+ NIFTY subclass for projection operators
+
+ Parameters
+ ----------
+ domain : space
+ The space wherein valid arguments live.
+ assign : ndarray, *optional*
+ Assignments of domain items to projection bands. An array
+ of integers, negative integers are associated with the
+ nullspace of the projection. (default: None)
+
+ Other Parameters
+ ----------------
+ log : bool, *optional*
+ Flag specifying if the spectral binning is performed on logarithmic
+ scale or not; if set, the number of used bins is set
+ automatically (if not given otherwise); by default no binning
+ is done (default: None).
+ nbin : integer, *optional*
+ Number of used spectral bins; if given `log` is set to ``False``;
+ integers below the minimum of 3 induce an automatic setting;
+ by default no binning is done (default: None).
+ binbounds : {list, array}, *optional*
+ User specific inner boundaries of the bins, which are preferred
+ over the above parameters; by default no binning is done
+ (default: None). vmin : {scalar, list, ndarray, field}, *optional*
+ Lower limit of the uniform distribution if ``random == "uni"``
+ (default: 0).
+
+ Notes
+ -----
+ The application of the projection operator features a ``band`` keyword
+ specifying a single projection band (see examples), a ``bandsup``
+ keyword specifying which projection bands to sum up, and a ``split``
+ keyword.
+
+ Examples
+ --------
+ >>> space = point_space(3)
+ >>> P = projection_operator(space, assign=[0, 1, 0])
+ >>> P.bands()
+ 2
+ >>> P([1, 2, 3], band=0) # equal to P.times(field(space,val=[1, 2, 3]))
+ <nifty.field>
+ >>> P([1, 2, 3], band=0).domain
+ <nifty.point_space>
+ >>> P([1, 2, 3], band=0).val # projection on band 0 (items 0 and 2)
+ array([ 1., 0., 3.])
+ >>> P([1, 2, 3], band=1).val # projection on band 1 (item 1)
+ array([ 0., 2., 0.])
+ >>> P([1, 2, 3])
+ <nifty.field>
+ >>> P([1, 2, 3]).domain
+ <nifty.nested_space>
+ >>> P([1, 2, 3]).val # projection on all bands
+ array([[ 1., 0., 3.],
+ [ 0., 2., 0.]])
+
+ Attributes
+ ----------
+ domain : space
+ The space wherein valid arguments live.
+ ind : ndarray
+ Assignments of domain items to projection bands.
+ sym : bool
+ Indicates whether the operator is self-adjoint or not
+ uni : bool
+ Indicates whether the operator is unitary or not
+ imp : bool
+ Indicates whether the incorporation of volume weights in
+ multiplications is already implemented in the `multiply`
+ instance methods or not
+ target : space
+ The space wherein the operator output lives
+
+ """
+ def __init__(self,domain,assign=None,**kwargs):
+ """
+ Sets the standard operator properties and `indexing`.
+
+ Parameters
+ ----------
+ domain : space
+ The space wherein valid arguments live.
+ assign : ndarray, *optional*
+ Assignments of domain items to projection bands. An array
+ of integers, negative integers are associated with the
+ nullspace of the projection. (default: None)
+
+ Returns
+ -------
+ None
+
+ Other Parameters
+ ----------------
+ log : bool, *optional*
+ Flag specifying if the spectral binning is performed on logarithmic
+ scale or not; if set, the number of used bins is set
+ automatically (if not given otherwise); by default no binning
+ is done (default: None).
+ nbin : integer, *optional*
+ Number of used spectral bins; if given `log` is set to ``False``;
+ integers below the minimum of 3 induce an automatic setting;
+ by default no binning is done (default: None).
+ binbounds : {list, array}, *optional*
+ User specific inner boundaries of the bins, which are preferred
+ over the above parameters; by default no binning is done
+ (default: None). vmin : {scalar, list, ndarray, field}, *optional*
+ Lower limit of the uniform distribution if ``random == "uni"``
+ (default: 0).
+
+ """
+ if(not isinstance(domain,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ self.domain = domain
+
+ ## check assignment(s)
+ if(assign is None):
+ try:
+ self.domain.set_power_indices(**kwargs)
+ except:
+ assign = np.arange(self.domain.dim(split=False),dtype=np.int)
+ else:
+ assign = self.domain.power_indices.get("pindex").flatten(order='C')
+ else:
+ assign = self.domain.enforce_shape(assign).astype(np.int).flatten(order='C')
+ ## build indexing
+ self.ind = [np.where(assign==ii)[0] for ii in xrange(np.max(assign,axis=None,out=None)+1) if ii in assign]
+
+ self.sym = True
+# about.infos.cprint("INFO: pseudo unitary projection operator.")
+ self.uni = False
+ self.imp = True
+
+ self.target = nested_space([point_space(len(self.ind),datatype=self.domain.datatype),self.domain])
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def bands(self):
+ """
+ Computes the number of projection bands
+
+ Returns
+ -------
+ bands : int
+ The number of projection bands
+ """
+ return len(self.ind)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def rho(self):
+ """
+ Computes the number of degrees of freedom per projection band.
+
+ Returns
+ -------
+ rho : ndarray
+ The number of degrees of freedom per projection band.
+ """
+ rho = np.empty(len(self.ind),dtype=np.int,order='C')
+ if(self.domain.dim(split=False)==self.domain.dof()): ## no hidden degrees of freedom
+ for ii in xrange(len(self.ind)):
+ rho[ii] = len(self.ind[ii])
+ else: ## hidden degrees of freedom
+ mof = np.round(np.real(self.domain.calc_weight(self.domain.get_meta_volume(total=False),power=-1).flatten(order='C')),decimals=0,out=None).astype(np.int) ## meta degrees of freedom
+ for ii in xrange(len(self.ind)):
+ rho[ii] = np.sum(mof[self.ind[ii]],axis=None,dtype=np.int,out=None)
+ return rho
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def _multiply(self,x,band=None,bandsup=None,**kwargs):
+ """
+ Applies the operator to a given field.
+
+ Parameters
+ ----------
+ x : field
+ Valid input field.
+ band : int, *optional*
+ Projection band whereon to project (default: None).
+ bandsup: {integer, list/array of integers}, *optional*
+ List of projection bands whereon to project and which to sum
+ up. The `band` keyword is prefered over `bandsup`
+ (default: None).
+
+ Returns
+ -------
+ Px : field
+ projected field(!)
+ """
+ if(band is not None):
+ band = int(band)
+ if(band>self.bands()-1)or(band<0):
+ raise TypeError(about._errors.cstring("ERROR: invalid band."))
+ Px = np.zeros(self.domain.dim(split=False),dtype=self.domain.datatype,order='C')
+ Px[self.ind[band]] += x.val.flatten(order='C')[self.ind[band]]
+ Px = field(self.domain,val=Px,target=x.target)
+ return Px
+
+ elif(bandsup is not None):
+ if(np.isscalar(bandsup)):
+ bandsup = np.arange(int(bandsup+1),dtype=np.int)
+ else:
+ bandsup = np.array(bandsup,dtype=np.int)
+ if(np.any(bandsup>self.bands()-1))or(np.any(bandsup<0)):
+ raise ValueError(about._errors.cstring("ERROR: invalid input."))
+ Px = np.zeros(self.domain.dim(split=False),dtype=self.domain.datatype,order='C')
+ x_ = x.val.flatten(order='C')
+ for bb in bandsup:
+ Px[self.ind[bb]] += x_[self.ind[bb]]
+ Px = field(self.domain,val=Px,target=x.target)
+ return Px
+
+ else:
+ Px = np.zeros((len(self.ind),self.domain.dim(split=False)),dtype=self.target.datatype,order='C')
+ x_ = x.val.flatten(order='C')
+ for bb in xrange(self.bands()):
+ Px[bb][self.ind[bb]] += x_[self.ind[bb]]
+ Px = field(self.target,val=Px,target=nested_space([point_space(len(self.ind),datatype=x.target.datatype),x.target]))
+ return Px
+
+ def _inverse_multiply(self,x,**kwargs):
+ raise AttributeError(about._errors.cstring("ERROR: singular operator.")) ## pseudo unitary
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def _debriefing(self,x,x_,target,inverse): ## > evaluates x and x_ after `multiply`
+ if(x_ is None):
+ return None
+ else:
+ ## weight if ...
+ if(not self.imp)and(not target.discrete)and(inverse):
+ x_ = x_.weight(power=-1,overwrite=False)
+ ## inspect x
+ if(isinstance(x,field)):
+ if(x_.domain==self.target):
+ ## repair ...
+ if(x_.domain.nest[-1]!=x.domain):
+ x_ = x_.transform(target=nested_space([point_space(len(self.ind),datatype=x.domain.datatype),x.domain]),overwrite=False) ## ... domain
+ if(x_.target.nest[-1]!=x.target):
+ x_.set_target(newtarget=nested_space([point_space(len(self.ind),datatype=x.target.datatype),x.target])) ## ... codomain
+ else:
+ ## repair ...
+ if(x_.domain!=x.domain):
+ x_ = x_.transform(target=x.domain,overwrite=False) ## ... domain
+ if(x_.target!=x.target):
+ x_.set_target(newtarget=x.target) ## ... codomain
+ return x_
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def times(self,x,**kwargs):
+ """
+ Applies the operator to a given object
+
+ Parameters
+ ----------
+ x : {scalar, ndarray, field}
+ Scalars are interpreted as constant arrays, and an array will
+ be interpreted as a field on the domain of the operator.
+
+ Returns
+ -------
+ Ox : {field, tuple of fields}
+ Mapped field on the target domain of the operator.
+
+ Other Parameters
+ ----------------
+ band : int, *optional*
+ Projection band whereon to project (default: None).
+ bandsup: {integer, list/array of integers}, *optional*
+ List of projection bands whereon to project and which to sum
+ up. The `band` keyword is prefered over `bandsup`
+ (default: None).
+ split: bool, *optional*
+ Whether to return a tuple of the projected and residual field;
+ applys only if `band` or `bandsup` is given
+ (default: False).
+
+ """
+ ## prepare
+ x_ = self._briefing(x,self.domain,False)
+ ## apply operator
+ x_ = self._multiply(x_,**kwargs)
+ ## evaluate
+ y = self._debriefing(x,x_,self.target,False)
+ ## split if ...
+ if(kwargs.get("split",False))and((kwargs.has_key("band"))or(kwargs.has_key("bandsup"))):
+ return y,x-y
+ else:
+ return y
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def pseudo_tr(self,x,**kwargs):
+ """
+ Computes the pseudo trace of a given object for all projection bands
+
+ Parameters
+ ----------
+ x : {field, operator}
+ The object whose pseudo-trace is to be computed. If the input is
+ a field, the pseudo trace equals the trace of
+ the projection operator mutliplied by a vector-vector operator
+ corresponding to the input field. This is also equal to the
+ pseudo inner product of the field with projected field itself.
+ If the input is a operator, the pseudo trace equals the trace of
+ the projection operator multiplied by the input operator.
+ 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)
+ save : bool, *optional*
+ whether all individual probing results are saved or not
+ (default: False)
+ path : string, *optional*
+ path wherein the results are saved (default: "tmp")
+ prefix : string, *optional*
+ prefix for all saved files (default: "")
+ loop : bool, *optional*
+ Indicates whether or not to perform a loop i.e., to
+ parallelise (default: False)
+
+ Returns
+ -------
+ tr : float
+ Pseudo trace for all projection bands
+ """
+ if(isinstance(x,operator)):
+ ## compute non-bare diagonal of the operator x
+ x = x.diag(bare=False,domain=self.domain,target=x.domain,var=False,**kwargs)
+ if(x is None):
+ raise TypeError(about._error.cstring("ERROR: 'NoneType' encountered."))
+
+ elif(isinstance(x,field)):
+ ## check domain
+ if(self.domain==x.domain):
+ x = x.val
+ else:
+ x = x.transform(target=self.domain,overwrite=False).val
+ ## compute non-bare diagonal of the vector-vector operator corresponding to the field x
+ x = x*np.conjugate(x)
+ ## weight
+ if(not self.domain.discrete):
+ x = self.domain.calc_weight(x,power=1)
+
+ else:
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+
+ x = np.real(x.flatten(order='C'))
+ if(not self.domain.dim(split=False)==self.domain.dof()):
+ x *= np.round(np.real(self.domain.calc_weight(self.domain.get_meta_volume(total=False),power=-1).flatten(order='C')),decimals=0,out=None).astype(np.int) ## meta degrees of freedom
+
+ tr = np.empty(self.bands(),dtype=x.dtype,order='C')
+ for bb in xrange(self.bands()):
+ tr[bb] = np.sum(x[self.ind[bb]],axis=None,dtype=None,out=None)
+ return tr
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def __repr__(self):
+ return "<nifty_core.projection_operator>"
+
+##-----------------------------------------------------------------------------
+
+
+
+##-----------------------------------------------------------------------------
+
+class vecvec_operator(operator):
+ """
+ .. __
+ .. __/ /__
+ .. __ __ _______ _______ __ __ _______ _______ /__ __/
+ .. | |/ / / __ / / ____/ | |/ / / __ / / ____/ /__/
+ .. | / / /____/ / /____ | / / /____/ / /____
+ .. |____/ \______/ \______/ |____/ \______/ \______/ operator class
+
+ NIFTY subclass for vector-vector operators
+
+ Parameters
+ ----------
+ domain : space, *optional*
+ The space wherein valid arguments live. If none is given, the
+ space of the field given in val is used. (default: None)
+ val : {scalar, ndarray, field}, *optional*
+ The field from which to construct the operator. For a scalar, a constant
+ field is defined having the value provided. If no domain
+ is given, val must be a field. (default: 1)
+
+ Attributes
+ ----------
+ domain : space
+ The space wherein valid arguments live.
+ val : ndarray
+ The field from which the operator is derived.
+ sym : bool
+ Indicates whether the operator is self-adjoint or not
+ uni : bool
+ Indicates whether the operator is unitary or not
+ imp : bool
+ Indicates whether the incorporation of volume weights in
+ multiplications is already implemented in the `multiply`
+ instance methods or not
+ target : space
+ The space wherein the operator output lives.
+ """
+ def __init__(self,domain=None,val=1):
+ """
+ Sets the standard operator properties and `values`.
+
+ Parameters
+ ----------
+ domain : space, *optional*
+ The space wherein valid arguments live. If none is given, the
+ space of the field given in val is used. (default: None)
+ val : {scalar, ndarray, field}, *optional*
+ The field from which to construct the operator. For a scalar, a constant
+ field is defined having the value provided. If no domain
+ is given, val must be a field. (default: 1)
+
+ Returns
+ -------
+ None
+ """
+ if(domain is None)and(isinstance(val,field)):
+ domain = val.domain
+ if(not isinstance(domain,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ self.domain = domain
+ self.val = self.domain.enforce_values(val,extend=True)
+ self.sym = True
+ self.uni = False
+ self.imp = False
+ self.target = self.domain
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def set_val(self,newval):
+ """
+ Sets the field values of the operator
+
+ Parameters
+ ----------
+ newval : {scalar, ndarray, field}
+ The new field values. For a scalar, a constant
+ diagonal is defined having the value provided. If no domain
+ is given, diag must be a field. (default: 1)
+
+ Returns
+ -------
+ None
+ """
+ self.val = self.domain.enforce_values(newval,extend=True)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def _multiply(self,x,**kwargs): ## > applies the operator to a given field
+ x_ = field(self.target,val=None,target=x.target)
+ x_.val = self.val*self.domain.calc_dot(self.val,x.val) ## bypasses self.domain.enforce_values
+ return x_
+
+ def _inverse_multiply(self,x,**kwargs):
+ raise AttributeError(about._errors.cstring("ERROR: singular 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(domain is None)or(domain==self.domain):
+ if(not self.domain.discrete):
+ return self.domain.calc_dot(self.val,self.domain.calc_weight(self.val,power=1))
+ else:
+ return self.domain.calc_dot(self.val,self.val)
+ else:
+ return super(vecvec_operator,self).tr(domain=domain,**kwargs) ## probing
+
+ def inverse_tr(self):
+ """
+ Inverse is ill-defined for this operator.
+ """
+ raise AttributeError(about._errors.cstring("ERROR: singular operator."))
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def diag(self,bare=False,domain=None,**kwargs):
+ """
+ Computes the diagonal of the operator.
+
+ Parameters
+ ----------
+ bare : bool, *optional*
+ Indicatese whether the diagonal entries are `bare` or not
+ (mandatory for the correct incorporation of volume weights)
+ (default: False)
+ 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).
+ save : bool, *optional*
+ whether all individual probing results are saved or not
+ (default: False)
+ path : string, *optional*
+ path wherein the results are saved (default: "tmp")
+ prefix : string, *optional*
+ prefix for all saved files (default: "")
+ loop : bool, *optional*
+ Indicates whether or not to perform a loop i.e., to
+ parallelise (default: False)
+
+ Returns
+ -------
+ diag : ndarray
+ The matrix diagonal
+ delta : float, *optional*
+ Probing variance of the trace. Returned if `var` is True in
+ of probing case.
+
+ Notes
+ -----
+ The ambiguity of `bare` or non-bare diagonal entries is based
+ on the choice of a matrix representation of the operator in
+ question. The naive choice of absorbing the volume weights
+ into the matrix leads to a matrix-vector calculus with the
+ non-bare entries which seems intuitive, though. The choice of
+ keeping matrix entries and volume weights separate deals with the
+ bare entries that allow for correct interpretation of the matrix
+ entries; e.g., as variance in case of an covariance operator.
+
+ """
+ if(domain is None)or(domain==self.domain):
+ diag = np.real(self.val*np.conjugate(self.val)) ## bare diagonal
+ ## weight if ...
+ if(not self.domain.discrete)and(not bare):
+ return self.domain.calc_weight(diag,power=1)
+ else:
+ return diag
+ else:
+ return super(vecvec_operator,self).diag(bare=bare,domain=domain,**kwargs) ## probing
+
+ def inverse_diag(self):
+ """
+ Inverse is ill-defined for this operator.
+
+ """
+ raise AttributeError(about._errors.cstring("ERROR: singular operator."))
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def det(self):
+ """
+ Computes the determinant of the operator.
+
+ Returns
+ -------
+ det : 0
+ The determinant
+
+ """
+ return 0
+
+ def inverse_det(self):
+ """
+ Inverse is ill-defined for this operator.
+
+ """
+ raise AttributeError(about._errors.cstring("ERROR: singular operator."))
+
+ def log_det(self):
+ """
+ Logarithm of the determinant is ill-defined for this singular operator.
+
+ """
+ raise AttributeError(about._errors.cstring("ERROR: singular operator."))
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def __repr__(self):
+ return "<nifty_core.vecvec_operator>"
+
+##-----------------------------------------------------------------------------
+
+
+
+##-----------------------------------------------------------------------------
+
+class response_operator(operator):
+ """
+ .. _____ _______ _______ ______ ______ __ ___ _______ _______
+ .. / __/ / __ / / _____/ / _ | / _ | / _ | / _____/ / __ /
+ .. / / / /____/ /_____ / / /_/ / / /_/ / / / / / /_____ / / /____/
+ .. /__/ \______/ /_______/ / ____/ \______/ /__/ /__/ /_______/ \______/ operator class
+ .. /__/
+
+ NIFTY subclass for response operators (of a certain family)
+
+ Any response operator handles Gaussian convolutions, itemwise masking,
+ and selective mappings.
+
+ Parameters
+ ----------
+ domain : space
+ The space wherein valid arguments live.
+ sigma : float, *optional*
+ The standard deviation of the Gaussian kernel. Zero indicates
+ no convolution. (default: 0)
+ mask : {scalar, ndarray}, *optional*
+ Masking values for arguments (default: 1)
+ assign : {list, ndarray}, *optional*
+ Assignments of codomain items to domain items. A list of
+ indices/ index tuples or a one/ two-dimensional array.
+ (default: None)
+ den : bool, *optional*
+ Whether to consider the arguments as densities or not.
+ Mandatory for the correct incorporation of volume weights.
+ (default: False)
+ target : space, *optional*
+ The space wherein the operator output lives (default: domain)
+
+ Attributes
+ ----------
+ domain : space
+ The space wherein valid arguments live.
+ sym : bool
+ Indicates whether the operator is self-adjoint or not.
+ uni : bool
+ Indicates whether the operator is unitary or not.
+ imp : bool
+ Indicates whether the incorporation of volume weights in
+ multiplications is already implemented in the `multiply`
+ instance methods or not.
+ target : space
+ The space wherein the operator output lives
+ sigma : float
+ The standard deviation of the Gaussian kernel. Zero indicates
+ no convolution.
+ mask : {scalar, ndarray}
+ Masking values for arguments
+ assign : {list, ndarray}
+ Assignments of codomain items to domain items. A list of
+ indices/ index tuples or a one/ two-dimensional array.
+ den : bool
+ Whether to consider the arguments as densities or not.
+ Mandatory for the correct incorporation of volume weights.
+ """
+ def __init__(self,domain,sigma=0,mask=1,assign=None,den=False,target=None):
+ """
+ Sets the standard properties and `density`, `sigma`, `mask` and `assignment(s)`.
+
+ Parameters
+ ----------
+ domain : space
+ The space wherein valid arguments live.
+ sigma : float, *optional*
+ The standard deviation of the Gaussian kernel. Zero indicates
+ no convolution. (default: 0)
+ mask : {scalar, ndarray}, *optional*
+ Masking values for arguments (default: 1)
+ assign : {list, ndarray}, *optional*
+ Assignments of codomain items to domain items. A list of
+ indices/ index tuples or a one/ two-dimensional array.
+ (default: None)
+ den : bool, *optional*
+ Whether to consider the arguments as densities or not.
+ Mandatory for the correct incorporation of volume weights.
+ (default: False)
+ target : space, *optional*
+ The space wherein the operator output lives (default: domain)
+
+ Returns
+ -------
+ None
+ """
+ if(not isinstance(domain,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ self.domain = domain
+
+ self.sym = False
+ self.uni = False
+ self.imp = False
+ self.den = bool(den)
+
+ self.mask = self.domain.enforce_values(mask,extend=False)
+
+ ## check sigma
+ if(sigma<0):
+ raise ValueError(about._errors.cstring("ERROR: invalid sigma."))
+ self.sigma = sigma
+
+ ## check assignment(s)
+ if(assign is None):
+ ## 1:1 assignment
+ assignments = self.domain.dim(split=False)
+ self.assign = None
+ elif(np.size(self.domain.dim(split=True))==1):
+ if(np.isscalar(assign)):
+ ## X:1 assignment
+ assignments = 1
+ if(assign[0]>=self.domain.dim(split=False))or(assign[0]<-self.domain.dim(split=False)):
+ raise IndexError(about._errors.cstring("ERROR: invalid bounds."))
+ self.assign = [int(assign)]
+ else:
+ assign = np.array(assign,dtype=np.int)
+ assignments = len(assign)
+ if(np.ndim(assign)!=1):
+ raise ValueError(about._errors.cstring("ERROR: invalid input."))
+ elif(np.any(assign>=self.domain.dim(split=False)))or(np.any(assign<-self.domain.dim(split=False))):
+ raise IndexError(about._errors.cstring("ERROR: invalid bounds."))
+ if(assignments==len(np.unique(assign,return_index=False,return_inverse=False))):
+ self.assign = assign.tolist()
+ else:
+ self.assign = assign
+ else:
+ if(np.isscalar(assign)):
+ raise ValueError(about._errors.cstring("ERROR: invalid input."))
+ else:
+ assign = np.array(assign,dtype=np.int)
+ assignments = np.size(assign,axis=0)
+ if(np.ndim(assign)!=2)or(np.size(assign,axis=1)!=np.size(self.domain.dim(split=True))):
+ raise ValueError(about._errors.cstring("ERROR: invalid input."))
+ 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])):
+ raise IndexError(about._errors.cstring("ERROR: invalid bounds."))
+ if(assignments==len(np.unique(np.ravel_multi_index(assign.T,self.domain.dim(split=True),mode="raise",order='C'),return_index=False,return_inverse=False))):
+ self.assign = assign.T.tolist()
+ else:
+ self.assign = assign
+
+ if(target is None):
+ ## set target
+ target = point_space(assignments,datatype=self.domain.datatype)
+ else:
+ ## check target
+ if(not isinstance(target,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ elif(not target.discrete):
+ raise ValueError(about._errors.cstring("ERROR: continuous codomain.")) ## discrete(!)
+ elif(np.size(target.dim(split=True))!=1):
+ raise ValueError(about._errors.cstring("ERROR: structured codomain.")) ## unstructured(!)
+ elif(assignments!=target.dim(split=False)):
+ raise ValueError(about._errors.cstring("ERROR: dimension mismatch ( "+str(assignments)+" <> "+str(target.dim(split=False))+" )."))
+ self.target = target
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def set_sigma(self,newsigma):
+ """
+ Sets the standard deviation of the response operator, indicating
+ the amount of convolution.
+
+ Parameters
+ ----------
+ sigma : float
+ The standard deviation of the Gaussian kernel. Zero indicates
+ no convolution.
+
+ Returns
+ -------
+ None
+ """
+ ## check sigma
+ if(newsigma<0):
+ raise ValueError(about._errors.cstring("ERROR: invalid sigma."))
+ self.sigma = newsigma
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def set_mask(self,newmask):
+ """
+ Sets the masking values of the response operator
+
+ Parameters
+ ----------
+ newmask : {scalar, ndarray}
+ masking values for arguments
+
+ Returns
+ -------
+ None
+ """
+ self.mask = self.domain.enforce_values(newmask,extend=False)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def _multiply(self,x,**kwargs): ## > applies the operator to a given field
+ ## smooth
+ x_ = self.domain.calc_smooth(x.val,sigma=self.sigma)
+ ## mask
+ x_ *= self.mask
+ ## assign and return
+ if(self.assign is None):
+ return field(self.target,val=x_,target=kwargs.get("target",None))
+ elif(isinstance(self.assign,list)):
+ return field(self.target,val=x_[self.assign],target=kwargs.get("target",None))
+ else:
+ 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
+ x_ = np.zeros(self.domain.dim(split=True),dtype=self.domain.datatype,order='C')
+ ## assign (transposed)
+ if(self.assign is None):
+ x_ = np.copy(x.val.flatten(order='C'))
+ elif(isinstance(self.assign,list)):
+ x_[self.assign] += x.val.flatten(order='C')
+ else:
+ for ii in xrange(np.size(self.assign,axis=0)):
+ x_[np.array([self.assign[ii]]).T.tolist()] += x[ii]
+ ## mask
+ x_ *= self.mask
+ ## smooth
+ x_ = self.domain.calc_smooth(x_,sigma=self.sigma)
+ #return x_
+ return field(self.domain,val=x_,target=kwargs.get("target",None))
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def _briefing(self,x,domain,inverse): ## > prepares x for `multiply`
+ ## inspect x
+ if(not isinstance(x,field)):
+ x_ = field(domain,val=x,target=None)
+ else:
+ ## check x.domain
+ if(domain==x.domain):
+ x_ = x
+ ## transform
+ else:
+ x_ = x.transform(target=domain,overwrite=False)
+ ## weight if ...
+ if(not self.imp)and(not domain.discrete)and(not inverse)and(self.den): ## respect density
+ x_ = x_.weight(power=1,overwrite=False)
+ return x_
+
+ def _debriefing(self,x,x_,target,inverse): ## > evaluates x and x_ after `multiply`
+ if(x_ is None):
+ return None
+ else:
+ ## inspect x_
+ if(not isinstance(x_,field)):
+ x_ = field(target,val=x_,target=None)
+ elif(x_.domain!=target):
+ raise ValueError(about._errors.cstring("ERROR: invalid output domain."))
+ ## weight if ...
+ if(not self.imp)and(not target.discrete)and(not self.den^inverse): ## respect density
+ x_ = x_.weight(power=-1,overwrite=False)
+ ## inspect x
+ if(isinstance(x,field)):
+ ## repair ...
+ if(self.domain==self.target!=x.domain):
+ x_ = x_.transform(target=x.domain,overwrite=False) ## ... domain
+ if(x_.domain==x.domain)and(x_.target!=x.target):
+ x_.set_target(newtarget=x.target) ## ... codomain
+ return x_
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def __repr__(self):
+ return "<nifty_core.response_operator>"
+
+##-----------------------------------------------------------------------------
+
+
+
+class invertible_operator(operator):
+ """
+ .. __ __ __ __ __
+ .. /__/ / /_ /__/ / / / /
+ .. __ __ ___ __ __ _______ _____ / _/ __ / /___ / / _______
+ .. / / / _ || |/ / / __ / / __/ / / / / / _ | / / / __ /
+ .. / / / / / / | / / /____/ / / / /_ / / / /_/ / / /_ / /____/
+ .. /__/ /__/ /__/ |____/ \______/ /__/ \___/ /__/ \______/ \___/ \______/ operator class
+
+ NIFTY subclass for invertible, self-adjoint (linear) operators
+
+ The invertible operator class is an abstract class for self-adjoint or
+ symmetric (linear) operators from which other more specific operator
+ subclassescan be derived. Such operators inherit an automated inversion
+ routine, namely conjugate gradient.
+
+ Parameters
+ ----------
+ domain : space
+ The space wherein valid arguments live.
+ uni : bool, *optional*
+ Indicates whether the operator is unitary or not.
+ (default: False)
+ imp : bool, *optional*
+ Indicates whether the incorporation of volume weights in
+ multiplications is already implemented in the `multiply`
+ instance methods or not (default: False).
+ para : {single object, tuple/list of objects}, *optional*
+ This is a freeform tuple/list of parameters that derivatives of
+ the operator class can use (default: None).
+
+ See Also
+ --------
+ operator
+
+ Notes
+ -----
+ This class is not meant to be instantiated. Operator classes derived
+ from this one only need a `_multiply` or `_inverse_multiply` instance
+ method to perform the other. However, one of them needs to be defined.
+
+ Attributes
+ ----------
+ domain : space
+ The space wherein valid arguments live.
+ sym : bool
+ Indicates whether the operator is self-adjoint or not.
+ uni : bool
+ Indicates whether the operator is unitary or not.
+ imp : bool
+ Indicates whether the incorporation of volume weights in
+ multiplications is already implemented in the `multiply`
+ instance methods or not.
+ target : space
+ The space wherein the operator output lives.
+ para : {single object, list of objects}
+ This is a freeform tuple/list of parameters that derivatives of
+ the operator class can use. Not used in the base operators.
+
+ """
+ def __init__(self,domain,uni=False,imp=False,para=None):
+ """
+ Sets the standard operator properties.
+
+ Parameters
+ ----------
+ domain : space
+ The space wherein valid arguments live.
+ uni : bool, *optional*
+ Indicates whether the operator is unitary or not.
+ (default: False)
+ imp : bool, *optional*
+ Indicates whether the incorporation of volume weights in
+ multiplications is already implemented in the `multiply`
+ instance methods or not (default: False).
+ para : {single object, tuple/list of objects}, *optional*
+ This is a freeform tuple/list of parameters that derivatives of
+ the operator class can use (default: None).
+
+ """
+ if(not isinstance(domain,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ self.domain = domain
+ self.sym = True
+ self.uni = bool(uni)
+
+ if(self.domain.discrete):
+ self.imp = True
+ else:
+ self.imp = bool(imp)
+
+ self.target = self.domain
+
+ if(para is not None):
+ self.para = para
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def _multiply(self,x,force=False,W=None,spam=None,reset=None,note=False,x0=None,tol=1E-4,clevel=1,limii=None,**kwargs):
+ """
+ Applies the invertible operator to a given field by invoking a
+ conjugate gradient.
+
+ Parameters
+ ----------
+ x : field
+ Valid input field.
+ force : bool
+ Indicates wheter to return a field instead of ``None`` in case
+ the conjugate gradient fails.
+
+ Returns
+ -------
+ OIIx : field
+ Mapped field with suitable domain.
+
+ See Also
+ --------
+ conjugate_gradient
+
+ Other Parameters
+ ----------------
+ W : {operator, function}, *optional*
+ Operator `W` that is a preconditioner on `A` and is applicable to a
+ field (default: None).
+ spam : function, *optional*
+ Callback function which is given the current `x` and iteration
+ counter each iteration (default: None).
+ reset : integer, *optional*
+ Number of iterations after which to restart; i.e., forget previous
+ conjugated directions (default: sqrt(b.dim())).
+ note : bool, *optional*
+ Indicates whether notes are printed or not (default: False).
+ x0 : field, *optional*
+ Starting guess for the minimization.
+ tol : scalar, *optional*
+ Tolerance specifying convergence; measured by current relative
+ residual (default: 1E-4).
+ clevel : integer, *optional*
+ Number of times the tolerance should be undershot before
+ exiting (default: 1).
+ limii : integer, *optional*
+ Maximum number of iterations performed (default: 10 * b.dim()).
+
+ """
+ x_,convergence = conjugate_gradient(self.inverse_times,x,W=W,spam=spam,reset=reset,note=note)(x0=x0,tol=tol,clevel=clevel,limii=limii)
+ ## check convergence
+ if(not convergence):
+ if(not force)or(x_ is None):
+ return None
+ about.warnings.cprint("WARNING: conjugate gradient failed.")
+ ## weight if ...
+ if(not self.imp): ## continiuos domain/target
+ x_.weight(power=-1,overwrite=True)
+ return x_
+
+ def _inverse_multiply(self,x,force=False,W=None,spam=None,reset=None,note=False,x0=None,tol=1E-4,clevel=1,limii=None,**kwargs):
+ """
+ Applies the inverse of the invertible operator to a given field by
+ invoking a conjugate gradient.
+
+ Parameters
+ ----------
+ x : field
+ Valid input field.
+ force : bool
+ Indicates wheter to return a field instead of ``None`` in case
+ the conjugate gradient fails.
+
+ Returns
+ -------
+ OIx : field
+ Mapped field with suitable domain.
+
+ See Also
+ --------
+ conjugate_gradient
+
+ Other Parameters
+ ----------------
+ W : {operator, function}, *optional*
+ Operator `W` that is a preconditioner on `A` and is applicable to a
+ field (default: None).
+ spam : function, *optional*
+ Callback function which is given the current `x` and iteration
+ counter each iteration (default: None).
+ reset : integer, *optional*
+ Number of iterations after which to restart; i.e., forget previous
+ conjugated directions (default: sqrt(b.dim())).
+ note : bool, *optional*
+ Indicates whether notes are printed or not (default: False).
+ x0 : field, *optional*
+ Starting guess for the minimization.
+ tol : scalar, *optional*
+ Tolerance specifying convergence; measured by current relative
+ residual (default: 1E-4).
+ clevel : integer, *optional*
+ Number of times the tolerance should be undershot before
+ exiting (default: 1).
+ limii : integer, *optional*
+ Maximum number of iterations performed (default: 10 * b.dim()).
+
+ """
+ x_,convergence = conjugate_gradient(self.times,x,W=W,spam=spam,reset=reset,note=note)(x0=x0,tol=tol,clevel=clevel,limii=limii)
+ ## check convergence
+ if(not convergence):
+ if(not force)or(x_ is None):
+ return None
+ about.warnings.cprint("WARNING: conjugate gradient failed.")
+ ## weight if ...
+ if(not self.imp): ## continiuos domain/target
+ x_.weight(power=1,overwrite=True)
+ return x_
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def __repr__(self):
+ return "<nifty_tools.invertible_operator>"
+
+##-----------------------------------------------------------------------------
+
+##-----------------------------------------------------------------------------
+
+class propagator_operator(operator):
+ """
+ .. __
+ .. / /_
+ .. _______ _____ ______ ______ ____ __ ____ __ ____ __ / _/ ______ _____
+ .. / _ / / __/ / _ | / _ | / _ / / _ / / _ / / / / _ | / __/
+ .. / /_/ / / / / /_/ / / /_/ / / /_/ / / /_/ / / /_/ / / /_ / /_/ / / /
+ .. / ____/ /__/ \______/ / ____/ \______| \___ / \______| \___/ \______/ /__/ operator class
+ .. /__/ /__/ /______/
+
+ NIFTY subclass for propagator operators (of a certain family)
+
+ The propagator operators :math:`D` implemented here have an inverse
+ formulation like :math:`(S^{-1} + M)`, :math:`(S^{-1} + N^{-1})`, or
+ :math:`(S^{-1} + R^\dagger N^{-1} R)` as appearing in Wiener filter
+ theory.
+
+ Parameters
+ ----------
+ S : operator
+ Covariance of the signal prior.
+ M : operator
+ Likelihood contribution.
+ R : operator
+ Response operator translating signal to (noiseless) data.
+ N : operator
+ Covariance of the noise prior or the likelihood, respectively.
+
+ See Also
+ --------
+ conjugate_gradient
+
+ Notes
+ -----
+ The propagator will puzzle the operators `S` and `M` or `R`, `N` or
+ only `N` together in the predefined from, a domain is set
+ automatically. The application of the inverse is done by invoking a
+ conjugate gradient.
+ Note that changes to `S`, `M`, `R` or `N` auto-update the propagator.
+
+ Examples
+ --------
+ >>> f = field(rg_space(4), val=[2, 4, 6, 8])
+ >>> S = power_operator(f.target, spec=1)
+ >>> N = diagonal_operator(f.domain, diag=1)
+ >>> D = propagator_operator(S=S, N=N) # D^{-1} = S^{-1} + N^{-1}
+ >>> D(f).val
+ array([ 1., 2., 3., 4.])
+
+ Attributes
+ ----------
+ domain : space
+ A space wherein valid arguments live.
+ codomain : space
+ An alternative space wherein valid arguments live; commonly the
+ codomain of the `domain` attribute.
+ sym : bool
+ Indicates that the operator is self-adjoint.
+ uni : bool
+ Indicates that the operator is not unitary.
+ imp : bool
+ Indicates that volume weights are implemented in the `multiply`
+ instance methods.
+ target : space
+ The space wherein the operator output lives.
+ _A1 : {operator, function}
+ Application of :math:`S^{-1}` to a field.
+ _A2 : {operator, function}
+ Application of all operations not included in `A1` to a field.
+ RN : {2-tuple of operators}, *optional*
+ Contains `R` and `N` if given.
+
+ """
+ def __init__(self,S=None,M=None,R=None,N=None):
+ """
+ Sets the standard operator properties and `codomain`, `_A1`, `_A2`,
+ and `RN` if required.
+
+ Parameters
+ ----------
+ S : operator
+ Covariance of the signal prior.
+ M : operator
+ Likelihood contribution.
+ R : operator
+ Response operator translating signal to (noiseless) data.
+ N : operator
+ Covariance of the noise prior or the likelihood, respectively.
+
+ """
+ ## check signal prior covariance
+ if(S is None):
+ raise Exception(about._errors.cstring("ERROR: insufficient input."))
+ elif(not isinstance(S,operator)):
+ raise ValueError(about._errors.cstring("ERROR: invalid input."))
+ space1 = S.domain
+
+ ## check likelihood (pseudo) covariance
+ if(M is None):
+ if(N is None):
+ raise Exception(about._errors.cstring("ERROR: insufficient input."))
+ elif(not isinstance(N,operator)):
+ raise ValueError(about._errors.cstring("ERROR: invalid input."))
+ if(R is None):
+ space2 = N.domain
+ elif(not isinstance(R,operator)):
+ raise ValueError(about._errors.cstring("ERROR: invalid input."))
+ else:
+ space2 = R.domain
+ elif(not isinstance(M,operator)):
+ raise ValueError(about._errors.cstring("ERROR: invalid input."))
+ else:
+ space2 = M.domain
+
+ ## set spaces
+ self.domain = space2
+ if(self.domain.check_codomain(space1)):
+ self.codomain = space1
+ else:
+ raise ValueError(about._errors.cstring("ERROR: invalid input."))
+ self.target = self.domain
+
+ ## define A1 == S_inverse
+ if(isinstance(S,diagonal_operator)):
+ self._A1 = S._inverse_multiply ## S.imp == True
+ else:
+ self._A1 = S.inverse_times
+
+ ## define A2 == M == R_adjoint N_inverse R == N_inverse
+ if(M is None):
+ if(R is not None):
+ self.RN = (R,N)
+ if(isinstance(N,diagonal_operator)):
+ self._A2 = self._standard_M_times_1
+ else:
+ self._A2 = self._standard_M_times_2
+ elif(isinstance(N,diagonal_operator)):
+ self._A2 = N._inverse_multiply ## N.imp == True
+ else:
+ self._A2 = N.inverse_times
+ elif(isinstance(M,diagonal_operator)):
+ self._A2 = M._multiply ## M.imp == True
+ else:
+ self._A2 = M.times
+
+ self.sym = True
+ self.uni = False
+ self.imp = True
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def _standard_M_times_1(self,x,**kwargs): ## applies > R_adjoint N_inverse R assuming N is diagonal
+ return self.RN[0].adjoint_times(self.RN[1]._inverse_multiply(self.RN[0].times(x))) ## N.imp = True
+
+ def _standard_M_times_2(self,x,**kwargs): ## applies > R_adjoint N_inverse R
+ return self.RN[0].adjoint_times(self.RN[1].inverse_times(self.RN[0].times(x)))
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def _inverse_multiply_1(self,x,**kwargs): ## > applies A1 + A2 in self.codomain
+ return self._A1(x,pseudo=True)+self._A2(x.transform(self.domain)).transform(self.codomain)
+
+ def _inverse_multiply_2(self,x,**kwargs): ## > applies A1 + A2 in self.domain
+ return self._A1(x.transform(self.codomain),pseudo=True).transform(self.domain)+self._A2(x)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def _briefing(self,x): ## > prepares x for `multiply`
+ ## inspect x
+ if(not isinstance(x,field)):
+ return field(self.domain,val=x,target=None),False
+ ## check x.domain
+ elif(x.domain==self.domain):
+ return x,False
+ elif(x.domain==self.codomain):
+ return x,True
+ ## transform
+ else:
+ return x.transform(target=self.codomain,overwrite=False),True
+
+ def _debriefing(self,x,x_,in_codomain): ## > evaluates x and x_ after `multiply`
+ if(x_ is None):
+ return None
+ ## inspect x
+ elif(isinstance(x,field)):
+ ## repair ...
+ if(in_codomain)and(x.domain!=self.codomain):
+ x_ = x_.transform(target=x.domain,overwrite=False) ## ... domain
+ if(x_.target!=x.target):
+ x_.set_target(newtarget=x.target) ## ... codomain
+ return x_
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def times(self,x,force=False,W=None,spam=None,reset=None,note=False,x0=None,tol=1E-4,clevel=1,limii=None,**kwargs):
+ """
+ Applies the propagator to a given object by invoking a
+ conjugate gradient.
+
+ Parameters
+ ----------
+ x : {scalar, list, array, field}
+ Scalars are interpreted as constant arrays, and an array will
+ be interpreted as a field on the domain of the operator.
+ force : bool
+ Indicates wheter to return a field instead of ``None`` in case
+ the conjugate gradient fails.
+
+ Returns
+ -------
+ Dx : field
+ Mapped field with suitable domain.
+
+ See Also
+ --------
+ conjugate_gradient
+
+ Other Parameters
+ ----------------
+ W : {operator, function}, *optional*
+ Operator `W` that is a preconditioner on `A` and is applicable to a
+ field (default: None).
+ spam : function, *optional*
+ Callback function which is given the current `x` and iteration
+ counter each iteration (default: None).
+ reset : integer, *optional*
+ Number of iterations after which to restart; i.e., forget previous
+ conjugated directions (default: sqrt(b.dim())).
+ note : bool, *optional*
+ Indicates whether notes are printed or not (default: False).
+ x0 : field, *optional*
+ Starting guess for the minimization.
+ tol : scalar, *optional*
+ Tolerance specifying convergence; measured by current relative
+ residual (default: 1E-4).
+ clevel : integer, *optional*
+ Number of times the tolerance should be undershot before
+ exiting (default: 1).
+ limii : integer, *optional*
+ Maximum number of iterations performed (default: 10 * b.dim()).
+
+ """
+ ## prepare
+ x_,in_codomain = self._briefing(x)
+ ## apply operator
+ if(in_codomain):
+ A = self._inverse_multiply_1
+ else:
+ A = self._inverse_multiply_2
+ x_,convergence = conjugate_gradient(A,x_,W=W,spam=spam,reset=reset,note=note)(x0=x0,tol=tol,clevel=clevel,limii=limii)
+ ## evaluate
+ if(not convergence):
+ if(not force):
+ return None
+ about.warnings.cprint("WARNING: conjugate gradient failed.")
+ return self._debriefing(x,x_,in_codomain)
+
+ def inverse_times(self,x,**kwargs):
+ """
+ Applies the inverse propagator to a given object.
+
+ Parameters
+ ----------
+ x : {scalar, list, array, field}
+ Scalars are interpreted as constant arrays, and an array will
+ be interpreted as a field on the domain of the operator.
+
+ Returns
+ -------
+ DIx : field
+ Mapped field with suitable domain.
+
+ """
+ ## prepare
+ x_,in_codomain = self._briefing(x)
+ ## apply operator
+ if(in_codomain):
+ x_ = self._inverse_multiply_1(x_)
+ else:
+ x_ = self._inverse_multiply_2(x_)
+ ## evaluate
+ return self._debriefing(x,x_,in_codomain)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def __repr__(self):
+ return "<nifty_tools.propagator_operator>"
+
+##-----------------------------------------------------------------------------
diff --git a/operators/nifty_probing.py b/operators/nifty_probing.py
new file mode 100644
index 0000000000000000000000000000000000000000..ed8875365ac4a308603cdda4dd6d69a1c3c9fcc0
--- /dev/null
+++ b/operators/nifty_probing.py
@@ -0,0 +1,1330 @@
+## NIFTY (Numerical Information Field Theory) has been developed at the
+## Max-Planck-Institute for Astrophysics.
+##
+## Copyright (C) 2015 Max-Planck-Society
+##
+## Author: Marco Selig
+## Project homepage: <http://www.mpa-garching.mpg.de/ift/nifty/>
+##
+## This program is free software: you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation, either version 3 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+## See the GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+from __future__ import division
+
+import os
+from sys import stdout as so
+import numpy as np
+from multiprocessing import Pool as mp
+from multiprocessing import Value as mv
+from multiprocessing import Array as ma
+
+from nifty.nifty_about import about
+from nifty.nifty_core import space, \
+ field
+
+
+
+##-----------------------------------------------------------------------------
+
+class _share(object):
+
+ __init__ = None
+
+ @staticmethod
+ def _init_share(_sum,_num,_var):
+ _share.sum = _sum
+ _share.num = _num
+ _share.var = _var
+
+##-----------------------------------------------------------------------------
+
+##=============================================================================
+
+class probing(object):
+ """
+ .. __ __
+ .. / / /__/
+ .. ______ _____ ______ / /___ __ __ ___ ____ __
+ .. / _ | / __/ / _ | / _ | / / / _ | / _ /
+ .. / /_/ / / / / /_/ / / /_/ / / / / / / / / /_/ /
+ .. / ____/ /__/ \______/ \______/ /__/ /__/ /__/ \____ / class
+ .. /__/ /______/
+
+ NIFTY class for probing (using multiprocessing)
+
+ This is the base NIFTY probing class from which other probing classes
+ (e.g. diagonal probing) are derived.
+
+ When called, a probing class instance evaluates an operator or a
+ function using random fields, whose components are random variables
+ with mean 0 and variance 1. When an instance is called it returns the
+ mean value of f(probe), where probe is a random field with mean 0 and
+ variance 1. The mean is calculated as 1/N Sum[ f(probe_i) ].
+
+ Parameters
+ ----------
+ op : operator
+ The operator specified by `op` is the operator to be probed.
+ If no operator is given, then probing will be done by applying
+ `function` to the probes. (default: None)
+ function : function, *optional*
+ If no operator has been specified as `op`, then specification of
+ `function` is non optional. This is the function, that is applied
+ to the probes. (default: `op.times`)
+ domain : space, *optional*
+ If no operator has been specified as `op`, then specification of
+ `domain` is non optional. This is the space that the probes live
+ in. (default: `op.domain`)
+ target : domain, *optional*
+ `target` is the codomain of `domain`
+ (default: `op.domain.get_codomain()`)
+ random : string, *optional*
+ the distribution from which the probes are drawn. `random` can be
+ either "pm1" or "gau". "pm1" is a uniform distribution over {+1,-1}
+ or {+1,+i,-1,-i}, respectively. "gau" is a normal distribution with
+ zero-mean and unit-variance (default: "pm1")
+ ncpu : int, *optional*
+ the number of cpus to be used from parallel probing. (default: 2)
+ nrun : int, *optional*
+ the number of probes to be evaluated. If `nrun<ncpu`, it will be
+ set to `ncpu`. (default: 8)
+ nper : int, *optional*
+ this number specifies how many probes will be evaluated by one
+ worker. Afterwards a new worker will be created to evaluate a chunk
+ of `nper` probes.
+ If for example `nper=nrun/ncpu`, then every worker will be created
+ for every cpu. This can lead to the case, that all workers but one
+ are already finished, but have to wait for the last worker that
+ might still have a considerable amount of evaluations left. This is
+ obviously not very effective.
+ If on the other hand `nper=1`, then for each evaluation a worker will
+ be created. In this case all cpus will work until nrun probes have
+ been evaluated.
+ It is recommended to leave `nper` as the default value. (default: 1)
+ var : bool, *optional*
+ If `var` is True, then the variance of the sampled function will
+ also be returned. The result is then a tuple with the mean in the
+ zeroth entry and the variance in the first entry. (default: False)
+
+
+ See Also
+ --------
+ diagonal_probing : A probing class to get the diagonal of an operator
+ trace_probing : A probing class to get the trace of an operator
+
+
+ Attributes
+ ----------
+ function : function
+ the function, that is applied to the probes
+ domain : space
+ the space, where the probes live in
+ target : space
+ the codomain of `domain`
+ random : string
+ the random number generator used to create the probes
+ (either "pm1" or "gau")
+ ncpu : int
+ the number of cpus used for probing
+ nrun : int
+ the number of probes to be evaluated, when the instance is called
+ nper : int
+ number of probes, that will be evaluated by one worker
+ var : bool
+ whether the variance will be additionally returned, when the
+ instance is called
+ quargs : dict
+ Keyword arguments passed to `function` in each call.
+
+ """
+ def __init__(self,op=None,function=None,domain=None,target=None,random="pm1",ncpu=2,nrun=8,nper=1,var=False,**quargs):
+ """
+ initializes a probing instance
+
+ Parameters
+ ----------
+ op : operator
+ The operator specified by `op` is the operator to be probed.
+ If no operator is given, then probing will be done by applying
+ `function` to the probes. (default: None)
+ function : function, *optional*
+ If no operator has been specified as `op`, then specification of
+ `function` is non optional. This is the function, that is applied
+ to the probes. (default: `op.times`)
+ domain : space, *optional*
+ If no operator has been specified as `op`, then specification of
+ `domain` is non optional. This is the space that the probes live
+ in. (default: `op.domain`)
+ target : domain, *optional*
+ `target` is the codomain of `domain`
+ (default: `op.domain.get_codomain()`)
+ random : string, *optional*
+ the distribution from which the probes are drawn. `random` can be
+ either "pm1" or "gau". "pm1" is a uniform distribution over {+1,-1}
+ or {+1,+i,-1,-i}, respectively. "gau" is a normal distribution with
+ zero-mean and unit-variance (default: "pm1")
+ ncpu : int, *optional*
+ the number of cpus to be used from parallel probing. (default: 2)
+ nrun : int, *optional*
+ the number of probes to be evaluated. If `nrun<ncpu`, it will be
+ set to `ncpu`. (default: 8)
+ nper : int, *optional*
+ this number specifies how many probes will be evaluated by one
+ worker. Afterwards a new worker will be created to evaluate a chunk
+ of `nper` probes.
+ If for example `nper=nrun/ncpu`, then every worker will be created
+ for every cpu. This can lead to the case, that all workers but one
+ are already finished, but have to wait for the last worker that
+ might still have a considerable amount of evaluations left. This is
+ obviously not very effective.
+ If on the other hand `nper=1`, then for each evaluation a worker will
+ be created. In this case all cpus will work until nrun probes have
+ been evaluated.
+ It is recommended to leave `nper` as the default value. (default: 1)
+ var : bool, *optional*
+ If `var` is True, then the variance of the sampled function will
+ also be returned. The result is then a tuple with the mean in the
+ zeroth entry and the variance in the first entry. (default: False)
+
+ """
+ if(op is None):
+ ## check whether callable
+ if(function is None)or(not hasattr(function,"__call__")):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ ## check given domain
+ if(domain is None)or(not isinstance(domain,space)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ else:
+ from nifty_operators import operator
+ if(not isinstance(op,operator)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ ## check whether callable
+ if(function is None)or(not hasattr(function,"__call__")):
+ function = op.times
+ elif(op==function):
+ function = op.times
+ ## check whether correctly bound
+ if(op!=function.im_self):
+ raise NameError(about._errors.cstring("ERROR: invalid input."))
+ ## check given shape and domain
+ if(domain is None)or(not isinstance(domain,space)):
+ if(function in [op.inverse_times,op.adjoint_times]):
+ domain = op.target
+ else:
+ domain = op.domain
+ else:
+ if(function in [op.inverse_times,op.adjoint_times]):
+ op.target.check_codomain(domain) ## a bit pointless
+ if(target is None)or(not isinstance(target,space)):
+ target = op.target
+ else:
+ op.domain.check_codomain(domain) ## a bit pointless
+ if(target is None)or(not isinstance(target,space)):
+ target = op.domain
+
+ self.function = function
+ self.domain = domain
+
+ ## check codomain
+ if(target is None):
+ target = self.domain.get_codomain()
+ else:
+ self.domain.check_codomain(target) ## a bit pointless
+ self.target = target
+
+ if(random not in ["pm1","gau"]):
+ raise ValueError(about._errors.cstring("ERROR: unsupported random key '"+str(random)+"'."))
+ self.random = random
+
+ self.ncpu = int(max(1,ncpu))
+ self.nrun = int(max(self.ncpu,nrun))
+ if(nper is None):
+ self.nper = None
+ else:
+ self.nper = int(max(1,min(self.nrun//self.ncpu,nper)))
+
+ self.var = bool(var)
+
+ self.quargs = quargs
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def configure(self,**kwargs):
+ """
+ changes the attributes of the instance
+
+ Parameters
+ ----------
+ random : string, *optional*
+ the random number generator used to create the probes (default: "pm1")
+ ncpu : int, *optional*
+ the number of cpus to be used for parallel probing. (default: 2)
+ nrun : int, *optional*
+ the number of probes to be evaluated. If `nrun<ncpu`, it will be
+ set to `ncpu`. (default: 8)
+ nper : int, *optional*
+ number of probes, that will be evaluated by one worker (default: 1)
+ var : bool, *optional*
+ whether the variance will be additionally returned (default: False)
+
+ """
+ if("random" in kwargs):
+ if(kwargs.get("random") not in ["pm1","gau"]):
+ raise ValueError(about._errors.cstring("ERROR: unsupported random key '"+str(kwargs.get("random"))+"'."))
+ self.random = kwargs.get("random")
+
+ if("ncpu" in kwargs):
+ self.ncpu = int(max(1,kwargs.get("ncpu")))
+ if("nrun" in kwargs):
+ self.nrun = int(max(self.ncpu,kwargs.get("nrun")))
+ if("nper" in kwargs):
+ if(kwargs.get("nper") is None):
+ self.nper = None
+ else:
+ self.nper = int(max(1,min(self.nrun//self.ncpu,kwargs.get("nper"))))
+
+ if("var" in kwargs):
+ self.var = bool(kwargs.get("var"))
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def gen_probe(self):
+ """
+ Generates a single probe
+
+ Returns
+ -------
+ probe : field
+ a random field living in `domain` with mean 0 and variance 1 in
+ each component
+
+ """
+ return field(self.domain,val=None,target=self.target,random=self.random)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def probing(self,idnum,probe):
+ """
+ Computes a single probing result given one probe
+
+ Parameters
+ ----------
+ probe : field
+ the field on which `function` will be applied
+ idnum : int
+ the identification number of the probing
+
+ Returns
+ -------
+ result : array-like
+ the result of applying `function` to `probe`. The exact type
+ depends on the function.
+
+ """
+ f = self.function(probe,**self.quargs)
+ if(isinstance(f,field)):
+ return f.val
+ else:
+ return f
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def evaluate(self,summa,num,var):
+ """
+ Evaluates the probing results.
+
+ Parameters
+ ----------
+ summa : numpy.array
+ Sum of all probing results.
+ num : int
+ Number of successful probings (not returning ``None``).
+ var : numpy.array
+ Sum of all squared probing results
+
+ Returns
+ -------
+ final : numpy.array
+ The final probing result; 1/N Sum[ probing(probe_i) ]
+ var : array-like
+ The variance of the final probing result;
+ 1/(N(N-1)) Sum[( probing(probe_i) - final )^2];
+ if the variance is returned, the return will be a tuple of
+ (`final`,`var`).
+
+ """
+ if(num<self.nrun):
+ about.infos.cflush(" ( %u probe(s) failed, effectiveness == %.1f%% )\n"%(self.nrun-num,100*num/self.nrun))
+ if(num==0):
+ about.warnings.cprint("WARNING: probing failed.")
+ return None
+ else:
+ about.infos.cflush("\n")
+
+ if(summa.size==1):
+ summa = summa.flatten(order='C')[0]
+ if(self.var):
+ var = var.flatten(order='C')[0]
+ if(np.iscomplexobj(summa))and(np.all(np.imag(summa)==0)):
+ summa = np.real(summa)
+
+ final = summa*(1/num)
+ if(self.var):
+ if(num==1):
+ about.warnings.cprint("WARNING: infinite variance.")
+ return final,None
+ else:
+ var = var*(1/(num*(num-1)))-np.real(np.conjugate(final)*final)*(1/(num-1))
+ return final,var
+ else:
+ return final
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def _progress(self,irun): ## > prints progress status by in upto 10 dots
+ tenths = 1+(10*irun//self.nrun)
+ about.infos.cflush(("\b")*10+('.')*tenths+(' ')*(10-tenths))
+
+ def _single_probing(self,zipped): ## > performs one probing operation
+ ## generate probe
+ np.random.seed(zipped[0])
+ probe = self.gen_probe()
+ ## do the actual probing
+ return self.probing(zipped[1],probe)
+
+ def _serial_probing(self,zipped): ## > performs the probing operation serially
+ try:
+ result = self._single_probing(zipped)
+ except Exception as exception:
+ raise exception
+ except BaseException: ## capture system-exiting exception including KeyboardInterrupt
+ raise Exception(about._errors.cstring("ERROR: unknown."))
+ else:
+ if(result is not None):
+ result = np.array(result).flatten(order='C')
+ if(np.iscomplexobj(result)):
+ _share.sum[0].acquire(block=True,timeout=None)
+ _share.sum[0][:] += np.real(result)
+ _share.sum[0].release()
+ _share.sum[1].acquire(block=True,timeout=None)
+ _share.sum[1][:] += np.imag(result)
+ _share.sum[1].release()
+ else:
+ _share.sum.acquire(block=True,timeout=None)
+ _share.sum[:] += result
+ _share.sum.release()
+ _share.num.acquire(block=True,timeout=None)
+ _share.num.value += 1
+ _share.num.release()
+ if(self.var):
+ _share.var.acquire(block=True,timeout=None)
+ _share.var[:] += np.real(np.conjugate(result)*result)
+ _share.var.release()
+ self._progress(_share.num.value)
+
+ def _parallel_probing(self): ## > performs the probing operations in parallel
+ ## define random seed
+ seed = np.random.randint(10**8,high=None,size=self.nrun)
+ ## get output shape
+ result = self.probing(0,field(self.domain,val=None,target=self.target))
+ if(np.isscalar(result))or(np.size(result)==1):
+ shape = np.ones(1,dtype=np.int,order='C')
+ else:
+ shape = np.array(np.array(result).shape)
+ ## define shared objects
+ if(np.iscomplexobj(result)):
+ _sum = (ma('d',np.zeros(np.prod(shape,axis=0,dtype=np.int,out=None),dtype=np.float64,order='C'),lock=True),ma('d',np.zeros(np.prod(shape,axis=0,dtype=np.int,out=None),dtype=np.float64,order='C'),lock=True)) ## tuple(real,imag)
+ else:
+ _sum = ma('d',np.zeros(np.prod(shape,axis=0,dtype=np.int,out=None),dtype=np.float64,order='C'),lock=True)
+ _num = mv('I',0,lock=True)
+ if(self.var):
+ _var = ma('d',np.zeros(np.prod(shape,axis=0,dtype=np.int,out=None),dtype=np.float64,order='C'),lock=True)
+ else:
+ _var = None
+ ## build pool
+ if(about.infos.status):
+ so.write(about.infos.cstring("INFO: multiprocessing "+(' ')*10))
+ so.flush()
+ pool = mp(processes=self.ncpu,initializer=_share._init_share,initargs=(_sum,_num,_var),maxtasksperchild=self.nper)
+ try:
+ ## pooling
+ pool.map_async(self._serial_probing,zip(seed,np.arange(self.nrun,dtype=np.int)),chunksize=None,callback=None).get(timeout=None)
+ ## close and join pool
+ about.infos.cflush(" done.")
+ pool.close()
+ pool.join()
+ except BaseException as exception:
+ ## terminate and join pool
+ about._errors.cprint("\nERROR: terminating pool.")
+ pool.terminate()
+ pool.join()
+ ## re-raise exception
+ raise exception ## traceback by looping
+ ## evaluate
+ if(np.iscomplexobj(result)):
+ _sum = (np.array(_sum[0][:])+np.array(_sum[1][:])*1j).reshape(shape) ## comlpex array
+ else:
+ _sum = np.array(_sum[:]).reshape(shape)
+ if(self.var):
+ _var = np.array(_var[:]).reshape(shape)
+ return self.evaluate(_sum,_num.value,_var)
+
+ def _nonparallel_probing(self): ## > performs the probing operations one after another
+ ## looping
+ if(about.infos.status):
+ so.write(about.infos.cstring("INFO: looping "+(' ')*10))
+ so.flush()
+ _sum = 0
+ _num = 0
+ _var = 0
+ for ii in xrange(self.nrun):
+ result = self._single_probing((np.random.randint(10**8,high=None,size=None),ii)) ## tuple(seed,idnum)
+ if(result is not None):
+ _sum += result ## result: {scalar, np.array}
+ _num += 1
+ if(self.var):
+ _var += np.real(np.conjugate(result)*result)
+ self._progress(_num)
+ about.infos.cflush(" done.")
+ ## evaluate
+ return self.evaluate(_sum,_num,_var)
+
+ def __call__(self,loop=False,**kwargs):
+ """
+
+ Starts the probing process.
+ All keyword arguments that can be given to `configure` can also be
+ given to `__call__` and have the same effect.
+
+ Parameters
+ ----------
+ loop : bool, *optional*
+ if `loop` is True, then multiprocessing will be disabled and
+ all probes are evaluated by a single worker (default: False)
+
+ Returns
+ -------
+ results : see **Returns** in `evaluate`
+
+ other parameters
+ ----------------
+ kwargs : see **Parameters** in `configure`
+
+ """
+ self.configure(**kwargs)
+ if(not about.multiprocessing.status)or(loop):
+ return self._nonparallel_probing()
+ else:
+ return self._parallel_probing()
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def __repr__(self):
+ return "<nifty_core.probing>"
+
+##=============================================================================
+
+
+
+##-----------------------------------------------------------------------------
+
+class trace_probing(probing):
+ """
+ .. __
+ .. / /_
+ .. / _/ _____ ____ __ _______ _______
+ .. / / / __/ / _ / / ____/ / __ /
+ .. / /_ / / / /_/ / / /____ / /____/
+ .. \___/ /__/ \______| \______/ \______/ probing class
+
+ NIFTY subclass for trace probing (using multiprocessing)
+
+ When called, a trace_probing class instance samples the trace of an
+ operator or a function using random fields, whose components are random
+ variables with mean 0 and variance 1. When an instance is called it
+ returns the mean value of the scalar product of probe and f(probe),
+ where probe is a random field with mean 0 and variance 1.
+ The mean is calculated as 1/N Sum[ probe_i.dot(f(probe_i)) ].
+
+ Parameters
+ ----------
+ op : operator
+ The operator specified by `op` is the operator to be probed.
+ If no operator is given, then probing will be done by applying
+ `function` to the probes. (default: None)
+ function : function, *optional*
+ If no operator has been specified as `op`, then specification of
+ `function` is non optional. This is the function, that is applied
+ to the probes. (default: `op.times`)
+ domain : space, *optional*
+ If no operator has been specified as `op`, then specification of
+ `domain` is non optional. This is the space that the probes live
+ in. (default: `op.domain`)
+ target : domain, *optional*
+ `target` is the codomain of `domain`
+ (default: `op.domain.get_codomain()`)
+ random : string, *optional*
+ the distribution from which the probes are drawn. `random` can be
+ either "pm1" or "gau". "pm1" is a uniform distribution over {+1,-1}
+ or {+1,+i,-1,-i}, respectively. "gau" is a normal distribution with
+ zero-mean and unit-variance (default: "pm1")
+ ncpu : int, *optional*
+ the number of cpus to be used from parallel probing. (default: 2)
+ nrun : int, *optional*
+ the number of probes to be evaluated. If `nrun<ncpu`, it will be
+ set to `ncpu`. (default: 8)
+ nper : int, *optional*
+ this number specifies how many probes will be evaluated by one
+ worker. Afterwards a new worker will be created to evaluate a chunk
+ of `nper` probes.
+ If for example `nper=nrun/ncpu`, then every worker will be created
+ for every cpu. This can lead to the case, that all workers but one
+ are already finished, but have to wait for the last worker that
+ might still have a considerable amount of evaluations left. This is
+ obviously not very effective.
+ If on the other hand `nper=1`, then for each evaluation a worker will
+ be created. In this case all cpus will work until nrun probes have
+ been evaluated.
+ It is recommended to leave `nper` as the default value. (default: 1)
+ var : bool, *optional*
+ If `var` is True, then the variance of the sampled function will
+ also be returned. The result is then a tuple with the mean in the
+ zeroth entry and the variance in the first entry. (default: False)
+
+
+ See Also
+ --------
+ probing : The base probing class
+ diagonal_probing : A probing class to get the diagonal of an operator
+ operator.tr : the trace function uses trace probing in the case of non
+ diagonal operators
+
+
+ Attributes
+ ----------
+ function : function
+ the function, that is applied to the probes
+ domain : space
+ the space, where the probes live in
+ target : space
+ the codomain of `domain`
+ random : string
+ the random number generator used to create the probes
+ (either "pm1" or "gau")
+ ncpu : int
+ the number of cpus used for probing
+ nrun : int
+ the number of probes to be evaluated, when the instance is called
+ nper : int
+ number of probes, that will be evaluated by one worker
+ var : bool
+ whether the variance will be additionally returned, when the
+ instance is called
+ quargs : dict
+ Keyword arguments passed to `function` in each call.
+
+ """
+ def __init__(self,op,function=None,domain=None,target=None,random="pm1",ncpu=2,nrun=8,nper=1,var=False,**quargs):
+ """
+ initializes a trace probing instance
+
+ Parameters
+ ----------
+ op : operator
+ The operator specified by `op` is the operator to be probed.
+ If no operator is given, then probing will be done by applying
+ `function` to the probes. (default: None)
+ function : function, *optional*
+ If no operator has been specified as `op`, then specification of
+ `function` is non optional. This is the function, that is applied
+ to the probes. (default: `op.times`)
+ domain : space, *optional*
+ If no operator has been specified as `op`, then specification of
+ `domain` is non optional. This is the space that the probes live
+ in. (default: `op.domain`)
+ target : domain, *optional*
+ `target` is the codomain of `domain`
+ (default: `op.domain.get_codomain()`)
+ random : string, *optional*
+ the distribution from which the probes are drawn. `random` can be
+ either "pm1" or "gau". "pm1" is a uniform distribution over {+1,-1}
+ or {+1,+i,-1,-i}, respectively. "gau" is a normal distribution with
+ zero-mean and unit-variance (default: "pm1")
+ ncpu : int, *optional*
+ the number of cpus to be used from parallel probing. (default: 2)
+ nrun : int, *optional*
+ the number of probes to be evaluated. If `nrun<ncpu`, it will be
+ set to `ncpu`. (default: 8)
+ nper : int, *optional*
+ this number specifies how many probes will be evaluated by one
+ worker. Afterwards a new worker will be created to evaluate a chunk
+ of `nper` probes.
+ If for example `nper=nrun/ncpu`, then every worker will be created
+ for every cpu. This can lead to the case, that all workers but one
+ are already finished, but have to wait for the last worker that
+ might still have a considerable amount of evaluations left. This is
+ obviously not very effective.
+ If on the other hand `nper=1`, then for each evaluation a worker will
+ be created. In this case all cpus will work until nrun probes have
+ been evaluated.
+ It is recommended to leave `nper` as the default value. (default: 1)
+ var : bool, *optional*
+ If `var` is True, then the variance of the sampled function will
+ also be returned. The result is then a tuple with the mean in the
+ zeroth entry and the variance in the first entry. (default: False)
+
+ """
+ from nifty_operators import operator
+ if(not isinstance(op,operator)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ elif(op.nrow()!=op.ncol()):
+ raise ValueError(about._errors.cstring("ERROR: trace ill-defined for "+str(op.nrow())+" x "+str(op.ncol())+" matrices."))
+
+ ## check whether callable
+ if(function is None)or(not hasattr(function,"__call__")):
+ function = op.times
+ elif(op==function):
+ function = op.times
+ ## check whether correctly bound
+ if(op!=function.im_self):
+ raise NameError(about._errors.cstring("ERROR: invalid input."))
+ self.function = function
+
+ ## check given domain
+ if(domain is None)or(not isinstance(domain,space)):
+ if(self.function in [op.inverse_times,op.adjoint_times]):
+ domain = op.target
+ else:
+ domain = op.domain
+ else:
+ 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:
+ target = op.domain
+ self.domain = domain
+
+ ## check codomain
+ if(target is None):
+ target = self.domain.get_codomain()
+ else:
+ self.domain.check_codomain(target) ## a bit pointless
+ self.target = target
+
+ ## check degrees of freedom
+ if(op.domain.dof()>self.domain.dof()):
+ about.infos.cprint("INFO: variant numbers of degrees of freedom ( "+str(op.domain.dof())+" / "+str(self.domain.dof())+" ).")
+
+ if(random not in ["pm1","gau"]):
+ raise ValueError(about._errors.cstring("ERROR: unsupported random key '"+str(random)+"'."))
+ self.random = random
+
+ self.ncpu = int(max(1,ncpu))
+ self.nrun = int(max(self.ncpu,nrun))
+ if(nper is None):
+ self.nper = None
+ else:
+ self.nper = int(max(1,min(self.nrun//self.ncpu,nper)))
+
+ self.var = bool(var)
+
+ self.quargs = quargs
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def probing(self,idnum,probe):
+ """
+ Computes a single probing result given one probe
+
+ Parameters
+ ----------
+ probe : field
+ the field on which `function` will be applied
+ idnum : int
+ the identification number of the probing
+
+ Returns
+ -------
+ result : float
+ the result of `probe.dot(function(probe))`
+ """
+ f = self.function(probe,**self.quargs) ## field
+ 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
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def evaluate(self,summa,num,var):
+ """
+ Evaluates the probing results.
+
+ Parameters
+ ----------
+ summa : scalar
+ Sum of all probing results.
+ num : int
+ Number of successful probings (not returning ``None``).
+ var : scalar
+ Sum of all squared probing results
+
+ Returns
+ -------
+ final : scalar
+ The final probing result; 1/N Sum[ probing(probe_i) ]
+ var : scalar
+ The variance of the final probing result;
+ 1/(N(N-1)) Sum[( probing(probe_i) - final )^2];
+ if the variance is returned, the return will be a tuple of
+ (`final`,`var`).
+
+ """
+ if(num<self.nrun):
+ about.infos.cflush(" ( %u probe(s) failed, effectiveness == %.1f%% )\n"%(self.nrun-num,100*num/self.nrun))
+ if(num==0):
+ about.warnings.cprint("WARNING: probing failed.")
+ return None
+ else:
+ about.infos.cflush("\n")
+
+ if(issubclass(self.domain.datatype,np.complexfloating)):
+ summa = np.real(summa)
+
+ final = summa/num
+ if(self.var):
+ if(num==1):
+ about.warnings.cprint("WARNING: infinite variance.")
+ return final,None
+ else:
+ var = var/(num*(num-1))-np.real(np.conjugate(final)*final)/(num-1)
+ return final,var
+ else:
+ return final
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def _serial_probing(self,zipped): ## > performs the probing operation serially
+ try:
+ result = self._single_probing(zipped)
+ except Exception as exception:
+ raise exception
+ except BaseException: ## capture system-exiting exception including KeyboardInterrupt
+ raise Exception(about._errors.cstring("ERROR: unknown."))
+ else:
+ if(result is not None):
+ if(isinstance(_share.sum,tuple)):
+ _share.sum[0].acquire(block=True,timeout=None)
+ _share.sum[0].value += np.real(result)
+ _share.sum[0].release()
+ _share.sum[1].acquire(block=True,timeout=None)
+ _share.sum[1].value += np.imag(result)
+ _share.sum[1].release()
+ else:
+ _share.sum.acquire(block=True,timeout=None)
+ _share.sum.value += result
+ _share.sum.release()
+ _share.num.acquire(block=True,timeout=None)
+ _share.num.value += 1
+ _share.num.release()
+ if(self.var):
+ _share.var.acquire(block=True,timeout=None)
+ _share.var.value += np.real(np.conjugate(result)*result)
+ _share.var.release()
+ self._progress(_share.num.value)
+
+ def _parallel_probing(self): ## > performs the probing operations in parallel
+ ## define random seed
+ seed = np.random.randint(10**8,high=None,size=self.nrun)
+ ## define shared objects
+ if(issubclass(self.domain.datatype,np.complexfloating)):
+ _sum = (mv('d',0,lock=True),mv('d',0,lock=True)) ## tuple(real,imag)
+ else:
+ _sum = mv('d',0,lock=True)
+ _num = mv('I',0,lock=True)
+ if(self.var):
+ _var = mv('d',0,lock=True)
+ else:
+ _var = None
+ ## build pool
+ if(about.infos.status):
+ so.write(about.infos.cstring("INFO: multiprocessing "+(' ')*10))
+ so.flush()
+ pool = mp(processes=self.ncpu,initializer=_share._init_share,initargs=(_sum,_num,_var),maxtasksperchild=self.nper)
+ try:
+ ## pooling
+ pool.map_async(self._serial_probing,zip(seed,np.arange(self.nrun,dtype=np.int)),chunksize=None,callback=None).get(timeout=None)
+ ## close and join pool
+ about.infos.cflush(" done.")
+ pool.close()
+ pool.join()
+ except BaseException as exception:
+ ## terminate and join pool
+ about._errors.cprint("\nERROR: terminating pool.")
+ pool.terminate()
+ pool.join()
+ ## re-raise exception
+ raise exception ## traceback by looping
+ ## evaluate
+ if(issubclass(self.domain.datatype,np.complexfloating)):
+ _sum = np.complex(_sum[0].value,_sum[1].value)
+ else:
+ _sum = _sum.value
+ if(self.var):
+ _var = _var.value
+ return self.evaluate(_sum,_num.value,_var)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def __repr__(self):
+ return "<nifty_core.trace_probing>"
+
+##-----------------------------------------------------------------------------
+
+
+
+##-----------------------------------------------------------------------------
+
+class diagonal_probing(probing):
+ """
+ .. __ __ __
+ .. / / /__/ / /
+ .. ____/ / __ ____ __ ____ __ ______ __ ___ ____ __ / /
+ .. / _ / / / / _ / / _ / / _ | / _ | / _ / / /
+ .. / /_/ / / / / /_/ / / /_/ / / /_/ / / / / / / /_/ / / /_
+ .. \______| /__/ \______| \___ / \______/ /__/ /__/ \______| \___/ probing class
+ .. /______/
+
+ NIFTY subclass for diagonal probing (using multiprocessing)
+
+ When called, a diagonal_probing class instance samples the diagonal of
+ an operator or a function using random fields, whose components are
+ random variables with mean 0 and variance 1. When an instance is called
+ it returns the mean value of probe*f(probe), where probe is a random
+ field with mean 0 and variance 1.
+ The mean is calculated as 1/N Sum[ probe_i*f(probe_i) ]
+ ('*' denoting component wise multiplication)
+
+ Parameters
+ ----------
+ op : operator
+ The operator specified by `op` is the operator to be probed.
+ If no operator is given, then probing will be done by applying
+ `function` to the probes. (default: None)
+ function : function, *optional*
+ If no operator has been specified as `op`, then specification of
+ `function` is non optional. This is the function, that is applied
+ to the probes. (default: `op.times`)
+ domain : space, *optional*
+ If no operator has been specified as `op`, then specification of
+ `domain` is non optional. This is the space that the probes live
+ in. (default: `op.domain`)
+ target : domain, *optional*
+ `target` is the codomain of `domain`
+ (default: `op.domain.get_codomain()`)
+ random : string, *optional*
+ the distribution from which the probes are drawn. `random` can be
+ either "pm1" or "gau". "pm1" is a uniform distribution over {+1,-1}
+ or {+1,+i,-1,-i}, respectively. "gau" is a normal distribution with
+ zero-mean and unit-variance (default: "pm1")
+ ncpu : int, *optional*
+ the number of cpus to be used for parallel probing. (default: 2)
+ nrun : int, *optional*
+ the number of probes to be evaluated. If `nrun<ncpu`, it will be
+ set to `ncpu`. (default: 8)
+ nper : int, *optional*
+ this number specifies how many probes will be evaluated by one
+ worker. Afterwards a new worker will be created to evaluate a chunk
+ of `nper` probes.
+ If for example `nper=nrun/ncpu`, then every worker will be created
+ for every cpu. This can lead to the case, that all workers but one
+ are already finished, but have to wait for the last worker that
+ might still have a considerable amount of evaluations left. This is
+ obviously not very effective.
+ If on the other hand `nper=1`, then for each evaluation a worker will
+ be created. In this case all cpus will work until nrun probes have
+ been evaluated.
+ It is recommended to leave `nper` as the default value. (default: 1)
+ var : bool, *optional*
+ If `var` is True, then the variance of the sampled function will
+ also be returned. The result is then a tuple with the mean in the
+ zeroth entry and the variance in the first entry. (default: False)
+ save : bool, *optional*
+ If `save` is True, then the probing results will be written to the
+ hard disk instead of being saved in the RAM. This is recommended
+ for high dimensional fields whose probes would otherwise fill up
+ the memory. (default: False)
+ path : string, *optional*
+ the path, where the probing results are saved, if `save` is True.
+ (default: "tmp")
+ prefix : string, *optional*
+ a prefix for the saved probing results. The saved results will be
+ named using that prefix and an 8-digit number
+ (e.g. "<prefix>00000001.npy"). (default: "")
+
+
+ See Also
+ --------
+ trace_probing : A probing class to get the trace of an operator
+ probing : The base probing class
+ operator.diag : The diag function uses diagonal probing in the case of
+ non diagonal operators
+ operator.hat : The hat function uses diagonal probing in the case of
+ non diagonal operators
+
+
+ Attributes
+ ----------
+ function : function
+ the function, that is applied to the probes
+ domain : space
+ the space, where the probes live in
+ target : space
+ the codomain of `domain`
+ random : string
+ the random number generator used to create the probes
+ (either "pm1" or "gau")
+ ncpu : int
+ the number of cpus used for probing
+ nrun : int
+ the number of probes to be evaluated, when the instance is called
+ nper : int
+ number of probes, that will be evaluated by one worker
+ var : bool
+ whether the variance will be additionally returned, when the
+ instance is called
+ save : {string, None}
+ the path and prefix for saved probe files. None in the case where
+ the probing results are stored in the RAM.
+ quargs : dict
+ Keyword arguments passed to `function` in each call.
+
+ """
+ def __init__(self,op,function=None,domain=None,target=None,random="pm1",ncpu=2,nrun=8,nper=1,var=False,save=False,path="tmp",prefix="",**quargs):
+ """
+ initializes a diagonal probing instance
+
+ Parameters
+ ----------
+ op : operator
+ The operator specified by `op` is the operator to be probed.
+ If no operator is given, then probing will be done by applying
+ `function` to the probes. (default: None)
+ function : function, *optional*
+ If no operator has been specified as `op`, then specification of
+ `function` is non optional. This is the function, that is applied
+ to the probes. (default: `op.times`)
+ domain : space, *optional*
+ If no operator has been specified as `op`, then specification of
+ `domain` is non optional. This is the space that the probes live
+ in. (default: `op.domain`)
+ target : domain, *optional*
+ `target` is the codomain of `domain`
+ (default: `op.domain.get_codomain()`)
+ random : string, *optional*
+ the distribution from which the probes are drawn. `random` can be
+ either "pm1" or "gau". "pm1" is a uniform distribution over {+1,-1}
+ or {+1,+i,-1,-i}, respectively. "gau" is a normal distribution with
+ zero-mean and unit-variance (default: "pm1")
+ ncpu : int, *optional*
+ the number of cpus to be used for parallel probing. (default: 2)
+ nrun : int, *optional*
+ the number of probes to be evaluated. If `nrun<ncpu`, it will be
+ set to `ncpu`. (default: 8)
+ nper : int, *optional*
+ this number specifies how many probes will be evaluated by one
+ worker. Afterwards a new worker will be created to evaluate a chunk
+ of `nper` probes.
+ If for example `nper=nrun/ncpu`, then every worker will be created
+ for every cpu. This can lead to the case, that all workers but one
+ are already finished, but have to wait for the last worker that
+ might still have a considerable amount of evaluations left. This is
+ obviously not very effective.
+ If on the other hand `nper=1`, then for each evaluation a worker will
+ be created. In this case all cpus will work until nrun probes have
+ been evaluated.
+ It is recommended to leave `nper` as the default value. (default: 8)
+ var : bool, *optional*
+ If `var` is True, then the variance of the sampled function will
+ also be returned. The result is then a tuple with the mean in the
+ zeroth entry and the variance in the first entry. (default: False)
+ save : bool, *optional*
+ If `save` is True, then the probing results will be written to the
+ hard disk instead of being saved in the RAM. This is recommended
+ for high dimensional fields whose probes would otherwise fill up
+ the memory. (default: False)
+ path : string, *optional*
+ the path, where the probing results are saved, if `save` is True.
+ (default: "tmp")
+ prefix : string, *optional*
+ a prefix for the saved probing results. The saved results will be
+ named using that prefix and an 8-digit number
+ (e.g. "<prefix>00000001.npy"). (default: "")
+
+ """
+ from nifty_operators import operator
+ if(not isinstance(op,operator)):
+ raise TypeError(about._errors.cstring("ERROR: invalid input."))
+ elif(op.nrow()!=op.ncol()):
+ raise ValueError(about._errors.cstring("ERROR: diagonal ill-defined for "+str(op.nrow())+" x "+str(op.ncol())+" matrices."))
+
+ ## check whether callable
+ if(function is None)or(not hasattr(function,"__call__")):
+ function = op.times
+ elif(op==function):
+ function = op.times
+ ## check whether correctly bound
+ if(op!=function.im_self):
+ raise NameError(about._errors.cstring("ERROR: invalid input."))
+ self.function = function
+
+ ## check given domain
+ if(domain is None)or(not isinstance(domain,space)):
+ if(self.function in [op.inverse_times,op.adjoint_times]):
+ domain = op.target
+ else:
+ domain = op.domain
+ else:
+ 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:
+ target = op.domain
+ self.domain = domain
+
+ ## check codomain
+ if(target is None):
+ target = self.domain.get_codomain()
+ else:
+ self.domain.check_codomain(target) ## a bit pointless
+ self.target = target
+
+ ## check degrees of freedom
+ if(self.domain.dof()>op.domain.dof()):
+ about.infos.cprint("INFO: variant numbers of degrees of freedom ( "+str(self.domain.dof())+" / "+str(op.domain.dof())+" ).")
+
+ if(random not in ["pm1","gau"]):
+ raise ValueError(about._errors.cstring("ERROR: unsupported random key '"+str(random)+"'."))
+ self.random = random
+
+ self.ncpu = int(max(1,ncpu))
+ self.nrun = int(max(self.ncpu,nrun))
+ if(nper is None):
+ self.nper = None
+ else:
+ self.nper = int(max(1,min(self.nrun//self.ncpu,nper)))
+
+ self.var = bool(var)
+
+ if(save):
+ path = os.path.expanduser(str(path))
+ if(not os.path.exists(path)):
+ os.makedirs(path)
+ self.save = os.path.join(path,str(prefix)) ## (back)slash inserted if needed
+ else:
+ self.save = None
+
+ self.quargs = quargs
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def configure(self,**kwargs):
+ """
+ changes the attributes of the instance
+
+ Parameters
+ ----------
+ random : string, *optional*
+ the random number generator used to create the probes
+ (default: "pm1")
+ ncpu : int, *optional*
+ the number of cpus to be used for parallel probing
+ (default: 2)
+ nrun : int, *optional*
+ the number of probes to be evaluated. If `nrun<ncpu`, it will
+ be set to `ncpu`. (default: 8)
+ nper : int, *optional*
+ number of probes, that will be evaluated by one worker
+ (default: 1)
+ var : bool, *optional*
+ whether the variance will be additionally returned
+ (default: False)
+ save : bool, *optional*
+ whether the individual probing results will be saved to the HDD
+ (default: False)
+ path : string, *optional*
+ the path, where the probing results are saved (default: "tmp")
+ prefix : string, *optional*
+ a prefix for the saved probing results (default: "")
+
+ """
+ if("random" in kwargs):
+ if(kwargs.get("random") not in ["pm1","gau"]):
+ raise ValueError(about._errors.cstring("ERROR: unsupported random key '"+str(kwargs.get("random"))+"'."))
+ self.random = kwargs.get("random")
+
+ if("ncpu" in kwargs):
+ self.ncpu = int(max(1,kwargs.get("ncpu")))
+ if("nrun" in kwargs):
+ self.nrun = int(max(self.ncpu,kwargs.get("nrun")))
+ if("nper" in kwargs):
+ if(kwargs.get("nper") is None):
+ self.nper = None
+ else:
+ self.nper = int(max(1,min(self.nrun//self.ncpu,kwargs.get("nper"))))
+
+ if("var" in kwargs):
+ self.var = bool(kwargs.get("var"))
+
+ if("save" in kwargs):
+ if(kwargs.get("save")):
+ if("path" in kwargs):
+ path = kwargs.get("path")
+ else:
+ if(self.save is not None):
+ about.warnings.cprint("WARNING: save path set to default.")
+ path = "tmp"
+ if("prefix" in kwargs):
+ prefix = kwargs.get("prefix")
+ else:
+ if(self.save is not None):
+ about.warnings.cprint("WARNING: save prefix set to default.")
+ prefix = ""
+ path = os.path.expanduser(str(path))
+ if(not os.path.exists(path)):
+ os.makedirs(path)
+ self.save = os.path.join(path,str(prefix)) ## (back)slash inserted if needed
+ else:
+ self.save = None
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def probing(self,idnum,probe):
+
+ """
+ Computes a single probing result given one probe
+
+ Parameters
+ ----------
+ probe : field
+ the field on which `function` will be applied
+ idnum : int
+ the identification number of the probing
+
+ Returns
+ -------
+ result : ndarray
+ the result of `probe*(function(probe))`
+ """
+ f = self.function(probe,**self.quargs) ## field
+ 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)
+ return result
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def _serial_probing(self,zipped): ## > performs the probing operation serially
+ try:
+ result = self._single_probing(zipped)
+ except Exception as exception:
+ raise exception
+ except BaseException: ## capture system-exiting exception including KeyboardInterrupt
+ raise Exception(about._errors.cstring("ERROR: unknown."))
+ else:
+ if(result is not None):
+ result = np.array(result).flatten(order='C')
+ if(isinstance(_share.sum,tuple)):
+ _share.sum[0].acquire(block=True,timeout=None)
+ _share.sum[0][:] += np.real(result)
+ _share.sum[0].release()
+ _share.sum[1].acquire(block=True,timeout=None)
+ _share.sum[1][:] += np.imag(result)
+ _share.sum[1].release()
+ else:
+ _share.sum.acquire(block=True,timeout=None)
+ _share.sum[:] += result
+ _share.sum.release()
+ _share.num.acquire(block=True,timeout=None)
+ _share.num.value += 1
+ _share.num.release()
+ if(self.var):
+ _share.var.acquire(block=True,timeout=None)
+ _share.var[:] += np.real(np.conjugate(result)*result)
+ _share.var.release()
+ self._progress(_share.num.value)
+
+ def _parallel_probing(self): ## > performs the probing operations in parallel
+ ## define random seed
+ seed = np.random.randint(10**8,high=None,size=self.nrun)
+ ## define shared objects
+ if(issubclass(self.domain.datatype,np.complexfloating)):
+ _sum = (ma('d',np.zeros(self.domain.dim(split=False),dtype=np.float64,order='C'),lock=True),ma('d',np.zeros(self.domain.dim(split=False),dtype=np.float64,order='C'),lock=True)) ## tuple(real,imag)
+ else:
+ _sum = ma('d',np.zeros(self.domain.dim(split=False),dtype=np.float64,order='C'),lock=True)
+ _num = mv('I',0,lock=True)
+ if(self.var):
+ _var = ma('d',np.zeros(self.domain.dim(split=False),dtype=np.float64,order='C'),lock=True)
+ else:
+ _var = None
+ ## build pool
+ if(about.infos.status):
+ so.write(about.infos.cstring("INFO: multiprocessing "+(' ')*10))
+ so.flush()
+ pool = mp(processes=self.ncpu,initializer=_share._init_share,initargs=(_sum,_num,_var),maxtasksperchild=self.nper)
+ try:
+ ## pooling
+ pool.map_async(self._serial_probing,zip(seed,np.arange(self.nrun,dtype=np.int)),chunksize=None,callback=None).get(timeout=None)
+ ## close and join pool
+ about.infos.cflush(" done.")
+ pool.close()
+ pool.join()
+ except BaseException as exception:
+ ## terminate and join pool
+ about._errors.cprint("\nERROR: terminating pool.")
+ pool.terminate()
+ pool.join()
+ ## re-raise exception
+ raise exception ## traceback by looping
+ ## evaluate
+ if(issubclass(self.domain.datatype,np.complexfloating)):
+ _sum = (np.array(_sum[0][:])+np.array(_sum[1][:])*1j).reshape(self.domain.dim(split=True)) ## comlpex array
+ else:
+ _sum = np.array(_sum[:]).reshape(self.domain.dim(split=True))
+ if(self.var):
+ _var = np.array(_var[:]).reshape(self.domain.dim(split=True))
+ return self.evaluate(_sum,_num.value,_var)
+
+ ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+ def __repr__(self):
+ return "<nifty_core.diagonal_probing>"
+
+##-----------------------------------------------------------------------------
+
+## IDEA: diagonal_inference