diff --git a/SimPEG/Solver.py b/SimPEG/Solver.py
deleted file mode 100644
index 2087ef5f..00000000
--- a/SimPEG/Solver.py
+++ /dev/null
@@ -1,75 +0,0 @@
-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
diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py
index 8107e847..f50e1cef 100644
--- a/SimPEG/__init__.py
+++ b/SimPEG/__init__.py
@@ -1,5 +1,11 @@
-import mesh
import utils
+from utils import Solver
+import mesh
import inverse
+<<<<<<< HEAD
from Solver import Solver
import visulize
+=======
+import forward
+import regularization
+>>>>>>> refs/remotes/origin/master
diff --git a/SimPEG/forward/DCProblem.py b/SimPEG/forward/DCProblem.py
new file mode 100644
index 00000000..132966a8
--- /dev/null
+++ b/SimPEG/forward/DCProblem.py
@@ -0,0 +1,249 @@
+from SimPEG.mesh import TensorMesh
+from SimPEG.forward import Problem, SyntheticProblem, ModelTransforms
+from SimPEG.tests import checkDerivative
+from SimPEG.utils import ModelBuilder, sdiag, mkvc
+from SimPEG import Solver
+import numpy as np
+import scipy.sparse as sp
+import scipy.sparse.linalg as linalg
+
+
+class DCProblem(ModelTransforms.LogModel, Problem):
+ """
+ **DCProblem**
+
+ Geophysical DC resistivity problem.
+
+ """
+ def __init__(self, mesh):
+ 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.
+
+ :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 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 mkvc(phi)
+
+ 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
+
+ .. 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 ) )
+ """
+ if u is None:
+ u = self.field(m)
+
+ u = self.reshapeFields(u)
+
+ 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 = np.empty_like(u)
+ for i, ui in enumerate(u.T): # loop over each column
+ dCdm[:, i] = D * ( sdiag( G * ui ) * ( Av_dm * ( mT_dm * v ) ) )
+
+ solve = Solver(dCdu)
+ Jv = - P * solve.solve(dCdm)
+ 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
+ A = self.createMatrix(m)
+ Av_dm = self.mesh.getFaceMassDeriv()
+ mT_dm = self.modelTransformDeriv(m)
+
+ dCdu = A.T
+ solve = Solver(dCdu)
+
+ w = solve.solve(P.T*v)
+
+ Jtv = 0
+ for i, ui in enumerate(u.T): # loop over each column
+ Jtv += sdiag( G * ui ) * ( D.T * w[:,i] )
+
+ Jtv = - mT_dm.T * ( Av_dm.T * Jtv )
+ return Jtv
+
+
+
+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
+
+
+
+if __name__ == '__main__':
+
+ from SimPEG.regularization import Regularization
+ from SimPEG import inverse
+ import matplotlib.pyplot as plt
+
+ # Create the mesh
+ h1 = np.ones(20)
+ h2 = np.ones(100)
+ mesh = TensorMesh([h1,h2])
+
+ # Create some parameters for the model
+ sig1 = np.log(1)
+ sig2 = np.log(0.01)
+
+ # Create a synthetic model from a block in a half-space
+ p0 = [5, 10]
+ p1 = [15, 50]
+ condVals = [sig1, sig2]
+ mSynth = ModelBuilder.defineBlockConductivity(p0,p1,mesh.gridCC,condVals)
+ plt.colorbar(mesh.plotImage(mSynth))
+ plt.show()
+
+ # 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 = 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)
+ u = synthetic.reshapeFields(u)
+ mesh.plotImage(u[:,10])
+ # plt.show()
+
+ # Now set up the problem to do some minimization
+ problem = DCProblem(mesh)
+ problem.P = P
+ problem.RHS = q
+ problem.dobs = dobs
+ problem.std = dobs*0 + 0.05
+ m0 = mesh.gridCC[:,0]*0+sig2
+
+ 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)
+
+ # Check Derivative
+ derChk = lambda m: [inv.dataObj(m), inv.dataObjDeriv(m)]
+ checkDerivative(derChk, mSynth)
+
+
+
+ print inv.dataObj(m0)
+ print inv.dataObj(mSynth)
+
+ m = inv.run(m0)
+
+ plt.colorbar(mesh.plotImage(m))
+ print m
+ plt.show()
+
+
+
+
+
+
diff --git a/SimPEG/forward/DCProblem/DCProblem.py b/SimPEG/forward/DCProblem/DCProblem.py
deleted file mode 100644
index 1df8897b..00000000
--- a/SimPEG/forward/DCProblem/DCProblem.py
+++ /dev/null
@@ -1,168 +0,0 @@
-from SimPEG.mesh 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))
diff --git a/SimPEG/forward/DCProblem/DCutils.py b/SimPEG/forward/DCProblem/DCutils.py
deleted file mode 100644
index f3445096..00000000
--- a/SimPEG/forward/DCProblem/DCutils.py
+++ /dev/null
@@ -1,29 +0,0 @@
-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
diff --git a/SimPEG/forward/DCProblem/__init__.py b/SimPEG/forward/DCProblem/__init__.py
deleted file mode 100644
index a868cf80..00000000
--- a/SimPEG/forward/DCProblem/__init__.py
+++ /dev/null
@@ -1,2 +0,0 @@
-from DCProblem import *
-from DCutils import *
diff --git a/SimPEG/forward/LinearProblem.py b/SimPEG/forward/LinearProblem.py
new file mode 100644
index 00000000..d30a5b4d
--- /dev/null
+++ b/SimPEG/forward/LinearProblem.py
@@ -0,0 +1,89 @@
+import numpy as np
+from SimPEG.mesh import TensorMesh
+from SimPEG.forward import Problem
+from SimPEG.regularization import Regularization
+from SimPEG.inverse import *
+import matplotlib.pyplot as plt
+
+
+class LinearProblem(Problem):
+ """docstring for LinearProblem"""
+
+ def dpred(self, m, u=None):
+ return self.G.dot(m)
+
+ def J(self, m, v, u=None):
+ return G.dot(v)
+
+ def Jt(self, m, v, u=None):
+ return G.T.dot(v)
+
+if __name__ == '__main__':
+ N = 100
+ h = np.ones(N)/N
+ M = TensorMesh([h])
+
+ nk = 20
+ jk = np.linspace(1.,20.,nk)
+ p = -0.25
+ q = 0.25
+
+
+
+ g = lambda k: np.exp(p*jk[k]*M.vectorCCx)*np.cos(2*np.pi*q*jk[k]*M.vectorCCx)
+
+ G = np.empty((nk, M.nC))
+
+ for i in range(nk):
+ G[i,:] = g(i)
+
+
+
+ plt.figure(1)
+ for i in range(nk):
+ plt.plot(G[i,:])
+
+
+ m_true = np.zeros(M.nC)
+ m_true[M.vectorCCx > 0.3] = 1.
+ m_true[M.vectorCCx > 0.45] = -0.5
+ m_true[M.vectorCCx > 0.6] = 0
+
+
+ d_true = G.dot(m_true)
+ noise = 0.1 * np.random.rand(d_true.size)
+
+ d_obs = d_true + noise
+
+ # plt.figure(3)
+ # plt.plot(d_true,'-o')
+ # plt.plot(d_obs,'r-o')
+
+
+
+
+
+ prob = LinearProblem(M)
+ prob.G = G
+ prob.dobs = d_obs
+ prob.std = np.ones_like(d_obs)*0.1
+
+ reg = Regularization(M)
+
+ opt = InexactGaussNewton(maxIter=20)
+
+ inv = Inversion(prob,reg,opt,beta0=1e-4)
+
+ m0 = np.zeros_like(m_true)
+
+ mrec = inv.run(m0)
+
+
+ plt.figure(2)
+
+ plt.plot(M.vectorCCx, m_true, 'b-')
+ plt.plot(M.vectorCCx, mrec, 'r-')
+
+
+
+ plt.show()
diff --git a/SimPEG/forward/ModelTransforms.py b/SimPEG/forward/ModelTransforms.py
new file mode 100644
index 00000000..ea89b974
--- /dev/null
+++ b/SimPEG/forward/ModelTransforms.py
@@ -0,0 +1,49 @@
+import numpy as np
+from SimPEG.utils import mkvc, sdiag
+
+class LogModel(object):
+ """docstring for LogModel"""
+ 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)))
diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py
index 5b716f1f..2e6831f7 100644
--- a/SimPEG/forward/Problem.py
+++ b/SimPEG/forward/Problem.py
@@ -1,5 +1,6 @@
import numpy as np
from SimPEG.utils import mkvc, sdiag
+import scipy.sparse as sp
norm = np.linalg.norm
@@ -49,16 +50,6 @@ class Problem(object):
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):
"""
@@ -72,6 +63,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):
@@ -83,16 +83,35 @@ class Problem(object):
def dobs(self, value):
self._dobs = value
- def evalFunction(self, m, doDerivative=True):
+ def dpred(self, m, u=None):
"""
- :param numpy.array m: model
- :param bool doDerivative: do you want to compute the derivative?
- :rtype: numpy.array
- :return: Jv
- """
- f = self.misfit(m)
+ Predicted data.
- return f, g, H
+ .. math::
+ d_\\text{pred} = Pu(m)
+ """
+ if u is None:
+ u = self.field(m)
+ return self.P*u
+
+ def dataResidual(self, m, u=None):
+ """
+ :param numpy.array m: geophysical model
+ :param numpy.array u: fields
+ :rtype: float
+ :return: data misfit
+
+ The data misfit:
+
+ .. math::
+
+ \mu_\\text{data} = \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.
+ """
+
+ return self.dpred(m, u=u) - self.dobs
def J(self, m, v, u=None):
"""
@@ -131,10 +150,38 @@ class Problem(object):
:rtype: numpy.array
:return: JTv
- Transpose of J
+ Effect of transpose of J on a vector v.
"""
pass
+
+ def J_approx(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
+
+ Approximate effect of J on a vector v
+
+ """
+ return self.J(m, v, u)
+
+ def Jt_approx(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
+
+ Approximate transpose of J*v
+
+ """
+ return self.Jt(m, v, u)
+
def field(self, m):
"""
The field given the model.
@@ -145,17 +192,6 @@ class Problem(object):
"""
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
@@ -168,13 +204,8 @@ class Problem(object):
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))
+ return m
def modelTransformDeriv(self, m):
"""
@@ -184,129 +215,10 @@ class Problem(object):
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)))
+ return sp.eye(m.size)
- 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):
@@ -337,3 +249,6 @@ class SyntheticProblem(object):
eps = np.linalg.norm(mkvc(dobs),2)*1e-5
Wd = 1/(abs(dobs)*std+eps)
return dobs, Wd
+
+
+
diff --git a/SimPEG/forward/__init__.py b/SimPEG/forward/__init__.py
index fe849d41..33c9a6b1 100644
--- a/SimPEG/forward/__init__.py
+++ b/SimPEG/forward/__init__.py
@@ -1,2 +1,4 @@
from Problem import *
import DCProblem
+from LinearProblem import LinearProblem
+import ModelTransforms
diff --git a/SimPEG/inverse/BetaSchedule.py b/SimPEG/inverse/BetaSchedule.py
new file mode 100644
index 00000000..fe197340
--- /dev/null
+++ b/SimPEG/inverse/BetaSchedule.py
@@ -0,0 +1,12 @@
+
+
+class Cooling(object):
+ """Simple Beta Schedule"""
+
+ beta0 = 1.e6
+ beta_coolingFactor = 5.
+
+ def getBeta(self):
+ if self._beta is None:
+ return beta0
+ return self._beta / beta_coolingFactor
diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py
new file mode 100644
index 00000000..3aa8ce58
--- /dev/null
+++ b/SimPEG/inverse/Inversion.py
@@ -0,0 +1,221 @@
+import numpy as np
+import scipy.sparse as sp
+from SimPEG.utils import sdiag, mkvc
+
+class Inversion(object):
+ """docstring for Inversion"""
+
+ maxIter = 10
+ name = 'SimPEG Inversion'
+
+ def __init__(self, prob, reg, opt, **kwargs):
+ self.prob = prob
+ self.reg = reg
+ self.opt = opt
+ self.opt.parent = self
+ self.setKwargs(**kwargs)
+
+ def setKwargs(self, **kwargs):
+ """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])
+ else:
+ raise Exception('%s attr is not recognized' % attr)
+
+ def printInit(self):
+ print "%s %s %s" % ('='*22, self.name, '='*22)
+ print " # beta phi_d phi_m f norm(dJ) #LS"
+ print "%s" % '-'*62
+
+ def printIter(self):
+ print "%3d %1.2e %1.2e %1.2e %1.2e %1.2e %3d" % (self.opt._iter, self._beta, self._phi_d_last, self._phi_m_last, self.opt.f, np.linalg.norm(self.opt.g), self.opt._iterLS)
+
+ @property
+ def Wd(self):
+ """
+ Standard deviation weighting matrix.
+ """
+ if getattr(self,'_Wd',None) is None:
+ eps = np.linalg.norm(mkvc(self.prob.dobs),2)*1e-5
+ self._Wd = 1/(abs(self.prob.dobs)*self.prob.std+eps)
+ return self._Wd
+
+ @property
+ def phi_d_target(self):
+ """
+ target for phi_d
+
+ By default this is the number of data.
+
+ Note that we do not set the target if it is None, but we return the default value.
+ """
+ if getattr(self, '_phi_d_target', None) is None:
+ return self.prob.dobs.size #
+ return self._phi_d_target
+ @phi_d_target.setter
+ def phi_d_target(self, value):
+ self._phi_d_target = value
+
+ def run(self, m0):
+ m = m0
+ self._iter = 0
+ self._beta = None
+ while True:
+ self._beta = self.getBeta()
+ m = self.opt.minimize(self.evalFunction,m)
+ if self.stoppingCriteria(): break
+ self._iter += 1
+ return m
+
+ beta0 = 1.e2
+ beta_coolingFactor = 5.
+
+ def getBeta(self):
+ if self._beta is None:
+ return self.beta0
+ return self._beta / self.beta_coolingFactor
+
+ def stoppingCriteria(self):
+ self._STOP = np.zeros(2,dtype=bool)
+ self._STOP[0] = self._iter >= self.maxIter
+ self._STOP[1] = self._phi_d_last <= self.phi_d_target
+ return np.any(self._STOP)
+
+
+ def evalFunction(self, m, return_g=True, return_H=True):
+
+ u = self.prob.field(m)
+ phi_d = self.dataObj(m, u)
+ phi_m = self.reg.modelObj(m)
+
+ self._phi_d_last = phi_d
+ self._phi_m_last = phi_m
+
+ f = phi_d + self._beta * phi_m
+
+ out = (f,)
+ if return_g:
+ phi_dDeriv = self.dataObjDeriv(m, u=u)
+ phi_mDeriv = self.reg.modelObjDeriv(m)
+
+ g = phi_dDeriv + self._beta * phi_mDeriv
+ out += (g,)
+
+ if return_H:
+ def H_fun(v):
+ phi_d2Deriv = self.dataObj2Deriv(m, v, u=u)
+ phi_m2Deriv = self.reg.modelObj2Deriv(m)*v
+
+ return phi_d2Deriv + self._beta * phi_m2Deriv
+
+ operator = sp.linalg.LinearOperator( (m.size, m.size), H_fun, dtype=float )
+ out += (operator,)
+ return out
+
+
+ def dataObj(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.
+ """
+ # TODO: ensure that this is a data is vector and Wd is a matrix.
+ R = self.Wd*self.prob.dataResidual(m, u=u)
+ R = mkvc(R)
+ return 0.5*np.vdot(R, R)
+
+ def dataObjDeriv(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.prob.field(m)
+
+ R = self.Wd*self.prob.dataResidual(m, u=u)
+
+ dmisfit = self.prob.Jt(m, self.Wd * R, u=u)
+
+ return dmisfit
+
+ def dataObj2Deriv(self, m, v, 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.prob.field(m)
+
+ R = self.Wd*self.prob.dataResidual(m, u=u)
+
+ # TODO: abstract to different norms a little cleaner.
+ # \/ it goes here. in l2 it is the identity.
+ dmisfit = self.prob.Jt_approx(m, self.Wd * self.Wd * self.prob.J_approx(m, v, u=u), u=u)
+
+ return dmisfit
+
diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/inverse/Optimize.py
index eee18b16..6221246f 100644
--- a/SimPEG/inverse/Optimize.py
+++ b/SimPEG/inverse/Optimize.py
@@ -2,54 +2,154 @@ import numpy as np
import matplotlib.pyplot as plt
from SimPEG.utils import mkvc, sdiag
norm = np.linalg.norm
+import scipy.sparse as sp
+from SimPEG import Solver
+
+try:
+ from pubsub import pub
+ doPub = True
+except Exception, e:
+ print 'Warning: you may not have the required pubsub installed, use pypubsub. You will not be able to listen to events.'
+ doPub = False
+
class Minimize(object):
- """docstring for Minimize"""
+ """
+
+ Minimize is a general class for derivative based optimization.
+
+
+ """
name = "GeneralOptimizationAlgorithm"
maxIter = 20
maxIterLS = 10
+ maxStep = np.inf
LSreduction = 1e-4
LSshorten = 0.5
- tolF = 1e-4
- tolX = 1e-4
- tolG = 1e-4
- eps = 1e-16
+ tolF = 1e-1
+ tolX = 1e-1
+ tolG = 1e-1
+ eps = 1e-5
- def __init__(self, problem, **kwargs):
- self.problem = problem
+ def __init__(self, **kwargs):
+ self._id = int(np.random.rand()*1e6) # create a unique identifier to this program to be used in pubsub
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])
else:
raise Exception('%s attr is not recognized' % attr)
- def minimize(self, x0):
+ def minimize(self, evalFunction, x0):
+ """
+ Minimizes the function (evalFunction) starting at the location x0.
+ :param def evalFunction: function handle that evaluates: f, g, H = F(x)
+ :param numpy.ndarray x0: starting location
+ :rtype: numpy.ndarray
+ :return: x, the last iterate of the optimization algorithm
+
+ evalFunction is a function handle::
+
+ (f[, g][, H]) = evalFunction(x, return_g=False, return_H=False )
+
+
+ Events are fired with the following inputs via pypubsub::
+
+ Minimize.printInit (minimize)
+ Minimize.evalFunction (minimize, f, g, H)
+ Minimize.printIter (minimize)
+ Minimize.searchDirection (minimize, p)
+ Minimize.scaleSearchDirection (minimize, p)
+ Minimize.modifySearchDirection (minimize, xt, passLS)
+ Minimize.endIteration (minimize, xt)
+ Minimize.printDone (minimize)
+
+ To hook into one of these events (must have pypubsub installed)::
+
+ from pubsub import pub
+ def listener(minimize,p):
+ print 'The search direction is: ', p
+ pub.subscribe(listener, 'Minimize.searchDirection')
+
+ You can use pubsub communication to debug your code, it is not used internally.
+
+
+ The algorithm for general minimization is as follows::
+
+ startup(x0)
+ printInit()
+
+ while True:
+ f, g, H = evalFunction(xc)
+ printIter()
+ if stoppingCriteria(): break
+ p = findSearchDirection()
+ p = scaleSearchDirection(p)
+ xt, passLS = modifySearchDirection(p)
+ if not passLS:
+ xt, caught = modifySearchDirectionBreak(p)
+ if not caught: return xc
+ doEndIteration(xt)
+
+ printDone()
+ return xc
+ """
+ self.evalFunction = evalFunction
self.startup(x0)
self.printInit()
while True:
- self.f, self.g, self.H = self.evalFunction(self.xc)
+ self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True)
+ if doPub: pub.sendMessage('Minimize.evalFunction', minimize=self, f=self.f, g=self.g, H=self.H)
self.printIter()
if self.stoppingCriteria(): break
p = self.findSearchDirection()
- xt, passLS = self.linesearch(p)
+ if doPub: pub.sendMessage('Minimize.searchDirection', minimize=self, p=p)
+ p = self.scaleSearchDirection(p)
+ if doPub: pub.sendMessage('Minimize.scaleSearchDirection', minimize=self, p=p)
+ xt, passLS = self.modifySearchDirection(p)
+ if doPub: pub.sendMessage('Minimize.modifySearchDirection', minimize=self, xt=xt, passLS=passLS)
if not passLS:
- xt = self.linesearchBreak(p)
+ xt, caught = self.modifySearchDirectionBreak(p)
+ if not caught: return self.xc
self.doEndIteration(xt)
+ if doPub: pub.sendMessage('Minimize.endIteration', minimize=self, xt=xt)
self.printDone()
return self.xc
+ @property
+ def parent(self):
+ """
+ This is the parent of the optimization routine.
+ """
+ return getattr(self, '_parent', None)
+ @parent.setter
+ def parent(self, value):
+ self._parent = value
+
def startup(self, x0):
+ """
+ **startup** is called at the start of any new minimize call.
+
+ This will set::
+
+ x0 = x0
+ xc = x0
+ _iter = _iterLS = 0
+
+ :param numpy.ndarray x0: initial x
+ :rtype: None
+ :return: None
+ """
self._iter = 0
self._iterLS = 0
self._STOP = np.zeros((5,1),dtype=bool)
@@ -59,29 +159,57 @@ class Minimize(object):
self.xOld = x0
def printInit(self):
+ """
+ **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.
+
+ """
+ if doPub: pub.sendMessage('Minimize.printInit', minimize=self)
+ if self.parent is not None and hasattr(self.parent, 'printInit'):
+ self.parent.printInit()
+ return
print "%s %s %s" % ('='*22, self.name, '='*22)
print "iter\tJc\t\tnorm(dJ)\tLS"
print "%s" % '-'*57
def printIter(self):
+ """
+ **printIter** is called directly after function evaluations.
+
+ If there is a parent object, printIter will check for a
+ parent.printIter function and call that.
+
+ """
+ if doPub: pub.sendMessage('Minimize.printIter', minimize=self)
+ if self.parent is not None and hasattr(self.parent, 'printIter'):
+ self.parent.printIter()
+ return
print "%3d\t%1.2e\t%1.2e\t%d" % (self._iter, self.f, norm(self.g), self._iterLS)
def printDone(self):
+ """
+ **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.
+
+ """
+ if doPub: pub.sendMessage('Minimize.printDone', minimize=self)
+ if self.parent is not None and hasattr(self.parent, 'printDone'):
+ self.parent.printDone()
+ return
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)))
+ # TODO: put controls on gradient value, min model update, and function value
+ if self._iter > 0:
+ 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
@@ -94,14 +222,87 @@ class Minimize(object):
self._STOP[4] = self._iter >= self.maxIter
return all(self._STOP[0:3]) | any(self._STOP[3:])
- def linesearch(self, p):
+ def projection(self, p):
+ """
+ projects the search direction.
+
+ by default, no projection is applied.
+
+ :param numpy.ndarray p: searchDirection
+ :rtype: numpy.ndarray
+ :return: p, projected search direction
+ """
+ return p
+
+ def findSearchDirection(self):
+ """
+ **findSearchDirection** should return an approximation of:
+
+ .. math::
+
+ H p = - g
+
+ Where you are solving for the search direction, p
+
+ The default is:
+
+ .. math::
+
+ H = I
+
+ p = - g
+
+ And corresponds to SteepestDescent.
+
+ The latest function evaluations are present in::
+
+ self.f, self.g, self.H
+
+ :rtype: numpy.ndarray
+ :return: p, Search Direction
+ """
+ return -self.g
+
+ def scaleSearchDirection(self, p):
+ """
+ **scaleSearchDirection** should scale the search direction if appropriate.
+
+ Set the parameter **maxStep** in the minimize object, to scale back the gradient to a maximum size.
+
+ :param numpy.ndarray p: searchDirection
+ :rtype: numpy.ndarray
+ :return: p, Scaled Search Direction
+ """
+
+ if self.maxStep < np.abs(p.max()):
+ p = self.maxStep*p/np.abs(p.max())
+ return p
+
+ def modifySearchDirection(self, p):
+ """
+ **modifySearchDirection** changes the search direction based on some sort of linesearch or trust-region criteria.
+
+ By default, an Armijo backtracking linesearch is preformed with the following parameters:
+
+ * maxIterLS, the maximum number of linesearch iterations
+ * LSreduction, the expected reduction expected, default: 1e-4
+ * LSshorten, how much the step is reduced, default: 0.5
+
+ If the linesearch is completed, and a descent direction is found, passLS is returned as True.
+
+ Else, a modifySearchDirectionBreak call is preformed.
+
+ :param numpy.ndarray p: searchDirection
+ :rtype: numpy.ndarray,bool
+ :return: (xt, passLS)
+ """
# 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)
+ 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
iterLS += 1
@@ -110,10 +311,37 @@ class Minimize(object):
self._iterLS = iterLS
return xt, iterLS < self.maxIterLS
- def linesearchBreak(self, p):
- raise Exception('The linesearch got broken. Boo.')
+ def modifySearchDirectionBreak(self, p):
+ """
+ Code is called if modifySearchDirection fails
+ to find a descent direction.
+
+ The search direction is passed as input and
+ this function must pass back both a new searchDirection,
+ and if the searchDirection break has been caught.
+
+ By default, no additional work is done, and the
+ evalFunction returns a False indicating the break was not caught.
+
+ :param numpy.ndarray p: searchDirection
+ :rtype: numpy.ndarray,bool
+ :return: (xt, breakCaught)
+ """
+ print 'The linesearch got broken. Boo.'
+ return p, False
def doEndIteration(self, xt):
+ """
+ **doEndIteration** is called at the end of each minimize iteration.
+
+ By default, function values and x locations are shuffled to store 1 past iteration in memory.
+
+ self.xc must be updated in this code.
+
+ :param numpy.ndarray xt: tested new iterate that ensures a descent direction.
+ :rtype: None
+ :return: None
+ """
# store old values
self.fOld = self.f
self.xOld, self.xc = self.xc, xt
@@ -123,7 +351,19 @@ class Minimize(object):
class GaussNewton(Minimize):
name = 'GaussNewton'
def findSearchDirection(self):
- return np.linalg.solve(self.H,-self.g)
+ return Solver(self.H).solve(-self.g)
+
+
+class InexactGaussNewton(Minimize):
+ name = 'InexactGaussNewton'
+
+ maxIterCG = 10
+ tolCG = 1e-5
+
+ def findSearchDirection(self):
+ # TODO: use BFGS as a preconditioner or gauss sidel of the WtW or solve WtW directly
+ p, info = sp.linalg.cg(self.H, -self.g, tol=self.tolCG, maxiter=self.maxIterCG)
+ return p
class SteepestDescent(Minimize):
@@ -133,18 +373,15 @@ class SteepestDescent(Minimize):
if __name__ == '__main__':
from SimPEG.tests import Rosenbrock, checkDerivative
+ import matplotlib.pyplot as plt
x0 = np.array([2.6, 3.7])
checkDerivative(Rosenbrock, x0, plotIt=False)
- xOpt = GaussNewton(Rosenbrock, maxIter=20).minimize(x0)
+
+ def listener1(minimize,p):
+ print 'hi: ', p
+ if doPub: pub.subscribe(listener1, 'Minimize.searchDirection')
+
+ xOpt = GaussNewton(maxIter=20,tolF=1e-10,tolX=1e-10,tolG=1e-10).minimize(Rosenbrock,x0)
print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1])
- xOpt = SteepestDescent(Rosenbrock, maxIter=20, maxIterLS=15).minimize(x0)
+ 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/inverse/__init__.py b/SimPEG/inverse/__init__.py
index b2a5e506..14bddce7 100644
--- a/SimPEG/inverse/__init__.py
+++ b/SimPEG/inverse/__init__.py
@@ -1 +1,3 @@
from Optimize import *
+from Inversion import *
+import BetaSchedule
diff --git a/SimPEG/mesh/BaseMesh.py b/SimPEG/mesh/BaseMesh.py
index 29c43dd1..6a9a8032 100644
--- a/SimPEG/mesh/BaseMesh.py
+++ b/SimPEG/mesh/BaseMesh.py
@@ -2,6 +2,7 @@ import numpy as np
from SimPEG.utils import mkvc
+
class BaseMesh(object):
"""
BaseMesh does all the counting you don't want to do.
@@ -216,6 +217,12 @@ class BaseMesh(object):
:rtype: int
:return: nC
+
+ .. plot::
+
+ from SimPEG.mesh import TensorMesh
+ import numpy as np
+ TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(centers=True,showIt=True)
"""
fget = lambda self: np.prod(self.n)
return locals()
@@ -270,6 +277,12 @@ class BaseMesh(object):
:rtype: int
:return: nN
+
+ .. plot::
+
+ from SimPEG.mesh import TensorMesh
+ import numpy as np
+ TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(nodes=True,showIt=True)
"""
fget = lambda self: np.prod(self.n + 1)
return locals()
@@ -324,6 +337,12 @@ class BaseMesh(object):
:rtype: numpy.array (dim, )
:return: [prod(nEx), prod(nEy), prod(nEz)]
+
+ .. plot::
+
+ from SimPEG.mesh import TensorMesh
+ import numpy as np
+ TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(edges=True,showIt=True)
"""
fget = lambda self: np.array([np.prod(x) for x in [self.nEx, self.nEy, self.nEz] if not x is None])
return locals()
@@ -378,6 +397,12 @@ class BaseMesh(object):
:rtype: numpy.array (dim, )
:return: [prod(nFx), prod(nFy), prod(nFz)]
+
+ .. plot::
+
+ from SimPEG.mesh import TensorMesh
+ import numpy as np
+ TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(faces=True,showIt=True)
"""
fget = lambda self: np.array([np.prod(x) for x in [self.nFx, self.nFy, self.nFz] if not x is None])
return locals()
diff --git a/SimPEG/mesh/Cyl1DMesh.py b/SimPEG/mesh/Cyl1DMesh.py
index 915bf0ef..93b82b25 100644
--- a/SimPEG/mesh/Cyl1DMesh.py
+++ b/SimPEG/mesh/Cyl1DMesh.py
@@ -5,8 +5,8 @@ from SimPEG.utils import mkvc, ndgrid, sdiag
class Cyl1DMesh(object):
"""
- Cyl1DMesh is a mesh class for cylindrically symmetric 1D problems
- """
+ Cyl1DMesh is a mesh class for cylindrically symmetric 1D problems
+ """
_meshType = 'CYL1D'
@@ -20,7 +20,7 @@ class Cyl1DMesh(object):
assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i)
# Ensure h contains 1D vectors
- self._h = [mkvc(x) for x in h]
+ self._h = [mkvc(x.astype(float)) for x in h]
if z0 is None:
z0 = 0
@@ -146,7 +146,7 @@ class Cyl1DMesh(object):
def vectorCCz():
doc = "Cell centered grid vector (1D) in the z direction"
- fget = lambda self: self.hz.cumsum() - self.hz/2 + self._z0
+ fget = lambda self: self.hz.cumsum() - self.hz/2 + self._z0
return locals()
vectorCCz = property(**vectorCCz())
@@ -177,7 +177,7 @@ class Cyl1DMesh(object):
self._gridFr = ndgrid([self.vectorNr, self.vectorCCz])
return self._gridFr
return locals()
- _gridFr = None
+ _gridFr = None
gridFr = property(**gridFr())
def gridFz():
@@ -187,7 +187,7 @@ class Cyl1DMesh(object):
self._gridFz = ndgrid([self.vectorCCr, self.vectorNz])
return self._gridFz
return locals()
- _gridFz = None
+ _gridFz = None
gridFz = property(**gridFz())
####################################################
@@ -350,23 +350,23 @@ class Cyl1DMesh(object):
np.all(loc[:,1]<=self.vectorNz.max()), \
"Points outside of mesh"
-
+
if locType=='fz':
Q = sp.lil_matrix((loc.shape[0], self.nF), dtype=float)
for i, iloc in enumerate(loc):
# Point is on a z-interface
- if np.any(np.abs(self.vectorNz-iloc[1])<0.001):
+ if np.any(np.abs(self.vectorNz-iloc[1])<0.001):
dFz = self.gridFz-iloc #Distance to z faces
dFz[dFz[:,0]>0,:] = np.inf #Looking for next face to the left...
indL = np.argmin(np.sum(dFz**2, axis=1)) #Closest one
if self.gridFz[indL,0] == self.vectorCCr.max(): #Point in outer half cell (linear extrapolation)
- zFL = self.gridFz[indL,:]
- zFLL = self.gridFz[indL-1,:]
+ zFL = self.gridFz[indL,:]
+ zFLL = self.gridFz[indL-1,:]
Q[i, indL+self.nFr] = (iloc[0] - zFLL[0])/(zFL[0] - zFLL[0])
Q[i, indL+self.nFr-1] = -(iloc[0] - zFL[0])/(zFL[0] - zFLL[0])
else:
- zFL = self.gridFz[indL,:]
+ zFL = self.gridFz[indL,:]
zFR = self.gridFz[indL+1,:]
Q[i,indL+self.nFr] = (zFR[0] - iloc[0])/(zFR[0] - zFL[0])
Q[i,indL+self.nFr+1] = (iloc[0] - zFL[0])/(zFR[0] - zFL[0])
@@ -400,7 +400,7 @@ class Cyl1DMesh(object):
Q[i, indAL+self.nFr-1] = -(dzB/DZ)*(drL/DR)
Q[i, indAL+self.nFr] = (dzB/DZ)*(drLL/DR)
else:
- indBR = indBL+1 # Face below and to the right
+ indBR = indBL+1 # Face below and to the right
indAR = indAL + 1 # Face above and to the right
zF_BR = self.gridFz[indBR,:]
diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/mesh/DiffOperators.py
index 110fe9bd..598b392a 100644
--- a/SimPEG/mesh/DiffOperators.py
+++ b/SimPEG/mesh/DiffOperators.py
@@ -161,6 +161,68 @@ class DiffOperators(object):
_cellGrad = None
cellGrad = property(**cellGrad())
+ def cellGradx():
+ doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
+
+ def fget(self):
+ if getattr(self, '_cellGradx', None) is None:
+ BC = ['neumann', 'neumann']
+ n = self.n
+ if(self.dim == 1):
+ G1 = ddxCellGrad(n[0], BC)
+ elif(self.dim == 2):
+ G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC))
+ elif(self.dim == 3):
+ G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC))
+ # Compute areas of cell faces & volumes
+ S = self.r(self.area, 'F','Fx', 'V')
+ V = self.vol
+ self._cellGradx = sdiag(S)*G1*sdiag(1/V)
+ return self._cellGradx
+ return locals()
+ cellGradx = property(**cellGradx())
+
+
+ def cellGrady():
+ doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
+ def fget(self):
+ if self.dim < 2:
+ return None
+ if getattr(self, '_cellGrady', None) is None:
+ BC = ['neumann', 'neumann']
+ n = self.n
+ if(self.dim == 2):
+ G2 = sp.kron(ddxCellGrad(n[1], BC), speye(n[0]))
+ elif(self.dim == 3):
+ G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0]))
+ # Compute areas of cell faces & volumes
+ S = self.r(self.area, 'F','Fy', 'V')
+ V = self.vol
+ self._cellGrady = sdiag(S)*G2*sdiag(1/V)
+ return self._cellGrady
+ return locals()
+ cellGrady = property(**cellGrady())
+
+
+
+ def cellGradz():
+ doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
+ def fget(self):
+ if self.dim < 3:
+ return None
+ if getattr(self, '_cellGradz', None) is None:
+ BC = ['neumann', 'neumann']
+ n = self.n
+ G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0]))
+ # Compute areas of cell faces & volumes
+ S = self.r(self.area, 'F','Fz', 'V')
+ V = self.vol
+ self._cellGradz = sdiag(S)*G3*sdiag(1/V)
+ return self._cellGradz
+ return locals()
+ cellGradz = property(**cellGradz())
+
+
def edgeCurl():
doc = "Construct the 3D curl operator."
diff --git a/SimPEG/mesh/InnerProducts.py b/SimPEG/mesh/InnerProducts.py
index 9fe84ac3..f0ac4ab0 100644
--- a/SimPEG/mesh/InnerProducts.py
+++ b/SimPEG/mesh/InnerProducts.py
@@ -81,9 +81,9 @@ class InnerProducts(object):
def getFaceInnerProduct(self, mu=None, returnP=False):
"""Wrapper function,
- :py:func:`SimPEG.InnerProducts.getEdgeInnerProduct`
+ :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getEdgeInnerProduct`
- :py:func:`SimPEG.InnerProducts.getEdgeInnerProduct2D`
+ :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getEdgeInnerProduct2D`
"""
if self.dim == 2:
return getFaceInnerProduct2D(self, mu, returnP)
@@ -93,9 +93,9 @@ class InnerProducts(object):
def getEdgeInnerProduct(self, sigma=None, returnP=False):
"""Wrapper function,
- :py:func:`SimPEG.InnerProducts.getFaceInnerProduct`
+ :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getFaceInnerProduct`
- :py:func:`SimPEG.InnerProducts.getFaceInnerProduct2D`
+ :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getFaceInnerProduct2D`
"""
if self.dim == 2:
return getEdgeInnerProduct2D(self, sigma, returnP)
diff --git a/SimPEG/mesh/LogicallyOrthogonalMesh.py b/SimPEG/mesh/LogicallyOrthogonalMesh.py
index 7b1e1bca..b510a754 100644
--- a/SimPEG/mesh/LogicallyOrthogonalMesh.py
+++ b/SimPEG/mesh/LogicallyOrthogonalMesh.py
@@ -38,7 +38,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
# Save nodes to private variable _gridN as vectors
self._gridN = np.ones((nodes[0].size, self.dim))
for i, node_i in enumerate(nodes):
- self._gridN[:, i] = mkvc(node_i)
+ self._gridN[:, i] = mkvc(node_i.astype(float))
def gridCC():
doc = "Cell-centered grid."
diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/mesh/TensorMesh.py
index a3329c14..6cc5d8a6 100644
--- a/SimPEG/mesh/TensorMesh.py
+++ b/SimPEG/mesh/TensorMesh.py
@@ -1,9 +1,10 @@
import numpy as np
+import scipy.sparse as sp
from BaseMesh import BaseMesh
from TensorView import TensorView
from DiffOperators import DiffOperators
from InnerProducts import InnerProducts
-from SimPEG.utils import ndgrid, mkvc
+from SimPEG.utils import ndgrid, mkvc, spzeros, interpmat
class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
@@ -38,7 +39,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i)
# Ensure h contains 1D vectors
- self._h = [mkvc(x) for x in h]
+ self._h = [mkvc(x.astype(float)) for x in h]
def __str__(self):
outStr = ' ---- {0:d}-D TensorMesh ---- '.format(self.dim)
@@ -156,7 +157,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridCC is None:
- self._gridCC = ndgrid([x for x in [self.vectorCCx, self.vectorCCy, self.vectorCCz] if not x is None])
+ self._gridCC = ndgrid(self.getTensor('CC'))
return self._gridCC
return locals()
_gridCC = None # Store grid by default
@@ -167,7 +168,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridN is None:
- self._gridN = ndgrid([x for x in [self.vectorNx, self.vectorNy, self.vectorNz] if not x is None])
+ self._gridN = ndgrid(self.getTensor('N'))
return self._gridN
return locals()
_gridN = None # Store grid by default
@@ -178,7 +179,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridFx is None:
- self._gridFx = ndgrid([x for x in [self.vectorNx, self.vectorCCy, self.vectorCCz] if not x is None])
+ self._gridFx = ndgrid(self.getTensor('Fx'))
return self._gridFx
return locals()
_gridFx = None # Store grid by default
@@ -189,7 +190,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridFy is None and self.dim > 1:
- self._gridFy = ndgrid([x for x in [self.vectorCCx, self.vectorNy, self.vectorCCz] if not x is None])
+ self._gridFy = ndgrid(self.getTensor('Fy'))
return self._gridFy
return locals()
_gridFy = None # Store grid by default
@@ -200,7 +201,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridFz is None and self.dim > 2:
- self._gridFz = ndgrid([x for x in [self.vectorCCx, self.vectorCCy, self.vectorNz] if not x is None])
+ self._gridFz = ndgrid(self.getTensor('Fz'))
return self._gridFz
return locals()
_gridFz = None # Store grid by default
@@ -211,7 +212,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridEx is None:
- self._gridEx = ndgrid([x for x in [self.vectorCCx, self.vectorNy, self.vectorNz] if not x is None])
+ self._gridEx = ndgrid(self.getTensor('Ex'))
return self._gridEx
return locals()
_gridEx = None # Store grid by default
@@ -222,7 +223,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridEy is None and self.dim > 1:
- self._gridEy = ndgrid([x for x in [self.vectorNx, self.vectorCCy, self.vectorNz] if not x is None])
+ self._gridEy = ndgrid(self.getTensor('Ey'))
return self._gridEy
return locals()
_gridEy = None # Store grid by default
@@ -233,7 +234,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridEz is None and self.dim > 2:
- self._gridEz = ndgrid([x for x in [self.vectorNx, self.vectorNy, self.vectorCCz] if not x is None])
+ self._gridEz = ndgrid(self.getTensor('Ez'))
return self._gridEz
return locals()
_gridEz = None # Store grid by default
@@ -312,6 +313,98 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
_edge = None
edge = property(**edge())
+ # --------------- Methods ---------------------
+
+ def getTensor(self, locType):
+ """ Returns a tensor list.
+
+ :param str locType: What tensor (see below)
+ :rtype: list
+ :return: list of the tensors that make up the mesh.
+
+ locType can be::
+
+ 'Ex' -> x-component of field defined on edges
+ 'Ey' -> y-component of field defined on edges
+ 'Ez' -> z-component of field defined on edges
+ 'Fx' -> x-component of field defined on faces
+ 'Fy' -> y-component of field defined on faces
+ 'Fz' -> z-component of field defined on faces
+ 'N' -> scalar field defined on nodes
+ 'CC' -> scalar field defined on cell centers
+ """
+
+ if locType is 'Fx':
+ ten = [self.vectorNx , self.vectorCCy, self.vectorCCz]
+ elif locType is 'Fy':
+ ten = [self.vectorCCx, self.vectorNy , self.vectorCCz]
+ elif locType is 'Fz':
+ ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ]
+ elif locType is 'Ex':
+ ten = [self.vectorCCx, self.vectorNy , self.vectorNz ]
+ elif locType is 'Ey':
+ ten = [self.vectorNx , self.vectorCCy, self.vectorNz ]
+ elif locType is 'Ez':
+ ten = [self.vectorNx , self.vectorNy , self.vectorCCz]
+ elif locType is 'CC':
+ ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz]
+ elif locType is 'N':
+ ten = [self.vectorNx , self.vectorNy , self.vectorNz ]
+
+ return [t for t in ten if t is not None]
+
+
+ def isInside(self, pts):
+ """
+ Determines if a set of points are inside a mesh.
+
+ :param numpy.ndarray pts: Location of points to test
+ :rtype numpy.ndarray
+ :return inside, numpy array of booleans
+ """
+
+ pts = np.atleast_2d(pts)
+ inside = (pts[:,0] >= self.vectorNx.min()) & (pts[:,0] <= self.vectorNx.max())
+ if self.dim > 1:
+ inside = inside & ((pts[:,1] >= self.vectorNy.min()) & (pts[:,1] <= self.vectorNy.max()))
+ if self.dim > 2:
+ inside = inside & ((pts[:,2] >= self.vectorNz.min()) & (pts[:,2] <= self.vectorNz.max()))
+ return inside
+
+ def getInterpolationMat(self, loc, locType):
+ """ Produces interpolation matrix
+
+ :param numpy.ndarray loc: Location of points to interpolate to
+ :param str locType: What to interpolate (see below)
+ :rtype: scipy.sparse.csr.csr_matrix
+ :return: M, the interpolation matrix
+
+ locType can be::
+
+ 'Ex' -> x-component of field defined on edges
+ 'Ey' -> y-component of field defined on edges
+ 'Ez' -> z-component of field defined on edges
+ 'Fx' -> x-component of field defined on faces
+ 'Fy' -> y-component of field defined on faces
+ 'Fz' -> z-component of field defined on faces
+ 'N' -> scalar field defined on nodes
+ 'CC' -> scalar field defined on cell centers
+ """
+
+ loc = np.atleast_2d(loc)
+ assert np.all(self.isInside(loc)), "Points outside of mesh"
+
+ ind = 0 if 'x' in locType else 1 if 'y' in locType else 2 if 'z' in locType else -1
+ if locType in ['Fx','Fy','Fz','Ex','Ey','Ez'] and self.dim >= ind:
+ nF_nE = self.nF if 'F' in locType else self.nE
+ components = [spzeros(loc.shape[0], n) for n in nF_nE]
+ components[ind] = interpmat(loc, *self.getTensor(locType))
+ Q = sp.hstack(components)
+ elif locType in ['CC', 'N']:
+ Q = interpmat(loc, *self.getTensor(locType))
+ else:
+ raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim))
+ return Q
if __name__ == '__main__':
print('Welcome to tensor mesh!')
@@ -327,5 +420,6 @@ if __name__ == '__main__':
h = h[:testDim]
M = TensorMesh(h)
+ print M
xn = M.plotGrid()
diff --git a/SimPEG/mesh/TensorView.py b/SimPEG/mesh/TensorView.py
index 687a520d..0b9ff7b3 100644
--- a/SimPEG/mesh/TensorView.py
+++ b/SimPEG/mesh/TensorView.py
@@ -267,6 +267,9 @@ class TensorView(object):
if faces:
ax.plot(xs1[:, 0], xs1[:, 1], 'g>')
ax.plot(xs2[:, 0], xs2[:, 1], 'g^')
+ if edges:
+ ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'c>')
+ ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'c^')
# Plot the grid lines
if lines:
diff --git a/SimPEG/regularization/Regularization.py b/SimPEG/regularization/Regularization.py
new file mode 100644
index 00000000..6f5970e6
--- /dev/null
+++ b/SimPEG/regularization/Regularization.py
@@ -0,0 +1,120 @@
+from SimPEG.utils import sdiag
+import numpy as np
+
+class Regularization(object):
+ """docstring for Regularization"""
+
+ @property
+ def mref(self):
+ if getattr(self, '_mref', None) is None:
+ self._mref = np.zeros(self.mesh.nC);
+ return self._mref
+ @mref.setter
+ def mref(self, value):
+ self._mref = value
+
+ @property
+ def Ws(self):
+ if getattr(self,'_Ws', None) is None:
+ self._Ws = sdiag(self.mesh.vol)
+ return self._Ws
+
+ @property
+ def Wx(self):
+ if getattr(self, '_Wx', None) is None:
+ a = self.mesh.r(self.mesh.area,'F','Fx','V')
+ self._Wx = sdiag(a)*self.mesh.cellGradx
+ return self._Wx
+
+ @property
+ def Wy(self):
+ if getattr(self, '_Wy', None) is None:
+ a = self.mesh.r(self.mesh.area,'F','Fy','V')
+ self._Wy = sdiag(a)*self.mesh.cellGrady
+ return self._Wy
+
+ @property
+ def Wz(self):
+ if getattr(self, '_Wz', None) is None:
+ a = self.mesh.r(self.mesh.area,'F','Fz','V')
+ self._Wz = sdiag(a)*self.mesh.cellGradz
+ return self._Wz
+
+
+
+ def __init__(self, mesh):
+ self.mesh = mesh
+ self._Wx = None
+ self._Wy = None
+ self._Wz = None
+ self.alpha_s = 1e-6
+ self.alpha_x = 1
+ self.alpha_y = 1
+ self.alpha_z = 1
+
+ def pnorm(self, r):
+ return 0.5*r.dot(r)
+
+ def modelObj(self, m):
+ mresid = m - self.mref
+
+ mobj = self.alpha_s * self.pnorm( self.Ws * mresid )
+
+ mobj += self.alpha_x * self.pnorm( self.Wx * mresid )
+
+ if self.mesh.dim > 1:
+ mobj += self.alpha_y * self.pnorm( self.Wy * mresid )
+ if self.mesh.dim > 2:
+ mobj += self.alpha_z * self.pnorm( self.Wz * mresid )
+
+ return mobj
+
+ def modelObjDeriv(self, m):
+ """
+
+ In 1D:
+
+ .. math::
+
+ m_{\\text{obj}} = {1 \over 2}\\alpha_s \left\| W_s (m- m_{\\text{ref}})\\right\|^2_2
+ + {1 \over 2}\\alpha_x \left\| W_x (m- m_{\\text{ref}})\\right\|^2_2
+
+ \\frac{ \partial m_{\\text{obj}} }{\partial m} =
+ \\alpha_s W_s^{\\top} W_s (m - m_{\\text{ref}}) +
+ \\alpha_x W_x^{\\top} W_x (m - m_{\\text{ref}})
+
+
+ \\frac{ \partial^2 m_{\\text{obj}} }{\partial m^2} =
+ \\alpha_s W_s^{\\top} W_s +
+ \\alpha_x W_x^{\\top} W_x
+
+ """
+
+ mresid = m - self.mref
+
+ mobjDeriv = self.alpha_s * self.Ws.T * ( self.Ws * mresid)
+
+ mobjDeriv = mobjDeriv + self.alpha_x * self.Wx.T * ( self.Wx * mresid)
+
+ if self.mesh.dim > 1:
+ mobjDeriv = mobjDeriv + self.alpha_y * self.Wy.T * ( self.Wy * mresid)
+ if self.mesh.dim > 2:
+ mobjDeriv = mobjDeriv + self.alpha_z * self.Wz.T * ( self.Wz * mresid)
+
+ return mobjDeriv
+
+
+ def modelObj2Deriv(self, m):
+ mresid = m - self.mref
+
+ mobj2Deriv = self.alpha_s * self.Ws.T * self.Ws
+
+ mobj2Deriv = mobj2Deriv + self.alpha_x * self.Wx.T * self.Wx
+
+ if self.mesh.dim > 1:
+ mobj2Deriv = mobj2Deriv + self.alpha_y * self.Wy.T * self.Wy
+ if self.mesh.dim > 2:
+ mobj2Deriv = mobj2Deriv + self.alpha_z * self.Wz.T * self.Wz
+
+ return mobj2Deriv
+
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
diff --git a/SimPEG/tests/HTMLTestRunner.py b/SimPEG/tests/HTMLTestRunner.py
new file mode 100644
index 00000000..af384971
--- /dev/null
+++ b/SimPEG/tests/HTMLTestRunner.py
@@ -0,0 +1,824 @@
+"""
+A TestRunner for use with the Python unit testing framework. It
+generates a HTML report to show the result at a glance.
+
+The simplest way to use this is to invoke its main method. E.g.
+
+ import unittest
+ import HTMLTestRunner
+
+ ... define your tests ...
+
+ if __name__ == '__main__':
+ HTMLTestRunner.main()
+
+
+For more customization options, instantiates a HTMLTestRunner object.
+HTMLTestRunner is a counterpart to unittest's TextTestRunner. E.g.
+
+ # output to a file
+ fp = file('my_report.html', 'wb')
+ runner = HTMLTestRunner.HTMLTestRunner(
+ stream=fp,
+ title='My unit test',
+ description='This demonstrates the report output by HTMLTestRunner.'
+ )
+
+ # Use an external stylesheet.
+ # See the Template_mixin class for more customizable options
+ runner.STYLESHEET_TMPL = ''
+
+ # run the test
+ runner.run(my_test_suite)
+
+
+------------------------------------------------------------------------
+Copyright (c) 2004-2007, Wai Yip Tung
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+* Redistributions of source code must retain the above copyright notice,
+ this list of conditions and the following disclaimer.
+* Redistributions in binary form must reproduce the above copyright
+ notice, this list of conditions and the following disclaimer in the
+ documentation and/or other materials provided with the distribution.
+* Neither the name Wai Yip Tung nor the names of its contributors may be
+ used to endorse or promote products derived from this software without
+ specific prior written permission.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
+IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
+TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
+PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
+OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
+PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+"""
+
+# URL: http://tungwaiyip.info/software/HTMLTestRunner.html
+
+__author__ = "Wai Yip Tung"
+__version__ = "0.8.2"
+
+
+"""
+Change History
+
+Version 0.8.2
+* Show output inline instead of popup window (Viorel Lupu).
+
+Version in 0.8.1
+* Validated XHTML (Wolfgang Borgert).
+* Added description of test classes and test cases.
+
+Version in 0.8.0
+* Define Template_mixin class for customization.
+* Workaround a IE 6 bug that it does not treat
+
+%(heading)s
+%(report)s
+%(ending)s
+
+
+