Commit 9c1ab786 authored by Ultimanet's avatar Ultimanet
Browse files

Removed nifty_explicit and nifty_power_conversion

parent 9a4c6546
## 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 nifty_core import *
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_core import about, \
space, \
field, \
operator,diagonal_operator,identity,vecvec_operator, \
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)