diff --git a/SimPEG/forward/DCProblem/DCProblem.py b/SimPEG/forward/DCProblem/DCProblem.py index 1df8897b..45b8a2c1 100644 --- a/SimPEG/forward/DCProblem/DCProblem.py +++ b/SimPEG/forward/DCProblem/DCProblem.py @@ -105,6 +105,10 @@ class DCProblem(Problem): if __name__ == '__main__': + + from SimPEG.regularization import Regularization + from SimPEG import inverse + # Create the mesh h1 = np.ones(100) h2 = np.ones(100) @@ -143,7 +147,7 @@ if __name__ == '__main__': dobs, Wd = synthetic.createData(mSynth, std=0.05) u = synthetic.field(mSynth) - mesh.plotImage(u[:,10], showIt=True) + mesh.plotImage(u[:,10], showIt=False) # Now set up the problem to do some minimization problem = DCProblem(mesh) @@ -153,8 +157,15 @@ if __name__ == '__main__': problem.dobs = dobs m0 = mesh.gridCC[:,0]*0+sig1 - print problem.misfit(m0) - print problem.misfit(mSynth) + # print problem.misfit(m0) + # print problem.misfit(mSynth) + + opt = inverse.InexactGaussNewton() + reg = Regularization(mesh) + + inv = inverse.Inversion(problem, reg, opt) + + inv.run(m0) # Check Derivative derChk = lambda m: [problem.misfit(m), problem.misfitDeriv(m)] @@ -166,3 +177,5 @@ if __name__ == '__main__': w = np.random.rand(dobs.shape[0]) print w.dot(problem.J(mSynth, v, u=u)) print v.dot(problem.Jt(mSynth, w, u=u)) + + diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py index 04f9771a..0333b8b4 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -62,6 +62,15 @@ class Problem(object): def P(self, value): self._P = value + @property + def std(self): + """ + Estimated Standard Deviations. + """ + return self._std + @std.setter + def std(self, value): + self._std = value @property def dobs(self): diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 99173af0..8d92e15e 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -1,4 +1,5 @@ import numpy as np +from SimPEG.utils import sdiag class Inversion(object): """docstring for Inversion""" @@ -11,20 +12,18 @@ class Inversion(object): self.opt = opt @property - def W(self): + def Wd(self): """ Standard deviation weighting matrix. """ - return self._W - @W.setter - def W(self, value): - self._W = value + return sdiag(1/self.prob.std) def run(self, m0): + m = m0 self._iter = 0 while True: self._beta = self.getBeta() - self.opt.minimize(self.evalFunction,m) + m = self.opt.minimize(self.evalFunction,m) if self.stoppingCriteria(): break self._iter += 1 diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index c3ad8dd2..43776d33 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -111,13 +111,16 @@ class Minimize(object): self._STOP[4] = self._iter >= self.maxIter return all(self._STOP[0:3]) | any(self._STOP[3:]) + def projection(self, p): + return p + def linesearch(self, p): # Armijo linesearch descent = np.inner(self.g, p) t = 1 iterLS = 0 while iterLS < self.maxIterLS: - xt = self.xc + t*p + xt = self.projection(self.xc + t*p) ft = self.evalFunction(xt, return_g=False, return_H=False) if ft < self.f + t*self.LSreduction*descent: break diff --git a/SimPEG/inverse/__init__.py b/SimPEG/inverse/__init__.py index b2a5e506..c4ca335b 100644 --- a/SimPEG/inverse/__init__.py +++ b/SimPEG/inverse/__init__.py @@ -1 +1,2 @@ from Optimize import * +from Inversion import * diff --git a/SimPEG/regularization/__init__.py b/SimPEG/regularization/__init__.py new file mode 100644 index 00000000..0230f8c3 --- /dev/null +++ b/SimPEG/regularization/__init__.py @@ -0,0 +1 @@ +from Regularization import Regularization