mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-13 17:45:30 +08:00
Rearrange Framework
This commit is contained in:
@@ -0,0 +1,140 @@
|
||||
import Utils, Survey, Problem, numpy as np, scipy.sparse as sp, gc
|
||||
|
||||
|
||||
def _splitForward(forward):
|
||||
assert forward.ispaired, 'The problem and survey must be paired.'
|
||||
if isinstance(forward, Survey.BaseSurvey):
|
||||
survey = forward
|
||||
prob = forward.prob
|
||||
elif isinstance(forward, Problem.BaseProblem):
|
||||
prob = forward
|
||||
survey = forward.survey
|
||||
else:
|
||||
raise Exception('The forward simulation must either be a problem or a survey.')
|
||||
|
||||
return prob, survey
|
||||
|
||||
|
||||
class BaseDataMisfit(object):
|
||||
"""BaseDataMisfit
|
||||
|
||||
.. note::
|
||||
|
||||
You should inherit from this class to create your own data misfit term.
|
||||
"""
|
||||
|
||||
__metaclass__ = Utils.SimPEGMetaClass
|
||||
|
||||
debug = False #: Print debugging information
|
||||
counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things
|
||||
|
||||
def __init__(self):
|
||||
pass
|
||||
|
||||
def splitForward(self, forward):
|
||||
"""splitForward(forward)
|
||||
|
||||
Split the forward simulation into a problem and a survey
|
||||
|
||||
:param Problem,Survey forward: forward simulation
|
||||
:rtype: Problem,Survey
|
||||
:return: (prob, survey)
|
||||
|
||||
"""
|
||||
prob, survey = _splitForward(forward)
|
||||
return prob, survey
|
||||
|
||||
@Utils.timeIt
|
||||
def dataObj(self, forward, m, u=None):
|
||||
"""dataObj(forward, m, u=None)
|
||||
|
||||
:param Problem,Survey forward: forward simulation
|
||||
:param numpy.array m: geophysical model
|
||||
:param numpy.array u: fields
|
||||
:rtype: float
|
||||
:return: data misfit
|
||||
|
||||
"""
|
||||
raise NotImplementedError('This method should be overwritten.')
|
||||
|
||||
@Utils.timeIt
|
||||
def dataObjDeriv(self, forward, m, u=None):
|
||||
"""dataObjDeriv(forward, m, u=None)
|
||||
|
||||
:param Problem,Survey forward: forward simulation
|
||||
:param numpy.array m: geophysical model
|
||||
:param numpy.array u: fields
|
||||
:rtype: numpy.array
|
||||
:return: data misfit derivative
|
||||
|
||||
"""
|
||||
raise NotImplementedError('This method should be overwritten.')
|
||||
|
||||
|
||||
@Utils.timeIt
|
||||
def dataObj2Deriv(self, forward, m, v, u=None):
|
||||
"""dataObj2Deriv(forward, m, v, u=None)
|
||||
|
||||
:param Problem,Survey forward: forward simulation
|
||||
: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
|
||||
|
||||
"""
|
||||
raise NotImplementedError('This method should be overwritten.')
|
||||
|
||||
|
||||
class l2_DataMisfit(object):
|
||||
"""
|
||||
|
||||
The data misfit with an l_2 norm:
|
||||
|
||||
.. math::
|
||||
|
||||
\mu_\\text{data} = {1\over 2}\left| \mathbf{W}_d (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self, **kwargs):
|
||||
pass
|
||||
|
||||
def getWd(self, survey):
|
||||
"""getWd(survey)
|
||||
|
||||
Get the data weighting matrix.
|
||||
|
||||
This is based on the norm of the data plus a noise floor.
|
||||
|
||||
:param Survey survey: geophysical survey
|
||||
:rtype: scipy.sparse.csr_matrix
|
||||
:return: Wd
|
||||
|
||||
"""
|
||||
eps = np.linalg.norm(Utils.mkvc(survey.dobs),2)*1e-5
|
||||
return Utils.sdiag(1/(abs(survey.dobs)*survey.std+eps))
|
||||
|
||||
@Utils.timeIt
|
||||
def dataObj(self, forward, m, u=None):
|
||||
"dataObj2Deriv(forward, m, u=None)"
|
||||
prob, survey = _splitForward(forward)
|
||||
Wd = self.getWd(survey)
|
||||
R = Wd * survey.residual(m, u=u)
|
||||
return 0.5*np.vdot(R, R)
|
||||
|
||||
@Utils.timeIt
|
||||
def dataObjDeriv(self, forward, m, u=None):
|
||||
"dataObj2Deriv(forward, m, u=None)"
|
||||
prob, survey = _splitForward(forward)
|
||||
if u is None: u = prob.fields(m)
|
||||
Wd = self.getWd(survey)
|
||||
return prob.Jtvec(m, Wd * (Wd * survey.residual(m, u=u)), u=u)
|
||||
|
||||
@Utils.timeIt
|
||||
def dataObj2Deriv(self, forward, m, v, u=None):
|
||||
"dataObj2Deriv(forward, m, v, u=None)"
|
||||
prob, survey = _splitForward(forward)
|
||||
if u is None: u = prob.fields(m)
|
||||
Wd = self.getWd(survey)
|
||||
return prob.Jtvec_approx(m, Wd * (Wd * prob.Jvec_approx(m, v, u=u)), u=u)
|
||||
+12
-10
@@ -21,15 +21,17 @@ class InversionDirective(object):
|
||||
self._inversion = i
|
||||
|
||||
@property
|
||||
def objFunc(self): return self.inversion.objFunc
|
||||
def invProb(self): return self.inversion.invProb
|
||||
@property
|
||||
def opt(self): return self.inversion.opt
|
||||
@property
|
||||
def reg(self): return self.inversion.objFunc.reg
|
||||
def reg(self): return self.invProb.reg
|
||||
@property
|
||||
def survey(self): return self.inversion.objFunc.survey
|
||||
def dmisfit(self): return self.invProb.dmisfit
|
||||
@property
|
||||
def prob(self): return self.inversion.objFunc.prob
|
||||
def survey(self): return self.invProb.survey
|
||||
@property
|
||||
def prob(self): return self.invProb.prob
|
||||
|
||||
def initialize(self):
|
||||
pass
|
||||
@@ -122,15 +124,15 @@ class BetaEstimate_ByEig(InversionDirective):
|
||||
|
||||
if self.debug: print 'Calculating the beta0 parameter.'
|
||||
|
||||
m = self.objFunc.m_current
|
||||
u = self.objFunc.u_current or self.prob.fields(m)
|
||||
m = self.invProb.m_current
|
||||
u = self.invProb.u_current or self.prob.fields(m)
|
||||
|
||||
x0 = np.random.rand(*m.shape)
|
||||
t = x0.dot(self.objFunc.dataObj2Deriv(m,x0,u=u))
|
||||
t = x0.dot(self.dmisfit.dataObj2Deriv(self.prob,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
|
||||
self.invProb.beta = self.beta0
|
||||
|
||||
|
||||
class BetaSchedule(InversionDirective):
|
||||
@@ -142,7 +144,7 @@ class BetaSchedule(InversionDirective):
|
||||
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
|
||||
self.invProb.beta /= self.coolingFactor
|
||||
|
||||
|
||||
# class UpdateReferenceModel(Parameter):
|
||||
@@ -154,5 +156,5 @@ class BetaSchedule(InversionDirective):
|
||||
# if mref is None:
|
||||
# if self.debug: print 'UpdateReferenceModel is using mref0'
|
||||
# mref = self.mref0
|
||||
# self.m_prev = self.objFunc.m_current
|
||||
# self.m_prev = self.invProb.m_current
|
||||
# return mref
|
||||
|
||||
@@ -51,12 +51,12 @@ def run(N, plotIt=True):
|
||||
M = prob.mesh
|
||||
|
||||
reg = Regularization.Tikhonov(mesh)
|
||||
objFunc = ObjFunction.BaseObjFunction(survey, reg)
|
||||
dmis = DataMisfit.l2_DataMisfit()
|
||||
opt = Optimization.InexactGaussNewton(maxIter=20)
|
||||
inv = Inversion.BaseInversion(objFunc, opt)
|
||||
invProb = InvProblem.BaseInvProblem(prob, reg, dmis, opt)
|
||||
beta = Directives.BetaSchedule()
|
||||
betaest = Directives.BetaEstimate_ByEig()
|
||||
inv.ruleList = Directives.DirectiveList(betaest, beta)
|
||||
inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest])
|
||||
m0 = np.zeros_like(survey.mtrue)
|
||||
|
||||
mrec = inv.run(m0)
|
||||
|
||||
@@ -0,0 +1,90 @@
|
||||
import Utils, Survey, Problem, numpy as np, scipy.sparse as sp, gc
|
||||
from Utils.SolverUtils import *
|
||||
from DataMisfit import _splitForward
|
||||
|
||||
class BaseInvProblem(object):
|
||||
"""BaseInvProblem(forward, reg, **kwargs)"""
|
||||
|
||||
__metaclass__ = Utils.SimPEGMetaClass
|
||||
|
||||
beta = 1.0 #: Trade-off parameter
|
||||
|
||||
debug = False #: Print debugging information
|
||||
counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things
|
||||
|
||||
reg = None #: Regularization
|
||||
dmisfit = None #: DataMisfit
|
||||
opt = None #: Optimization program
|
||||
|
||||
u_current = None #: The most current evaluated field
|
||||
m_current = None #: The most current model
|
||||
|
||||
|
||||
def __init__(self, forward, reg, dmisfit, opt, **kwargs):
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
self.prob, self.survey = _splitForward(forward)
|
||||
self.reg = reg
|
||||
self.dmisfit = dmisfit
|
||||
self.opt = opt
|
||||
|
||||
@Utils.callHooks('startup')
|
||||
def startup(self, m0):
|
||||
"""startup(m0)
|
||||
|
||||
Called when inversion is first starting.
|
||||
"""
|
||||
if self.debug: print 'Calling InvProblem.startup'
|
||||
|
||||
if self.reg.mref is None:
|
||||
print 'Regularization has not set mref. SimPEG.InvProblem will set it to m0.'
|
||||
self.reg.mref = m0
|
||||
|
||||
self.phi_d = np.nan
|
||||
self.phi_m = np.nan
|
||||
|
||||
self.m_current = m0
|
||||
|
||||
print 'Setting bfgsH0 to the inverse of the modelObj2Deriv. Done using direct methods.'
|
||||
self.opt.bfgsH0 = Solver(self.reg.modelObj2Deriv(self.m_current))
|
||||
|
||||
@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
|
||||
forward = self.prob
|
||||
gc.collect()
|
||||
|
||||
u = self.prob.fields(m)
|
||||
self.u_current = u
|
||||
|
||||
phi_d = self.dmisfit.dataObj(forward, m, u=u)
|
||||
phi_m = self.reg.modelObj(m)
|
||||
|
||||
self.dpred = self.survey.dpred(m, u=u) # This is a cheap matrix vector calculation.
|
||||
|
||||
self.phi_d, self.phi_d_last = phi_d, self.phi_d
|
||||
self.phi_m, self.phi_m_last = phi_m, self.phi_m
|
||||
|
||||
f = phi_d + self.beta * phi_m
|
||||
|
||||
out = (f,)
|
||||
if return_g:
|
||||
phi_dDeriv = self.dmisfit.dataObjDeriv(forward, 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.dmisfit.dataObj2Deriv(forward, m, v, u=u)
|
||||
phi_m2Deriv = self.reg.modelObj2Deriv(m, v=v)
|
||||
|
||||
return phi_d2Deriv + self.beta * phi_m2Deriv
|
||||
|
||||
H = sp.linalg.LinearOperator( (m.size, m.size), H_fun, dtype=m.dtype )
|
||||
out += (H,)
|
||||
return out if len(out) > 1 else out[0]
|
||||
+12
-14
@@ -5,7 +5,7 @@ import Directives
|
||||
|
||||
|
||||
class BaseInversion(object):
|
||||
"""BaseInversion(objFunc, opt, **kwargs)
|
||||
"""BaseInversion(invProb, opt, **kwargs)
|
||||
"""
|
||||
|
||||
__metaclass__ = Utils.SimPEGMetaClass
|
||||
@@ -30,23 +30,21 @@ class BaseInversion(object):
|
||||
self._directiveList = value
|
||||
self._directiveList.inversion = self
|
||||
|
||||
def __init__(self, objFunc, opt, **kwargs):
|
||||
def __init__(self, invProb, **kwargs):
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
|
||||
self.objFunc = objFunc
|
||||
self.objFunc.parent = self
|
||||
self.invProb = invProb
|
||||
|
||||
self.opt = opt
|
||||
opt.callback = self._optCallback
|
||||
self.opt.parent = self
|
||||
self.opt = invProb.opt
|
||||
self.opt.callback = self._optCallback
|
||||
|
||||
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)
|
||||
# # 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)
|
||||
|
||||
@Utils.timeIt
|
||||
def run(self, m0):
|
||||
@@ -55,9 +53,9 @@ class BaseInversion(object):
|
||||
Runs the inversion!
|
||||
|
||||
"""
|
||||
self.objFunc.startup(m0)
|
||||
self.invProb.startup(m0)
|
||||
self.directiveList.call('initialize')
|
||||
self.m = self.opt.minimize(self.objFunc.evalFunction, m0)
|
||||
self.m = self.opt.minimize(self.invProb.evalFunction, m0)
|
||||
self.directiveList.call('finish')
|
||||
|
||||
return self.m
|
||||
|
||||
@@ -1,220 +0,0 @@
|
||||
import Utils, Survey, Problem, numpy as np, scipy.sparse as sp, gc
|
||||
|
||||
class BaseObjFunction(object):
|
||||
"""BaseObjFunction(forward, reg, **kwargs)"""
|
||||
|
||||
__metaclass__ = Utils.SimPEGMetaClass
|
||||
|
||||
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
|
||||
|
||||
surveyPair = Survey.BaseSurvey
|
||||
problemPair = Problem.BaseProblem
|
||||
|
||||
name = 'Base Objective Function' #: Name of the objective function
|
||||
|
||||
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
|
||||
|
||||
@property
|
||||
def inv(self): return self.parent
|
||||
@property
|
||||
def objFunc(self): return self
|
||||
@property
|
||||
def opt(self): return getattr(self.parent,'opt',None)
|
||||
|
||||
|
||||
def __init__(self, forward, reg, **kwargs):
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
|
||||
assert forward.ispaired, 'The forward problem and survey must be paired.'
|
||||
if isinstance(forward, self.surveyPair):
|
||||
self.survey = forward
|
||||
self.prob = forward.prob
|
||||
elif isinstance(forward, self.problemPair):
|
||||
self.prob = forward
|
||||
self.survey = forward.survey
|
||||
|
||||
|
||||
self.reg = reg
|
||||
self.reg.parent = self
|
||||
|
||||
|
||||
@Utils.callHooks('startup')
|
||||
def startup(self, m0):
|
||||
"""startup(m0)
|
||||
|
||||
Called when inversion is first starting.
|
||||
"""
|
||||
if self.debug: print 'Calling ObjFunction.startup'
|
||||
|
||||
if self.reg.mref is None:
|
||||
print 'Regularization has not set mref. SimPEG.ObjFunction will set it to m0.'
|
||||
self.reg.mref = m0
|
||||
|
||||
self.phi_d = np.nan
|
||||
self.phi_m = np.nan
|
||||
|
||||
self.m_current = m0
|
||||
|
||||
@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
|
||||
gc.collect()
|
||||
|
||||
u = self.prob.fields(m)
|
||||
self.u_current = u
|
||||
|
||||
phi_d = self.dataObj(m, u=u)
|
||||
phi_m = self.reg.modelObj(m)
|
||||
|
||||
self.dpred = self.survey.dpred(m, u=u) # This is a cheap matrix vector calculation.
|
||||
|
||||
self.phi_d, self.phi_d_last = phi_d, self.phi_d
|
||||
self.phi_m, self.phi_m_last = phi_m, self.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(m, v=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.survey.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.prob.fields(m)
|
||||
|
||||
R = self.survey.residualWeighted(m, u=u)
|
||||
|
||||
dmisfit = self.prob.Jtvec(m, self.survey.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.fields(m)
|
||||
|
||||
R = self.survey.residualWeighted(m, u=u)
|
||||
|
||||
# TODO: abstract to different norms a little cleaner.
|
||||
# \/ it goes here. in l2 it is the identity.
|
||||
dmisfit = self.prob.Jtvec_approx(m, self.survey.Wd * self.survey.Wd * self.prob.Jvec_approx(m, v, u=u), u=u)
|
||||
|
||||
return dmisfit
|
||||
+2
-28
@@ -1,5 +1,5 @@
|
||||
import Utils, numpy as np, scipy.sparse as sp
|
||||
from Utils.SolverUtils import Solver, SolverCG
|
||||
from Utils.SolverUtils import *
|
||||
norm = np.linalg.norm
|
||||
|
||||
|
||||
@@ -663,13 +663,7 @@ class BFGS(Minimize, Remember):
|
||||
Must be a SimPEG.Solver
|
||||
"""
|
||||
if getattr(self,'_bfgsH0',None) is None:
|
||||
# Check if it has been set by the user and the default is not being used.
|
||||
if self.parent is None:
|
||||
self._bfgsH0 = Solver(sp.identity(self.xc.size).tocsc(), flag='D')
|
||||
else:
|
||||
print 'Setting bfgsH0 to the inverse of the modelObj2Deriv. Done using direct methods.'
|
||||
objFunc = self.parent.objFunc
|
||||
self._bfgsH0 = Solver(objFunc.reg.modelObj2Deriv(objFunc.m_current))
|
||||
self._bfgsH0 = SolverDiag(sp.identity(self.xc.size))
|
||||
return self._bfgsH0
|
||||
|
||||
@bfgsH0.setter
|
||||
@@ -878,23 +872,3 @@ class NewtonRoot(object):
|
||||
break
|
||||
|
||||
return x
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
from SimPEG.Tests import Rosenbrock, checkDerivative
|
||||
import matplotlib.pyplot as plt
|
||||
x0 = np.array([2.6, 3.7])
|
||||
checkDerivative(Rosenbrock, x0, plotIt=False)
|
||||
|
||||
xOpt = GaussNewton(maxIter=20,tolF=1e-10,tolX=1e-10,tolG=1e-10).minimize(Rosenbrock,x0)
|
||||
print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1])
|
||||
xOpt = SteepestDescent(maxIter=30, maxIterLS=15,tolF=1e-10,tolX=1e-10,tolG=1e-10).minimize(Rosenbrock, x0)
|
||||
print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1])
|
||||
|
||||
|
||||
print 'test the newtonRoot finding.'
|
||||
fun = lambda x, return_g=True: np.sin(x) if not return_g else ( np.sin(x), Utils.sdiag( np.cos(x) ) )
|
||||
x = np.array([np.pi-0.3, np.pi+0.1, 0])
|
||||
pnt = NewtonRoot(comments=True).root(fun,x)
|
||||
print pnt
|
||||
|
||||
@@ -597,45 +597,6 @@ class BaseSurvey(object):
|
||||
"""
|
||||
return Utils.mkvc(self.dpred(m, u=u) - self.dobs)
|
||||
|
||||
|
||||
@property
|
||||
def Wd(self):
|
||||
"""
|
||||
Data weighting matrix. This is a covariance matrix used in::
|
||||
|
||||
def residualWeighted(m,u=None):
|
||||
return self.Wd*self.residual(m, u=u)
|
||||
|
||||
By default, this is based on the norm of the data plus a noise floor.
|
||||
|
||||
"""
|
||||
if getattr(self,'_Wd',None) is None:
|
||||
print 'SimPEG is making Survey.Wd to be norm of the data plus a floor.'
|
||||
eps = np.linalg.norm(Utils.mkvc(self.dobs),2)*1e-5
|
||||
self._Wd = 1/(abs(self.dobs)*self.std+eps)
|
||||
return self._Wd
|
||||
@Wd.setter
|
||||
def Wd(self, value):
|
||||
self._Wd = value
|
||||
|
||||
def residualWeighted(self, m, u=None):
|
||||
"""residualWeighted(m, u=None)
|
||||
|
||||
:param numpy.array m: geophysical model
|
||||
:param numpy.array u: fields
|
||||
:rtype: numpy.array
|
||||
:return: weighted data residual
|
||||
|
||||
The weighted data residual:
|
||||
|
||||
.. math::
|
||||
|
||||
\mu_\\text{data}^{\\text{weighted}} = \mathbf{W}_d(\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs})
|
||||
|
||||
Where \\\\(W_d\\\\) is a covariance matrix that weights the data residual.
|
||||
"""
|
||||
return Utils.mkvc(self.Wd*self.residual(m, u=u))
|
||||
|
||||
@property
|
||||
def isSynthetic(self):
|
||||
"Check if the data is synthetic."
|
||||
|
||||
@@ -114,3 +114,36 @@ def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5):
|
||||
Solver = SolverWrapD(sp.linalg.spsolve, factorize=False)
|
||||
SolverLU = SolverWrapD(sp.linalg.splu, factorize=True)
|
||||
SolverCG = SolverWrapI(sp.linalg.cg)
|
||||
|
||||
|
||||
class SolverDiag(object):
|
||||
"""docstring for SolverDiag"""
|
||||
def __init__(self, A):
|
||||
self.A = A
|
||||
self._diagonal = A.diagonal()
|
||||
|
||||
def __mul__(self, rhs):
|
||||
n = self.A.shape[0]
|
||||
assert rhs.size % n == 0, 'Incorrect shape of rhs.'
|
||||
nrhs = rhs.size // n
|
||||
|
||||
if len(rhs.shape) == 1 or rhs.shape[1] == 1:
|
||||
x = self._solve1(rhs)
|
||||
else:
|
||||
x = self._solveM(rhs)
|
||||
|
||||
if nrhs == 1:
|
||||
return x.flatten()
|
||||
elif nrhs > 1:
|
||||
return x.reshape((n,nrhs), order='F')
|
||||
|
||||
def _solve1(self, rhs):
|
||||
return rhs.flatten()/self._diagonal
|
||||
|
||||
def _solveM(self, rhs):
|
||||
n = self.A.shape[0]
|
||||
nrhs = rhs.size // n
|
||||
return rhs/self._diagonal.repeat(nrhs).reshape((n,nrhs))
|
||||
|
||||
def clean(self):
|
||||
pass
|
||||
|
||||
+2
-1
@@ -7,7 +7,8 @@ import Maps
|
||||
import Problem
|
||||
import Survey
|
||||
import Regularization
|
||||
import ObjFunction
|
||||
import DataMisfit
|
||||
import InvProblem
|
||||
import Optimization
|
||||
import Directives
|
||||
import Inversion
|
||||
|
||||
@@ -0,0 +1,52 @@
|
||||
.. _api_DataMisfit:
|
||||
|
||||
|
||||
Data Misfit
|
||||
***********
|
||||
|
||||
The data misfit using an l_2 norm is:
|
||||
|
||||
.. math::
|
||||
|
||||
\mu_\text{data} = {1\over 2}\left| \mathbf{W}_d (\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}_d (\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 \\\(\\mathbf{W}_d\\\) 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}_d \mathbf{R}
|
||||
|
||||
The second derivative is:
|
||||
|
||||
.. math::
|
||||
|
||||
\frac{\partial^2 \mu_\text{data}}{\partial^2 \mathbf{m}} = \mathbf{J}^\top \mathbf{W}_d \mathbf{W}_d \mathbf{J}
|
||||
|
||||
|
||||
The API
|
||||
=======
|
||||
|
||||
.. autoclass:: SimPEG.DataMisfit.BaseDataMisfit
|
||||
:members:
|
||||
:undoc-members:
|
||||
|
||||
Common Data Misfits
|
||||
===================
|
||||
|
||||
l2 norm
|
||||
-------
|
||||
|
||||
.. autoclass:: SimPEG.DataMisfit.l2_DataMisfit
|
||||
:members:
|
||||
:undoc-members:
|
||||
@@ -20,6 +20,7 @@ import sys, os
|
||||
|
||||
sys.path.append('../')
|
||||
|
||||
|
||||
# -- General configuration -----------------------------------------------------
|
||||
|
||||
# If your documentation needs a minimal Sphinx version, state it here.
|
||||
@@ -242,3 +243,5 @@ texinfo_documents = [
|
||||
|
||||
# How to display URL addresses: 'footnote', 'no', or 'inline'.
|
||||
#texinfo_show_urls = 'footnote'
|
||||
|
||||
autodoc_member_order = 'bysource'
|
||||
|
||||
+2
-1
@@ -47,8 +47,9 @@ Inversion
|
||||
*********
|
||||
|
||||
.. toctree::
|
||||
:maxdepth: 2
|
||||
:maxdepth: 3
|
||||
|
||||
api_DataMisfit
|
||||
api_Inverse
|
||||
api_Parameters
|
||||
|
||||
|
||||
Binary file not shown.
|
Before Width: | Height: | Size: 19 KiB After Width: | Height: | Size: 45 KiB |
Reference in New Issue
Block a user