Commit c65a30a8 authored by theos's avatar theos

First steps towards a PropagatorOperator

parent 280e2cde
......@@ -48,8 +48,9 @@ from random import Random
from nifty_simple_math import *
from nifty_utilities import *
from field_types import FieldType,\
FieldArray
from field_types import *
from minimization import *
from spaces import *
......
# -*- coding: utf-8 -*-
from conjugate_gradient import ConjugateGradient
# -*- coding: utf-8 -*-
import numpy as np
from nifty.config import notification, about
class ConjugateGradient(object):
"""
.. _______ ____ __
.. / _____/ / _ /
.. / /____ __ / /_/ / __
.. \______//__/ \____ //__/ class
.. /______/
NIFTY tool class for conjugate gradient
This tool minimizes :math:`A x = b` with respect to `x` given `A` and
`b` using a conjugate gradient; i.e., a step-by-step minimization
relying on conjugated gradient directions. Further, `A` is assumed to
be a positive definite and self-adjoint operator. The use of a
preconditioner `W` that is roughly the inverse of `A` is optional.
For details on the methodology refer to [#]_, for details on usage and
output, see the notes below.
Parameters
----------
A : {operator, function}
Operator `A` applicable to a field.
b : field
Resulting field of the operation `A(x)`.
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).
See Also
--------
scipy.sparse.linalg.cg
Notes
-----
After initialization by `__init__`, the minimizer is started by calling
it using `__call__`, which takes additional parameters. Notifications,
if enabled, will state the iteration number, current step widths
`alpha` and `beta`, the current relative residual `delta` that is
compared to the tolerance, and the convergence level if changed.
The minimizer will exit in three states: DEAD if alpha becomes
infinite, QUIT if the maximum number of iterations is reached, or DONE
if convergence is achieved. Returned will be the latest `x` and the
latest convergence level, which can evaluate ``True`` for the exit
states QUIT and DONE.
References
----------
.. [#] J. R. Shewchuk, 1994, `"An Introduction to the Conjugate
Gradient Method Without the Agonizing Pain"
<http://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf>`_
Examples
--------
>>> b = field(point_space(2), val=[1, 9])
>>> A = diagonal_operator(b.domain, diag=[4, 3])
>>> x,convergence = conjugate_gradient(A, b, note=True)(tol=1E-4, clevel=3)
iteration : 00000001 alpha = 3.3E-01 beta = 1.3E-03 delta = 3.6E-02
iteration : 00000002 alpha = 2.5E-01 beta = 7.6E-04 delta = 1.0E-03
iteration : 00000003 alpha = 3.3E-01 beta = 2.5E-04 delta = 1.6E-05 convergence level : 1
iteration : 00000004 alpha = 2.5E-01 beta = 1.8E-06 delta = 2.1E-08 convergence level : 2
iteration : 00000005 alpha = 2.5E-01 beta = 2.2E-03 delta = 1.0E-09 convergence level : 3
... done.
>>> bool(convergence)
True
>>> x.val # yields 1/4 and 9/3
array([ 0.25, 3. ])
Attributes
----------
A : {operator, function}
Operator `A` applicable to a field.
x : field
Current field.
b : field
Resulting field of the operation `A(x)`.
W : {operator, function}
Operator `W` that is a preconditioner on `A` and is applicable to a
field; can be ``None``.
spam : function
Callback function which is given the current `x` and iteration
counter each iteration; can be ``None``.
reset : integer
Number of iterations after which to restart; i.e., forget previous
conjugated directions (default: sqrt(b.dim)).
note : notification
Notification instance.
"""
def __init__(self):
pass
def __call__(self, A, b, x0=None, W=None, spam=None, reset=None,
note=False, **kwargs):
"""
Runs the conjugate gradient minimization.
Parameters
----------
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).
Returns
-------
x : field
Latest `x` of the minimization.
convergence : integer
Latest convergence level indicating whether the minimization
has converged or not.
"""
if hasattr(A, "__call__"):
self.A = A
else:
raise AttributeError(about._errors.cstring(
"ERROR: A must be callable!"))
self.b = b
if W is None or hasattr(W, "__call__"):
self.W = W
else:
raise AttributeError(about._errors.cstring(
"ERROR: W must be None or callable!"))
self.spam = spam
if reset is None:
self.reset = max(2, int(np.sqrt(b.domain.dim)))
else:
self.reset = max(2, int(reset))
self.note = notification(default=bool(note))
self.x = self.b.copy_empty()
self.x.set_val(new_val=0)
self.x.set_val(new_val=x0)
if self.W is None:
return self._calc_without(**kwargs)
else:
return self._calc_with(**kwargs)
def _calc_without(self, tol=1E-4, clevel=1, limii=None):
clevel = int(clevel)
if limii is None:
limii = 10*self.b.domain.dim
else:
limii = int(limii)
r = self.b-self.A(self.x)
print ('r', r.val)
d = self.b.copy_empty()
d.set_val(new_val = r.get_val())
gamma = r.dot(d)
if gamma==0:
return self.x, clevel+1
delta_ = np.absolute(gamma)**(-0.5)
convergence = 0
ii = 1
while(True):
# print ('gamma', gamma)
q = self.A(d)
# print ('q', q.val)
alpha = gamma/d.dot(q) ## positive definite
if np.isfinite(alpha) == False:
self.note.cprint(
"\niteration : %08u alpha = NAN\n... dead."%ii)
return self.x, 0
self.x += d * alpha
# print ('x', self.x.val)
if np.signbit(np.real(alpha)) == True:
about.warnings.cprint(
"WARNING: positive definiteness of A violated.")
r = self.b-self.A(self.x)
elif (ii%self.reset) == 0:
r = self.b-self.A(self.x)
else:
r -= q * alpha
# print ('r', r.val)
gamma_ = gamma
gamma = r.dot(r)
# print ('gamma', gamma)
beta = max(0, gamma/gamma_) ## positive definite
# print ('d*beta', beta, (d*beta).val)
d = r + d*beta
# print ('d', d.val)
delta = delta_*np.absolute(gamma)**0.5
self.note.cflush(
"\niteration : %08u alpha = %3.1E beta = %3.1E delta = %3.1E"\
%(ii,np.real(alpha),np.real(beta),np.real(delta)))
if gamma == 0:
convergence = clevel+1
self.note.cprint(" convergence level : INF\n... done.")
break
elif np.absolute(delta)<tol:
convergence += 1
self.note.cflush(" convergence level : %u"%convergence)
if(convergence==clevel):
self.note.cprint("\n... done.")
break
else:
convergence = max(0, convergence-1)
if ii==limii:
self.note.cprint("\n... quit.")
break
if (self.spam is not None):
self.spam(self.x, ii)
ii += 1
if (self.spam is not None):
self.spam(self.x,ii)
return self.x, convergence
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def _calc_with(self, tol=1E-4, clevel=1, limii=None): ## > runs cg with preconditioner
clevel = int(clevel)
if(limii is None):
limii = 10*self.b.domain.dim
else:
limii = int(limii)
r = self.b-self.A(self.x)
d = self.W(r)
gamma = r.dot(d)
if gamma==0:
return self.x, clevel+1
delta_ = np.absolute(gamma)**(-0.5)
convergence = 0
ii = 1
while(True):
q = self.A(d)
alpha = gamma/d.dot(q) ## positive definite
if np.isfinite(alpha) == False:
self.note.cprint(
"\niteration : %08u alpha = NAN\n... dead."%ii)
return self.x, 0
self.x += d * alpha ## update
if np.signbit(np.real(alpha)) == True:
about.warnings.cprint(
"WARNING: positive definiteness of A violated.")
r = self.b-self.A(self.x)
elif (ii%self.reset) == 0:
r = self.b-self.A(self.x)
else:
r -= q * alpha
s = self.W(r)
gamma_ = gamma
gamma = r.dot(s)
if np.signbit(np.real(gamma)) == True:
about.warnings.cprint(
"WARNING: positive definiteness of W violated.")
beta = max(0, gamma/gamma_) ## positive definite
d = s + d*beta ## conjugated gradient
delta = delta_*np.absolute(gamma)**0.5
self.note.cflush(
"\niteration : %08u alpha = %3.1E beta = %3.1E delta = %3.1E"\
%(ii,np.real(alpha),np.real(beta),np.real(delta)))
if gamma==0:
convergence = clevel+1
self.note.cprint(" convergence level : INF\n... done.")
break
elif np.absolute(delta)<tol:
convergence += 1
self.note.cflush(" convergence level : %u"%convergence)
if convergence==clevel:
self.note.cprint("\n... done.")
break
else:
convergence = max(0, convergence-1)
if ii==limii:
self.note.cprint("\n... quit.")
break
if (self.spam is not None):
self.spam(self.x,ii)
ii += 1
if (self.spam is not None):
self.spam(self.x,ii)
return self.x, convergence
##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def __repr__(self):
return "<nifty_tools.conjugate_gradient>"
##=============================================================================
......@@ -29,6 +29,8 @@ from endomorphic_operator import EndomorphicOperator
from fft_operator import *
from propagator_operator import PropagatorOperator
from nifty_operators import operator,\
diagonal_operator,\
power_operator,\
......
......@@ -163,15 +163,17 @@ class DiagonalOperator(EndomorphicOperator):
# the one of x, reshape the local data of self and apply it directly
active_axes = []
if spaces is None:
for axes in x.domain_axes:
active_axes += axes
if self.domain != ():
for axes in x.domain_axes:
active_axes += axes
else:
for space_index in spaces:
active_axes += x.domain_axes[space_index]
if types is None:
for axes in x.field_type_axes:
active_axes += axes
if self.field_type != ():
for axes in x.field_type_axes:
active_axes += axes
else:
for type_index in types:
active_axes += x.field_type_axes[type_index]
......
......@@ -186,7 +186,7 @@ class LinearOperator(object):
"match."))
else:
for i, field_type_index in enumerate(types):
if x.field_types[field_type_index] != self_field_type[i]:
if x.field_type[field_type_index] != self_field_type[i]:
raise ValueError(about._errors.cstring(
"ERROR: The operator's and and field's field_type "
"don't match."))
......
# -*- coding: utf-8 -*-
from propagator_operator import PropagatorOperator
# -*- coding: utf-8 -*-
from nifty.config import about
from nifty.minimization import ConjugateGradient
from nifty.operators.linear_operator import LinearOperator
class PropagatorOperator(LinearOperator):
"""
.. __
.. / /_
.. _______ _____ ______ ______ ____ __ ____ __ ____ __ / _/ ______ _____
.. / _ / / __/ / _ | / _ | / _ / / _ / / _ / / / / _ | / __/
.. / /_/ / / / / /_/ / / /_/ / / /_/ / / /_/ / / /_/ / / /_ / /_/ / / /
.. / ____/ /__/ \______/ / ____/ \______| \___ / \______| \___/ \______/ /__/ 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.
"""
# ---Overwritten properties and methods---
def __init__(self, S=None, M=None, R=None, N=None, inverter=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.
"""
self.S = S
self.S_inverse_times = self.S.inverse_times
# build up the likelihood contribution
(self.M_times,
M_domain,
M_field_type,
M_target,
M_field_type_target) = self._build_likelihood_contribution(M, R, N)
# assert that S and M have matching domains
if not (self.domain == M_domain and
self.field_type == M_target and
self.target == M_target and
self.field_type_target == M_field_type_target):
raise ValueError(about._errors.cstring(
"ERROR: The domains and targets of the prior " +
"signal covariance and the likelihood contribution must be " +
"the same in the sense of '=='."))
if inverter is not None:
self.inverter = inverter
else:
self.inverter = conjugate_gradient()
# ---Mandatory properties and methods---
@property
def domain(self):
return self.S.domain
@property
def field_type(self):
return self.S.field_type
@property
def target(self):
return self.S.target
@property
def field_type_target(self):
return self.S.field_type_target
@property
def implemented(self):
return True
@property
def symmetric(self):
return True
@property
def unitary(self):
return False
# ---Added properties and methods---
def _build_likelihood_contribution(self, M, R, N):
# if a M is given, return its times method and its domains
# supplier and discard R and N
if M is not None:
return (M.times, M.domain, M.field_type, M.target, M.cotarget)
if N is not None:
if R is not None:
return (lambda z: R.adjoint_times(N.inverse_times(R.times(z))),
R.domain, R.field_type, R.domain, R.field_type)
else:
return (N.inverse_times,
N.domain, N.field_type, N.target, N.field_type_target)
else:
raise ValueError(about._errors.cstring(
"ERROR: At least M or N must be given."))
def _multiply(self, x, W=None, spam=None, reset=None, note=False,
x0=None, tol=1E-4, clevel=1, limii=None, **kwargs):
if W is None:
W = self.S
(result, convergence) = self.inverter(A=self._inverse_multiply,
b=x,
W=W,
spam=spam,
reset=reset,
note=note,
x0=x0,
tol=tol,
clevel=clevel,
limii=limii)
# evaluate
if not convergence:
about.warnings.cprint("WARNING: conjugate gradient failed.")
return result
def _inverse_multiply(self, x, **kwargs):
result = self.S_inverse_times(x)
result += self.M_times(x)
return result
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment