diff --git a/SimPEG/Inversion.py b/SimPEG/Inversion.py index 0bef2d92..5b2dfa4d 100644 --- a/SimPEG/Inversion.py +++ b/SimPEG/Inversion.py @@ -4,7 +4,7 @@ from Optimization import Remember, IterationPrinters, StoppingCriteria class BaseInversion(object): - """BaseInversion(prob, reg, opt, data, **kwargs) + """BaseInversion(objFunc, opt, **kwargs) """ __metaclass__ = Utils.Save.Savable diff --git a/SimPEG/ObjFunction.py b/SimPEG/ObjFunction.py index d88d5740..fb2fe2f4 100644 --- a/SimPEG/ObjFunction.py +++ b/SimPEG/ObjFunction.py @@ -1,16 +1,16 @@ from SimPEG import Utils, np, sp class BaseObjFunction(object): - """docstring for BaseObjFunction""" + """BaseObjFunction(data, reg, **kwargs)""" __metaclass__ = Utils.Save.Savable - beta = Utils.ParameterProperty('beta', default=None, doc='Regularization trade-off parameter') + beta = Utils.ParameterProperty('beta', default=1, doc='Regularization trade-off parameter') debug = False #: Print debugging information counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things - name = 'BaseObjFunction' #: Name of the objective function + name = 'Base Objective Function' #: Name of the objective function u_current = None #: The most current evaluated field m_current = None #: The most current model @@ -25,12 +25,27 @@ class BaseObjFunction(object): print 'Objective function has switched to a new parent!' self._parent = p + @property + def inv(self): return self.parent + @property + def objFunc(self): return self + @property + def opt(self): return getattr(self.parent,'opt',None) + @property + def prob(self): return self.data.prob + @property + def mesh(self): return self.data.prob.mesh + @property + def model(self): return self.data.prob.model + def __init__(self, data, reg, **kwargs): Utils.setKwargs(self, **kwargs) self.data = data + self.reg = reg + self.reg.parent = self @Utils.callHooks('startup') @@ -41,7 +56,7 @@ class BaseObjFunction(object): """ if self.debug: print 'Calling ObjFunction.startup' - if not hasattr(self.reg, '_mref'): + if self.reg.mref is None: print 'Regularization has not set mref. SimPEG will set it to m0.' self.reg.mref = m0 @@ -226,8 +241,8 @@ class BetaSchedule(Utils.Parameter): self.beta = self.estimateBeta0() opt = self.parent.parent.opt - if opt._iter > 0 and opt._iter % self.coolingRate == 0: - if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % opt._iter + if opt.iter > 0 and opt.iter % self.coolingRate == 0: + if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % opt.iter self.beta /= self.coolingFactor return self.beta diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py index 35b71160..2362dd20 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -9,11 +9,11 @@ class StoppingCriteria(object): """docstring for StoppingCriteria""" iteration = { "str": "%d : maxIter = %3d <= iter = %3d", - "left": lambda M: M.maxIter, "right": lambda M: M._iter, + "left": lambda M: M.maxIter, "right": lambda M: M.iter, "stopType": "critical"} iterationLS = { "str": "%d : maxIterLS = %3d <= iterLS = %3d", - "left": lambda M: M.maxIterLS, "right": lambda M: M._iterLS, + "left": lambda M: M.maxIterLS, "right": lambda M: M.iterLS, "stopType": "critical"} armijoGoldstein = { "str": "%d : ft = %1.4e <= alp*descent = %1.4e", @@ -21,11 +21,11 @@ class StoppingCriteria(object): "stopType": "optimal"} tolerance_f = { "str": "%d : |fc-fOld| = %1.4e <= tolF*(1+|f0|) = %1.4e", - "left": lambda M: 1 if M._iter==0 else abs(M.f-M.f_last), "right": lambda M: 0 if M._iter==0 else M.tolF*(1+abs(M.f0)), + "left": lambda M: 1 if M.iter==0 else abs(M.f-M.f_last), "right": lambda M: 0 if M.iter==0 else M.tolF*(1+abs(M.f0)), "stopType": "optimal"} moving_x = { "str": "%d : |xc-x_last| = %1.4e <= tolX*(1+|x0|) = %1.4e", - "left": lambda M: 1 if M._iter==0 else norm(M.xc-M.x_last), "right": lambda M: 0 if M._iter==0 else M.tolX*(1+norm(M.x0)), + "left": lambda M: 1 if M.iter==0 else norm(M.xc-M.x_last), "right": lambda M: 0 if M.iter==0 else M.tolX*(1+norm(M.x0)), "stopType": "optimal"} tolerance_g = { "str": "%d : |proj(x-g)-x| = %1.4e <= tolG = %1.4e", @@ -56,12 +56,12 @@ class StoppingCriteria(object): class IterationPrinters(object): """docstring for IterationPrinters""" - iteration = {"title": "#", "value": lambda M: M._iter, "width": 5, "format": "%3d"} + iteration = {"title": "#", "value": lambda M: M.iter, "width": 5, "format": "%3d"} f = {"title": "f", "value": lambda M: M.f, "width": 10, "format": "%1.2e"} norm_g = {"title": "|proj(x-g)-x|", "value": lambda M: norm(M.projection(M.xc - M.g) - M.xc), "width": 15, "format": "%1.2e"} - totalLS = {"title": "LS", "value": lambda M: M._iterLS, "width": 5, "format": "%d"} + totalLS = {"title": "LS", "value": lambda M: M.iterLS, "width": 5, "format": "%d"} - iterationLS = {"title": "#", "value": lambda M: (M._iter, M._iterLS), "width": 5, "format": "%3d.%d"} + iterationLS = {"title": "#", "value": lambda M: (M.iter, M.iterLS), "width": 5, "format": "%3d.%d"} LS_ft = {"title": "ft", "value": lambda M: M._LS_ft, "width": 10, "format": "%1.2e"} LS_t = {"title": "t", "value": lambda M: M._LS_t, "width": 10, "format": "%0.5f"} LS_armijoGoldstein = {"title": "f + alp*g.T*p", "value": lambda M: M.f + M.LSreduction*M._LS_descent, "width": 16, "format": "%1.2e"} @@ -188,15 +188,15 @@ class Minimize(object): x0 = x0 xc = x0 - _iter = _iterLS = 0 + iter = iterLS = 0 :param numpy.ndarray x0: initial x :rtype: None :return: None """ - self._iter = 0 - self._iterLS = 0 + self.iter = 0 + self.iterLS = 0 x0 = self.projection(x0) # ensure that we start of feasible. self.x0 = x0 @@ -268,7 +268,7 @@ class Minimize(object): pass def stoppingCriteria(self, inLS=False): - if self._iter == 0: + if self.iter == 0: self.f0 = self.f self.g0 = self.g return Utils.checkStoppers(self, self.stoppers if not inLS else self.stoppersLS) @@ -360,21 +360,21 @@ class Minimize(object): """ # Projected Armijo linesearch self._LS_t = 1 - self._iterLS = 0 - while self._iterLS < self.maxIterLS: + self.iterLS = 0 + while self.iterLS < self.maxIterLS: self._LS_xt = self.projection(self.xc + self._LS_t*p) self._LS_ft = self.evalFunction(self._LS_xt, return_g=False, return_H=False) self._LS_descent = np.inner(self.g, self._LS_xt - self.xc) # this takes into account multiplying by t, but is important for projection. if self.stoppingCriteria(inLS=True): break - self._iterLS += 1 + self.iterLS += 1 self._LS_t = self.LSshorten*self._LS_t if self.debugLS: - if self._iterLS == 1: self.printInit(inLS=True) + if self.iterLS == 1: self.printInit(inLS=True) self.printIter(inLS=True) - if self.debugLS and self._iterLS > 0: self.printDone(inLS=True) + if self.debugLS and self.iterLS > 0: self.printDone(inLS=True) - return self._LS_xt, self._iterLS < self.maxIterLS + return self._LS_xt, self.iterLS < self.maxIterLS @Utils.count def modifySearchDirectionBreak(self, p): @@ -416,7 +416,7 @@ class Minimize(object): # store old values self.f_last = self.f self.x_last, self.xc = self.xc, xt - self._iter += 1 + self.iter += 1 if self.debug: self.printDone() @@ -613,7 +613,7 @@ class ProjectedGradient(Minimize, Remember): f_current_decrease = self.f_last - self.f self.comment = '' - if self._iter < 1: + if self.iter < 1: # Note that this is reset on every CG iteration. self.f_decrease_max = -np.inf else: @@ -684,7 +684,7 @@ class BFGS(Minimize, Remember): return self.bfgs(-self.g) def _doEndIteration_BFGS(self, xt): - if self._iter is 0: + if self.iter is 0: self.g_last = self.g return @@ -817,7 +817,7 @@ class NewtonRoot(object): """ if self.comments: print 'Newton Method:\n' - self._iter = 0 + self.iter = 0 while True: r, J = fun(x, return_g=True) @@ -851,10 +851,10 @@ class NewtonRoot(object): rt = fun(xt, return_g=False) x = xt - self._iter += 1 + self.iter += 1 if norm(rt) < self.tol: break - if self._iter > self.maxIter: + if self.iter > self.maxIter: print 'NewtonRoot stopped by maxIters. norm: %4.4e' % norm(rt) break diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index ea6cd753..d483c2e5 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -25,14 +25,30 @@ class BaseRegularization(object): self.mesh = mesh self.model = model + mref = Utils.ParameterProperty('mref', default=None, doc='Reference model.') + @property - def mref(self): - if getattr(self, '_mref', None) is None: - return np.zeros(self.model.nP); - return self._mref - @mref.setter - def mref(self, value): - self._mref = value + def parent(self): + """This is the parent of the regularization.""" + return getattr(self,'_parent',None) + @parent.setter + def parent(self, p): + if getattr(self,'_parent',None) is not None: + print 'Regularization has switched to a new parent!' + self._parent = p + + @property + def inv(self): return self.parent.inv + @property + def objFunc(self): return self.parent + @property + def reg(self): return self + @property + def opt(self): return self.parent.opt + @property + def prob(self): return self.parent.prob + @property + def data(self): return self.parent.data @property diff --git a/SimPEG/Utils/Save.py b/SimPEG/Utils/Save.py index 55de4b4b..a9c77191 100644 --- a/SimPEG/Utils/Save.py +++ b/SimPEG/Utils/Save.py @@ -57,7 +57,7 @@ class SimPEGTable: # At the start of every iteration we will create a inversion iteration node. def _doStartIteration_hdf5_inv(invObj): - invObj._invNodeIt = invObj._invNode.addGroup('%d'%(invObj._iter+1)) + invObj._invNodeIt = invObj._invNode.addGroup('%d'%(invObj.iter+1)) preIteration(invObj._invNodeIt) invObj.hook(_doStartIteration_hdf5_inv, overwrite=True) @@ -78,7 +78,7 @@ class SimPEGTable: invObj.hook(_finish_hdf5_inv, overwrite=True) def _doStartIteration_hdf5_opt(optObj): - optObj._optNodeIt = optObj.parent._invNode.addGroup('%d.%d'%(optObj.parent._iter, optObj._iter)) + optObj._optNodeIt = optObj.parent._invNode.addGroup('%d.%d'%(optObj.parent.iter, optObj.iter)) preIteration(optObj._optNodeIt) invObj.opt.hook(_doStartIteration_hdf5_opt, overwrite=True) diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 8c0f899b..8853500e 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -299,25 +299,30 @@ class Parameter(object): hook(self._parent, _startup_paramProperty, name=startupName, overwrite=True) @property - def opt(self): - return self.parent.parent.opt - + def inv(self): return self.parent.inv @property - def objFunc(self): - return self.parent - + def objFunc(self): return self.parent.objFunc @property - def reg(self): - return self.parent.reg + def opt(self): return self.parent.opt + @property + def reg(self): return self.parent.reg + @property + def data(self): return self.parent.data + @property + def prob(self): return self.parent.prob + @property + def model(self): return self.parent.model + @property + def mesh(self): return self.parent.mesh def initialize(self): pass def get(self): if (self.current is None or - not self.opt._iter == self.currentIter): + not self.opt.iter == self.currentIter): self.current = self.nextIter() - self.currentIter = self.opt._iter + self.currentIter = self.opt.iter return self.current def nextIter(self):