From 67b067d93897620b47a329e483424eb07d0fa69a Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 17 Jan 2014 14:03:08 -0800 Subject: [PATCH] Major Reorganization (Things are likely still broken...) Added Parameter to Utils, which hints at where we are going with functions as parameters. --- SimPEG/Data.py | 39 +- SimPEG/Examples/DC.py | 4 +- SimPEG/Examples/Linear.py | 2 +- SimPEG/Inverse/BetaSchedule.py | 12 - SimPEG/Inverse/Inversion.py | 383 ------------------ SimPEG/Inverse/__init__.py | 4 - SimPEG/Inversion.py | 188 +++++++++ SimPEG/Model.py | 52 ++- SimPEG/ObjFunction.py | 260 ++++++++++++ .../{Inverse/Optimize.py => Optimization.py} | 14 +- SimPEG/{Inverse => }/Regularization.py | 2 +- SimPEG/Tests/test_model.py | 4 +- SimPEG/Tests/test_optimizers.py | 12 +- SimPEG/Tests/test_problem.py | 2 +- SimPEG/Utils/__init__.py | 62 +++ SimPEG/Utils/interputils.py | 4 +- SimPEG/__init__.py | 4 +- 17 files changed, 575 insertions(+), 473 deletions(-) delete mode 100644 SimPEG/Inverse/BetaSchedule.py delete mode 100644 SimPEG/Inverse/Inversion.py delete mode 100644 SimPEG/Inverse/__init__.py create mode 100644 SimPEG/Inversion.py create mode 100644 SimPEG/ObjFunction.py rename SimPEG/{Inverse/Optimize.py => Optimization.py} (98%) rename SimPEG/{Inverse => }/Regularization.py (99%) diff --git a/SimPEG/Data.py b/SimPEG/Data.py index dc9239e1..17eebaeb 100644 --- a/SimPEG/Data.py +++ b/SimPEG/Data.py @@ -1,37 +1,6 @@ import Utils import numpy as np -def requiresProblem(f): - """ - Use this to wrap a funciton:: - - @requiresProblem - def dpred(self): - pass - - This wrapper will ensure that a problem has been bound to the data. - If a problem is not bound an Exception will be raised, and an nice error message printed. - """ - extra = """ - To use data.%s(), SimPEG requires that a problem be bound to the data. - If a problem has not been bound, an Exception will be raised. - To bind a problem to the Data object:: - - data.setProblem(myProblem) - - """ % f.__name__ - from functools import wraps - @wraps(f) - def requiresProblemWrapper(self,*args,**kwargs): - if getattr(self, 'prob', None) is None: - raise Exception(extra) - return f(self,*args,**kwargs) - - doc = requiresProblemWrapper.__doc__ - requiresProblemWrapper.__doc__ = ('' if doc is None else doc) + extra - - return requiresProblemWrapper - class BaseData(object): """Data holds the observed data, and the standard deviations.""" @@ -55,7 +24,7 @@ class BaseData(object): prob.data = self @Utils.count - @requiresProblem + @Utils.requires('prob') def dpred(self, m, u=None): """ Create the projected data from a model. @@ -68,7 +37,7 @@ class BaseData(object): Where P is a projection of the fields onto the data space. """ if u is None: u = self.prob.field(m) - return self.projectField(u) + return Utils.mkvc(self.projectField(u)) @Utils.count @@ -99,7 +68,7 @@ class BaseData(object): \mu_\\text{data} = \mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs} """ - return self.dpred(m, u=u) - self.dobs + return Utils.mkvc(self.dpred(m, u=u) - self.dobs) @property @@ -136,7 +105,7 @@ class BaseData(object): Where W_d is a covariance matrix that weights the data residual. """ - return self.Wd*self.residual(m, u=u) + return Utils.mkvc(self.Wd*self.residual(m, u=u)) @property def RHS(self): diff --git a/SimPEG/Examples/DC.py b/SimPEG/Examples/DC.py index a05e3e3c..28b1466f 100644 --- a/SimPEG/Examples/DC.py +++ b/SimPEG/Examples/DC.py @@ -193,7 +193,7 @@ if __name__ == '__main__': p0 = [5, 10] p1 = [15, 50] condVals = [sig1, sig2] - mSynth = Utils.ModelBuilder.defineBlockConductivity(p0,p1,M.gridCC,condVals) + mSynth = Utils.ModelBuilder.defineBlockConductivity(M.gridCC,p0,p1,condVals) plt.colorbar(M.plotImage(mSynth)) # plt.show() @@ -217,7 +217,7 @@ if __name__ == '__main__': u = prob.field(mSynth) u = data.reshapeFields(u) M.plotImage(u[:,10]) - # plt.show() + plt.show() # Now set up the prob to do some minimization # prob.dobs = dobs diff --git a/SimPEG/Examples/Linear.py b/SimPEG/Examples/Linear.py index 618ff260..acaaf38b 100644 --- a/SimPEG/Examples/Linear.py +++ b/SimPEG/Examples/Linear.py @@ -1,4 +1,4 @@ -from SimPEG import Mesh, Model, Problem, Data, Inverse, np +from SimPEG import Mesh, Model, Problem, Data, Inversion, np import matplotlib.pyplot as plt diff --git a/SimPEG/Inverse/BetaSchedule.py b/SimPEG/Inverse/BetaSchedule.py deleted file mode 100644 index c225c0c0..00000000 --- a/SimPEG/Inverse/BetaSchedule.py +++ /dev/null @@ -1,12 +0,0 @@ - - -class Cooling(object): - """Simple Beta Schedule""" - - beta0 = None #: The initial beta value, set to none means that it will be approximated in the first iteration. - beta_coolingFactor = 2. - - def getBeta(self): - if self._beta is None: - return self.beta0 - return self._beta / self.beta_coolingFactor diff --git a/SimPEG/Inverse/Inversion.py b/SimPEG/Inverse/Inversion.py deleted file mode 100644 index 2e7d7f50..00000000 --- a/SimPEG/Inverse/Inversion.py +++ /dev/null @@ -1,383 +0,0 @@ -import SimPEG -from SimPEG import Utils, sp, np -from Optimize import Remember -from BetaSchedule import Cooling -from SimPEG.Inverse import IterationPrinters, StoppingCriteria - -class BaseInversion(object): - """BaseInversion(prob, reg, opt, data, **kwargs) - """ - - __metaclass__ = Utils.Save.Savable - - maxIter = 1 #: Maximum number of iterations - name = 'BaseInversion' - - debug = False #: Print debugging information - - comment = '' #: Used by some functions to indicate what is going on in the algorithm - counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things - - beta0 = None #: The initial Beta (regularization parameter) - beta0_ratio = 0.1 #: When beta0 is set to None, estimateBeta0 is used with this ratio - - def __init__(self, prob, reg, opt, data, **kwargs): - Utils.setKwargs(self, **kwargs) - self.prob = prob - self.reg = reg - self.opt = opt - self.data = data - self.opt.parent = self - - self.stoppers = [StoppingCriteria.iteration] - - # Check if we have inserted printers into the optimization - if IterationPrinters.phi_d not in self.opt.printers: - self.opt.printers.insert(1,IterationPrinters.beta) - self.opt.printers.insert(2,IterationPrinters.phi_d) - self.opt.printers.insert(3,IterationPrinters.phi_m) - - if not hasattr(opt, '_bfgsH0') and hasattr(opt, 'bfgsH0'): # Check if it has been set by the user and the default is not being used. - print 'Setting bfgsH0 to the inverse of the modelObj2Deriv. Done using direct methods.' - opt.bfgsH0 = SimPEG.Solver(reg.modelObj2Deriv()) - - - @property - def phi_d_target(self): - """ - target for phi_d - - By default this is the number of data. - - Note that we do not set the target if it is None, but we return the default value. - """ - if getattr(self, '_phi_d_target', None) is None: - return self.data.dobs.size # - return self._phi_d_target - - @phi_d_target.setter - def phi_d_target(self, value): - self._phi_d_target = value - - @Utils.timeIt - def run(self, m0): - """run(m0) - - Runs the inversion! - - """ - self.startup(m0) - while True: - self.doStartIteration() - self.m = self.opt.minimize(self.evalFunction, self.m) - self.doEndIteration() - if self.stoppingCriteria(): break - - self.printDone() - self.finish() - - return self.m - - @Utils.callHooks('startup') - def startup(self, m0): - """ - **startup** is called at the start of any new run call. - - :param numpy.ndarray x0: initial x - :rtype: None - :return: None - """ - - if not hasattr(self.reg, '_mref'): - print 'Regularization has not set mref. SimPEG will set it to m0.' - self.reg.mref = m0 - - self.m = m0 - self._iter = 0 - self._beta = None - self.phi_d_last = np.nan - self.phi_m_last = np.nan - - @Utils.callHooks('doStartIteration') - def doStartIteration(self): - """ - **doStartIteration** is called at the end of each run iteration. - - :rtype: None - :return: None - """ - self._beta = self.getBeta() - - - @Utils.callHooks('doEndIteration') - def doEndIteration(self): - """ - **doEndIteration** is called at the end of each run iteration. - - :rtype: None - :return: None - """ - # store old values - self.phi_d_last = self.phi_d - self.phi_m_last = self.phi_m - self._iter += 1 - - def getBeta(self): - return self.beta0 - - def estimateBeta0(self, u=None, ratio=0.1): - """estimateBeta0(u=None, ratio=0.1) - - 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}} - - - :param numpy.array u: fields - :param float ratio: desired ratio of the eigenvalues, default is 0.1 - :rtype: float - :return: beta0 - """ - if u is None: - u = self.prob.field(self.m) - - x0 = np.random.rand(*self.m.shape) - t = x0.dot(self.dataObj2Deriv(self.m,x0,u=u)) - b = x0.dot(self.reg.modelObj2Deriv()*x0) - return ratio*(t/b) - - def stoppingCriteria(self): - if self.debug: print 'checking stoppingCriteria' - return Utils.checkStoppers(self, self.stoppers) - - - def printDone(self): - """ - **printDone** is called at the end of the inversion routine. - - """ - Utils.printStoppers(self, self.stoppers) - - @Utils.callHooks('finish') - def finish(self): - """finish() - - **finish** is called at the end of the optimization. - """ - pass - - @Utils.timeIt - def evalFunction(self, m, return_g=True, return_H=True): - """evalFunction(m, return_g=True, return_H=True) - - - """ - - u = self.prob.field(m) - - if self._iter is 0 and self._beta is None: - self._beta = self.beta0 = self.estimateBeta0(u=u,ratio=self.beta0_ratio) - - phi_d = self.dataObj(m, u) - phi_m = self.reg.modelObj(m) - - self.dpred = self.data.dpred(m, u=u) # This is a cheap matrix vector calculation. - self.phi_d = phi_d - self.phi_m = phi_m - - f = phi_d + self._beta * phi_m - - out = (f,) - if return_g: - phi_dDeriv = self.dataObjDeriv(m, u=u) - phi_mDeriv = self.reg.modelObjDeriv(m) - - g = phi_dDeriv + self._beta * phi_mDeriv - out += (g,) - - if return_H: - def H_fun(v): - phi_d2Deriv = self.dataObj2Deriv(m, v, u=u) - phi_m2Deriv = self.reg.modelObj2Deriv()*v - - return phi_d2Deriv + self._beta * phi_m2Deriv - - operator = sp.linalg.LinearOperator( (m.size, m.size), H_fun, dtype=m.dtype ) - out += (operator,) - return out if len(out) > 1 else out[0] - - @Utils.timeIt - def dataObj(self, m, u=None): - """dataObj(m, u=None) - - :param numpy.array m: geophysical model - :param numpy.array u: fields - :rtype: float - :return: data misfit - - The data misfit using an l_2 norm is: - - .. math:: - - \mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2 - - Where P is a projection matrix that brings the field on the full domain to the data measurement locations; - u is the field of interest; d_obs is the observed data; and W is the weighting matrix. - """ - # TODO: ensure that this is a data is vector and Wd is a matrix. - R = self.data.residualWeighted(m, u=u) - R = Utils.mkvc(R) - return 0.5*np.vdot(R, R) - - @Utils.timeIt - def dataObjDeriv(self, m, u=None): - """dataObjDeriv(m, u=None) - - :param numpy.array m: geophysical model - :param numpy.array u: fields - :rtype: numpy.array - :return: data misfit derivative - - The data misfit using an l_2 norm is: - - .. math:: - - \mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2 - - If the field, u, is provided, the calculation of the data is fast: - - .. math:: - - \mathbf{d}_\\text{pred} = \mathbf{Pu(m)} - - \mathbf{R} = \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) - - Where P is a projection matrix that brings the field on the full domain to the data measurement locations; - u is the field of interest; d_obs is the observed data; and W is the weighting matrix. - - The derivative of this, with respect to the model, is: - - .. math:: - - \\frac{\partial \mu_\\text{data}}{\partial \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ R} - - """ - if u is None: - u = self.prob.field(m) - - R = self.data.residualWeighted(m, u=u) - - dmisfit = self.prob.Jt(m, self.data.Wd * R, u=u) - - return dmisfit - - @Utils.timeIt - def dataObj2Deriv(self, m, v, u=None): - """dataObj2Deriv(m, v, u=None) - - :param numpy.array m: geophysical model - :param numpy.array v: vector to multiply - :param numpy.array u: fields - :rtype: numpy.array - :return: data misfit derivative - - The data misfit using an l_2 norm is: - - .. math:: - - \mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2 - - If the field, u, is provided, the calculation of the data is fast: - - .. math:: - - \mathbf{d}_\\text{pred} = \mathbf{Pu(m)} - - \mathbf{R} = \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) - - Where P is a projection matrix that brings the field on the full domain to the data measurement locations; - u is the field of interest; d_obs is the observed data; and W is the weighting matrix. - - The derivative of this, with respect to the model, is: - - .. math:: - - \\frac{\partial \mu_\\text{data}}{\partial \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ R} - - \\frac{\partial^2 \mu_\\text{data}}{\partial^2 \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ W J} - - """ - if u is None: - u = self.prob.field(m) - - R = self.data.residualWeighted(m, u=u) - - # TODO: abstract to different norms a little cleaner. - # \/ it goes here. in l2 it is the identity. - dmisfit = self.prob.Jt_approx(m, self.data.Wd * self.data.Wd * self.prob.J_approx(m, v, u=u), u=u) - - return dmisfit - - def save(self, group): - group.attrs['phi_d'] = self.phi_d - group.attrs['phi_m'] = self.phi_m - group.setArray('m', self.m) - group.setArray('dpred', self.dpred) - -class Inversion(Cooling, Remember, BaseInversion): - - maxIter = 10 - name = "SimPEG Inversion" - - def __init__(self, prob, reg, opt, data, **kwargs): - BaseInversion.__init__(self, prob, reg, opt, data, **kwargs) - - self.stoppers.append(StoppingCriteria.phi_d_target_Inversion) - - if StoppingCriteria.phi_d_target_Minimize not in self.opt.stoppers: - self.opt.stoppers.append(StoppingCriteria.phi_d_target_Minimize) - -class TimeSteppingInversion(Remember, BaseInversion): - """ - A slightly different view on regularization parameters, - let Beta be viewed as 1/dt, and timestep by updating the - reference model every optimization iteration. - """ - maxIter = 1 - name = "Time-Stepping SimPEG Inversion" - - def __init__(self, prob, reg, opt, data, **kwargs): - BaseInversion.__init__(self, prob, reg, opt, data, **kwargs) - - self.stoppers.append(StoppingCriteria.phi_d_target_Inversion) - - if StoppingCriteria.phi_d_target_Minimize not in self.opt.stoppers: - self.opt.stoppers.append(StoppingCriteria.phi_d_target_Minimize) - - def _startup_TimeSteppingInversion(self, m0): - - def _doEndIteration_updateMref(self, xt): - if self.debug: 'Updating the reference model.' - self.parent.reg.mref = self.xc - - self.opt.hook(_doEndIteration_updateMref, overwrite=True) diff --git a/SimPEG/Inverse/__init__.py b/SimPEG/Inverse/__init__.py deleted file mode 100644 index d83a8269..00000000 --- a/SimPEG/Inverse/__init__.py +++ /dev/null @@ -1,4 +0,0 @@ -from Optimize import * -from Inversion import * -from Regularization import Regularization -import BetaSchedule diff --git a/SimPEG/Inversion.py b/SimPEG/Inversion.py new file mode 100644 index 00000000..63a42bb3 --- /dev/null +++ b/SimPEG/Inversion.py @@ -0,0 +1,188 @@ +import SimPEG +from SimPEG import Utils, sp, np +from Optimization import Remember, IterationPrinters, StoppingCriteria + +class BaseInversion(object): + """BaseInversion(prob, reg, opt, data, **kwargs) + """ + + __metaclass__ = Utils.Save.Savable + + maxIter = 1 #: Maximum number of iterations + name = 'BaseInversion' + + debug = False #: Print debugging information + + comment = '' #: Used by some functions to indicate what is going on in the algorithm + counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things + + def __init__(self, dataObj, opt, **kwargs): + Utils.setKwargs(self, **kwargs) + + self.dataObj = dataObj + self.dataObj.parent = self + + self.opt = opt + self.opt.parent = self + + self.stoppers = [StoppingCriteria.iteration] + + # Check if we have inserted printers into the optimization + if IterationPrinters.phi_d not in self.opt.printers: + self.opt.printers.insert(1,IterationPrinters.beta) + self.opt.printers.insert(2,IterationPrinters.phi_d) + self.opt.printers.insert(3,IterationPrinters.phi_m) + + if not hasattr(opt, '_bfgsH0') and hasattr(opt, 'bfgsH0'): # Check if it has been set by the user and the default is not being used. + #TODO: I don't think that this if statement is working... + print 'Setting bfgsH0 to the inverse of the modelObj2Deriv. Done using direct methods.' + opt.bfgsH0 = SimPEG.Solver(reg.modelObj2Deriv()) + + + @property + def phi_d_target(self): + """ + target for phi_d + + By default this is the number of data. + + Note that we do not set the target if it is None, but we return the default value. + """ + if getattr(self, '_phi_d_target', None) is None: + return self.data.dobs.size # + return self._phi_d_target + + @phi_d_target.setter + def phi_d_target(self, value): + self._phi_d_target = value + + @Utils.timeIt + def run(self, m0): + """run(m0) + + Runs the inversion! + + """ + self.startup(m0) + while True: + self.doStartIteration() + self.m = self.opt.minimize(self.evalFunction, self.m) + self.doEndIteration() + if self.stoppingCriteria(): break + + self.printDone() + self.finish() + + return self.m + + @Utils.callHooks('startup') + def startup(self, m0): + """ + **startup** is called at the start of any new run call. + + :param numpy.ndarray x0: initial x + :rtype: None + :return: None + """ + + if not hasattr(self.reg, '_mref'): + print 'Regularization has not set mref. SimPEG will set it to m0.' + self.reg.mref = m0 + + self.m = m0 + self._iter = 0 + self._beta = None + self.phi_d_last = np.nan + self.phi_m_last = np.nan + + @Utils.callHooks('doStartIteration') + def doStartIteration(self): + """ + **doStartIteration** is called at the end of each run iteration. + + :rtype: None + :return: None + """ + self._beta = self.getBeta() + + + @Utils.callHooks('doEndIteration') + def doEndIteration(self): + """ + **doEndIteration** is called at the end of each run iteration. + + :rtype: None + :return: None + """ + # store old values + self.phi_d_last = self.phi_d + self.phi_m_last = self.phi_m + self._iter += 1 + + + def stoppingCriteria(self): + if self.debug: print 'checking stoppingCriteria' + return Utils.checkStoppers(self, self.stoppers) + + + def printDone(self): + """ + **printDone** is called at the end of the inversion routine. + + """ + Utils.printStoppers(self, self.stoppers) + + @Utils.callHooks('finish') + def finish(self): + """finish() + + **finish** is called at the end of the optimization. + """ + pass + + + def save(self, group): + group.attrs['phi_d'] = self.phi_d + group.attrs['phi_m'] = self.phi_m + group.setArray('m', self.m) + group.setArray('dpred', self.dpred) + + + +# class Inversion(Cooling, Remember, BaseInversion): + +# maxIter = 10 +# name = "SimPEG Inversion" + +# def __init__(self, prob, reg, opt, data, **kwargs): +# BaseInversion.__init__(self, prob, reg, opt, data, **kwargs) + +# self.stoppers.append(StoppingCriteria.phi_d_target_Inversion) + +# if StoppingCriteria.phi_d_target_Minimize not in self.opt.stoppers: +# self.opt.stoppers.append(StoppingCriteria.phi_d_target_Minimize) + +# class TimeSteppingInversion(Remember, BaseInversion): +# """ +# A slightly different view on regularization parameters, +# let Beta be viewed as 1/dt, and timestep by updating the +# reference model every optimization iteration. +# """ +# maxIter = 1 +# name = "Time-Stepping SimPEG Inversion" + +# def __init__(self, prob, reg, opt, data, **kwargs): +# BaseInversion.__init__(self, prob, reg, opt, data, **kwargs) + +# self.stoppers.append(StoppingCriteria.phi_d_target_Inversion) + +# if StoppingCriteria.phi_d_target_Minimize not in self.opt.stoppers: +# self.opt.stoppers.append(StoppingCriteria.phi_d_target_Minimize) + +# def _startup_TimeSteppingInversion(self, m0): + +# def _doEndIteration_updateMref(self, xt): +# if self.debug: 'Updating the reference model.' +# self.parent.reg.mref = self.xc + +# self.opt.hook(_doEndIteration_updateMref, overwrite=True) diff --git a/SimPEG/Model.py b/SimPEG/Model.py index 75912deb..6eb3c297 100644 --- a/SimPEG/Model.py +++ b/SimPEG/Model.py @@ -2,14 +2,18 @@ from SimPEG import Utils, np, sp class BaseModel(object): - """SimPEG Model""" + """ + SimPEG Model + + """ __metaclass__ = Utils.Save.Savable counter = None #: A SimPEG.Utils.Counter object + mesh = None #: A SimPEG Mesh - def __init__(self): - pass + def __init__(self, mesh): + self.mesh = mesh def transform(self, m): """ @@ -19,13 +23,22 @@ class BaseModel(object): The *transform* changes the model into the physical property. - A common example of this is to invert for electrical conductivity - in log space. In this case, your model will be log(sigma) and to - get back to sigma, you can take the exponential: - """ return m + def transformInverse(self, D): + """ + :param numpy.array D: physical property + :rtype: numpy.array + :return: model + + The *transformInverse* changes the physical property into the model. + + .. note:: The *transformInverse* may not be easy to create in general. + + """ + raise NotImplementedError('The transformInverse is not implemented.') + def transformDeriv(self, m): """ :param numpy.array m: model @@ -37,16 +50,16 @@ class BaseModel(object): """ return sp.identity(m.size) - def example(self, mesh, type=None): - return np.random.rand(mesh.nC) + def example(self, modelType=None): + return np.random.rand(self.mesh.nC) class LogModel(BaseModel): """SimPEG LogModel""" - def __init__(self, **kwargs): - BaseModel.__init__(self, **kwargs) + def __init__(self, mesh, **kwargs): + BaseModel.__init__(self, mesh, **kwargs) def transform(self, m): """ @@ -68,6 +81,23 @@ class LogModel(BaseModel): """ return np.exp(Utils.mkvc(m)) + + def transformInverse(self, D): + """ + :param numpy.array D: physical property + :rtype: numpy.array + :return: model + + The *transformInverse* changes the physical property into the model. + + .. math:: + + m = \log{\sigma} + + """ + return np.log(Utils.mkvc(D)) + + def transformDeriv(self, m): """ :param numpy.array m: model diff --git a/SimPEG/ObjFunction.py b/SimPEG/ObjFunction.py new file mode 100644 index 00000000..38323cd3 --- /dev/null +++ b/SimPEG/ObjFunction.py @@ -0,0 +1,260 @@ +from SimPEG import Utils + +class BaseObjFunction(object): + """docstring for BaseObjFunction""" + + __metaclass__ = Utils.Save.Savable + + beta = None #: Regularization trade-off parameter + + name = 'BaseObjFunction' #: Name of the objective function + + counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things + + u_current = None #: The most current evaluated field + m_current = None #: The most current model + + @property + def parent(self): + """This is the parent of the objective function.""" + return getattr(self,'_parent',None) + @parent.setter + def parent(self, p): + if getattr(self,'_parent',None) is not None: + print 'Objective function has switched to a new parent!' + self._parent = p + + + def __init__(self, data, reg, **kwargs): + Utils.setKwargs(self, **kwargs) + + self.data = data + self.reg = reg + + @Utils.timeIt + def evalFunction(self, m, return_g=True, return_H=True): + """evalFunction(m, return_g=True, return_H=True) + """ + + self.u_current = None + self.m_current = m + + u = self.data.prob.field(m) + self.u_current = u + + if self._iter is 0 and self._beta is None: + self._beta = self.beta0 = self.estimateBeta0(u=u,ratio=self.beta0_ratio) + + phi_d = self.dataObj(m, u=u) + phi_m = self.reg.modelObj(m) + + self.dpred = self.data.dpred(m, u=u) # This is a cheap matrix vector calculation. + self.phi_d = phi_d + self.phi_m = phi_m + + f = phi_d + self._beta * phi_m + + out = (f,) + if return_g: + phi_dDeriv = self.dataObjDeriv(m, u=u) + phi_mDeriv = self.reg.modelObjDeriv(m) + + g = phi_dDeriv + self._beta * phi_mDeriv + out += (g,) + + if return_H: + def H_fun(v): + phi_d2Deriv = self.dataObj2Deriv(m, v, u=u) + phi_m2Deriv = self.reg.modelObj2Deriv()*v + + return phi_d2Deriv + self._beta * phi_m2Deriv + + operator = sp.linalg.LinearOperator( (m.size, m.size), H_fun, dtype=m.dtype ) + out += (operator,) + return out if len(out) > 1 else out[0] + + @Utils.timeIt + def dataObj(self, m, u=None): + """dataObj(m, u=None) + + :param numpy.array m: geophysical model + :param numpy.array u: fields + :rtype: float + :return: data misfit + + The data misfit using an l_2 norm is: + + .. math:: + + \mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2 + + Where P is a projection matrix that brings the field on the full domain to the data measurement locations; + u is the field of interest; d_obs is the observed data; and W is the weighting matrix. + """ + # TODO: ensure that this is a data is vector and Wd is a matrix. + R = self.data.residualWeighted(m, u=u) + return 0.5*np.vdot(R, R) + + @Utils.timeIt + def dataObjDeriv(self, m, u=None): + """dataObjDeriv(m, u=None) + + :param numpy.array m: geophysical model + :param numpy.array u: fields + :rtype: numpy.array + :return: data misfit derivative + + The data misfit using an l_2 norm is: + + .. math:: + + \mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2 + + If the field, u, is provided, the calculation of the data is fast: + + .. math:: + + \mathbf{d}_\\text{pred} = \mathbf{Pu(m)} + + \mathbf{R} = \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) + + Where P is a projection matrix that brings the field on the full domain to the data measurement locations; + u is the field of interest; d_obs is the observed data; and W is the weighting matrix. + + The derivative of this, with respect to the model, is: + + .. math:: + + \\frac{\partial \mu_\\text{data}}{\partial \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ R} + + """ + if u is None: u = self.data.prob.field(m) + + R = self.data.residualWeighted(m, u=u) + + dmisfit = self.data.prob.Jt(m, self.data.Wd * R, u=u) + + return dmisfit + + @Utils.timeIt + def dataObj2Deriv(self, m, v, u=None): + """dataObj2Deriv(m, v, u=None) + + :param numpy.array m: geophysical model + :param numpy.array v: vector to multiply + :param numpy.array u: fields + :rtype: numpy.array + :return: data misfit derivative + + The data misfit using an l_2 norm is: + + .. math:: + + \mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2 + + If the field, u, is provided, the calculation of the data is fast: + + .. math:: + + \mathbf{d}_\\text{pred} = \mathbf{Pu(m)} + + \mathbf{R} = \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) + + Where P is a projection matrix that brings the field on the full domain to the data measurement locations; + u is the field of interest; d_obs is the observed data; and W is the weighting matrix. + + The derivative of this, with respect to the model, is: + + .. math:: + + \\frac{\partial \mu_\\text{data}}{\partial \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ R} + + \\frac{\partial^2 \mu_\\text{data}}{\partial^2 \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ W J} + + """ + if u is None: u = self.data.prob.field(m) + + R = self.data.residualWeighted(m, u=u) + + # TODO: abstract to different norms a little cleaner. + # \/ it goes here. in l2 it is the identity. + dmisfit = self.data.prob.Jt_approx(m, self.data.Wd * self.data.Wd * self.data.prob.J_approx(m, v, u=u), u=u) + + return dmisfit + + + +class BetaSchedule(Utils.Parameter): + """BetaSchedule""" + + beta0 = 'guess' #: The initial Beta (regularization parameter) + beta0_ratio = 0.1 #: When beta0 is set to 'guess', estimateBeta0 is used with this ratio + + coolingFactor = 2. + coolingRate = 3 + + beta = None #: Beta parameter + + def __init__(self, **kwargs): + Utils.setKwargs(self, **kwargs) + + def initialize(self): + self.beta = self.beta0 + + @Utils.requires('parent') + def get(self): + if self.beta is 'guess': + self.beta = self.estimateBeta0() + + invesion = self.parent.parent + if inversion._iter > 0 and inversion._iter % self.coolingRate == 0: + self.beta /= self.coolingFactor + + return self.beta + + @Utils.requires('parent') + def estimateBeta0(self, u=None): + """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}} + + + :param numpy.array u: fields + :param float ratio: desired ratio of the eigenvalues, default is 0.1 + :rtype: float + :return: beta0 + """ + objFunc = + data = objFunc.data + + m = invesion.m + + if u is None: + u = data.prob.field(m) + + x0 = np.random.rand(*m.shape) + t = x0.dot(objFunc.dataObj2Deriv(m,x0,u=u)) + b = x0.dot(objFunc.reg.modelObj2Deriv()*x0) + return self.beta0_ratio*(t/b) diff --git a/SimPEG/Inverse/Optimize.py b/SimPEG/Optimization.py similarity index 98% rename from SimPEG/Inverse/Optimize.py rename to SimPEG/Optimization.py index f11cd7b1..3157b5c1 100644 --- a/SimPEG/Inverse/Optimize.py +++ b/SimPEG/Optimization.py @@ -1,5 +1,4 @@ from SimPEG import Solver, Utils, sp, np -import matplotlib.pyplot as plt norm = np.linalg.norm @@ -101,6 +100,7 @@ class Minimize(object): comment = '' #: Used by some functions to indicate what is going on in the algorithm counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things + parent = None #: This is the parent of the optimization routine. def __init__(self, **kwargs): self.stoppers = [StoppingCriteria.tolerance_f, StoppingCriteria.moving_x, StoppingCriteria.tolerance_g, StoppingCriteria.norm_g, StoppingCriteria.iteration] @@ -179,16 +179,6 @@ class Minimize(object): return self.xc - @property - def parent(self): - """ - This is the parent of the optimization routine. - """ - return getattr(self, '_parent', None) - @parent.setter - def parent(self, value): - self._parent = value - @Utils.callHooks('startup') def startup(self, x0): """ @@ -873,7 +863,7 @@ class NewtonRoot(object): if __name__ == '__main__': - from SimPEG.tests import Rosenbrock, checkDerivative + from SimPEG.Tests import Rosenbrock, checkDerivative import matplotlib.pyplot as plt x0 = np.array([2.6, 3.7]) checkDerivative(Rosenbrock, x0, plotIt=False) diff --git a/SimPEG/Inverse/Regularization.py b/SimPEG/Regularization.py similarity index 99% rename from SimPEG/Inverse/Regularization.py rename to SimPEG/Regularization.py index 1569d767..e592d643 100644 --- a/SimPEG/Inverse/Regularization.py +++ b/SimPEG/Regularization.py @@ -1,6 +1,6 @@ from SimPEG import Utils, np, sp -class Regularization(object): +class BaseRegularization(object): """**Regularization** Here we will define regularization of a model, m, in general however, this should be thought of as (m-m_ref) but otherwise it is exactly the same: diff --git a/SimPEG/Tests/test_model.py b/SimPEG/Tests/test_model.py index 946224f1..0c3b2c54 100644 --- a/SimPEG/Tests/test_model.py +++ b/SimPEG/Tests/test_model.py @@ -18,8 +18,8 @@ class ModelTests(unittest.TestCase): print 'SimPEG.Model.BaseModel: Testing Model Transform' for M in dir(Model): if 'Model' not in M: continue - model = getattr(Model, M)() - m = model.example(self.mesh2) + model = getattr(Model, M)(self.mesh2) + m = model.example() passed = checkDerivative(lambda m : [model.transform(m), model.transformDeriv(m)], m, plotIt=False) self.assertTrue(passed) diff --git a/SimPEG/Tests/test_optimizers.py b/SimPEG/Tests/test_optimizers.py index 12c8c9b5..3132beb9 100644 --- a/SimPEG/Tests/test_optimizers.py +++ b/SimPEG/Tests/test_optimizers.py @@ -4,7 +4,7 @@ from SimPEG.Mesh import TensorMesh from SimPEG.Utils import sdiag import numpy as np import scipy.sparse as sp -from SimPEG import Inverse +from SimPEG import Optimization from SimPEG.Tests import getQuadratic, Rosenbrock TOL = 1e-2 @@ -16,7 +16,7 @@ class TestOptimizers(unittest.TestCase): self.b = np.array([-5,-5]) def test_GN_Rosenbrock(self): - GN = Inverse.GaussNewton() + GN = Optimization.GaussNewton() xopt = GN.minimize(Rosenbrock,np.array([0,0])) x_true = np.array([1.,1.]) print 'xopt: ', xopt @@ -24,7 +24,7 @@ class TestOptimizers(unittest.TestCase): self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True) def test_GN_quadratic(self): - GN = Inverse.GaussNewton() + GN = Optimization.GaussNewton() xopt = GN.minimize(getQuadratic(self.A,self.b),np.array([0,0])) x_true = np.array([5.,5.]) print 'xopt: ', xopt @@ -32,7 +32,7 @@ class TestOptimizers(unittest.TestCase): self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True) def test_ProjGradient_quadraticBounded(self): - PG = Inverse.ProjectedGradient(debug=True) + PG = Optimization.ProjectedGradient(debug=True) PG.lower, PG.upper = -2, 2 xopt = PG.minimize(getQuadratic(self.A,self.b),np.array([0,0])) x_true = np.array([2.,2.]) @@ -42,7 +42,7 @@ class TestOptimizers(unittest.TestCase): def test_ProjGradient_quadratic1Bound(self): myB = np.array([-5,1]) - PG = Inverse.ProjectedGradient() + PG = Optimization.ProjectedGradient() PG.lower, PG.upper = -2, 2 xopt = PG.minimize(getQuadratic(self.A,myB),np.array([0,0])) x_true = np.array([2.,-1.]) @@ -53,7 +53,7 @@ class TestOptimizers(unittest.TestCase): def test_NewtonRoot(self): fun = lambda x, return_g=True: np.sin(x) if not return_g else ( np.sin(x), sdiag( np.cos(x) ) ) x = np.array([np.pi-0.3, np.pi+0.1, 0]) - xopt = Inverse.NewtonRoot(comments=False).root(fun,x) + xopt = Optimization.NewtonRoot(comments=False).root(fun,x) x_true = np.array([np.pi,np.pi,0]) print 'Newton Root Finding' print 'xopt: ', xopt diff --git a/SimPEG/Tests/test_problem.py b/SimPEG/Tests/test_problem.py index fa94e6b8..e1d30337 100644 --- a/SimPEG/Tests/test_problem.py +++ b/SimPEG/Tests/test_problem.py @@ -14,7 +14,7 @@ class ProblemTests(unittest.TestCase): c = np.array([1, 4]) self.mesh2 = Mesh.TensorMesh([a, b], np.array([3, 5])) self.p2 = Problem.BaseProblem(self.mesh2, None) - self.reg = Inverse.Regularization(self.mesh2) + self.reg = Regularization.BaseRegularization(self.mesh2) def test_regularization(self): derChk = lambda m: [self.reg.modelObj(m), self.reg.modelObjDeriv(m)] diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index b48a856e..b9094a42 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -139,6 +139,42 @@ def dependentProperty(name, value, children, doc): setattr(self, name, val) return property(fget=fget, fset=fset, doc=doc) +def requires(var): + """ + Use this to wrap a funciton:: + + @requires('prob') + def dpred(self): + pass + + This wrapper will ensure that a problem has been bound to the data. + If a problem is not bound an Exception will be raised, and an nice error message printed. + """ + def requiresVar(f): + if var is 'prob': + extra = """ + To use data.%s(), SimPEG requires that a problem be bound to the data. + If a problem has not been bound, an Exception will be raised. + To bind a problem to the Data object:: + + data.setProblem(myProblem) + + """ % f.__name__ + else: + extra = """ + To use *%s* method, SimPEG requires that the %s be specified. + """ % (f.__name__, var) + @wraps(f) + def requiresVarWrapper(self,*args,**kwargs): + if getattr(self, var, None) is None: + raise Exception(extra) + return f(self,*args,**kwargs) + + doc = requiresVarWrapper.__doc__ + requiresVarWrapper.__doc__ = ('' if doc is None else doc) + extra + + return requiresVarWrapper + return requiresVar class Counter(object): """ @@ -232,6 +268,32 @@ def timeIt(f): return wrapper +class Parameter(object): + """Parameter""" + + debug = False #: Print debugging information + + def __init__(self, *args, **kwargs): + pass + + @property + def parent(self): + """This is the parent of the Parameter instance.""" + return getattr(self,'_parent',None) + @parent.setter + def parent(self, p): + if getattr(self,'_parent',None) is not None: + print 'Parameter has switched to a new parent!' + self._parent = p + + def initialize(self): + pass + + def get(self): + raise NotImplementedError('Getting the Parameter is not yet implemented.') + + + if __name__ == '__main__': class MyClass(object): def __init__(self, url): diff --git a/SimPEG/Utils/interputils.py b/SimPEG/Utils/interputils.py index 43cbff6b..38308492 100644 --- a/SimPEG/Utils/interputils.py +++ b/SimPEG/Utils/interputils.py @@ -14,10 +14,10 @@ def _interp_point_1D(x, xr_i): """ # TODO: This fails if the point is on the outside of the mesh. We may want to replace this by extrapolation? im = np.argmin(abs(x-xr_i)) - if xr_i - x[im] >= 0: # Point on the left + if xr_i - x[im] >= 0: # Point on the left ind_x1 = im ind_x2 = im+1 - elif xr_i - x[im] < 0: # Point on the right + elif xr_i - x[im] < 0: # Point on the right ind_x1 = im-1 ind_x2 = im dx1 = xr_i - x[ind_x1] diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 904766d2..01607c80 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -6,7 +6,9 @@ import Mesh import Model import Problem import Data -import Inverse +import Inversion +import Optimization +import Regularization import Examples import Tests