diff --git a/SimPEG/DataMisfit.py b/SimPEG/DataMisfit.py index 417b94e1..c0bb38ad 100644 --- a/SimPEG/DataMisfit.py +++ b/SimPEG/DataMisfit.py @@ -1,20 +1,6 @@ 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 @@ -28,25 +14,16 @@ class BaseDataMisfit(object): 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 + def __init__(self, survey, **kwargs): + assert survey.ispaired, 'The survey must be paired to a problem.' + if isinstance(survey, Survey.BaseSurvey): + self.survey = survey + self.prob = survey.prob + Utils.setKwargs(self,**kwargs) @Utils.timeIt - def dataObj(self, forward, m, u=None): - """dataObj(forward, m, u=None) + def dataObj(self, m, u=None): + """dataObj(m, u=None) :param Problem,Survey forward: forward simulation :param numpy.array m: geophysical model @@ -58,8 +35,8 @@ class BaseDataMisfit(object): raise NotImplementedError('This method should be overwritten.') @Utils.timeIt - def dataObjDeriv(self, forward, m, u=None): - """dataObjDeriv(forward, m, u=None) + def dataObjDeriv(self, m, u=None): + """dataObjDeriv(m, u=None) :param Problem,Survey forward: forward simulation :param numpy.array m: geophysical model @@ -72,8 +49,8 @@ class BaseDataMisfit(object): @Utils.timeIt - def dataObj2Deriv(self, forward, m, v, u=None): - """dataObj2Deriv(forward, m, v, u=None) + def dataObj2Deriv(self, m, v, u=None): + """dataObj2Deriv(m, v, u=None) :param Problem,Survey forward: forward simulation :param numpy.array m: geophysical model @@ -100,7 +77,7 @@ class BaseDataMisfit(object): return survey.nD -class l2_DataMisfit(object): +class l2_DataMisfit(BaseDataMisfit): """ The data misfit with an l_2 norm: @@ -111,44 +88,52 @@ class l2_DataMisfit(object): """ - def __init__(self, **kwargs): - pass + def __init__(self, survey, **kwargs): + BaseDataMisfit.__init__(self, survey, **kwargs) - def getWd(self, survey): + @property + def Wd(self): """getWd(survey) - Get the data weighting matrix. + The data weighting matrix. - This is based on the norm of the data plus a noise floor. + The default 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)) + + if getattr(self, '_Wd', None) is None: + print 'SimPEG.l2_DataMisfit is creating default weightings for Wd.' + survey = self.survey + eps = np.linalg.norm(Utils.mkvc(survey.dobs),2)*1e-5 + self._Wd = Utils.sdiag(1/(abs(survey.dobs)*survey.std+eps)) + return self._Wd + + @Wd.setter + def Wd(self, value): + self._Wd = value @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) + def dataObj(self, m, u=None): + "dataObj2Deriv(m, u=None)" + prob = self.prob + survey = self.survey + R = self.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) + def dataObjDeriv(self, m, u=None): + "dataObj2Deriv(m, u=None)" + prob = self.prob + survey = self.survey 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) + return prob.Jtvec(m, self.Wd * (self.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) + def dataObj2Deriv(self, m, v, u=None): + "dataObj2Deriv(m, v, u=None)" + prob = self.prob 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) + return prob.Jtvec_approx(m, self.Wd * (self.Wd * prob.Jvec_approx(m, v, u=u)), u=u) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 2e2de55d..14637cf5 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -27,9 +27,9 @@ class InversionDirective(object): @property def dmisfit(self): return self.invProb.dmisfit @property - def survey(self): return self.invProb.survey + def survey(self): return self.dmisfit.survey @property - def prob(self): return self.invProb.prob + def prob(self): return self.dmisfit.prob def initialize(self): pass @@ -126,7 +126,7 @@ class BetaEstimate_ByEig(InversionDirective): u = self.invProb.u_current or self.prob.fields(m) x0 = np.random.rand(*m.shape) - t = x0.dot(self.dmisfit.dataObj2Deriv(self.prob,m,x0,u=u)) + t = x0.dot(self.dmisfit.dataObj2Deriv(m,x0,u=u)) b = x0.dot(self.reg.modelObj2Deriv(m, v=x0)) self.beta0 = self.beta0_ratio*(t/b) diff --git a/SimPEG/Examples/Linear.py b/SimPEG/Examples/Linear.py index 1a36e541..4307ceec 100644 --- a/SimPEG/Examples/Linear.py +++ b/SimPEG/Examples/Linear.py @@ -51,9 +51,9 @@ def run(N, plotIt=True): M = prob.mesh reg = Regularization.Tikhonov(mesh) - dmis = DataMisfit.l2_DataMisfit() + dmis = DataMisfit.l2_DataMisfit(survey) opt = Optimization.InexactGaussNewton(maxIter=20) - invProb = InvProblem.BaseInvProblem(prob, reg, dmis, opt) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt) beta = Directives.BetaSchedule() betaest = Directives.BetaEstimate_ByEig() inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest]) diff --git a/SimPEG/InvProblem.py b/SimPEG/InvProblem.py index b982c7c7..5f3d1abb 100644 --- a/SimPEG/InvProblem.py +++ b/SimPEG/InvProblem.py @@ -1,9 +1,11 @@ import Utils, Survey, Problem, numpy as np, scipy.sparse as sp, gc from Utils.SolverUtils import * -from DataMisfit import _splitForward +import DataMisfit +import Regularization + class BaseInvProblem(object): - """BaseInvProblem(forward, reg, **kwargs)""" + """BaseInvProblem(dmisfit, reg, opt)""" __metaclass__ = Utils.SimPEGMetaClass @@ -20,12 +22,16 @@ class BaseInvProblem(object): m_current = None #: The most current model - def __init__(self, forward, reg, dmisfit, opt, **kwargs): + def __init__(self, dmisfit, reg, opt, **kwargs): Utils.setKwargs(self, **kwargs) - self.prob, self.survey = _splitForward(forward) - self.reg = reg + assert isinstance(dmisfit, DataMisfit.BaseDataMisfit), 'dmisfit must be a DataMisfit class.' + assert isinstance(reg, Regularization.BaseRegularization), 'reg must be a Regularization class.' self.dmisfit = dmisfit + self.reg = reg self.opt = opt + self.prob, self.survey = dmisfit.prob, dmisfit.survey + #TODO: Remove: (and make iteration printers better!) + self.opt.parent = self @Utils.callHooks('startup') def startup(self, m0): @@ -60,7 +66,7 @@ class BaseInvProblem(object): u = self.prob.fields(m) self.u_current = u - phi_d = self.dmisfit.dataObj(forward, m, u=u) + phi_d = self.dmisfit.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. @@ -72,7 +78,7 @@ class BaseInvProblem(object): out = (f,) if return_g: - phi_dDeriv = self.dmisfit.dataObjDeriv(forward, m, u=u) + phi_dDeriv = self.dmisfit.dataObjDeriv(m, u=u) phi_mDeriv = self.reg.modelObjDeriv(m) g = phi_dDeriv + self.beta * phi_mDeriv @@ -80,7 +86,7 @@ class BaseInvProblem(object): if return_H: def H_fun(v): - phi_d2Deriv = self.dmisfit.dataObj2Deriv(forward, m, v, u=u) + phi_d2Deriv = self.dmisfit.dataObj2Deriv(m, v, u=u) phi_m2Deriv = self.reg.modelObj2Deriv(m, v=v) return phi_d2Deriv + self.beta * phi_m2Deriv diff --git a/SimPEG/Inversion.py b/SimPEG/Inversion.py index cd571afe..691f143e 100644 --- a/SimPEG/Inversion.py +++ b/SimPEG/Inversion.py @@ -33,7 +33,8 @@ class BaseInversion(object): self._directiveList = value self._directiveList.inversion = self - def __init__(self, invProb, **kwargs): + def __init__(self, invProb, directiveList=[], **kwargs): + self.directiveList = directiveList Utils.setKwargs(self, **kwargs) self.invProb = invProb @@ -43,11 +44,11 @@ class BaseInversion(object): 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): diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index ec481145..59572bbd 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -312,29 +312,6 @@ class Mesh2Mesh(IdentityMap): """ Takes a model on one mesh are translates it to another mesh. - .. plot:: - - from SimPEG import * - M = Mesh.TensorMesh([100,100]) - h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)]) - h1 = h1/h1.sum() - M2 = Mesh.TensorMesh([h1,h1]) - V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50) - v = Utils.mkvc(V) - modh = Maps.Mesh2Mesh([M,M2]) - modH = Maps.Mesh2Mesh([M2,M]) - H = modH.transform(v) - h = modh.transform(H) - ax = plt.subplot(131) - M.plotImage(v, ax=ax) - ax.set_title('Fine Mesh (Original)') - ax = plt.subplot(132) - M2.plotImage(H,clim=[0,1],ax=ax) - ax.set_title('Course Mesh') - ax = plt.subplot(133) - M.plotImage(h,clim=[0,1],ax=ax) - ax.set_title('Fine Mesh (Interpolated)') - """ def __init__(self, meshes, **kwargs): diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py index e63a5fe2..fc7941d8 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -5,6 +5,7 @@ norm = np.linalg.norm __all__ = ['Minimize', 'Remember', 'SteepestDescent', 'BFGS', 'GaussNewton', 'InexactGaussNewton', 'ProjectedGradient', 'NewtonRoot', 'StoppingCriteria', 'IterationPrinters'] +SolverICG = SolverWrapI(sp.linalg.cg, checkAccuracy=False) class StoppingCriteria(object): """docstring for StoppingCriteria""" @@ -72,9 +73,9 @@ class IterationPrinters(object): bSet = {"title": "bSet", "value": lambda M: np.sum(M.bindingSet(M.xc)), "width": 8, "format": "%d"} comment = {"title": "Comment", "value": lambda M: M.comment, "width": 12, "format": "%s"} - beta = {"title": "beta", "value": lambda M: M.parent.objFunc.beta, "width": 10, "format": "%1.2e"} - phi_d = {"title": "phi_d", "value": lambda M: M.parent.objFunc.phi_d, "width": 10, "format": "%1.2e"} - phi_m = {"title": "phi_m", "value": lambda M: M.parent.objFunc.phi_m, "width": 10, "format": "%1.2e"} + beta = {"title": "beta", "value": lambda M: M.parent.beta, "width": 10, "format": "%1.2e"} + phi_d = {"title": "phi_d", "value": lambda M: M.parent.phi_d, "width": 10, "format": "%1.2e"} + phi_m = {"title": "phi_m", "value": lambda M: M.parent.phi_m, "width": 10, "format": "%1.2e"} class Minimize(object): @@ -772,7 +773,7 @@ class InexactGaussNewton(BFGS, Minimize, Remember): @Utils.timeIt def findSearchDirection(self): - Hinv = SolverCG(self.H, M=self.approxHinv, tol=self.tolCG, maxiter=self.maxIterCG) + Hinv = SolverICG(self.H, M=self.approxHinv, tol=self.tolCG, maxiter=self.maxIterCG) p = Hinv * (-self.g) return p diff --git a/docs/api_Maps.rst b/docs/api_Maps.rst index 19c5e89b..606f1282 100644 --- a/docs/api_Maps.rst +++ b/docs/api_Maps.rst @@ -157,6 +157,33 @@ Vertical 1D Map Mesh to Mesh Map ---------------- + +.. plot:: + + from SimPEG import * + import matplotlib.pyplot as plt + M = Mesh.TensorMesh([100,100]) + h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)]) + h1 = h1/h1.sum() + M2 = Mesh.TensorMesh([h1,h1]) + V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50) + v = Utils.mkvc(V) + modh = Maps.Mesh2Mesh([M,M2]) + modH = Maps.Mesh2Mesh([M2,M]) + H = modH.transform(v) + h = modh.transform(H) + ax = plt.subplot(131) + M.plotImage(v, ax=ax) + ax.set_title('Fine Mesh (Original)') + ax = plt.subplot(132) + M2.plotImage(H,clim=[0,1],ax=ax) + ax.set_title('Course Mesh') + ax = plt.subplot(133) + M.plotImage(h,clim=[0,1],ax=ax) + ax.set_title('Fine Mesh (Interpolated)') + plt.show() + + .. autoclass:: SimPEG.Maps.Mesh2Mesh :members: :undoc-members: