## 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: ## ## 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 . """ .. __ ____ __ .. /__/ / _/ / /_ .. __ ___ __ / /_ / _/ __ __ .. / _ | / / / _/ / / / / / / .. / / / / / / / / / /_ / /_/ / .. /__/ /__/ /__/ /__/ \___/ \___ / tools .. /______/ A nifty set of tools. ## TODO: *DESCRIPTION* """ from __future__ import division #import numpy as np from 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 = self.domain 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 elif(id(self.inverse_times)==id(invertible_operator.inverse_times)): ## avoid infinite recursion A = super(invertible_operator,self).inverse_times else: A = self.inverse_times 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_,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 elif(id(self.times)==id(invertible_operator.times)): ## avoid infinite recursion A = super(invertible_operator,self).times else: A = self.times 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_,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_,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) ##----------------------------------------------------------------------------- ##============================================================================= class conjugate_gradient(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). 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 two states: 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 all exit states. References ---------- .. [#] J. R. Shewchuk, 1994, `"An Introduction to the Conjugate Gradient Method Without the Agonizing Pain" ``_ 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,A,b,W=None,spam=None,reset=None,note=False): """ Initializes the conjugate_gradient and sets the attributes (except for `x`). 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). """ if(hasattr(A,"__call__")): self.A = A ## applies A else: raise AttributeError(about._errors.cstring("ERROR: invalid input.")) self.b = b if(W is None)or(hasattr(W,"__call__")): self.W = W ## applies W ~ A_inverse else: raise AttributeError(about._errors.cstring("ERROR: invalid input.")) self.spam = spam ## serves as callback given x and iteration number if(reset is None): ## 2 < reset ~ sqrt(dim) self.reset = max(2,int(np.sqrt(b.domain.dim(split=False)))) else: self.reset = max(2,int(reset)) self.note = notification(default=bool(note)) ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def __call__(self,x0=None,**kwargs): ## > runs cg with/without preconditioner """ 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. """ self.x = field(self.b.domain,val=x0,target=self.b.target) 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): ## > runs cg without preconditioner if(limii is None): limii = 10*self.b.domain.dim(split=False) d = r = self.b-self.A(self.x) gamma = r.dot(d) delta_ = np.absolute(gamma)**(-0.5) convergence = 0 ii = 1 while(True): q = self.A(d) alpha = gamma/d.dot(q) ## positive definite self.x += alpha*d if(ii%self.reset==0)or(np.signbit(np.real(alpha))): r = self.b-self.A(self.x) else: r -= alpha*q gamma_ = gamma gamma = r.dot(r) beta = max(0,gamma/gamma_) ## positive definite d = r+beta*d delta = delta_*np.absolute(gamma)**0.5 self.note.cflush("\niteration : %08u alpha = %3.1E beta = %3.1E delta = %3.1E"%(ii,alpha,beta,delta)) if(ii==limii): self.note.cprint("\n... quit.") break if(gamma==0): convergence = clevel+1 self.note.cprint(" convergence level : INF\n... done.") break elif(np.absolute(delta) runs cg with preconditioner if(limii is None): limii = 10*self.b.domain.dim(split=False) r = self.b-self.A(self.x) d = self.W(r) gamma = r.dot(d) delta_ = np.absolute(gamma)**(-0.5) convergence = 0 ii = 1 while(True): q = self.A(d) alpha = gamma/d.dot(q) ## positive definite self.x += alpha*d ## update if(ii%self.reset==0)or(np.signbit(np.real(alpha))): r = self.b-self.A(self.x) else: r -= alpha*q s = self.W(r) gamma_ = gamma gamma = r.dot(s) beta = max(0,gamma/gamma_) ## positive definite d = s+beta*d ## conjugated gradient delta = delta_*np.absolute(gamma)**0.5 self.note.cflush("\niteration : %08u alpha = %3.1E beta = %3.1E delta = %3.1E"%(ii,alpha,beta,delta)) if(ii==limii): self.note.cprint("\n... quit.") break if(gamma==0): convergence = clevel self.note.cprint(" convergence level : INF\n... done.") break elif(np.absolute(delta)`_ Examples -------- >>> def egg(x): ... E = 0.5*x.dot(x) # energy E(x) -- a two-dimensional parabola ... g = x # gradient ... return E,g >>> x = field(point_space(2), val=[1, 3]) >>> x,convergence = steepest_descent(egg, note=True)(x0=x, tol=1E-4, clevel=3) iteration : 00000001 alpha = 1.0E+00 delta = 6.5E-01 iteration : 00000002 alpha = 2.0E+00 delta = 1.4E-01 iteration : 00000003 alpha = 1.6E-01 delta = 2.1E-03 iteration : 00000004 alpha = 2.6E-03 delta = 3.0E-04 iteration : 00000005 alpha = 2.0E-04 delta = 5.3E-05 convergence level : 1 iteration : 00000006 alpha = 8.2E-05 delta = 4.4E-06 convergence level : 2 iteration : 00000007 alpha = 6.6E-06 delta = 3.1E-06 convergence level : 3 ... done. >>> bool(convergence) True >>> x.val # approximately zero array([ -6.87299426e-07 -2.06189828e-06]) Attributes ---------- x : field Current field. eggs : function Given the current `x` it returns the tuple of energy and gradient. spam : function Callback function which is given the current `x` and iteration counter each iteration; can be ``None``. a : {4-tuple} Numbers obeying 0 < a1 ~ a2 < 1 ~ a3 < a4 that modify the step widths (default: (0.2,0.5,1,2)). c : {2-tuple} Numbers obeying 0 < c1 < c2 < 1 that specify the Wolfe-conditions (default: (1E-4,0.9)). note : notification Notification instance. """ def __init__(self,eggs,spam=None,a=(0.2,0.5,1,2),c=(1E-4,0.9),note=False): """ Initializes the steepest_descent and sets the attributes (except for `x`). Parameters ---------- eggs : function Given the current `x` it returns the tuple of energy and gradient. spam : function, *optional* Callback function which is given the current `x` and iteration counter each iteration (default: None). a : {4-tuple}, *optional* Numbers obeying 0 < a1 ~ a2 < 1 ~ a3 < a4 that modify the step widths (default: (0.2,0.5,1,2)). c : {2-tuple}, *optional* Numbers obeying 0 < c1 < c2 < 1 that specify the Wolfe-conditions (default: (1E-4,0.9)). note : bool, *optional* Indicates whether notes are printed or not (default: False). """ self.eggs = eggs ## returns energy and gradient self.spam = spam ## serves as callback given x and iteration number self.a = a ## 0 < a1 ~ a2 < 1 ~ a3 < a4 self.c = c ## 0 < c1 < c2 < 1 self.note = notification(default=bool(note)) ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def __call__(self,x0,alpha=1,tol=1E-4,clevel=8,limii=100000): """ Runs the steepest descent minimization. Parameters ---------- x0 : field Starting guess for the minimization. alpha : scalar, *optional* Starting step width to be multiplied with normalized gradient (default: 1). tol : scalar, *optional* Tolerance specifying convergence; measured by maximal change in `x` (default: 1E-4). clevel : integer, *optional* Number of times the tolerance should be undershot before exiting (default: 8). limii : integer, *optional* Maximum number of iterations performed (default: 100,000). Returns ------- x : field Latest `x` of the minimization. convergence : integer Latest convergence level indicating whether the minimization has converged or not. """ if(not isinstance(x0,field)): raise TypeError(about._errors.cstring("ERROR: invalid input.")) self.x = x0 E,g = self.eggs(self.x) ## energy and gradient norm = g.norm() ## gradient norm convergence = 0 ii = 1 while(True): x_,E,g,alpha,a = self._get_alpha(E,g,norm,alpha) ## "news",alpha,a if(alpha is None): self.note.cprint("\niteration : %08u alpha < 1.0E-13\n... dead."%ii) break else: delta = np.absolute(g.val).max()*(alpha/norm) self.note.cflush("\niteration : %08u alpha = %3.1E delta = %3.1E"%(ii,alpha,delta)) ## update self.x = x_ alpha *= a norm = g.norm() ## gradient norm if(ii==limii): self.note.cprint("\n... quit.") break elif(delta determines the new alpha while(True): ## Wolfe conditions wolfe,x_,E_,g_,a = self._check_wolfe(E,g,norm,alpha) # wolfe,x_,E_,g_,a = self._check_strong_wolfe(E,g,norm,alpha) if(wolfe): return x_,E_,g_,alpha,a else: alpha *= a if(alpha<1E-13): return None,None,None,None,None ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def _check_wolfe(self,E,g,norm,alpha): ## > checks the Wolfe conditions x_ = self._get_x(g,norm,alpha) pg = norm E_,g_ = self.eggs(x_) if(E_>E+self.c[0]*alpha*pg): if(E_ checks the strong Wolfe conditions # # x_ = self._get_x(g,norm,alpha) # pg = norm # E_,g_ = self.eggs(x_) # if(E_>E+self.c[0]*alpha*pg): # if(E_self.c[1]*np.absolute(pg)): # return True,x_,E_,g_,self.a[3] # return True,x_,E_,g_,self.a[2] ##+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ def _get_x(self,g,norm,alpha): ## > updates x return self.x-g*(alpha/norm) ##=============================================================================