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

bugfix in hp.get_plot; improvements in probing._parallel_probing;...

bugfix in hp.get_plot; improvements in probing._parallel_probing; nifty_tools.py added; version number updated.
parent d876195c
......@@ -484,7 +484,7 @@ class _about(object): ## nifty support class for global settings
"""
## version
self._version = "0.4.2"
self._version = "0.5.0"
## switches and notifications
self._errors = notification(default=True,ccode=notification._code)
......@@ -2151,7 +2151,7 @@ class rg_space(space):
.. / __/ / _ /
.. / / / /_/ /
.. /__/ \____ / space class
.. /_____/
.. /______/
NIFTY subclass for spaces of regular Cartesian grids.
......@@ -5427,6 +5427,7 @@ class hp_space(space):
cmap = pl.cm.jet ## default
cmap.set_under(color='k',alpha=0.0) ## transparent box
hp.mollview(x,fig=None,rot=None,coord=None,unit=unit,xsize=800,title=title,nest=False,min=vmin,max=vmax,flip="astro",remove_dip=False,remove_mono=False,gal_cut=0,format="%g",format2="%g",cbar=cbar,cmap=cmap,notext=False,norm=norm,hold=False,margins=None,sub=None)
fig = pl.gcf()
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)
......@@ -7437,8 +7438,6 @@ def arctanh(x):
else:
return np.arctanh(np.array(x))
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def sqrt(x):
"""
Returns the square root of a given object.
......@@ -7466,8 +7465,6 @@ def sqrt(x):
else:
return np.sqrt(np.array(x))
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def exp(x):
"""
Returns the exponential of a given object.
......@@ -7541,8 +7538,6 @@ def log(x,base=None):
else:
raise ValueError(about._errors.cstring("ERROR: invalid input basis."))
## - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
def conjugate(x):
"""
Computes the complex conjugate of a given object.
......@@ -7565,9 +7560,6 @@ def conjugate(x):
##-----------------------------------------------------------------------------
##-----------------------------------------------------------------------------
......@@ -7679,7 +7671,7 @@ class operator(object):
self.uni = bool(uni)
if(self.domain.discrete):
self.imp = False
self.imp = True
else:
self.imp = bool(imp)
......@@ -9245,7 +9237,7 @@ class power_operator(diagonal_operator):
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def set_power(self,newspec,bare=None,pindex=None,**kwargs):
def set_power(self,newspec,bare=True,pindex=None,**kwargs):
"""
Sets the power spectrum of the diagonal operator
......@@ -9280,14 +9272,12 @@ class power_operator(diagonal_operator):
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).
(default: None).
"""
if(bare is None):
about.warnings.cprint("WARNING: bare keyword set to default.")
bare = True
# if(bare is None):
# about.warnings.cprint("WARNING: bare keyword set to default.")
# bare = True
## check implicit pindex
if(pindex is None):
try:
......@@ -10615,13 +10605,15 @@ class probing(object):
pool = mp(processes=self.ncpu,initializer=None,initargs=(),maxtasksperchild=self.nper)
try:
## retrieve results
results = pool.map_async(self._serial_probing,zip(seed,np.arange(self.nrun,dtype=np.int)),chunksize=None,callback=None).get(timeout=None)
## close pool
results = pool.map(self._serial_probing,zip(seed,np.arange(self.nrun,dtype=np.int)),chunksize=None)#,callback=None).get(timeout=None) ## map_async replaced
## close and join pool
about.infos.cflush(" done.")
pool.close()
pool.join()
except:
## terminate pool
## terminate and join pool
pool.terminate()
pool.join()
raise Exception(about._errors.cstring("ERROR: unknown. NOTE: pool terminated.")) ## traceback by looping
## cleanup
results = [rr for rr in results if(rr is not None)]
......
## 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/>.
"""
.. __ ____ __
.. /__/ / _/ / /_
.. __ ___ __ / /_ / _/ __ __
.. / _ | / / / _/ / / / / / /
.. / / / / / / / / / /_ / /_/ /
.. /__/ /__/ /__/ /__/ \___/ \___ / tools
.. /______/
A nifty set of tools.
## TODO: *DESCRIPTION*
"""
from __future__ import division
#import numpy as np
from nifty.nifty_core import *
##-----------------------------------------------------------------------------
class invertible_operator(operator):
"""
.. __ __ __ __ __
.. /__/ / /_ /__/ / / / /
.. __ __ ___ __ __ _______ _____ / _/ __ / /___ / / _______
.. / / / _ || |/ / / __ / / __/ / / / / / _ | / / / __ /
.. / / / / / / | / / /____/ / / / /_ / / / /_/ / / /_ / /____/
.. /__/ /__/ /__/ |____/ \______/ /__/ \___/ /__/ \______/ \___/ \______/ operator class
NIFTY subclass for invertible, self-adjoint (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.
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
-----
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 = target
if(para is not None):
self.para = para
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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.
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`` is
forced incase the conjugate gradient fails.
Returns
-------
Ox : 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_ = self._briefing(x,self.domain,False)
## apply operator
if(self.imp):
A = self._inverse_multiply
else:
A = self.inverse_times
x_,converged = conjugate_gradient(A,x_,W=W,spam=spam,reset=reset,note=note)(x0=x0,tol=tol,clevel=clevel,limii=limii)
## evaluate
if(not converged):
if(not force):
return None
about.warnings.cprint("WARNING: conjugate gradient failed.")
return self._debriefing(x,x_,self.target,False)
def inverse_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.
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`` is
forced incase 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()).
"""
## 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
if(self.imp):
A = self._multiply
else:
A = self.times
x_,converged = conjugate_gradient(A,x_,W=W,spam=spam,reset=reset,note=note)(x0=x0,tol=tol,clevel=clevel,limii=limii)
## evaluate
if(not converged):
if(not force):
return None
about.warnings.cprint("WARNING: conjugate gradient failed.")
return self._debriefing(x,x_,self.domain,True)
##-----------------------------------------------------------------------------
##-----------------------------------------------------------------------------
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): ## 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): ## 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)+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)).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.
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`` is
forced incase 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_,converged = conjugate_gradient(A,x_,W=W,spam=spam,reset=reset,note=note)(x0=x0,tol=tol,clevel=clevel,limii=limii)