diff --git a/SimPEG/GaussNewton.py b/SimPEG/GaussNewton.py deleted file mode 100644 index bd6d37c3..00000000 --- a/SimPEG/GaussNewton.py +++ /dev/null @@ -1,141 +0,0 @@ -import numpy as np -import matplotlib.pyplot as plt -from pylab import norm - -def GaussNewton(fctn, x0,maxIter=20, maxIterLS=10, LSreduction=1e-4, tolJ=1e-3, tolX=1e-3, - tolG=1e-3, eps=1e-16, xStop=[]): - """ - GaussNewton Optimization - - Input: - ------ - fctn - objective Function (lambda function) - x0 - starting guess - - Output: - ------- - xOpt - numerical optimizer - """ - # initial output - print "%s GaussNewton %s" % ('='*22,'='*22) - print "iter\tJc\t\tnorm(dJ)\tLS" - print "%s" % '-'*57 - - # evaluate stopping criteria - if xStop==[]: - xStop=x0 - Jstop = fctn(xStop) - print "%3d\t%1.2e" % (-1, Jstop[0]) - - # initialize - xc = x0 - STOP = np.zeros((5,1),dtype=bool) - iterLS=0; iter=0 - - Jold = Jstop - xOld=xc - while 1: - # evaluate objective function - Jc,dJ,H = fctn(xc) - print "%3d\t%1.2e\t%1.2e\t%d" % (iter, Jc[0],norm(dJ),iterLS) - - # check stopping rules - STOP[0] = (iter>0) & (abs(Jc[0]-Jold[0]) <= tolJ*(1+abs(Jstop[0]))) - STOP[1] = (iter>0) & (norm(xc-xOld) <= tolX*(1+norm(x0))) - STOP[2] = norm(dJ) <= tolG*(1+abs(Jstop[0])) - STOP[3] = norm(dJ) <= 1e3*eps - STOP[4] = (iter >= maxIter) - if all(STOP[0:3]) | any(STOP[3:]): - break - - # get search direction - dx = np.linalg.solve(H,-dJ) - - # Armijo linesearch - descent = np.dot(dJ.T,dx) - LS =0; t = 1; iterLS=1 - while (iterLS eps] - order1 = order1[E1[1:] > eps] - belowTol = order1.size == 0 and order0.size > 0 - correctOrder = order1.size > 0 and np.mean(order1) > tolerance * expectedOrder - - passTest = belowTol or correctOrder - - if passTest: - print "%s PASS! %s\n" % ('='*25, '='*25) - else: - print "%s\n%s FAIL! %s\n%s" % ('*'*57, '<'*25, '>'*25, '*'*57) - - if plotIt: - plt.figure() - plt.clf() - plt.loglog(t, E0, 'b') - plt.loglog(t, E1, 'g--') - plt.title('checkDerivative') - plt.xlabel('h') - plt.ylabel('error of Taylor approximation') - plt.legend(['0th order', '1st order'], loc='upper left') - plt.show() - - return passTest - if __name__ == '__main__': x0 = np.array([2.6, 3.7]) checkDerivative(Rosenbrock, x0, plotIt=False) diff --git a/SimPEG/tests/OrderTest.py b/SimPEG/tests/TestUtils.py similarity index 56% rename from SimPEG/tests/OrderTest.py rename to SimPEG/tests/TestUtils.py index fb305d32..79340da7 100644 --- a/SimPEG/tests/OrderTest.py +++ b/SimPEG/tests/TestUtils.py @@ -1,5 +1,7 @@ -import sys -sys.path.append('../../') +import numpy as np +import matplotlib.pyplot as plt +from pylab import norm +from SimPEG.utils import mkvc from SimPEG import TensorMesh, utils, LogicallyOrthogonalMesh import numpy as np import unittest @@ -16,46 +18,47 @@ class OrderTest(unittest.TestCase): Given are an operator A and its discretization A[h]. For a given test function f and h --> 0 we compare: - error(h) = \| A[h](f) - A(f) \|_{\infty} + .. math:: + error(h) = \| A[h](f) - A(f) \|_{\infty} Note that you can provide any norm. Test is passed when estimated rate order of convergence is at least within the specified tolerance of the estimated rate supplied by the user. - Minimal example for a curl operator: + Minimal example for a curl operator:: - class TestCURL(OrderTest): - name = "Curl" + class TestCURL(OrderTest): + name = "Curl" - def getError(self): - # For given Mesh, generate A[h], f and A(f) and return norm of error. + def getError(self): + # For given Mesh, generate A[h], f and A(f) and return norm of error. - fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x)) - sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y)) + fun = lambda x: np.cos(x) # i (cos(y)) + j (cos(z)) + k (cos(x)) + sol = lambda x: np.sin(x) # i (sin(z)) + j (sin(x)) + k (sin(y)) - Ex = fun(self.M.gridEx[:, 1]) - Ey = fun(self.M.gridEy[:, 2]) - Ez = fun(self.M.gridEz[:, 0]) - f = np.concatenate((Ex, Ey, Ez)) + Ex = fun(self.M.gridEx[:, 1]) + Ey = fun(self.M.gridEy[:, 2]) + Ez = fun(self.M.gridEz[:, 0]) + f = np.concatenate((Ex, Ey, Ez)) - Fx = sol(self.M.gridFx[:, 2]) - Fy = sol(self.M.gridFy[:, 0]) - Fz = sol(self.M.gridFz[:, 1]) - Af = np.concatenate((Fx, Fy, Fz)) + Fx = sol(self.M.gridFx[:, 2]) + Fy = sol(self.M.gridFy[:, 0]) + Fz = sol(self.M.gridFz[:, 1]) + Af = np.concatenate((Fx, Fy, Fz)) - # Generate DIV matrix - Ah = self.M.edgeCurl + # Generate DIV matrix + Ah = self.M.edgeCurl - curlE = Ah*E - err = np.linalg.norm((Ah*f -Af), np.inf) - return err + curlE = Ah*E + err = np.linalg.norm((Ah*f -Af), np.inf) + return err - def test_order(self): - # runs the test - self.orderTest() + def test_order(self): + # runs the test + self.orderTest() See also: test_operatorOrder.py @@ -159,5 +162,78 @@ class OrderTest(unittest.TestCase): print '' self.assertTrue(passTest) -if __name__ == '__main__': - unittest.main() +def Rosenbrock(x): + """Rosenbrock function for testing GaussNewton scheme""" + + f = 100*(x[1]-x[0]**2)**2+(1-x[0])**2 + g = np.array([2*(200*x[0]**3-200*x[0]*x[1]+x[0]-1), 200*(x[1]-x[0]**2)]) + H = np.array([[-400*x[1]+1200*x[0]**2+2, -400*x[0]], [-400*x[0], 200]]) + return f, g, H + +def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None): + """ + Basic derivative check + + Compares error decay of 0th and 1st order Taylor approximation at point + x0 for a randomized search direction. + + :param lambda fctn: function handle + :param numpy.array x0: point at which to check derivative + :param int num: number of times to reduce step length, h + :param bool plotIt: if you would like to plot + :param numpy.array dx: step direction + :rtype: bool + :return: did you pass the test?! + + """ + + print "%s checkDerivative %s" % ('='*20, '='*20) + print "iter\th\t\t|J0-Jt|\t\t|J0+h*dJ'*dx-Jt|\tOrder\n%s" % ('-'*57) + + Jc = fctn(x0) + + x0 = mkvc(x0) + + if dx is None: + dx = np.random.randn(len(x0)) + + t = np.logspace(-1, -num, num) + E0 = np.ones(t.shape) + E1 = np.ones(t.shape) + + l2norm = lambda x: np.sqrt(np.inner(x, x)) # because np.norm breaks if they are scalars? + 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 + 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]) + + tolerance = 0.9 + expectedOrder = 2 + eps = 1e-10 + order0 = order0[E0[1:] > eps] + order1 = order1[E1[1:] > eps] + belowTol = order1.size == 0 and order0.size > 0 + correctOrder = order1.size > 0 and np.mean(order1) > tolerance * expectedOrder + + passTest = belowTol or correctOrder + + if passTest: + print "%s PASS! %s\n" % ('='*25, '='*25) + else: + print "%s\n%s FAIL! %s\n%s" % ('*'*57, '<'*25, '>'*25, '*'*57) + + if plotIt: + plt.figure() + plt.clf() + plt.loglog(t, E0, 'b') + plt.loglog(t, E1, 'g--') + plt.title('checkDerivative') + plt.xlabel('h') + plt.ylabel('error of Taylor approximation') + plt.legend(['0th order', '1st order'], loc='upper left') + plt.show() + + return passTest diff --git a/SimPEG/tests/__init__.py b/SimPEG/tests/__init__.py new file mode 100644 index 00000000..f896cfd0 --- /dev/null +++ b/SimPEG/tests/__init__.py @@ -0,0 +1,2 @@ +import TestUtils +from TestUtils import checkDerivative, Rosenbrock, OrderTest diff --git a/SimPEG/tests/test_massMatrices.py b/SimPEG/tests/test_massMatrices.py index 13037156..79626ca1 100644 --- a/SimPEG/tests/test_massMatrices.py +++ b/SimPEG/tests/test_massMatrices.py @@ -1,6 +1,6 @@ import numpy as np import unittest -from OrderTest import OrderTest +from TestUtils import OrderTest # MATLAB code: diff --git a/SimPEG/tests/test_operators.py b/SimPEG/tests/test_operators.py index ab377089..6906c580 100644 --- a/SimPEG/tests/test_operators.py +++ b/SimPEG/tests/test_operators.py @@ -1,6 +1,6 @@ import numpy as np import unittest -from OrderTest import OrderTest +from TestUtils import OrderTest MESHTYPES = ['uniformTensorMesh', 'uniformLOM', 'rotateLOM'] call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1]) diff --git a/SimPEG/tests/test_tensorMesh.py b/SimPEG/tests/test_tensorMesh.py index 32d47353..300e6634 100644 --- a/SimPEG/tests/test_tensorMesh.py +++ b/SimPEG/tests/test_tensorMesh.py @@ -1,9 +1,7 @@ import numpy as np import unittest -import sys -sys.path.append('../') -from TensorMesh import TensorMesh -from OrderTest import OrderTest +from SimPEG import TensorMesh +from TestUtils import OrderTest from scipy.sparse.linalg import dsolve diff --git a/SimPEG/tests/test_utils.py b/SimPEG/tests/test_utils.py index 0dcdd362..f0b867ec 100644 --- a/SimPEG/tests/test_utils.py +++ b/SimPEG/tests/test_utils.py @@ -1,8 +1,6 @@ import numpy as np import unittest -import sys -sys.path.append('../') -from utils import mkvc, ndgrid, indexCube, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal +from SimPEG.utils import mkvc, ndgrid, indexCube, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal class TestSequenceFunctions(unittest.TestCase): diff --git a/docs/api_GaussNewton.rst b/docs/api_GaussNewton.rst deleted file mode 100644 index 65012347..00000000 --- a/docs/api_GaussNewton.rst +++ /dev/null @@ -1,8 +0,0 @@ -.. _api_GaussNewton: - -Gauss Newton -************ - -.. automodule:: SimPEG.GaussNewton - :members: - :undoc-members: diff --git a/docs/api_Optimize.rst b/docs/api_Optimize.rst new file mode 100644 index 00000000..9765bd3b --- /dev/null +++ b/docs/api_Optimize.rst @@ -0,0 +1,8 @@ +.. _api_Inverse: + +Optimize +******** + +.. automodule:: SimPEG.inverse.Optimize + :members: + :undoc-members: diff --git a/docs/api_Tests.rst b/docs/api_Tests.rst new file mode 100644 index 00000000..164f5b5d --- /dev/null +++ b/docs/api_Tests.rst @@ -0,0 +1,8 @@ +.. _api_Tests: + +Testing SimPEG +************** + +.. automodule:: SimPEG.tests.TestUtils + :members: + :undoc-members: diff --git a/docs/index.rst b/docs/index.rst index 48ebe940..b08624d5 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -30,14 +30,6 @@ Meshing & Operators api_DiffOperators api_InnerProducts -Inversion -========= - -.. toctree:: - :maxdepth: 2 - - api_GaussNewton - Forward Problems ================ @@ -46,6 +38,22 @@ Forward Problems api_Problem +Inversion +========= + +.. toctree:: + :maxdepth: 2 + + api_Optimize + +Testing SimPEG +============== + +.. toctree:: + :maxdepth: 2 + + api_Tests + Project Index & Search ======================