mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-14 11:18:18 +08:00
Moved testing code to the test directory. added __init__.py and some documentation.
This commit is contained in:
Rowan Cockett
@@ -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<maxIterLS):
|
||||
xt = xc + t*dx
|
||||
Jt = fctn(xt)
|
||||
LS = Jt[0]<Jc[0]+t*LSreduction*descent
|
||||
if LS:
|
||||
break
|
||||
iterLS = iterLS+1
|
||||
t = .5*t
|
||||
|
||||
# store old values
|
||||
Jold = Jc; xOld = xc
|
||||
# update
|
||||
xc = xt
|
||||
iter = iter +1
|
||||
|
||||
print "%s STOP! %s" % ('-'*25,'-'*25)
|
||||
print "%d : |Jc-Jold| = %1.4e <= tolJ*(1+|Jstop|) = %1.4e" % (STOP[0],abs(Jc[0]-Jold[0]),tolJ*(1+abs(Jstop[0])))
|
||||
print "%d : |xc-xOld| = %1.4e <= tolX*(1+|x0|) = %1.4e" % (STOP[1],norm(xc-xOld),tolX*(1+norm(x0)))
|
||||
print "%d : |dJ| = %1.4e <= tolG*(1+|Jstop|) = %1.4e" % (STOP[2],norm(dJ),tolG*(1+abs(Jstop[0])))
|
||||
print "%d : |dJ| = %1.4e <= 1e3*eps = %1.4e" % (STOP[3],norm(dJ),1e3*eps)
|
||||
print "%d : iter = %3d\t <= maxIter\t = %3d" % (STOP[4],iter,maxIter)
|
||||
print "%s DONE! %s\n" % ('='*25,'='*25)
|
||||
|
||||
return xc
|
||||
|
||||
def Rosenbrock(x):
|
||||
"""
|
||||
Rosenbrock function for testing GaussNewton scheme
|
||||
"""
|
||||
J = 100*(x[1]-x[0]**2)**2+(1-x[0])**2
|
||||
dJ = 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]],dtype=float);
|
||||
|
||||
return J,dJ,H
|
||||
|
||||
def checkDerivative(fctn,x0):
|
||||
"""
|
||||
Basic derivative check
|
||||
|
||||
Compares error decay of 0th and 1st order Taylor approximation at point
|
||||
x0 for a randomized search direction.
|
||||
|
||||
Input:
|
||||
------
|
||||
fctn - function handle
|
||||
x0 - point at which to check derivative
|
||||
"""
|
||||
|
||||
print "%s checkDerivative %s" % ('='*20,'='*20)
|
||||
print "iter\th\t\t|J0-Jt|\t\t|J0+h*dJ'*dx-Jt|"
|
||||
|
||||
Jc,dJ,H = fctn(x0)
|
||||
|
||||
dx = np.random.randn(len(x0),1)
|
||||
|
||||
t = np.logspace(-1,-10,10)
|
||||
E0 = np.zeros(t.shape)
|
||||
E1 = np.zeros(t.shape)
|
||||
|
||||
for i in range(0,10):
|
||||
Jt = fctn(x0+t[i]*dx)
|
||||
E0[i] = norm(Jt[0]-Jc[0]) # 0th order Taylor
|
||||
E1[i] = norm(Jt[0]-Jc[0]-t[i]*np.dot(dJ.T,dx)) # 1st order Taylor
|
||||
|
||||
print "%d\t%1.2e\t%1.3e\t%1.3e" % (i,t[i],E0[i],E1[i])
|
||||
|
||||
print "%s DONE! %s\n" % ('='*25,'='*25)
|
||||
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
|
||||
|
||||
if __name__ == '__main__':
|
||||
x0 = np.array([[2.6],[3.7]])
|
||||
fctn = lambda x:Rosenbrock(x)
|
||||
checkDerivative(fctn,x0)
|
||||
xOpt = GaussNewton(fctn,x0,maxIter=20)
|
||||
print "xOpt=[%f,%f]" % (xOpt[0],xOpt[1])
|
||||
|
||||
@@ -1,6 +1,6 @@
|
||||
from SimPEG import TensorMesh
|
||||
from SimPEG.forward import Problem, SyntheticProblem
|
||||
from SimPEG.inverse import checkDerivative
|
||||
from SimPEG.tests import checkDerivative
|
||||
from SimPEG.utils import ModelBuilder, sdiag
|
||||
import numpy as np
|
||||
import scipy.sparse.linalg as linalg
|
||||
|
||||
@@ -191,10 +191,6 @@ class Problem(object):
|
||||
"""
|
||||
return sdiag(np.exp(mkvc(m)))
|
||||
|
||||
def _test_modelTransformDeriv(self):
|
||||
m = np.random.rand(5)
|
||||
return checkDerivative(lambda m : [self.modelTransform(m), self.modelTransformDeriv(m)], m)
|
||||
|
||||
def misfit(self, m, u=None):
|
||||
"""
|
||||
:param numpy.array m: geophysical model
|
||||
|
||||
@@ -129,81 +129,6 @@ class SteepestDescent(Minimize):
|
||||
def findSearchDirection(self):
|
||||
return -self.g
|
||||
|
||||
|
||||
|
||||
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.
|
||||
|
||||
Input:
|
||||
------
|
||||
fctn - function handle
|
||||
x0 - point at which to check derivative
|
||||
"""
|
||||
|
||||
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
|
||||
|
||||
if __name__ == '__main__':
|
||||
x0 = np.array([2.6, 3.7])
|
||||
checkDerivative(Rosenbrock, x0, plotIt=False)
|
||||
|
||||
@@ -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
|
||||
@@ -0,0 +1,2 @@
|
||||
import TestUtils
|
||||
from TestUtils import checkDerivative, Rosenbrock, OrderTest
|
||||
@@ -1,6 +1,6 @@
|
||||
import numpy as np
|
||||
import unittest
|
||||
from OrderTest import OrderTest
|
||||
from TestUtils import OrderTest
|
||||
|
||||
|
||||
# MATLAB code:
|
||||
|
||||
@@ -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])
|
||||
|
||||
@@ -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
|
||||
|
||||
|
||||
|
||||
@@ -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):
|
||||
|
||||
@@ -1,8 +0,0 @@
|
||||
.. _api_GaussNewton:
|
||||
|
||||
Gauss Newton
|
||||
************
|
||||
|
||||
.. automodule:: SimPEG.GaussNewton
|
||||
:members:
|
||||
:undoc-members:
|
||||
@@ -0,0 +1,8 @@
|
||||
.. _api_Inverse:
|
||||
|
||||
Optimize
|
||||
********
|
||||
|
||||
.. automodule:: SimPEG.inverse.Optimize
|
||||
:members:
|
||||
:undoc-members:
|
||||
@@ -0,0 +1,8 @@
|
||||
.. _api_Tests:
|
||||
|
||||
Testing SimPEG
|
||||
**************
|
||||
|
||||
.. automodule:: SimPEG.tests.TestUtils
|
||||
:members:
|
||||
:undoc-members:
|
||||
+16
-8
@@ -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
|
||||
======================
|
||||
|
||||
Reference in New Issue
Block a user