wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-08 04:50:44 +08:00
Code Issues Packages Projects Releases Wiki Activity

Parameters --> Rules

Browse Source 
See the Linear example for updates on how to migrate to this version.
This commit is contained in:
rowanc1
2014-05-14 14:30:33 -07:00
parent e4670e129d
commit ab249d31b3
9 changed files with 189 additions and 199 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Inversion.py
+18 -8
View File
@@ -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')
SimPEG/Maps.py
+1 -1
View File
@@ -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):
SimPEG/ObjFunction.py
+2 -2
View File
@@ -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
SimPEG/Optimization.py
+3 -2
View File
@@ -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)
SimPEG/Parameters.py
-177
View File
@@ -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
SimPEG/Regularization.py
+2 -3
View File
@@ -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."""
SimPEG/Rules.py
+158
View File
@@ -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
SimPEG/__init__.py
+1 -1
View File
@@ -9,8 +9,8 @@ import Survey
import Regularization
import ObjFunction
import Optimization
import Rules
import Inversion
import Parameters
import Tests
Tutorials/Linear.py
+4 -5
View File
@@ -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-')