diff --git a/README.md b/README.md index c2eafa26..37e0a743 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ ![SimPEG](https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png) -Simulation and Parameter Estimation in Geoscience - A python package for simulation and gradient based parameter estimation in the context of geophysical applications. +Simulation and Parameter Estimation in Geophysics - A python package for simulation and gradient based parameter estimation in the context of geophysical applications. The vision is to create a package for finite volume simulation with applications to geophysical imaging and subsurface flow. To enable the understanding of the many different components, this package has the following features: @@ -10,4 +10,20 @@ The vision is to create a package for finite volume simulation with applications * supports 1D, 2D and 3D problems * designed for large-scale inversions +Documentation: +[http://simpeg.readthedocs.org/en/latest/](http://simpeg.readthedocs.org/en/latest/) + +Code: +[https://github.com/simpeg/simpeg](https://github.com/simpeg/simpeg) + +Tests: +[https://travis-ci.org/simpeg/simpeg](https://travis-ci.org/simpeg/simpeg) + +Build Status: [![Build Status](https://travis-ci.org/simpeg/simpeg.png)](https://travis-ci.org/simpeg/simpeg) + +Bugs & Issues: +[https://github.com/simpeg/simpeg/issues](https://github.com/simpeg/simpeg/issues) + +Code Snippets & Tutorials: +[http://www.row1.ca/simpeg](http://www.row1.ca/simpeg) diff --git a/SimPEG/Data.py b/SimPEG/Data.py index 6c94811a..948b1c83 100644 --- a/SimPEG/Data.py +++ b/SimPEG/Data.py @@ -1,5 +1,4 @@ -import Utils -import numpy as np +import Utils, numpy as np class BaseData(object): diff --git a/SimPEG/Model.py b/SimPEG/Model.py index 8482edde..33cfeef8 100644 --- a/SimPEG/Model.py +++ b/SimPEG/Model.py @@ -1,4 +1,4 @@ -from SimPEG import Utils, np, sp +import Utils, Parameters, numpy as np, scipy.sparse as sp class BaseModel(object): diff --git a/SimPEG/ObjFunction.py b/SimPEG/ObjFunction.py index fb2fe2f4..ea2fe831 100644 --- a/SimPEG/ObjFunction.py +++ b/SimPEG/ObjFunction.py @@ -1,11 +1,11 @@ -from SimPEG import Utils, np, sp +import Utils, Parameters, numpy as np, scipy.sparse as sp class BaseObjFunction(object): """BaseObjFunction(data, reg, **kwargs)""" __metaclass__ = Utils.Save.Savable - beta = Utils.ParameterProperty('beta', default=1, doc='Regularization trade-off parameter') + beta = Parameters.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 @@ -213,84 +213,3 @@ class BaseObjFunction(object): 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, *args, **kwargs): - Utils.Parameter.__init__(self, *args, **kwargs) - Utils.setKwargs(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() - - 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 - self.beta /= self.coolingFactor - - 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}} - - - :param numpy.array u: fields - :param float ratio: desired ratio of the eigenvalues, default is 0.1 - :rtype: float - :return: beta0 - """ - objFunc = self.parent - data = objFunc.data - - m = objFunc.m_current - u = objFunc.u_current - - 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/Optimization.py b/SimPEG/Optimization.py index 2362dd20..4e5963f6 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -1,4 +1,5 @@ -from SimPEG import Solver, Utils, sp, np +import Utils, numpy as np, scipy.sparse as sp +from Solver import Solver norm = np.linalg.norm diff --git a/SimPEG/Parameters.py b/SimPEG/Parameters.py new file mode 100644 index 00000000..2a07aa69 --- /dev/null +++ b/SimPEG/Parameters.py @@ -0,0 +1,164 @@ +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 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): + self.current = self.nextIter() + self.currentIter = self.opt.iter + 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 + data = objFunc.data + + m = objFunc.m_current + u = objFunc.u_current + + 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) + + +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 diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index 501b1c45..e53c4675 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -1,6 +1,4 @@ -import Utils, Data -import scipy.sparse as sp -import numpy as np +import Utils, Data, numpy as np, scipy.sparse as sp class BaseProblem(object): diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index bc754592..aa62449a 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -1,4 +1,4 @@ -from SimPEG import Utils, Model, np, sp +import Utils, Model, Parameters, numpy as np, scipy.sparse as sp class BaseRegularization(object): """ @@ -23,7 +23,7 @@ class BaseRegularization(object): assert isinstance(model, self.modelPair), "Incorrect model for this regularization" self.model = model - mref = Utils.ParameterProperty('mref', default=None, doc='Reference model.') + mref = Parameters.ParameterProperty('mref', default=None, doc='Reference model.') @property def parent(self): diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 8853500e..437f1387 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -267,83 +267,6 @@ def timeIt(f): return wrapper -class Parameter(object): - """Parameter""" - - debug = False #: Print debugging information - - current = None #: This hold - currentIter = 0 - - 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): - 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() - - 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 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): - self.current = self.nextIter() - self.currentIter = self.opt.iter - 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) - - if __name__ == '__main__': class MyClass(object): def __init__(self, url): diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 10d6d4e8..4967608e 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -10,6 +10,7 @@ import Regularization import ObjFunction import Optimization import Inversion +import Parameters import Examples import Tests diff --git a/docs/api_Optimize.rst b/docs/api_Inverse.rst similarity index 99% rename from docs/api_Optimize.rst rename to docs/api_Inverse.rst index 462ae0be..21f4386e 100644 --- a/docs/api_Optimize.rst +++ b/docs/api_Inverse.rst @@ -1,6 +1,15 @@ .. _api_Inverse: +Regularization +************** + +.. automodule:: SimPEG.Regularization + :show-inheritance: + :members: + :undoc-members: + + Objective Function ****************** @@ -27,10 +36,3 @@ Inversion :undoc-members: -Regularization -************** - -.. automodule:: SimPEG.Regularization - :show-inheritance: - :members: - :undoc-members: diff --git a/docs/api_Parameters.rst b/docs/api_Parameters.rst new file mode 100644 index 00000000..e90b58a5 --- /dev/null +++ b/docs/api_Parameters.rst @@ -0,0 +1,11 @@ +.. _api_Parameters: + + +Parameters +********** + +.. automodule:: SimPEG.Parameters + :show-inheritance: + :members: + :undoc-members: + :inherited-members: diff --git a/docs/index.rst b/docs/index.rst index b2882cf0..5ba56947 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -3,7 +3,7 @@ :alt: SimPEG :align: center -SimPEG (Simulation and Parameter Estimation in Geoscience) is a python +SimPEG (Simulation and Parameter Estimation in Geophysics) is a python package for simulation and gradient based parameter estimation in the context of geoscience applications. @@ -16,11 +16,10 @@ these goals, this package has the following features: * provides a framework for geophysical and hydrogeologic problems * supports 1D, 2D and 3D problems - -.. raw:: html - - - +.. image:: simpeg-framework.png + :width: 400 px + :alt: Framework + :align: center Meshing & Operators =================== @@ -49,7 +48,8 @@ Inversion .. toctree:: :maxdepth: 2 - api_Optimize + api_Inverse + api_Parameters Testing SimPEG ============== diff --git a/docs/simpeg-framework.png b/docs/simpeg-framework.png new file mode 100644 index 00000000..ad8e5c1c Binary files /dev/null and b/docs/simpeg-framework.png differ