diff --git a/SimPEG/forward/DCProblem.py b/SimPEG/forward/DCProblem.py index c2a615ed..132966a8 100644 --- a/SimPEG/forward/DCProblem.py +++ b/SimPEG/forward/DCProblem.py @@ -7,6 +7,7 @@ import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as linalg + class DCProblem(ModelTransforms.LogModel, Problem): """ **DCProblem** @@ -18,6 +19,11 @@ class DCProblem(ModelTransforms.LogModel, Problem): super(DCProblem, self).__init__(mesh) self.mesh.setCellGradBC('neumann') + def reshapeFields(self, u): + if len(u.shape) == 1: + u = u.reshape([-1, self.RHS.shape[1]], order='F') + return u + def createMatrix(self, m): """ Makes the matrix A(m) for the DC resistivity problem. @@ -38,11 +44,25 @@ class DCProblem(ModelTransforms.LogModel, Problem): A = D*Msig*G return A.tocsc() + def dpred(self, m, u=None): + """ + Predicted data. + + .. math:: + d_\\text{pred} = Pu(m) + """ + if u is None: + u = self.field(m) + + u = self.reshapeFields(u) + + return mkvc(self.P*u) + def field(self, m): A = self.createMatrix(m) solve = Solver(A) phi = solve.solve(self.RHS) - return phi + return mkvc(phi) def J(self, m, v, u=None): """ @@ -69,6 +89,8 @@ class DCProblem(ModelTransforms.LogModel, Problem): if u is None: u = self.field(m) + u = self.reshapeFields(u) + P = self.P D = self.mesh.faceDiv G = self.mesh.cellGrad @@ -83,13 +105,18 @@ class DCProblem(ModelTransforms.LogModel, Problem): dCdm[:, i] = D * ( sdiag( G * ui ) * ( Av_dm * ( mT_dm * v ) ) ) solve = Solver(dCdu) - # solve = linalg.factorized(dCdu) - Jv = - P * solve.solve(dCdm) - return Jv + return mkvc(Jv) def Jt(self, m, v, u=None): """Takes data, turns it into a model..ish""" + + if u is None: + u = self.field(m) + + u = self.reshapeFields(u) + v = self.reshapeFields(v) + P = self.P D = self.mesh.faceDiv G = self.mesh.cellGrad @@ -147,7 +174,7 @@ if __name__ == '__main__': import matplotlib.pyplot as plt # Create the mesh - h1 = np.ones(100) + h1 = np.ones(20) h2 = np.ones(100) mesh = TensorMesh([h1,h2]) @@ -156,12 +183,12 @@ if __name__ == '__main__': sig2 = np.log(0.01) # Create a synthetic model from a block in a half-space - p0 = [20, 20] - p1 = [50, 50] + p0 = [5, 10] + p1 = [15, 50] condVals = [sig1, sig2] mSynth = ModelBuilder.defineBlockConductivity(p0,p1,mesh.gridCC,condVals) plt.colorbar(mesh.plotImage(mSynth)) - # plt.show() + plt.show() # Set up the projection nelec = 50 @@ -184,7 +211,9 @@ if __name__ == '__main__': dobs, Wd = synthetic.createData(mSynth, std=0.05) u = synthetic.field(mSynth) - # mesh.plotImage(u[:,10], showIt=False) + u = synthetic.reshapeFields(u) + mesh.plotImage(u[:,10]) + # plt.show() # Now set up the problem to do some minimization problem = DCProblem(mesh) @@ -194,39 +223,16 @@ if __name__ == '__main__': problem.std = dobs*0 + 0.05 m0 = mesh.gridCC[:,0]*0+sig2 - - - # Adjoint Test - u = np.random.rand(mesh.nC, problem.RHS.shape[1]) - v = np.random.rand(mesh.nC) - w = np.random.rand(*dobs.shape) - Jv = mkvc(problem.J(mSynth, v, u=u)) - print mkvc(w).dot(Jv) - print v.dot(problem.Jt(mSynth, w, u=u)) - - # Check Derivative - dm = np.random.randn(*m0.shape) - for alp in np.logspace(-2,-6, 5): - a = problem.dpred(m0) - b = problem.dpred(m0 + alp*dm) - c = problem.J(m0, alp*dm) - print np.linalg.norm(a-b), np.linalg.norm(a-b+c) - - - # derChk = lambda m: [problem.dpred(m), problem.J(mSynth,m)] - # checkDerivative(derChk, mSynth) - - - opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=3) + opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) reg = Regularization(mesh) - - inv = inverse.Inversion(problem, reg, opt) + inv = inverse.Inversion(problem, reg, opt, beta0=1e4) # Check Derivative derChk = lambda m: [inv.dataObj(m), inv.dataObjDeriv(m)] checkDerivative(derChk, mSynth) + print inv.dataObj(m0) print inv.dataObj(mSynth) diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 827e12f8..3aa8ce58 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -16,7 +16,7 @@ class Inversion(object): self.setKwargs(**kwargs) def setKwargs(self, **kwargs): - # Set the variables, throw an error if they don't exist. + """Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist.""" for attr in kwargs: if hasattr(self, attr): setattr(self, attr, kwargs[attr]) diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py index d4829b7b..6221246f 100644 --- a/SimPEG/inverse/Optimize.py +++ b/SimPEG/inverse/Optimize.py @@ -39,8 +39,7 @@ class Minimize(object): self.setKwargs(**kwargs) def setKwargs(self, **kwargs): - """Sets key word arguments (kwargs) that are present in the object.""" - # Set the variables, throw an error if they don't exist. + """Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist.""" for attr in kwargs: if hasattr(self, attr): setattr(self, attr, kwargs[attr]) @@ -161,7 +160,7 @@ class Minimize(object): def printInit(self): """ - printInit is called at the beginning of the optimization routine. + **printInit** is called at the beginning of the optimization routine. If there is a parent object, printInit will check for a parent.printInit function and call that. @@ -177,7 +176,7 @@ class Minimize(object): def printIter(self): """ - printIter is called directly after function evaluations. + **printIter** is called directly after function evaluations. If there is a parent object, printIter will check for a parent.printIter function and call that. @@ -191,7 +190,7 @@ class Minimize(object): def printDone(self): """ - printDone is called at the end of the optimization routine. + **printDone** is called at the end of the optimization routine. If there is a parent object, printDone will check for a parent.printDone function and call that. @@ -386,12 +385,3 @@ if __name__ == '__main__': 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]) - - def simplePass(x): - return np.sin(x), sdiag(np.cos(x)) - - def simpleFail(x): - return np.sin(x), -sdiag(np.cos(x)) - - checkDerivative(simplePass, np.random.randn(5), plotIt=False) - checkDerivative(simpleFail, np.random.randn(5), plotIt=False) diff --git a/SimPEG/tests/TestUtils.py b/SimPEG/tests/TestUtils.py index ff7ee668..29f9086f 100644 --- a/SimPEG/tests/TestUtils.py +++ b/SimPEG/tests/TestUtils.py @@ -1,12 +1,15 @@ import numpy as np import matplotlib.pyplot as plt from pylab import norm -from SimPEG.utils import mkvc +from SimPEG.utils import mkvc, sdiag from SimPEG import utils from SimPEG.mesh import TensorMesh, LogicallyOrthogonalMesh import numpy as np import unittest +import inspect +happiness = ['The test be workin!', 'You get a gold star!', 'Yay passed!', 'Happy little convergence test!', 'That was easy!', 'Testing is important.', 'You are awesome.', 'Go Test Go!', 'Once upon a time, a happy little test passed.', 'And then everyone was happy.'] +sadness = ['No gold star for you.','Try again soon.','Thankfully, persistence is a great substitute for talent.','It might be easier to call this a feature...','Coffee break?', 'Boooooooo :(', 'Testing is important. Do it again.'] class OrderTest(unittest.TestCase): """ @@ -157,9 +160,10 @@ class OrderTest(unittest.TestCase): print '---------------------------------------------' passTest = np.mean(np.array(order)) > self.tolerance*self._expectedOrder if passTest: - print ['The test be workin!', 'You get a gold star!', 'Yay passed!', 'Happy little convergence test!', 'That was easy!'][np.random.randint(5)] + print happiness[np.random.randint(len(happiness))] else: print 'Failed to pass test on ' + self._meshType + '.' + print sadness[np.random.randint(len(sadness))] print '' self.assertTrue(passTest) @@ -222,7 +226,11 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): for i in range(num): Jt = fctn(x0+t[i]*dx) E0[i] = l2norm(Jt[0]-Jc[0]) # 0th order Taylor - E1[i] = l2norm(Jt[0]-Jc[0]-t[i]*Jc[1].dot(dx)) # 1st order Taylor + if inspect.isfunction(Jc[1]): + E1[i] = l2norm(Jt[0]-Jc[0]-t[i]*Jc[1](dx)) # 1st order Taylor + else: + # We assume it is a numpy.ndarray + E1[i] = l2norm(Jt[0]-Jc[0]-t[i]*Jc[1].dot(dx)) # 1st order Taylor order0 = np.log10(E0[:-1]/E0[1:]) order1 = np.log10(E1[:-1]/E1[1:]) print "%d\t%1.2e\t%1.3e\t\t%1.3e\t\t%1.3f" % (i, t[i], E0[i], E1[i], np.nan if i == 0 else order1[i-1]) @@ -238,9 +246,12 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): passTest = belowTol or correctOrder if passTest: - print "%s PASS! %s\n" % ('='*25, '='*25) + print "%s PASS! %s" % ('='*25, '='*25) + print happiness[np.random.randint(len(happiness))]+'\n' else: print "%s\n%s FAIL! %s\n%s" % ('*'*57, '<'*25, '>'*25, '*'*57) + print sadness[np.random.randint(len(sadness))]+'\n' + if plotIt: plt.figure() @@ -254,3 +265,19 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): plt.show() return passTest + + +if __name__ == '__main__': + + def simplePass(x): + return np.sin(x), sdiag(np.cos(x)) + + def simpleFunction(x): + return np.sin(x), lambda xi: sdiag(np.cos(x))*xi + + def simpleFail(x): + return np.sin(x), -sdiag(np.cos(x)) + + checkDerivative(simplePass, np.random.randn(5), plotIt=False) + checkDerivative(simpleFunction, np.random.randn(5), plotIt=False) + checkDerivative(simpleFail, np.random.randn(5), plotIt=False) diff --git a/SimPEG/tests/test_forward_DCproblem.py b/SimPEG/tests/test_forward_DCproblem.py index 6ffa4650..c6e6f9c2 100644 --- a/SimPEG/tests/test_forward_DCproblem.py +++ b/SimPEG/tests/test_forward_DCproblem.py @@ -3,9 +3,11 @@ import unittest from SimPEG.mesh import TensorMesh from SimPEG.utils import ModelBuilder, sdiag from SimPEG.forward import Problem, SyntheticProblem -from SimPEG.forward.DCProblem import DCProblem, DCutils +from SimPEG.forward.DCProblem import * from TestUtils import checkDerivative from scipy.sparse.linalg import dsolve +from SimPEG.regularization import Regularization +from SimPEG import inverse class DCProblemTests(unittest.TestCase): @@ -34,7 +36,7 @@ class DCProblemTests(unittest.TestCase): elecend = 0.5+spacelec*(nelec-1) elecLocR = np.linspace(elecini, elecend, nelec) rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5 - q, Q, rxmidloc = DCutils.genTxRxmat(nelec, spacelec, surfloc, elecini, mesh) + q, Q, rxmidloc = genTxRxmat(nelec, spacelec, surfloc, elecini, mesh) P = Q.T # Create some data @@ -52,22 +54,27 @@ class DCProblemTests(unittest.TestCase): problem.RHS = q problem.W = Wd problem.dobs = dobs + problem.std = dobs*0 + 0.05 + opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6) + reg = Regularization(mesh) + inv = inverse.Inversion(problem, reg, opt, beta0=1e4) + + self.inv = inv + self.reg = reg self.p = problem self.mesh = mesh self.m0 = mSynth self.dobs = dobs - def test_misfit(self): - print 'SimPEG.forward.DCProblem: Testing Misfit' - derChk = lambda m: [self.p.misfit(m), self.p.misfitDeriv(m)] + derChk = lambda m: [self.p.dpred(m), lambda mx: self.p.J(self.m0, mx)] passed = checkDerivative(derChk, self.m0, plotIt=False) self.assertTrue(passed) def test_adjoint(self): # Adjoint Test - u = np.random.rand(self.mesh.nC) + u = np.random.rand(self.mesh.nC*self.p.RHS.shape[1]) v = np.random.rand(self.mesh.nC) w = np.random.rand(self.dobs.shape[0]) wtJv = w.dot(self.p.J(self.m0, v, u=u)) @@ -75,6 +82,13 @@ class DCProblemTests(unittest.TestCase): passed = (wtJv - vtJtw) < 1e-10 self.assertTrue(passed) + def test_dataObj(self): + derChk = lambda m: [self.inv.dataObj(m), self.inv.dataObjDeriv(m)] + checkDerivative(derChk, self.m0, plotIt=False) + + def test_modelObj(self): + derChk = lambda m: [self.reg.modelObj(m), self.reg.modelObjDeriv(m)] + checkDerivative(derChk, self.m0, plotIt=False) if __name__ == '__main__': diff --git a/SimPEG/tests/test_forward_problem.py b/SimPEG/tests/test_forward_problem.py index d913a697..ba49efd8 100644 --- a/SimPEG/tests/test_forward_problem.py +++ b/SimPEG/tests/test_forward_problem.py @@ -2,6 +2,7 @@ import numpy as np import unittest from SimPEG.mesh import TensorMesh from SimPEG.forward import Problem +from SimPEG.regularization import Regularization from TestUtils import checkDerivative from scipy.sparse.linalg import dsolve @@ -15,7 +16,7 @@ class ProblemTests(unittest.TestCase): c = np.array([1, 4]) self.mesh2 = TensorMesh([a, b], np.array([3, 5])) self.p2 = Problem(self.mesh2) - + self.reg = Regularization(self.mesh2) def test_modelTransform(self): print 'SimPEG.forward.Problem: Testing Model Transform' @@ -23,6 +24,13 @@ class ProblemTests(unittest.TestCase): passed = checkDerivative(lambda m : [self.p2.modelTransform(m), self.p2.modelTransformDeriv(m)], m, plotIt=False) self.assertTrue(passed) + def test_regularization(self): + derChk = lambda m: [self.reg.modelObj(m), self.reg.modelObjDeriv(m)] + mSynth = np.random.randn(self.mesh2.nC) + checkDerivative(derChk, mSynth, plotIt=False) + + + if __name__ == '__main__': unittest.main() diff --git a/SimPEG/tests/test_utils.py b/SimPEG/tests/test_utils.py index f0b867ec..058dba56 100644 --- a/SimPEG/tests/test_utils.py +++ b/SimPEG/tests/test_utils.py @@ -1,6 +1,28 @@ import numpy as np import unittest from SimPEG.utils import mkvc, ndgrid, indexCube, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal +from SimPEG.tests import checkDerivative + + +class TestCheckDerivative(unittest.TestCase): + + def test_simplePass(self): + def simplePass(x): + return np.sin(x), sdiag(np.cos(x)) + passed = checkDerivative(simplePass, np.random.randn(5), plotIt=False) + self.assertTrue(passed, True) + + def test_simpleFunction(self): + def simpleFunction(x): + return np.sin(x), lambda xi: sdiag(np.cos(x))*xi + passed = checkDerivative(simpleFunction, np.random.randn(5), plotIt=False) + self.assertTrue(passed, True) + + def test_simpleFail(self): + def simpleFail(x): + return np.sin(x), -sdiag(np.cos(x)) + passed = checkDerivative(simpleFail, np.random.randn(5), plotIt=False) + self.assertTrue(not passed, True) class TestSequenceFunctions(unittest.TestCase): @@ -85,5 +107,6 @@ class TestSequenceFunctions(unittest.TestCase): self.assertTrue(np.linalg.norm(Z3.todense().ravel(), 2) < 1e-12) + if __name__ == '__main__': unittest.main()