wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-10 10:13:17 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merged in inverseProblem (pull request #13)

Browse Source 
Documentation, Utilities, Testing, DC Problem
This commit is contained in:
rowanc1
2013-10-16 09:20:04 -07:00
parent 9cf76966c2 4851295738
commit 83291451d3
30 changed files with 1140 additions and 244 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/DiffOperators.py
+8 -4
View File
@@ -1,6 +1,6 @@
import numpy as np
from scipy import sparse as sp
from utils import mkvc, sdiag, speye, kron3, spzeros
from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros
def ddx(n):
@@ -287,15 +287,19 @@ class DiffOperators(object):
nodalVectorAve = property(**nodalVectorAve())
def getEdgeMass(self, materialProp=None):
"""mass matix for products of edge functions w'*M(materialProp)*e"""
"""mass matrix for products of edge functions w'*M(materialProp)*e"""
if(materialProp is None):
materialProp = np.ones(self.nC)
Av = self.edgeAve
return sdiag(Av.T * (self.vol * mkvc(materialProp)))
def getFaceMass(self, materialProp=None):
"""mass matix for products of edge functions w'*M(materialProp)*e"""
"""mass matrix for products of face functions w'*M(materialProp)*f"""
if(materialProp is None):
materialProp = np.ones(self.nC)
Av = self.faceAve
return sdiag(Av.T*(self.vol*mkvc(materialProp)))
return sdiag(Av.T * (self.vol * mkvc(materialProp)))
def getFaceMassDeriv(self):
Av = self.faceAve
return Av.T * sdiag(self.vol)
SimPEG/GaussNewton.py
-141
View File
@@ -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])
SimPEG/InnerProducts.py
+1 -1
View File
@@ -1,5 +1,5 @@
from scipy import sparse as sp
from utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
from SimPEG.utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
import numpy as np
SimPEG/Solver.py
+75
View File
@@ -0,0 +1,75 @@
import numpy as np
import scipy.sparse.linalg as linalg
class Solver(object):
"""docstring for Solver"""
def __init__(self, A, doDirect=True, flag=None, options={}):
assert type(doDirect) is bool, 'doDirect must be a boolean'
assert flag in [None, 'L', 'U', 'D'], "flag must be set to None, 'L', 'U', or 'D'"
self.A = A
self.dsolve = None
self.doDirect = doDirect
self.flag = flag
self.options = options
def solve(self, b):
if self.flag is None and self.doDirect:
return self.solveDirect(b, **self.options)
elif self.flag is None and not self.doDirect:
return self.solveIter(b, **self.options)
elif self.flag == 'U':
return self.solveBackward(b)
elif self.flag == 'L':
return self.solveForward(b)
elif self.flag == 'D':
return self.solveDiagonal(b)
else:
raise Exception('Unknown flag.')
pass
def clean(self):
"""Cleans up the memory"""
del self.dsolve
self.dsolve = None
def solveDirect(self, b, backend='scipy'):
assert np.shape(self.A)[1] == np.shape(b)[0], 'Dimension mismatch'
if self.dsolve is None:
self.A = self.A.tocsc() # for efficiency
self.dsolve = linalg.factorized(self.A)
if len(b.shape) == 1 or b.shape[1] == 1:
# Just one RHS
return self.dsolve(b)
# Multiple RHSs
X = np.empty_like(b)
for i in range(b.shape[1]):
X[:,i] = self.dsolve(b[:,i])
return X
def solveIter(self, b, M=None, iterSolver='CG'):
pass
def solveBackward(self, b):
pass
def solveForward(self, b):
pass
def solveDiagonal(self, b):
diagA = self.A.diagonal()
if len(b.shape) == 1 or b.shape[1] == 1:
# Just one RHS
return b/diagA
# Multiple RHSs
X = np.empty_like(b)
for i in range(b.shape[1]):
X[:,i] = b[:,i]/diagA
return X
SimPEG/__init__.py
+2
View File
@@ -1,3 +1,5 @@
from TensorMesh import TensorMesh
from LogicallyOrthogonalMesh import LogicallyOrthogonalMesh
import utils
import inverse
from Solver import Solver
SimPEG/forward/DCProblem/DCProblem.py
+168
View File
@@ -0,0 +1,168 @@
from SimPEG import TensorMesh
from SimPEG.forward import Problem, SyntheticProblem
from SimPEG.tests import checkDerivative
from SimPEG.utils import ModelBuilder, sdiag
import numpy as np
import scipy.sparse.linalg as linalg
import DCutils
class DCProblem(Problem):
"""
**DCProblem**
Geophysical DC resistivity problem.
"""
def __init__(self, mesh):
super(DCProblem, self).__init__(mesh)
self.mesh.setCellGradBC('neumann')
def createMatrix(self, m):
"""
Makes the matrix A(m) for the DC resistivity problem.
:param numpy.array m: model
:rtype: scipy.csc_matrix
:return: A(m)
.. math::
c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0
Where M() is the mass matrix and mT is the model transform.
"""
D = self.mesh.faceDiv
G = self.mesh.cellGrad
sigma = self.modelTransform(m)
Msig = self.mesh.getFaceMass(sigma)
A = D*Msig*G
return A.tocsc()
def field(self, m):
A = self.createMatrix(m)
solve = linalg.factorized(A)
nRHSs = self.RHS.shape[1] # Number of RHSs
phi = np.zeros((self.mesh.nC, nRHSs)) + np.nan
for ii in range(nRHSs):
phi[:,ii] = solve(self.RHS[:,ii])
return phi
def J(self, m, v, u=None, solve=None):
"""
:param numpy.array m: model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:rtype: numpy.array
:return: Jv
.. math::
c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0
\\nabla_u (A(m)u - q) = A(m)
\\nabla_m (A(m)u - q) = G\\text{sdiag}(Du)\\nabla_m(M(mT(m)))
Where M() is the mass matrix and mT is the model transform.
.. math::
J = - P \left( \\nabla_u c(m, u) \\right)^{-1} \\nabla_m c(m, u)
J(v) = - P ( A(m)^{-1} ( G\\text{sdiag}(Du)\\nabla_m(M(mT(m))) v ) )
"""
P = self.P
D = self.mesh.faceDiv
G = self.mesh.cellGrad
A = self.createMatrix(m)
Av_dm = self.mesh.getFaceMassDeriv()
mT_dm = self.modelTransformDeriv(m)
dCdu = A
dCdm = D * ( sdiag( G * u ) * ( Av_dm * ( mT_dm * v ) ) )
if solve is None:
solve = linalg.factorized(dCdu)
Jv = - P * solve(dCdm)
return Jv
def Jt(self, m, v, u=None, solve=None):
P = self.P
D = self.mesh.faceDiv
G = self.mesh.cellGrad
A = self.createMatrix(m)
Av_dm = self.mesh.getFaceMassDeriv()
mT_dm = self.modelTransformDeriv(m)
dCdu = A.T
if solve is None:
solve = linalg.factorized(dCdu.tocsc())
w = solve(P.T*v)
Jtv = - mT_dm.T * ( Av_dm.T * ( sdiag( G * u ) * ( D.T * w ) ) )
return Jtv
if __name__ == '__main__':
# Create the mesh
h1 = np.ones(100)
h2 = np.ones(100)
mesh = TensorMesh([h1,h2])
# Create some parameters for the model
sig1 = 1
sig2 = 0.01
# Create a synthetic model from a block in a half-space
p0 = [20, 20]
p1 = [50, 50]
condVals = [sig1, sig2]
mSynth = ModelBuilder.defineBlockConductivity(p0,p1,mesh.gridCC,condVals)
mesh.plotImage(mSynth, showIt=False)
# Set up the projection
nelec = 50
spacelec = 2
surfloc = 0.5
elecini = 0.5
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)
P = Q.T
# Create some data
class syntheticDCProblem(DCProblem, SyntheticProblem):
pass
synthetic = syntheticDCProblem(mesh);
synthetic.P = P
synthetic.RHS = q
dobs, Wd = synthetic.createData(mSynth, std=0.05)
u = synthetic.field(mSynth)
mesh.plotImage(u[:,10], showIt=True)
# Now set up the problem to do some minimization
problem = DCProblem(mesh)
problem.P = P
problem.RHS = q
problem.W = Wd
problem.dobs = dobs
m0 = mesh.gridCC[:,0]*0+sig1
print problem.misfit(m0)
print problem.misfit(mSynth)
# Check Derivative
derChk = lambda m: [problem.misfit(m), problem.misfitDeriv(m)]
checkDerivative(derChk, mSynth)
# Adjoint Test
u = np.random.rand(mesh.nC)
v = np.random.rand(mesh.nC)
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))
SimPEG/forward/DCProblem/DCutils.py
+29
View File
@@ -0,0 +1,29 @@
import numpy as np
import scipy.sparse as sp
def genTxRxmat(nelec, spacelec, surfloc, elecini, mesh):
""" Generate projection matrix (Q) and """
elecend = 0.5+spacelec*(nelec-1)
elecLocR = np.linspace(elecini, elecend, nelec)
elecLocT = elecLocR+1
nrx = nelec-1
ntx = nelec-1
q = np.zeros((mesh.nC, ntx))
Q = np.zeros((mesh.nC, nrx))
for i in range(nrx):
rxind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i]))
rxind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i+1]))
txind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i]))
txind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i+1]))
q[txind1,i] = 1
q[txind2,i] = -1
Q[rxind1,i] = 1
Q[rxind2,i] = -1
Q = sp.csr_matrix(Q)
rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5
return q, Q, rxmidLoc
SimPEG/forward/DCProblem/__init__.py
+2
View File
@@ -0,0 +1,2 @@
from DCProblem import *
from DCutils import *
SimPEG/forward/Problem.py
+339
View File
@@ -0,0 +1,339 @@
import numpy as np
from SimPEG.utils import mkvc, sdiag
norm = np.linalg.norm
class Problem(object):
"""
Problem is the base class for all geophysical forward problems in SimPEG.
The problem is a partial differential equation of the form:
.. math::
c(m, u) = 0
Here, m is the model and u is the field (or fields).
Given the model, m, we can calculate the fields u(m),
however, the data we collect is a subset of the fields,
and can be defined by a linear projection, P.
.. math::
d_\\text{pred} = Pu(m)
We are interested in how changing the model transforms the data,
as such we can take write the Taylor expansion:
.. math::
Pu(m + hv) = Pu(m) + hP\\frac{\partial u(m)}{\partial m} v + \mathcal{O}(h^2 \left\| v \\right\| )
We can linearize and define the sensitivity matrix as:
.. math::
J = P\\frac{\partial u}{\partial m}
The sensitivity matrix, and it's transpose will be used in the inverse problem
to (locally) find how model parameters change the data, and optimize!
"""
def __init__(self, mesh):
self.mesh = mesh
@property
def RHS(self):
"""
Source matrix.
"""
return self._RHS
@RHS.setter
def RHS(self, value):
self._RHS = value
@property
def W(self):
"""
Standard deviation weighting matrix.
"""
return self._W
@W.setter
def W(self, value):
self._W = value
@property
def P(self):
"""
Projection matrix.
.. math::
d_\\text{pred} = Pu(m)
"""
return self._P
@P.setter
def P(self, value):
self._P = value
@property
def dobs(self):
"""
Observed data.
"""
return self._dobs
@dobs.setter
def dobs(self, value):
self._dobs = value
def evalFunction(self, m, doDerivative=True):
"""
:param numpy.array m: model
:param bool doDerivative: do you want to compute the derivative?
:rtype: numpy.array
:return: Jv
"""
f = self.misfit(m)
return f, g, H
def J(self, m, v, u=None):
"""
:param numpy.array m: model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:rtype: numpy.array
:return: Jv
Working with the general PDE, c(m, u) = 0, where m is the model and u is the field,
the sensitivity is defined as:
.. math::
J = P\\frac{\partial u}{\partial m}
We can take the derivative of the PDE:
.. math::
\\nabla_m c(m, u) \delta m + \\nabla_u c(m, u) \delta u = 0
If the forward problem is invertible, then we can rearrange for du/dm:
.. math::
J = - P \left( \\nabla_u c(m, u) \\right)^{-1} \\nabla_m c(m, u)
This can often be computed given a vector (i.e. J(v)) rather than stored, as J is a large dense matrix.
"""
pass
def Jt(self, m, v, u=None):
"""
:param numpy.array m: model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:rtype: numpy.array
:return: JTv
Transpose of J
"""
pass
def field(self, m):
"""
The field given the model.
.. math::
u(m)
"""
pass
def dpred(self, m, u=None):
"""
Predicted data.
.. math::
d_\\text{pred} = Pu(m)
"""
if u is None:
u = self.field(m)
return self.P*u
def modelTransform(self, m):
"""
:param numpy.array m: model
:rtype: numpy.array
:return: transformed model
The modelTransform changes the model into the physical property.
A common example of this is to invert for electrical conductivity
in log space. In this case, your model will be log(sigma) and to
get back to sigma, you can take the exponential:
.. math::
m = \log{\sigma}
\exp{m} = \exp{\log{\sigma}} = \sigma
"""
return np.exp(mkvc(m))
def modelTransformDeriv(self, m):
"""
:param numpy.array m: model
:rtype: scipy.csr_matrix
:return: derivative of transformed model
The modelTransform changes the model into the physical property.
The modelTransformDeriv provides the derivative of the modelTransform.
If the model transform is:
.. math::
m = \log{\sigma}
\exp{m} = \exp{\log{\sigma}} = \sigma
Then the derivative is:
.. math::
\\frac{\partial \exp{m}}{\partial m} = \\text{sdiag}(\exp{m})
"""
return sdiag(np.exp(mkvc(m)))
def misfit(self, m, u=None):
"""
:param numpy.array m: geophysical model
:param numpy.array u: fields
:rtype: float
:return: data misfit
The data misfit using an l_2 norm is:
.. math::
\mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2
Where P is a projection matrix that brings the field on the full domain to the data measurement locations;
u is the field of interest; d_obs is the observed data; and W is the weighting matrix.
"""
R = self.W*(self.dpred(m, u=u) - self.dobs)
R = mkvc(R)
return 0.5*R.dot(R)
def misfitDeriv(self, m, u=None):
"""
:param numpy.array m: geophysical model
:param numpy.array u: fields
:rtype: numpy.array
:return: data misfit derivative
The data misfit using an l_2 norm is:
.. math::
\mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2
If the field, u, is provided, the calculation of the data is fast:
.. math::
\mathbf{d}_\\text{pred} = \mathbf{Pu(m)}
\mathbf{R} = \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs})
Where P is a projection matrix that brings the field on the full domain to the data measurement locations;
u is the field of interest; d_obs is the observed data; and W is the weighting matrix.
The derivative of this, with respect to the model, is:
.. math::
\\frac{\partial \mu_\\text{data}}{\partial \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ R}
"""
if u is None:
u = self.field(m)
R = self.W*(self.dpred(m, u=u) - self.dobs)
dmisfit = 0
for i in range(self.RHS.shape[1]): # Loop over each right hand side
dmisfit += self.Jt(m, self.W[:,i]*R[:,i], u=u[:,i])
return dmisfit
def misfitDerivDeriv(self, m, u=None):
"""
:param numpy.array m: geophysical model
:param numpy.array u: fields
:rtype: numpy.array
:return: data misfit derivative
The data misfit using an l_2 norm is:
.. math::
\mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2
If the field, u, is provided, the calculation of the data is fast:
.. math::
\mathbf{d}_\\text{pred} = \mathbf{Pu(m)}
\mathbf{R} = \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs})
Where P is a projection matrix that brings the field on the full domain to the data measurement locations;
u is the field of interest; d_obs is the observed data; and W is the weighting matrix.
The derivative of this, with respect to the model, is:
.. math::
\\frac{\partial \mu_\\text{data}}{\partial \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ R}
\\frac{\partial^2 \mu_\\text{data}}{\partial^2 \mathbf{m}} = \mathbf{J}^\\top \mathbf{W \circ W J}
"""
if u is None:
u = self.field(m)
R = self.W*(self.dpred(m, u=u) - self.dobs)
dmisfit = 0
for i in range(self.RHS.shape[1]): # Loop over each right hand side
dmisfit += self.Jt(m, self.W[:,i]*R[:,i], u=u[:,i])
return dmisfit
class SyntheticProblem(object):
"""
Has helpful functions when dealing with synthetic problems
To use this class, inherit to your problem::
class mySyntheticExample(Problem, SyntheticProblem):
pass
"""
def createData(self, m, std=0.05):
"""
:param numpy.array m: geophysical model
:param numpy.array std: standard deviation
:rtype: numpy.array, numpy.array
:return: dobs, Wd
Create synthetic data given a model, and a standard deviation.
Returns the observed data with random Gaussian noise
and Wd which is the same size as data, and can be used to weight the inversion.
"""
dobs = self.dpred(m)
dobs = dobs
noise = std*abs(dobs)*np.random.randn(*dobs.shape)
dobs = dobs+noise
eps = np.linalg.norm(mkvc(dobs),2)*1e-5
Wd = 1/(abs(dobs)*std+eps)
return dobs, Wd
SimPEG/forward/__init__.py
+2
View File
@@ -0,0 +1,2 @@
from Problem import *
import DCProblem
SimPEG/inverse/Optimize.py
+150
View File
@@ -0,0 +1,150 @@
import numpy as np
import matplotlib.pyplot as plt
from SimPEG.utils import mkvc, sdiag
norm = np.linalg.norm
class Minimize(object):
"""docstring for Minimize"""
name = "GeneralOptimizationAlgorithm"
maxIter = 20
maxIterLS = 10
LSreduction = 1e-4
LSshorten = 0.5
tolF = 1e-4
tolX = 1e-4
tolG = 1e-4
eps = 1e-16
def __init__(self, problem, **kwargs):
self.problem = problem
self.setKwargs(**kwargs)
def setKwargs(self, **kwargs):
# Set the variables, throw an error if they don't exist.
for attr in kwargs:
if hasattr(self, attr):
setattr(self, attr, kwargs[attr])
else:
raise Exception('%s attr is not recognized' % attr)
def minimize(self, x0):
self.startup(x0)
self.printInit()
while True:
self.f, self.g, self.H = self.evalFunction(self.xc)
self.printIter()
if self.stoppingCriteria(): break
p = self.findSearchDirection()
xt, passLS = self.linesearch(p)
if not passLS:
xt = self.linesearchBreak(p)
self.doEndIteration(xt)
self.printDone()
return self.xc
def startup(self, x0):
self._iter = 0
self._iterLS = 0
self._STOP = np.zeros((5,1),dtype=bool)
self.x0 = x0
self.xc = x0
self.xOld = x0
def printInit(self):
print "%s %s %s" % ('='*22, self.name, '='*22)
print "iter\tJc\t\tnorm(dJ)\tLS"
print "%s" % '-'*57
def printIter(self):
print "%3d\t%1.2e\t%1.2e\t%d" % (self._iter, self.f, norm(self.g), self._iterLS)
def printDone(self):
print "%s STOP! %s" % ('-'*25,'-'*25)
print "%d : |fc-fOld| = %1.4e <= tolF*(1+|fStop|) = %1.4e" % (self._STOP[0], abs(self.f-self.fOld), self.tolF*(1+abs(self.fStop)))
print "%d : |xc-xOld| = %1.4e <= tolX*(1+|x0|) = %1.4e" % (self._STOP[1], norm(self.xc-self.xOld), self.tolX*(1+norm(self.x0)))
print "%d : |g| = %1.4e <= tolG*(1+|fStop|) = %1.4e" % (self._STOP[2], norm(self.g), self.tolG*(1+abs(self.fStop)))
print "%d : |g| = %1.4e <= 1e3*eps = %1.4e" % (self._STOP[3], norm(self.g), 1e3*self.eps)
print "%d : iter = %3d\t <= maxIter\t = %3d" % (self._STOP[4], self._iter, self.maxIter)
print "%s DONE! %s\n" % ('='*25,'='*25)
def evalFunction(self, x, doDerivative=True):
f, g, H = self.problem(x)
return f, g, H
def findSearchDirection(self):
return -self.g
def stoppingCriteria(self):
if self._iter == 0:
self.fStop = self.f # Save this for stopping criteria
# check stopping rules
self._STOP[0] = self._iter > 0 and (abs(self.f-self.fOld) <= self.tolF*(1+abs(self.fStop)))
self._STOP[1] = self._iter > 0 and (norm(self.xc-self.xOld) <= self.tolX*(1+norm(self.x0)))
self._STOP[2] = norm(self.g) <= self.tolG*(1+abs(self.fStop))
self._STOP[3] = norm(self.g) <= 1e3*self.eps
self._STOP[4] = self._iter >= self.maxIter
return all(self._STOP[0:3]) | any(self._STOP[3:])
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
ft, temp, temp = self.evalFunction(xt, doDerivative=False)
if ft < self.f + t*self.LSreduction*descent:
break
iterLS += 1
t = self.LSshorten*t
self._iterLS = iterLS
return xt, iterLS < self.maxIterLS
def linesearchBreak(self, p):
raise Exception('The linesearch got broken. Boo.')
def doEndIteration(self, xt):
# store old values
self.fOld = self.f
self.xOld, self.xc = self.xc, xt
self._iter += 1
class GaussNewton(Minimize):
name = 'GaussNewton'
def findSearchDirection(self):
return np.linalg.solve(self.H,-self.g)
class SteepestDescent(Minimize):
name = 'SteepestDescent'
def findSearchDirection(self):
return -self.g
if __name__ == '__main__':
from SimPEG.tests import Rosenbrock, checkDerivative
x0 = np.array([2.6, 3.7])
checkDerivative(Rosenbrock, x0, plotIt=False)
xOpt = GaussNewton(Rosenbrock, maxIter=20).minimize(x0)
print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1])
xOpt = SteepestDescent(Rosenbrock, maxIter=20, maxIterLS=15).minimize(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)
SimPEG/inverse/__init__.py
+1
View File
@@ -0,0 +1 @@
from Optimize import *
SimPEG/tests/OrderTest.py → SimPEG/tests/TestUtils.py
+104 -28
View File
@@ -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
SimPEG/tests/__init__.py
+2
View File
@@ -0,0 +1,2 @@
import TestUtils
from TestUtils import checkDerivative, Rosenbrock, OrderTest
SimPEG/tests/test_forward_DCproblem.py
+81
View File
@@ -0,0 +1,81 @@
import numpy as np
import unittest
from SimPEG import TensorMesh
from SimPEG.utils import ModelBuilder, sdiag
from SimPEG.forward import Problem, SyntheticProblem
from SimPEG.forward.DCProblem import DCProblem, DCutils
from TestUtils import checkDerivative
from scipy.sparse.linalg import dsolve
class DCProblemTests(unittest.TestCase):
def setUp(self):
# Create the mesh
h1 = np.ones(20)
h2 = np.ones(20)
mesh = TensorMesh([h1,h2])
# Create some parameters for the model
sig1 = 1
sig2 = 0.01
# Create a synthetic model from a block in a half-space
p0 = [2, 2]
p1 = [5, 5]
condVals = [sig1, sig2]
mSynth = ModelBuilder.defineBlockConductivity(p0,p1,mesh.gridCC,condVals)
# Set up the projection
nelec = 10
spacelec = 2
surfloc = 0.5
elecini = 0.5
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)
P = Q.T
# Create some data
class syntheticDCProblem(DCProblem, SyntheticProblem):
pass
synthetic = syntheticDCProblem(mesh);
synthetic.P = P
synthetic.RHS = q
dobs, Wd = synthetic.createData(mSynth, std=0.05)
# Now set up the problem to do some minimization
problem = DCProblem(mesh)
problem.P = P
problem.RHS = q
problem.W = Wd
problem.dobs = dobs
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)]
passed = checkDerivative(derChk, self.m0, plotIt=False)
self.assertTrue(passed)
def test_adjoint(self):
# Adjoint Test
u = np.random.rand(self.mesh.nC)
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))
vtJtw = v.dot(self.p.Jt(self.m0, w, u=u))
passed = (wtJv - vtJtw) < 1e-10
self.assertTrue(passed)
if __name__ == '__main__':
unittest.main()
SimPEG/tests/test_forward_problem.py
+28
View File
@@ -0,0 +1,28 @@
import numpy as np
import unittest
from SimPEG import TensorMesh
from SimPEG.forward import Problem
from TestUtils import checkDerivative
from scipy.sparse.linalg import dsolve
class ProblemTests(unittest.TestCase):
def setUp(self):
a = np.array([1, 1, 1])
b = np.array([1, 2])
c = np.array([1, 4])
self.mesh2 = TensorMesh([a, b], np.array([3, 5]))
self.p2 = Problem(self.mesh2)
def test_modelTransform(self):
print 'SimPEG.forward.Problem: Testing Model Transform'
m = np.random.rand(self.mesh2.nC)
passed = checkDerivative(lambda m : [self.p2.modelTransform(m), self.p2.modelTransformDeriv(m)], m, plotIt=False)
self.assertTrue(passed)
if __name__ == '__main__':
unittest.main()
SimPEG/tests/test_massMatrices.py
+1 -1
View File
@@ -1,6 +1,6 @@
import numpy as np
import unittest
from OrderTest import OrderTest
from TestUtils import OrderTest
# MATLAB code:
SimPEG/tests/test_operators.py
+1 -1
View File
@@ -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])
SimPEG/tests/test_tensorMesh.py
+2 -4
View File
@@ -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
SimPEG/tests/test_utils.py
+1 -3
View File
@@ -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):
SimPEG/utils/ModelBuilder.py
+42 -31
View File
@@ -3,22 +3,20 @@ import numpy as np
def getIndecesBlock(p0,p1,ccMesh):
"""
Creates a vector containing the block indexes in the cell centerd mesh.
Returns a tuple
Creates a vector containing the block indexes in the cell centerd mesh.
Returns a tuple
The block is defined by the points
p0 : describe the position of the left upper front corner, and
p1 : describe the position of the right bottom back corner.
The block is defined by the points
ccMesh represents the cell-centered mesh
p0, describe the position of the left upper front corner, and
The points p0 and p1 must live in the the same dimensional space as the mesh.
p1, describe the position of the right bottom back corner.
ccMesh represents the cell-centered mesh
The points p0 and p1 must live in the the same dimensional space as the mesh.
"""
# Validation of the input
assert type(p0) == np.ndarray, "Vector must be a numpy array"
assert type(p1) == np.ndarray, "Vector must be a numpy array"
# Validation: p0 and p1 live in the same dimensional space
assert len(p0) == len(p1), "Dimension mismatch. len(p0) != len(p1)"
@@ -47,7 +45,7 @@ def getIndecesBlock(p0,p1,ccMesh):
ind = np.where(indX & indY)
else:
elif dimMesh == 3:
# Define the points
x1 = p0[0]
y1 = p0[1]
@@ -68,9 +66,9 @@ def getIndecesBlock(p0,p1,ccMesh):
def defineBlockConductivity(p0,p1,ccMesh,condVals):
"""
Build a block with the conductivity specified by condVal. Returns an array.
condVals[0] conductivity of the block
condVals[1] conductivity of the ground
Build a block with the conductivity specified by condVal. Returns an array.
condVals[0] conductivity of the block
condVals[1] conductivity of the ground
"""
sigma = np.zeros(ccMesh.shape[0]) + condVals[1]
ind = getIndecesBlock(p0,p1,ccMesh)
@@ -84,7 +82,8 @@ def defineTwoLayeredConductivity(depth,ccMesh,condVals):
Define a two layered model. Depth of the first layer must be specified.
CondVals vector with the conductivity values of the layers. Eg:
Convention to number the layers:
Convention to number the layers::
<----------------------------|------------------------------------>
0 depth zf
1st layer 2nd layer
@@ -98,13 +97,16 @@ def defineTwoLayeredConductivity(depth,ccMesh,condVals):
# Identify 1st cell centered reference point
p0[0] = ccMesh[0,0]
p0[1] = ccMesh[0,1]
p0[2] = ccMesh[0,2]
if dim>1: p0[1] = ccMesh[0,1]
if dim>2: p0[2] = ccMesh[0,2]
# Identify the last cell-centered reference point
p1[0] = ccMesh[-1,0]
p1[1] = ccMesh[-1,1]
p1[2] = ccMesh[-1,2] - depth;
if dim>1: p1[1] = ccMesh[-1,1]
if dim>2: p1[2] = ccMesh[-1,2]
# The depth is always defined on the last one.
p1[len(p1)-1] -= depth
ind = getIndecesBlock(p0,p1,ccMesh)
@@ -117,23 +119,24 @@ def scalarConductivity(ccMesh,pFunction):
Define the distribution conductivity in the mesh according to the
analytical expression given in pFunction
"""
xCC = ccMesh[:,0]
yCC = ccMesh[:,1]
zCC = ccMesh[:,2]
dim = np.size(ccMesh[0,:])
CC = [ccMesh[:,0]]
if dim>1: CC.append(ccMesh[:,1])
if dim>2: CC.append(ccMesh[:,2])
sigma = pFunction(xCC,yCC,zCC)
sigma = pFunction(*CC)
return sigma
if __name__ == '__main__':
import sys
sys.path.append('../')
from TensorMesh import TensorMesh
from SimPEG import TensorMesh
from matplotlib import pyplot as plt
# Define the mesh
testDim = 3
testDim = 2
h1 = 0.3*np.ones(7)
h1[0] = 0.5
h1[-1] = 0.6
@@ -157,8 +160,8 @@ if __name__ == '__main__':
# ------------------- Test conductivities! --------------------------
print('Testing 1 block conductivity')
p0 = np.array([0.5,0.5,0.5])
p1 = np.array([1.0,1.0,1.0])
p0 = np.array([0.5,0.5,0.5])[:testDim]
p1 = np.array([1.0,1.0,1.0])[:testDim]
condVals = np.array([100,1e-6])
sigma = defineBlockConductivity(p0,p1,ccMesh,condVals)
@@ -167,6 +170,7 @@ if __name__ == '__main__':
print sigma.shape
M.plotImage(sigma)
print 'Done with block! :)'
plt.show()
# -----------------------------------------
print('Testing the two layered model')
@@ -178,11 +182,17 @@ if __name__ == '__main__':
M.plotImage(sigma)
print sigma
print 'layer model!'
plt.show()
# -----------------------------------------
print('Testing scalar conductivity')
pFunction = lambda x,y,z: np.exp(x+y+z)
if testDim == 1:
pFunction = lambda x: np.exp(x)
elif testDim == 2:
pFunction = lambda x,y: np.exp(x+y)
elif testDim == 3:
pFunction = lambda x,y,z: np.exp(x+y+z)
sigma = scalarConductivity(ccMesh,pFunction)
@@ -190,5 +200,6 @@ if __name__ == '__main__':
M.plotImage(sigma)
print sigma
print 'Scalar conductivity defined!'
plt.show()
# -----------------------------------------
SimPEG/utils/__init__.py
+4 -1
View File
@@ -1,4 +1,7 @@
import matutils
import sputils
import lomutils
import ModelBuilder
from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind
from sputils import spzeros, kron3, speye, sdiag
from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube, exampleLomGird
import ModelBuilder
SimPEG/utils/lomutils.py
+11 -16
View File
@@ -31,19 +31,18 @@ def volTetra(xyz, A, B, C, D):
"""
Returns the volume for tetrahedras volume specified by the indexes A to D.
:param numpy.array xyz: X,Y,Z vertex vector
:param numpy.array A,B,C,D: vert index of the tetrahedra
:rtype: numpy.array
:return: V, volume of the tetrahedra
Input:
xyz - X,Y,Z vertex vector
A,B,C,D - vert index of the tetrahedra
Algorithm http://en.wikipedia.org/wiki/Tetrahedron#Volume
Output:
V - volume
.. math::
Algorithm: http://en.wikipedia.org/wiki/Tetrahedron#Volume
V = {1 \over 3} A h
V = 1/3 A * h
V = 1/6 | ( a - d ) o ( ( b - d ) X ( c - d ) ) |
V = {1 \over 6} | ( a - d ) \cdot ( ( b - d ) ( c - d ) ) |
"""
@@ -69,7 +68,7 @@ def indexCube(nodes, gridSize, n=None):
Output:
index - index in the order asked e.g. 'ABCD' --> (A,B,C,D)
TWO DIMENSIONS:
TWO DIMENSIONS::
node(i,j) node(i,j+1)
A -------------- B
@@ -81,7 +80,7 @@ def indexCube(nodes, gridSize, n=None):
node(i+1,j) node(i+1,j+1)
THREE DIMENSIONS:
THREE DIMENSIONS::
node(i,j,k+1) node(i,j+1,k+1)
E --------------- F
@@ -97,10 +96,6 @@ def indexCube(nodes, gridSize, n=None):
D -------------- C
node(i+1,j,k) node(i+1,j+1,k)
@author Rowan Cockett
Last modified on: 2013/07/26
"""
assert type(nodes) == str, "Nodes must be a str variable: e.g. 'ABCD'"
@@ -211,7 +206,7 @@ def faceInfo(xyz, A, B, C, D, average=True, normalizeNormals=True):
#
# So also could be viewed as the average parallelogram.
#
# WARNING: This does not compute correctly for concave quadrilaterals
# TODO: This does not compute correctly for concave quadrilaterals
area = (length(nA)+length(nB)+length(nC)+length(nD))/4
return N, area
SimPEG/utils/matutils.py
+2 -2
View File
@@ -4,7 +4,7 @@ import numpy as np
def mkvc(x, numDims=1):
"""Creates a vector with the number of dimension specified
e.g.:
e.g.::
a = np.array([1, 2, 3])
@@ -43,7 +43,7 @@ def ndgrid(*args, **kwargs):
The inputs can be a list or separate arguments.
e.g.
e.g.::
a = np.array([1, 2, 3])
b = np.array([1, 2])
docs/api_GaussNewton.rst
-8
View File
@@ -1,8 +0,0 @@
.. _api_GaussNewton:
Gauss Newton
************
.. automodule:: SimPEG.GaussNewton
:members:
:undoc-members:
docs/api_Optimize.rst
+8
View File
@@ -0,0 +1,8 @@
.. _api_Inverse:
Optimize
********
.. automodule:: SimPEG.inverse.Optimize
:members:
:undoc-members:
docs/api_Problem.rst
+26
View File
@@ -0,0 +1,26 @@
.. _api_Problem:
Problem
*******
.. automodule:: SimPEG.forward.Problem
:members:
:undoc-members:
DCProblem
*********
.. automodule:: SimPEG.forward.DCProblem.DCProblem
:members:
:undoc-members:
DCutils
*******
.. automodule:: SimPEG.forward.DCProblem.DCutils
:members:
:undoc-members:
docs/api_Tests.rst
+8
View File
@@ -0,0 +1,8 @@
.. _api_Tests:
Testing SimPEG
**************
.. automodule:: SimPEG.tests.TestUtils
:members:
:undoc-members:
docs/api_Utils.rst
+21
View File
@@ -0,0 +1,21 @@
.. _api_Utils:
Utilities
*********
.. automodule:: SimPEG.utils.matutils
:members:
:undoc-members:
.. automodule:: SimPEG.utils.sputils
:members:
:undoc-members:
.. automodule:: SimPEG.utils.lomutils
:members:
:undoc-members:
.. automodule:: SimPEG.utils.ModelBuilder
:members:
:undoc-members:
docs/index.rst
+21 -3
View File
@@ -30,20 +30,38 @@ Meshing & Operators
api_DiffOperators
api_InnerProducts
Forward Problems
================
.. toctree::
:maxdepth: 2
api_Problem
Inversion
=========
.. toctree::
:maxdepth: 2
api_GaussNewton
api_Optimize
Example Problems
================
Testing SimPEG
==============
.. toctree::
:maxdepth: 2
api_Tests
Utility Codes
=============
.. toctree::
:maxdepth: 2
api_Utils
Project Index & Search