diff --git a/SimPEG/DataMisfit.py b/SimPEG/DataMisfit.py new file mode 100644 index 00000000..870b5880 --- /dev/null +++ b/SimPEG/DataMisfit.py @@ -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) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index cff073ac..0334e0e4 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -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 diff --git a/SimPEG/Examples/Linear.py b/SimPEG/Examples/Linear.py index 3fa02c11..1a36e541 100644 --- a/SimPEG/Examples/Linear.py +++ b/SimPEG/Examples/Linear.py @@ -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) diff --git a/SimPEG/InvProblem.py b/SimPEG/InvProblem.py new file mode 100644 index 00000000..b982c7c7 --- /dev/null +++ b/SimPEG/InvProblem.py @@ -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] diff --git a/SimPEG/Inversion.py b/SimPEG/Inversion.py index 62772476..b380b3df 100644 --- a/SimPEG/Inversion.py +++ b/SimPEG/Inversion.py @@ -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 diff --git a/SimPEG/ObjFunction.py b/SimPEG/ObjFunction.py deleted file mode 100644 index a5332921..00000000 --- a/SimPEG/ObjFunction.py +++ /dev/null @@ -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 diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py index c8314af5..e63a5fe2 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -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 diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 206c8be3..5d0da5c4 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -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." diff --git a/SimPEG/Utils/SolverUtils.py b/SimPEG/Utils/SolverUtils.py index 50317b63..386a3f98 100644 --- a/SimPEG/Utils/SolverUtils.py +++ b/SimPEG/Utils/SolverUtils.py @@ -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 diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 7d894e65..5ba6ea55 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -7,7 +7,8 @@ import Maps import Problem import Survey import Regularization -import ObjFunction +import DataMisfit +import InvProblem import Optimization import Directives import Inversion diff --git a/docs/api_DataMisfit.rst b/docs/api_DataMisfit.rst new file mode 100644 index 00000000..3033acce --- /dev/null +++ b/docs/api_DataMisfit.rst @@ -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: diff --git a/docs/conf.py b/docs/conf.py index 86aacce4..f2ab0db1 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -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' diff --git a/docs/index.rst b/docs/index.rst index 4a45f25f..7b5f91e1 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -47,8 +47,9 @@ Inversion ********* .. toctree:: - :maxdepth: 2 + :maxdepth: 3 + api_DataMisfit api_Inverse api_Parameters diff --git a/docs/simpeg-framework.png b/docs/simpeg-framework.png index 0d020d63..158ea885 100644 Binary files a/docs/simpeg-framework.png and b/docs/simpeg-framework.png differ