diff --git a/SimPEG/Inversion.py b/SimPEG/Inversion.py index 8dbee17d..fc541e59 100644 --- a/SimPEG/Inversion.py +++ b/SimPEG/Inversion.py @@ -1,6 +1,7 @@ import SimPEG from SimPEG import Utils, sp, np from Optimization import Remember, IterationPrinters, StoppingCriteria +import Rules class BaseInversion(object): @@ -15,6 +16,18 @@ class BaseInversion(object): counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things + @property + def ruleList(self): + if getattr(self,'_ruleList', None) is None: + self._ruleList = Rules.RuleList(inversion=self) + return self._ruleList + + @ruleList.setter + def ruleList(self, value): + assert isinstance(value, Rules.RuleList), 'Must be a RuleList' + self._ruleList = value + self._ruleList.inversion = self + def __init__(self, objFunc, opt, **kwargs): Utils.setKwargs(self, **kwargs) @@ -22,6 +35,7 @@ class BaseInversion(object): self.objFunc.parent = self self.opt = opt + opt.callback = self._optCallback self.opt.parent = self self.stoppers = [StoppingCriteria.iteration] @@ -40,15 +54,11 @@ class BaseInversion(object): """ self.objFunc.startup(m0) + self.ruleList.call('initialize') self.m = self.opt.minimize(self.objFunc.evalFunction, m0) - self.finish() + self.ruleList.call('finish') return self.m - @Utils.callHooks('finish') - def finish(self): - """finish() - - **finish** is called at the end of the optimization. - """ - pass + def _optCallback(self, xt): + self.ruleList.call('endIter') diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index c828d2ea..882af375 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -1,4 +1,4 @@ -import Utils, Parameters, numpy as np, scipy.sparse as sp +import Utils, numpy as np, scipy.sparse as sp from Tests import checkDerivative class IdentityMap(object): diff --git a/SimPEG/ObjFunction.py b/SimPEG/ObjFunction.py index 8f8a3dc9..a5332921 100644 --- a/SimPEG/ObjFunction.py +++ b/SimPEG/ObjFunction.py @@ -1,11 +1,11 @@ -import Utils, Parameters, Survey, Problem, numpy as np, scipy.sparse as sp, gc +import Utils, Survey, Problem, numpy as np, scipy.sparse as sp, gc class BaseObjFunction(object): """BaseObjFunction(forward, reg, **kwargs)""" __metaclass__ = Utils.SimPEGMetaClass - beta = Parameters.ParameterProperty('beta', default=1, doc='Regularization trade-off parameter') + beta = 1.0 #: Regularization trade-off parameter debug = False #: Print debugging information counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py index 77b6772e..f9979730 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -424,14 +424,15 @@ class Minimize(object): :rtype: None :return: None """ - if self.callback is not None: - self.callback(xt) # store old values self.f_last = self.f self.x_last, self.xc = self.xc, xt self.iter += 1 if self.debug: self.printDone() + if self.callback is not None: + self.callback(xt) + def save(self, group): group.setArray('searchDirection', self.searchDirection) diff --git a/SimPEG/Parameters.py b/SimPEG/Parameters.py deleted file mode 100644 index e138ee17..00000000 --- a/SimPEG/Parameters.py +++ /dev/null @@ -1,177 +0,0 @@ -import Utils, numpy as np - - -class Parameter(object): - """Parameter""" - - debug = False #: Print debugging information - - current = None #: This hold - currentIter = 0 - - def __init__(self, **kwargs): - Utils.setKwargs(self, **kwargs) - - @property - def parent(self): - """This is the parent of the Parameter instance.""" - return getattr(self,'_parent',None) - @parent.setter - def parent(self, p): - startupName = '_startup_paramProperty_'+self._propertyName - if getattr(self,'_parent',None) is not None: - delattr(self._parent,startupName) - print 'Warning: Parameter %s has switched to a new parent.' % self._propertyName - if self.debug: print '%s function has been deleted' % startupName - self._parent = p - - prop = self - def _startup_paramProperty(self, *args): - if prop.debug: print 'initializing %s' % prop._propertyName - prop.initialize() - - Utils.hook(self._parent, _startup_paramProperty, name=startupName, overwrite=True) - - @property - def inv(self): return self.parent.inv - @property - def objFunc(self): return self.parent.objFunc - @property - def opt(self): return self.parent.opt - @property - def reg(self): return self.parent.reg - @property - def survey(self): return self.parent.survey - @property - def prob(self): return self.parent.prob - @property - def mapping(self): return self.parent.mapping - @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): - self.current = self.nextIter() - self.currentIter = getattr(self.opt, 'iter', 0) - return self.current - - def nextIter(self): - raise NotImplementedError('Getting the Parameter is not yet implemented.') - - -def ParameterProperty(name, default=None, doc=""): - def getter(self): - out = getattr(self,'_'+name,default) - if isinstance(out, Parameter): - out = out.get() - return out - def setter(self, value): - if isinstance(value, Parameter): - value._propertyName = name - value.parent = self - setattr(self, '_'+name, value) - - return property(fget=getter, fset=setter, doc=doc) - - -class BetaEstimate(Parameter): - """BetaEstimate""" - - beta0 = 'guess' #: The initial Beta (regularization parameter) - beta0_ratio = 0.1 #: When beta0 is set to 'guess', estimateBeta0 is used with this ratio - - beta = None #: Beta parameter - - def __init__(self, **kwargs): - Parameter.__init__(self, **kwargs) - - def initialize(self): - self.beta = self.beta0 - - @Utils.requires('parent') - def nextIter(self): - if self.beta is 'guess': - if self.debug: print 'BetaSchedule is estimating Beta0.' - self.beta = self.estimateBeta0() - return self.beta - - @Utils.requires('parent') - def estimateBeta0(self): - """estimateBeta0(u=None) - - The initial beta is calculated by comparing the estimated - eigenvalues of JtJ and WtW. - - To estimate the eigenvector of **A**, we will use one iteration - of the *Power Method*: - - .. math:: - - \mathbf{x_1 = A x_0} - - Given this (very course) approximation of the eigenvector, - we can use the *Rayleigh quotient* to approximate the largest eigenvalue. - - .. math:: - - \lambda_0 = \\frac{\mathbf{x^\\top A x}}{\mathbf{x^\\top x}} - - We will approximate the largest eigenvalue for both JtJ and WtW, and - use some ratio of the quotient to estimate beta0. - - .. math:: - - \\beta_0 = \gamma \\frac{\mathbf{x^\\top J^\\top J x}}{\mathbf{x^\\top W^\\top W x}} - - :rtype: float - :return: beta0 - """ - objFunc = self.parent - survey = objFunc.survey - - m = objFunc.m_current - u = objFunc.u_current - - if u is None: - u = survey.prob.fields(m) - - x0 = np.random.rand(*m.shape) - t = x0.dot(objFunc.dataObj2Deriv(m,x0,u=u)) - b = x0.dot(objFunc.reg.modelObj2Deriv(m, v=x0)) - return self.beta0_ratio*(t/b) - - -class BetaSchedule(BetaEstimate): - """BetaSchedule""" - - coolingFactor = 2. - coolingRate = 3 - - @Utils.requires('parent') - def nextIter(self): - if self.beta is 'guess': - if self.debug: print 'BetaSchedule is estimating Beta0.' - self.beta = self.estimateBeta0() - - if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0: - if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter - self.beta /= self.coolingFactor - - return self.beta - - -class UpdateReferenceModel(Parameter): - - mref0 = None - - def nextIter(self): - mref = getattr(self, 'm_prev', None) - if mref is None: - if self.debug: print 'UpdateReferenceModel is using mref0' - mref = self.mref0 - self.m_prev = self.objFunc.m_current - return mref diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 5318b7da..90dfe4d5 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -1,4 +1,4 @@ -import Utils, Maps, Mesh, Parameters, numpy as np, scipy.sparse as sp +import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp class BaseRegularization(object): """ @@ -20,6 +20,7 @@ class BaseRegularization(object): mapping = None #: A SimPEG.Map instance. mesh = None #: A SimPEG.Mesh instance. + mref = None #: Reference model. def __init__(self, mesh, mapping=None, **kwargs): Utils.setKwargs(self, **kwargs) @@ -28,8 +29,6 @@ class BaseRegularization(object): self.mapping = mapping or Maps.IdentityMap(mesh) self.mapping._assertMatchesPair(self.mapPair) - mref = Parameters.ParameterProperty('mref', default=None, doc='Reference model.') - @property def parent(self): """This is the parent of the regularization.""" diff --git a/SimPEG/Rules.py b/SimPEG/Rules.py new file mode 100644 index 00000000..aed5e9d2 --- /dev/null +++ b/SimPEG/Rules.py @@ -0,0 +1,158 @@ +import Utils, numpy as np + +class InversionRule(object): + """InversionRule""" + + debug = False #: Print debugging information + + current = None #: This hold + + def __init__(self, **kwargs): + Utils.setKwargs(self, **kwargs) + + @property + def inversion(self): + """This is the inversion of the InversionRule instance.""" + return getattr(self,'_inversion',None) + @inversion.setter + def inversion(self, i): + if getattr(self,'_inversion',None) is not None: + print 'Warning: InversionRule %s has switched to a new inversion.' % self.__name__ + self._inversion = i + + @property + def objFunc(self): return self.inversion.objFunc + @property + def opt(self): return self.inversion.opt + @property + def reg(self): return self.inversion.objFunc.reg + @property + def survey(self): return self.inversion.objFunc.survey + @property + def prob(self): return self.inversion.objFunc.prob + + def initialize(self): + pass + + def endIter(self): + pass + + def finish(self): + pass + +class RuleList(object): + + rList = None #: The list of Rules + + def __init__(self, *rules, **kwargs): + self.rList = [] + for r in rules: + assert isinstance(r, InversionRule), 'All rules must be InversionRules not %s' % r.__name__ + self.rList.append(r) + Utils.setKwargs(self, **kwargs) + + @property + def debug(self): + return getattr(self, '_debug', False) + @debug.setter + def debug(self, value): + for r in self.rList: + r.debug = value + self._debug = value + + @property + def inversion(self): + """This is the inversion of the InversionRule instance.""" + return getattr(self,'_inversion',None) + @inversion.setter + def inversion(self, i): + if self.inversion is i: return + if getattr(self,'_inversion',None) is not None: + print 'Warning: %s has switched to a new inversion.' % self.__name__ + for r in self.rList: + r.inversion = i + self._inversion = i + + def call(self, ruleType): + if self.rList is None: + if self.debug: 'RuleList is None, no rules to call!' + return + + rules = ['initialize', 'endIter', 'finish'] + assert ruleType in rules, 'Rule type must be in ["%s"]' % '", "'.join(rules) + for r in self.rList: + getattr(r, ruleType)() + + +class BetaEstimate_ByEig(InversionRule): + """BetaEstimate""" + + beta0 = None #: The initial Beta (regularization parameter) + beta0_ratio = 0.1 #: estimateBeta0 is used with this ratio + + def initialize(self): + """ + The initial beta is calculated by comparing the estimated + eigenvalues of JtJ and WtW. + + To estimate the eigenvector of **A**, we will use one iteration + of the *Power Method*: + + .. math:: + + \mathbf{x_1 = A x_0} + + Given this (very course) approximation of the eigenvector, + we can use the *Rayleigh quotient* to approximate the largest eigenvalue. + + .. math:: + + \lambda_0 = \\frac{\mathbf{x^\\top A x}}{\mathbf{x^\\top x}} + + We will approximate the largest eigenvalue for both JtJ and WtW, and + use some ratio of the quotient to estimate beta0. + + .. math:: + + \\beta_0 = \gamma \\frac{\mathbf{x^\\top J^\\top J x}}{\mathbf{x^\\top W^\\top W x}} + + :rtype: float + :return: beta0 + """ + + if self.debug: print 'Calculating the beta0 parameter.' + + m = self.objFunc.m_current + u = self.objFunc.u_current or self.prob.fields(m) + + x0 = np.random.rand(*m.shape) + t = x0.dot(self.objFunc.dataObj2Deriv(m,x0,u=u)) + b = x0.dot(self.reg.modelObj2Deriv(m, v=x0)) + self.beta0 = self.beta0_ratio*(t/b) + + self.objFunc.beta = self.beta0 + + +class BetaSchedule(InversionRule): + """BetaSchedule""" + + coolingFactor = 2. + coolingRate = 3 + + def endIter(self): + if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0: + if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter + self.objFunc.beta /= self.coolingFactor + + +# class UpdateReferenceModel(Parameter): + +# mref0 = None + +# def nextIter(self): +# mref = getattr(self, 'm_prev', None) +# if mref is None: +# if self.debug: print 'UpdateReferenceModel is using mref0' +# mref = self.mref0 +# self.m_prev = self.objFunc.m_current +# return mref diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 49061e09..aac0c7f8 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -9,8 +9,8 @@ import Survey import Regularization import ObjFunction import Optimization +import Rules import Inversion -import Parameters import Tests diff --git a/Tutorials/Linear.py b/Tutorials/Linear.py index 1f4bf94e..59f42cc0 100644 --- a/Tutorials/Linear.py +++ b/Tutorials/Linear.py @@ -43,8 +43,6 @@ def example(N): mtrue[mesh.vectorCCx > 0.45] = -0.5 mtrue[mesh.vectorCCx > 0.6] = 0 - - prob = LinearProblem(mesh, G) survey = prob.createSyntheticSurvey(mtrue, std=0.01) @@ -59,10 +57,12 @@ if __name__ == '__main__': M = prob.mesh reg = Regularization.Tikhonov(mesh) - beta = Parameters.BetaSchedule() - objFunc = ObjFunction.BaseObjFunction(survey, reg, beta=beta) + objFunc = ObjFunction.BaseObjFunction(survey, reg) opt = Optimization.InexactGaussNewton(maxIter=20) inv = Inversion.BaseInversion(objFunc, opt) + beta = Rules.BetaSchedule() + betaest = Rules.BetaEstimate_ByEig() + inv.ruleList = Rules.RuleList(betaest, beta) m0 = np.zeros_like(survey.mtrue) mrec = inv.run(m0) @@ -72,7 +72,6 @@ if __name__ == '__main__': plt.plot(prob.G[i,:]) plt.figure(2) - plt.plot(M.vectorCCx, survey.mtrue, 'b-') plt.plot(M.vectorCCx, mrec, 'r-')