diff --git a/.gitignore b/.gitignore
index 5c39733b..7678b15c 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,4 +3,3 @@
*.sublime-project
*.sublime-workspace
docs/_build/
-myNotebooks/*
diff --git a/SimPEG/GCEtools/gceStartup.txt b/GCEtools/gceStartup.txt
similarity index 100%
rename from SimPEG/GCEtools/gceStartup.txt
rename to GCEtools/gceStartup.txt
diff --git a/SimPEG/GCEtools/startup.sh b/GCEtools/startup.sh
similarity index 100%
rename from SimPEG/GCEtools/startup.sh
rename to GCEtools/startup.sh
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 00000000..f94a23fd
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,20 @@
+The MIT License (MIT)
+
+Copyright (c) 2013-2014 SimPEG Developers
+
+Permission is hereby granted, free of charge, to any person obtaining a copy of
+this software and associated documentation files (the "Software"), to deal in
+the Software without restriction, including without limitation the rights to
+use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
+the Software, and to permit persons to whom the Software is furnished to do so,
+subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
+FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
+COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
+IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
diff --git a/README.md b/README.md
index c2eafa26..37e0a743 100644
--- a/README.md
+++ b/README.md
@@ -1,6 +1,6 @@

-Simulation and Parameter Estimation in Geoscience - A python package for simulation and gradient based parameter estimation in the context of geophysical applications.
+Simulation and Parameter Estimation in Geophysics - A python package for simulation and gradient based parameter estimation in the context of geophysical applications.
The vision is to create a package for finite volume simulation with applications to geophysical imaging and subsurface flow. To enable the understanding of the many different components, this package has the following features:
@@ -10,4 +10,20 @@ The vision is to create a package for finite volume simulation with applications
* supports 1D, 2D and 3D problems
* designed for large-scale inversions
+Documentation:
+[http://simpeg.readthedocs.org/en/latest/](http://simpeg.readthedocs.org/en/latest/)
+
+Code:
+[https://github.com/simpeg/simpeg](https://github.com/simpeg/simpeg)
+
+Tests:
+[https://travis-ci.org/simpeg/simpeg](https://travis-ci.org/simpeg/simpeg)
+
+Build Status:
[](https://travis-ci.org/simpeg/simpeg)
+
+Bugs & Issues:
+[https://github.com/simpeg/simpeg/issues](https://github.com/simpeg/simpeg/issues)
+
+Code Snippets & Tutorials:
+[http://www.row1.ca/simpeg](http://www.row1.ca/simpeg)
diff --git a/SimPEG/Data.py b/SimPEG/Data.py
new file mode 100644
index 00000000..fa37b09a
--- /dev/null
+++ b/SimPEG/Data.py
@@ -0,0 +1,147 @@
+import Utils, numpy as np
+
+
+class BaseData(object):
+ """Data holds the observed data, and the standard deviations."""
+
+ __metaclass__ = Utils.Save.Savable
+
+ std = None #: Estimated Standard Deviations
+ dobs = None #: Observed data
+ dtrue = None #: True data, if data is synthetic
+ mtrue = None #: True model, if data is synthetic
+
+ counter = None #: A SimPEG.Utils.Counter object
+
+ def __init__(self, **kwargs):
+ Utils.setKwargs(self, **kwargs)
+
+ @property
+ def prob(self):
+ """
+ The geophysical problem that explains this data, use::
+
+ data.pair(prob)
+ """
+ return getattr(self, '_prob', None)
+
+ def pair(self, p):
+ """Bind a problem to this data instance using pointers"""
+ assert hasattr(p, 'dataPair'), "Problem must have an attribute 'dataPair'."
+ assert isinstance(self, p.dataPair), "Problem requires data object must be an instance of a %s class."%(p.dataPair.__name__)
+ if p.ispaired:
+ raise Exception("The problem object is already paired to a data. Use prob.unpair()")
+ self._prob = p
+ p._data = self
+
+ def unpair(self):
+ """Unbind a problem from this data instance"""
+ if not self.ispaired: return
+ self.prob._data = None
+ self._prob = None
+
+ @property
+ def ispaired(self): return self.prob is not None
+
+ @Utils.count
+ @Utils.requires('prob')
+ def dpred(self, m, u=None):
+ """
+ Create the projected data from a model.
+ The field, u, (if provided) will be used for the predicted data
+ instead of recalculating the fields (which may be expensive!).
+
+ .. math::
+ d_\\text{pred} = P(u(m))
+
+ Where P is a projection of the fields onto the data space.
+ """
+ if u is None: u = self.prob.field(m)
+ return Utils.mkvc(self.projectField(u))
+
+
+ @Utils.count
+ def projectField(self, u):
+ """
+ This function projects the fields onto the data space.
+
+
+ .. math::
+ d_\\text{pred} = P(u(m))
+ """
+ return u
+
+ #TODO: def projectFieldDeriv(self, u): Does this need to be made??!
+
+ @Utils.count
+ def residual(self, m, u=None):
+ """
+ :param numpy.array m: geophysical model
+ :param numpy.array u: fields
+ :rtype: float
+ :return: data residual
+
+ The data residual:
+
+ .. math::
+
+ \mu_\\text{data} = \mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}
+
+ """
+ return Utils.mkvc(self.dpred(m, u=u) - self.dobs)
+
+
+ @property
+ def Wd(self):
+ """
+ Data weighting matrix. This is a covariance matrix used in::
+
+ def data.residualWeighted(m,u=None):
+ return self.Wd*self.residual(m, u=u)
+
+ By default, this is based on the norm of the data plus a noise floor.
+
+ """
+ if getattr(self,'_Wd',None) is None:
+ eps = np.linalg.norm(Utils.mkvc(self.dobs),2)*1e-5
+ self._Wd = 1/(abs(self.dobs)*self.std+eps)
+ return self._Wd
+ @Wd.setter
+ def Wd(self, value):
+ self._Wd = value
+
+ def residualWeighted(self, m, u=None):
+ """
+ :param numpy.array m: geophysical model
+ :param numpy.array u: fields
+ :rtype: float
+ :return: data residual
+
+ The weighted data residual:
+
+ .. math::
+
+ \mu_\\text{data}^{\\text{weighted}} = \mathbf{W}_d(\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs})
+
+ Where W_d is a covariance matrix that weights the data residual.
+ """
+ return Utils.mkvc(self.Wd*self.residual(m, u=u))
+
+ @property
+ def RHS(self):
+ """
+ Source matrix.
+ """
+ return getattr(self, '_RHS', None)
+ @RHS.setter
+ def RHS(self, value):
+ self._RHS = value
+
+ @property
+ def isSynthetic(self):
+ "Check if the data is synthetic."
+ return (self.mtrue is not None)
+
+if __name__ == '__main__':
+ d = BaseData()
+ d.dpred()
diff --git a/SimPEG/examples/DC.py b/SimPEG/Examples/DC.py
similarity index 68%
rename from SimPEG/examples/DC.py
rename to SimPEG/Examples/DC.py
index ad98d6c2..074f5e3d 100644
--- a/SimPEG/examples/DC.py
+++ b/SimPEG/Examples/DC.py
@@ -1,27 +1,53 @@
-from SimPEG.mesh import TensorMesh
-from SimPEG.forward import Problem, 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
+from SimPEG import *
-class DCProblem(ModelTransforms.LogModel, Problem):
+
+class DCData(Data.BaseData):
+ """
+ **DCData**
+
+ Geophysical DC resistivity data.
+
+ """
+
+ P = None #: projection
+
+ def __init__(self, **kwargs):
+ Data.BaseData.__init__(self, **kwargs)
+ Utils.setKwargs(self, **kwargs)
+
+ def reshapeFields(self, u):
+ if len(u.shape) == 1:
+ u = u.reshape([-1, self.RHS.shape[1]], order='F')
+ return u
+
+ def projectField(self, u):
+ """
+ Predicted data.
+
+ .. math::
+ d_\\text{pred} = Pu(m)
+ """
+ u = self.reshapeFields(u)
+ return Utils.mkvc(self.P*u)
+
+
+
+class DCProblem(Problem.BaseProblem):
"""
**DCProblem**
Geophysical DC resistivity problem.
"""
- def __init__(self, mesh):
- Problem.__init__(self, 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
+ dataPair = DCData
+
+ def __init__(self, mesh, model, **kwargs):
+ Problem.BaseProblem.__init__(self, mesh, model)
+ self.mesh.setCellGradBC('neumann')
+ Utils.setKwargs(self, **kwargs)
+
def createMatrix(self, m):
"""
@@ -38,30 +64,16 @@ class DCProblem(ModelTransforms.LogModel, Problem):
"""
D = self.mesh.faceDiv
G = self.mesh.cellGrad
- sigma = self.modelTransform(m)
+ sigma = self.model.transform(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)
+ phi = solve.solve(self.data.RHS)
+ return Utils.mkvc(phi)
def J(self, m, v, u=None):
"""
@@ -88,24 +100,24 @@ class DCProblem(ModelTransforms.LogModel, Problem):
if u is None:
u = self.field(m)
- u = self.reshapeFields(u)
+ u = self.data.reshapeFields(u)
- P = self.P
+ P = self.data.P
D = self.mesh.faceDiv
G = self.mesh.cellGrad
A = self.createMatrix(m)
Av_dm = self.mesh.getFaceMassDeriv()
- mT_dm = self.modelTransformDeriv(m)
+ mT_dm = self.model.transformDeriv(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 ) ) )
+ dCdm[:, i] = D * ( Utils.sdiag( G * ui ) * ( Av_dm * ( mT_dm * v ) ) )
solve = Solver(dCdu)
Jv = - P * solve.solve(dCdm)
- return mkvc(Jv)
+ return Utils.mkvc(Jv)
def Jt(self, m, v, u=None):
"""Takes data, turns it into a model..ish"""
@@ -113,15 +125,15 @@ class DCProblem(ModelTransforms.LogModel, Problem):
if u is None:
u = self.field(m)
- u = self.reshapeFields(u)
- v = self.reshapeFields(v)
+ u = self.data.reshapeFields(u)
+ v = self.data.reshapeFields(v)
- P = self.P
+ P = self.data.P
D = self.mesh.faceDiv
G = self.mesh.cellGrad
A = self.createMatrix(m)
Av_dm = self.mesh.getFaceMassDeriv()
- mT_dm = self.modelTransformDeriv(m)
+ mT_dm = self.model.transformDeriv(m)
dCdu = A.T
solve = Solver(dCdu)
@@ -130,7 +142,7 @@ class DCProblem(ModelTransforms.LogModel, Problem):
Jtv = 0
for i, ui in enumerate(u.T): # loop over each column
- Jtv += sdiag( G * ui ) * ( D.T * w[:,i] )
+ Jtv += Utils.sdiag( G * ui ) * ( D.T * w[:,i] )
Jtv = - mT_dm.T * ( Av_dm.T * Jtv )
return Jtv
@@ -165,16 +177,13 @@ def genTxRxmat(nelec, spacelec, surfloc, elecini, mesh):
return q, Q, rxmidLoc
-
if __name__ == '__main__':
-
- from SimPEG import inverse
import matplotlib.pyplot as plt
# Create the mesh
h1 = np.ones(20)
h2 = np.ones(100)
- mesh = TensorMesh([h1,h2])
+ M = Mesh.TensorMesh([h1,h2])
# Create some parameters for the model
sig1 = np.log(1)
@@ -184,9 +193,9 @@ if __name__ == '__main__':
p0 = [5, 10]
p1 = [15, 50]
condVals = [sig1, sig2]
- mSynth = ModelBuilder.defineBlockConductivity(p0,p1,mesh.gridCC,condVals)
- plt.colorbar(mesh.plotImage(mSynth))
- plt.show()
+ mSynth = Utils.ModelBuilder.defineBlockConductivity(M.gridCC,p0,p1,condVals)
+ plt.colorbar(M.plotImage(mSynth))
+ # plt.show()
# Set up the projection
nelec = 50
@@ -196,41 +205,40 @@ if __name__ == '__main__':
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)
+ q, Q, rxmidloc = genTxRxmat(nelec, spacelec, surfloc, elecini, M)
P = Q.T
+ model = Model.LogModel(M)
+ prob = DCProblem(M, model)
+
# Create some data
- problem = DCProblem(mesh)
- problem.P = P
- problem.RHS = q
- dobs, Wd = problem.createSyntheticData(mSynth, std=0.05)
+ data = prob.createSyntheticData(mSynth, std=0.05, P=P, RHS=q)
- u = problem.field(mSynth)
- u = problem.reshapeFields(u)
- mesh.plotImage(u[:,10])
- # plt.show()
+ u = prob.field(mSynth)
+ u = data.reshapeFields(u)
+ M.plotImage(u[:,10])
+ plt.show()
- # Now set up the problem to do some minimization
- problem.dobs = dobs
- problem.std = dobs*0 + 0.05
- m0 = mesh.gridCC[:,0]*0+sig2
+ # Now set up the prob to do some minimization
+ # prob.dobs = dobs
+ # prob.std = dobs*0 + 0.05
+ m0 = M.gridCC[:,0]*0+sig2
- opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
- reg = inverse.Regularization(mesh)
- inv = inverse.Inversion(problem, reg, opt, beta0=1e4)
+ reg = Regularization.Tikhonov(model)
+ objFunc = ObjFunction.BaseObjFunction(data, reg)
+ opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=3, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
+ inv = Inversion.BaseInversion(objFunc, opt)
# Check Derivative
- derChk = lambda m: [inv.dataObj(m), inv.dataObjDeriv(m)]
- checkDerivative(derChk, mSynth)
+ derChk = lambda m: [objFunc.dataObj(m), objFunc.dataObjDeriv(m)]
+ # Tests.checkDerivative(derChk, mSynth)
-
-
- print inv.dataObj(m0)
- print inv.dataObj(mSynth)
+ print objFunc.dataObj(m0)
+ print objFunc.dataObj(mSynth)
m = inv.run(m0)
- plt.colorbar(mesh.plotImage(m))
+ plt.colorbar(M.plotImage(m))
print m
plt.show()
diff --git a/SimPEG/examples/Linear.py b/SimPEG/Examples/Linear.py
similarity index 62%
rename from SimPEG/examples/Linear.py
rename to SimPEG/Examples/Linear.py
index dbee7dff..a83b26f1 100644
--- a/SimPEG/examples/Linear.py
+++ b/SimPEG/Examples/Linear.py
@@ -1,14 +1,15 @@
-from SimPEG import mesh, forward, inverse, np
+from SimPEG import *
import matplotlib.pyplot as plt
-class LinearProblem(forward.Problem):
+class LinearProblem(Problem.BaseProblem):
"""docstring for LinearProblem"""
- def __init__(self, *args, **kwargs):
- forward.Problem.__init__(self, *args, **kwargs)
+ def __init__(self, mesh, model, G, **kwargs):
+ Problem.BaseProblem.__init__(self, mesh, model, **kwargs)
+ self.G = G
- def dpred(self, m, u=None):
+ def field(self, m, u=None):
return self.G.dot(m)
def J(self, m, v, u=None):
@@ -20,7 +21,7 @@ class LinearProblem(forward.Problem):
def example(N):
h = np.ones(N)/N
- M = mesh.TensorMesh([h])
+ M = Mesh.TensorMesh([h])
nk = 20
jk = np.linspace(1.,20.,nk)
@@ -39,21 +40,24 @@ def example(N):
mtrue[M.vectorCCx > 0.45] = -0.5
mtrue[M.vectorCCx > 0.6] = 0
- prob = LinearProblem(M)
- prob.G = G
+
+
+ model = Model.BaseModel(M)
+ prob = LinearProblem(M, model, G)
data = prob.createSyntheticData(mtrue, std=0.01)
- return prob, data
+ return prob, data, model
if __name__ == '__main__':
- prob, data = example(100)
+ prob, data, model = example(100)
M = prob.mesh
- reg = inverse.Regularization(M)
- opt = inverse.InexactGaussNewton(maxIter=20)
- inv = inverse.Inversion(prob,reg,opt,data)
+ reg = Regularization.Tikhonov(model)
+ objFunc = ObjFunction.BaseObjFunction(data, reg)
+ opt = Optimization.InexactGaussNewton(maxIter=20)
+ inv = Inversion.BaseInversion(objFunc, opt)
m0 = np.zeros_like(data.mtrue)
mrec = inv.run(m0)
diff --git a/SimPEG/examples/__init__.py b/SimPEG/Examples/__init__.py
similarity index 100%
rename from SimPEG/examples/__init__.py
rename to SimPEG/Examples/__init__.py
diff --git a/SimPEG/Inversion.py b/SimPEG/Inversion.py
new file mode 100644
index 00000000..813d201b
--- /dev/null
+++ b/SimPEG/Inversion.py
@@ -0,0 +1,125 @@
+import SimPEG
+from SimPEG import Utils, sp, np
+from Optimization import Remember, IterationPrinters
+
+
+class BaseInversion(object):
+ """BaseInversion(objFunc, opt, **kwargs)
+ """
+
+ __metaclass__ = Utils.Save.Savable
+
+ name = 'BaseInversion'
+
+ debug = False #: Print debugging information
+
+ counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things
+
+ def __init__(self, objFunc, opt, **kwargs):
+ Utils.setKwargs(self, **kwargs)
+
+ self.objFunc = objFunc
+ self.objFunc.parent = self
+
+ self.opt = opt
+ self.opt.parent = self
+
+ self.stoppers = [StoppingCriteria.iteration]
+
+ # Check if we have inserted printers into the optimization
+ if IterationPrinters.phi_d not in self.opt.printers:
+ self.opt.printers.insert(1,IterationPrinters.beta)
+ self.opt.printers.insert(2,IterationPrinters.phi_d)
+ self.opt.printers.insert(3,IterationPrinters.phi_m)
+
+ if not hasattr(opt, '_bfgsH0') and hasattr(opt, 'bfgsH0'): # Check if it has been set by the user and the default is not being used.
+ #TODO: I don't think that this if statement is working...
+ print 'Setting bfgsH0 to the inverse of the modelObj2Deriv. Done using direct methods.'
+ opt.bfgsH0 = SimPEG.Solver(objFunc.reg.modelObj2Deriv())
+
+
+ #TODO: Move this to the data class?
+ @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.data.dobs.size #
+ return self._phi_d_target
+
+ @phi_d_target.setter
+ def phi_d_target(self, value):
+ self._phi_d_target = value
+
+ @Utils.timeIt
+ def run(self, m0):
+ """run(m0)
+
+ Runs the inversion!
+
+ """
+ self.objFunc.startup(m0)
+ self.m = self.opt.minimize(self.objFunc.evalFunction, m0)
+ self.finish()
+
+ return self.m
+
+ @Utils.callHooks('finish')
+ def finish(self):
+ """finish()
+
+ **finish** is called at the end of the optimization.
+ """
+ pass
+
+
+ def save(self, group):
+ group.attrs['phi_d'] = self.phi_d
+ group.attrs['phi_m'] = self.phi_m
+ group.setArray('m', self.m)
+ group.setArray('dpred', self.dpred)
+
+
+
+# class Inversion(Cooling, Remember, BaseInversion):
+
+# maxIter = 10
+# name = "SimPEG Inversion"
+
+# def __init__(self, prob, reg, opt, data, **kwargs):
+# BaseInversion.__init__(self, prob, reg, opt, data, **kwargs)
+
+# self.stoppers.append(StoppingCriteria.phi_d_target_Inversion)
+
+# if StoppingCriteria.phi_d_target_Minimize not in self.opt.stoppers:
+# self.opt.stoppers.append(StoppingCriteria.phi_d_target_Minimize)
+
+# class TimeSteppingInversion(Remember, BaseInversion):
+# """
+# A slightly different view on regularization parameters,
+# let Beta be viewed as 1/dt, and timestep by updating the
+# reference model every optimization iteration.
+# """
+# maxIter = 1
+# name = "Time-Stepping SimPEG Inversion"
+
+# def __init__(self, prob, reg, opt, data, **kwargs):
+# BaseInversion.__init__(self, prob, reg, opt, data, **kwargs)
+
+# self.stoppers.append(StoppingCriteria.phi_d_target_Inversion)
+
+# if StoppingCriteria.phi_d_target_Minimize not in self.opt.stoppers:
+# self.opt.stoppers.append(StoppingCriteria.phi_d_target_Minimize)
+
+# def _startup_TimeSteppingInversion(self, m0):
+
+# def _doEndIteration_updateMref(self, xt):
+# if self.debug: 'Updating the reference model.'
+# self.parent.reg.mref = self.xc
+
+# self.opt.hook(_doEndIteration_updateMref, overwrite=True)
diff --git a/SimPEG/mesh/BaseMesh.py b/SimPEG/Mesh/BaseMesh.py
similarity index 91%
rename from SimPEG/mesh/BaseMesh.py
rename to SimPEG/Mesh/BaseMesh.py
index 6151d20b..672b5c02 100644
--- a/SimPEG/mesh/BaseMesh.py
+++ b/SimPEG/Mesh/BaseMesh.py
@@ -1,5 +1,5 @@
import numpy as np
-from SimPEG import utils
+from SimPEG import Utils
class BaseMesh(object):
@@ -78,7 +78,7 @@ class BaseMesh(object):
x_array = np.ones((x.size, len(x)))
# Unwrap it and put it in a np array
for i, xi in enumerate(x):
- x_array[:, i] = utils.mkvc(xi)
+ x_array[:, i] = Utils.mkvc(xi)
x = x_array
assert type(x) == np.ndarray, "x must be a numpy array"
@@ -91,17 +91,17 @@ class BaseMesh(object):
if format == 'M':
return xx.reshape(nn, order='F')
elif format == 'V':
- return utils.mkvc(xx)
+ return Utils.mkvc(xx)
def switchKernal(xx):
"""Switches over the different options."""
if xType in ['CC', 'N']:
- nn = (self.n) if xType == 'CC' else (self.n+1)
+ nn = (self._n) if xType == 'CC' else (self._n+1)
assert xx.size == np.prod(nn), "Number of elements must not change."
return outKernal(xx, nn)
elif xType in ['F', 'E']:
# This will only deal with components of fields, not full 'F' or 'E'
- xx = utils.mkvc(xx) # unwrap it in case it is a matrix
+ xx = Utils.mkvc(xx) # unwrap it in case it is a matrix
nn = self.nFv if xType == 'F' else self.nEv
nn = np.r_[0, nn]
@@ -147,16 +147,6 @@ class BaseMesh(object):
else:
return switchKernal(x)
- def n():
- doc = """
- Number of Cells in each dimension (array of integers)
-
- :rtype: numpy.array
- :return: n
- """
- fget = lambda self: self._n
- return locals()
- n = property(**n())
def dim():
doc = """
@@ -176,7 +166,7 @@ class BaseMesh(object):
:rtype: int
:return: nCx
"""
- fget = lambda self: self.n[0]
+ fget = lambda self: self._n[0]
return locals()
nCx = property(**nCx())
@@ -190,7 +180,7 @@ class BaseMesh(object):
def fget(self):
if self.dim > 1:
- return self.n[1]
+ return self._n[1]
else:
return None
return locals()
@@ -205,7 +195,7 @@ class BaseMesh(object):
def fget(self):
if self.dim > 2:
- return self.n[2]
+ return self._n[2]
else:
return None
return locals()
@@ -219,12 +209,12 @@ class BaseMesh(object):
:return: nC
.. plot::
+ :include-source:
- from SimPEG.mesh import TensorMesh
- import numpy as np
- TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(centers=True,showIt=True)
+ from SimPEG import Mesh, np
+ Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(centers=True,showIt=True)
"""
- fget = lambda self: np.prod(self.n)
+ fget = lambda self: np.prod(self._n)
return locals()
nC = property(**nC())
@@ -260,7 +250,7 @@ class BaseMesh(object):
def fget(self):
if self.dim > 1:
- return self.n[1] + 1
+ return self._n[1] + 1
else:
return None
return locals()
@@ -276,7 +266,7 @@ class BaseMesh(object):
def fget(self):
if self.dim > 2:
- return self.n[2] + 1
+ return self._n[2] + 1
else:
return None
return locals()
@@ -290,12 +280,12 @@ class BaseMesh(object):
:return: nN
.. plot::
+ :include-source:
- from SimPEG.mesh import TensorMesh
- import numpy as np
- TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(nodes=True,showIt=True)
+ from SimPEG import Mesh, np
+ Mesh.TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(nodes=True,showIt=True)
"""
- fget = lambda self: np.prod(self.n + 1)
+ fget = lambda self: np.prod(self.nCv + 1)
return locals()
nN = property(**nN())
@@ -361,10 +351,10 @@ class BaseMesh(object):
:return: [prod(nEx), prod(nEy), prod(nEz)]
.. plot::
+ :include-source:
- from SimPEG.mesh import TensorMesh
- import numpy as np
- TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(edges=True,showIt=True)
+ from SimPEG import Mesh, np
+ Mesh.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()
@@ -433,10 +423,10 @@ class BaseMesh(object):
:return: [prod(nFx), prod(nFy), prod(nFz)]
.. plot::
+ :include-source:
- from SimPEG.mesh import TensorMesh
- import numpy as np
- TensorMesh([np.ones(n) for n in [2,3]]).plotGrid(faces=True,showIt=True)
+ from SimPEG import Mesh, np
+ Mesh.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
similarity index 99%
rename from SimPEG/mesh/Cyl1DMesh.py
rename to SimPEG/Mesh/Cyl1DMesh.py
index e22e12b9..2b291c54 100644
--- a/SimPEG/mesh/Cyl1DMesh.py
+++ b/SimPEG/Mesh/Cyl1DMesh.py
@@ -1,7 +1,7 @@
import numpy as np
import scipy.sparse as sp
from scipy.constants import pi
-from SimPEG.utils import mkvc, ndgrid, sdiag
+from SimPEG.Utils import mkvc, ndgrid, sdiag
class Cyl1DMesh(object):
"""
@@ -84,7 +84,7 @@ class Cyl1DMesh(object):
doc = "Total number of cells in each direction"
fget = lambda self: np.array([self.nCx, self.nCz])
return locals()
- nCv = property(**nCv())
+ nCv = property(**nCv())
def nNr():
doc = "Number of nodes in the radial direction"
diff --git a/SimPEG/mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py
similarity index 97%
rename from SimPEG/mesh/DiffOperators.py
rename to SimPEG/Mesh/DiffOperators.py
index d384ca17..978f964b 100644
--- a/SimPEG/mesh/DiffOperators.py
+++ b/SimPEG/Mesh/DiffOperators.py
@@ -1,6 +1,6 @@
import numpy as np
from scipy import sparse as sp
-from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros, ddx, av, avExtrap
+from SimPEG.Utils import mkvc, sdiag, speye, kron3, spzeros, ddx, av, avExtrap
def checkBC(bc):
@@ -129,7 +129,7 @@ class DiffOperators(object):
def fget(self):
if(self._faceDiv is None):
# The number of cell centers in each direction
- n = self.n
+ n = self.nCv
# Compute faceDivergence operator on faces
if(self.dim == 1):
D = ddx(n[0])
@@ -158,7 +158,7 @@ class DiffOperators(object):
def fget(self):
if(self._faceDivx is None):
# The number of cell centers in each direction
- n = self.n
+ n = self.nCv
# Compute faceDivergence operator on faces
if(self.dim == 1):
D1 = ddx(n[0])
@@ -183,7 +183,7 @@ class DiffOperators(object):
if(self.dim < 2): return None
if(self._faceDivy is None):
# The number of cell centers in each direction
- n = self.n
+ n = self.nCv
# Compute faceDivergence operator on faces
if(self.dim == 2):
D2 = sp.kron(ddx(n[1]), speye(n[0]))
@@ -225,7 +225,7 @@ class DiffOperators(object):
def fget(self):
if(self._nodalGrad is None):
# The number of cell centers in each direction
- n = self.n
+ n = self.nCv
# Compute divergence operator on faces
if(self.dim == 1):
G = ddx(n[0])
@@ -253,7 +253,7 @@ class DiffOperators(object):
if(self._nodalLaplacian is None):
print 'Warning: Laplacian has not been tested rigorously.'
# The number of cell centers in each direction
- n = self.n
+ n = self.nCv
# Compute divergence operator on faces
if(self.dim == 1):
D1 = sdiag(1./self.hx) * ddx(mesh.nCx)
@@ -291,7 +291,7 @@ class DiffOperators(object):
"""
if(type(BC) is str):
- BC = [BC for _ in self.n] # Repeat the str self.dim times
+ BC = [BC for _ in self.nCv] # Repeat the str self.dim times
elif(type(BC) is list):
assert len(BC) == self.dim, 'BC list must be the size of your mesh'
else:
@@ -313,7 +313,7 @@ class DiffOperators(object):
def fget(self):
if(self._cellGrad is None):
BC = self.setCellGradBC(self._cellGradBC_list)
- n = self.n
+ n = self.nCv
if(self.dim == 1):
G = ddxCellGrad(n[0], BC[0])
elif(self.dim == 2):
@@ -340,7 +340,7 @@ class DiffOperators(object):
def fget(self):
if(self._cellGradBC is None):
BC = self.setCellGradBC(self._cellGradBC_list)
- n = self.n
+ n = self.nCv
if(self.dim == 1):
G = ddxCellGradBC(n[0], BC[0])
elif(self.dim == 2):
@@ -367,7 +367,7 @@ class DiffOperators(object):
def fget(self):
if getattr(self, '_cellGradx', None) is None:
BC = ['neumann', 'neumann']
- n = self.n
+ n = self.nCv
if(self.dim == 1):
G1 = ddxCellGrad(n[0], BC)
elif(self.dim == 2):
@@ -388,7 +388,7 @@ class DiffOperators(object):
if self.dim < 2: return None
if getattr(self, '_cellGrady', None) is None:
BC = ['neumann', 'neumann']
- n = self.n
+ n = self.nCv
if(self.dim == 2):
G2 = sp.kron(ddxCellGrad(n[1], BC), speye(n[0]))
elif(self.dim == 3):
@@ -466,7 +466,7 @@ class DiffOperators(object):
def fget(self):
if(self._aveF2CC is None):
- n = self.n
+ n = self.nCv
if(self.dim == 1):
self._aveF2CC = av(n[0])
elif(self.dim == 2):
@@ -486,7 +486,7 @@ class DiffOperators(object):
def fget(self):
if(self._aveCC2F is None):
- n = self.n
+ n = self.nCv
if(self.dim == 1):
self._aveCC2F = avExtrap(n[0])
elif(self.dim == 2):
@@ -507,7 +507,7 @@ class DiffOperators(object):
def fget(self):
if(self._aveE2CC is None):
# The number of cell centers in each direction
- n = self.n
+ n = self.nCv
if(self.dim == 1):
raise Exception('Edge Averaging does not make sense in 1D: Use Identity?')
elif(self.dim == 2):
@@ -528,7 +528,7 @@ class DiffOperators(object):
def fget(self):
if(self._aveN2CC is None):
# The number of cell centers in each direction
- n = self.n
+ n = self.nCv
if(self.dim == 1):
self._aveN2CC = av(n[0])
elif(self.dim == 2):
@@ -546,7 +546,7 @@ class DiffOperators(object):
def fget(self):
if(self._aveN2E is None):
# The number of cell centers in each direction
- n = self.n
+ n = self.nCv
if(self.dim == 1):
self._aveN2E = av(n[0])
elif(self.dim == 2):
@@ -567,7 +567,7 @@ class DiffOperators(object):
def fget(self):
if(self._aveN2F is None):
# The number of cell centers in each direction
- n = self.n
+ n = self.nCv
if(self.dim == 1):
self._aveN2F = av(n[0])
elif(self.dim == 2):
diff --git a/SimPEG/mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py
similarity index 99%
rename from SimPEG/mesh/InnerProducts.py
rename to SimPEG/Mesh/InnerProducts.py
index 45853071..9e7c1d93 100644
--- a/SimPEG/mesh/InnerProducts.py
+++ b/SimPEG/Mesh/InnerProducts.py
@@ -1,5 +1,5 @@
from scipy import sparse as sp
-from SimPEG.utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
+from SimPEG.Utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
import numpy as np
diff --git a/SimPEG/mesh/LogicallyOrthogonalMesh.py b/SimPEG/Mesh/LogicallyOrthogonalMesh.py
similarity index 74%
rename from SimPEG/mesh/LogicallyOrthogonalMesh.py
rename to SimPEG/Mesh/LogicallyOrthogonalMesh.py
index b3dcd095..e39c0f54 100644
--- a/SimPEG/mesh/LogicallyOrthogonalMesh.py
+++ b/SimPEG/Mesh/LogicallyOrthogonalMesh.py
@@ -1,4 +1,4 @@
-from SimPEG import utils, np
+from SimPEG import Utils, np
from BaseMesh import BaseMesh
from DiffOperators import DiffOperators
from InnerProducts import InnerProducts
@@ -7,8 +7,8 @@ from LomView import LomView
# Some helper functions.
length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5
length3D = lambda x: (x[:, 0]**2 + x[:, 1]**2 + x[:, 2]**2)**0.5
-normalize2D = lambda x: x/np.kron(np.ones((1, 2)), utils.mkvc(length2D(x), 2))
-normalize3D = lambda x: x/np.kron(np.ones((1, 3)), utils.mkvc(length3D(x), 2))
+normalize2D = lambda x: x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2))
+normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2))
class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
@@ -17,11 +17,16 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
Example of a logically orthogonal mesh:
- .. plot:: examples/mesh/plot_LogicallyOrthogonalMesh.py
+ .. plot::
+ :include-source:
+ from SimPEG import Mesh, Utils
+ X, Y = Utils.exampleLomGird([3,3],'rotate')
+ M = Mesh.LogicallyOrthogonalMesh([X, Y])
+ M.plotGrid(showIt=True)
"""
- __metaclass__ = utils.Save.Savable
+ __metaclass__ = Utils.Save.Savable
_meshType = 'LOM'
@@ -40,7 +45,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] = utils.mkvc(node_i.astype(float))
+ self._gridN[:, i] = Utils.mkvc(node_i.astype(float))
def gridCC():
doc = "Cell-centered grid."
@@ -71,10 +76,10 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
if self._gridFx is None:
N = self.r(self.gridN, 'N', 'N', 'M')
if self.dim == 2:
- XY = [utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
+ XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
self._gridFx = np.c_[XY[0], XY[1]]
elif self.dim == 3:
- XYZ = [utils.mkvc(0.25 * (n[:, :-1, :-1] + n[:, :-1, 1:] + n[:, 1:, :-1] + n[:, 1:, 1:])) for n in N]
+ XYZ = [Utils.mkvc(0.25 * (n[:, :-1, :-1] + n[:, :-1, 1:] + n[:, 1:, :-1] + n[:, 1:, 1:])) for n in N]
self._gridFx = np.c_[XYZ[0], XYZ[1], XYZ[2]]
return self._gridFx
return locals()
@@ -88,10 +93,10 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
if self._gridFy is None:
N = self.r(self.gridN, 'N', 'N', 'M')
if self.dim == 2:
- XY = [utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N]
+ XY = [Utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N]
self._gridFy = np.c_[XY[0], XY[1]]
elif self.dim == 3:
- XYZ = [utils.mkvc(0.25 * (n[:-1, :, :-1] + n[:-1, :, 1:] + n[1:, :, :-1] + n[1:, :, 1:])) for n in N]
+ XYZ = [Utils.mkvc(0.25 * (n[:-1, :, :-1] + n[:-1, :, 1:] + n[1:, :, :-1] + n[1:, :, 1:])) for n in N]
self._gridFy = np.c_[XYZ[0], XYZ[1], XYZ[2]]
return self._gridFy
return locals()
@@ -104,7 +109,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
def fget(self):
if self._gridFz is None and self.dim == 3:
N = self.r(self.gridN, 'N', 'N', 'M')
- XYZ = [utils.mkvc(0.25 * (n[:-1, :-1, :] + n[:-1, 1:, :] + n[1:, :-1, :] + n[1:, 1:, :])) for n in N]
+ XYZ = [Utils.mkvc(0.25 * (n[:-1, :-1, :] + n[:-1, 1:, :] + n[1:, :-1, :] + n[1:, 1:, :])) for n in N]
self._gridFz = np.c_[XYZ[0], XYZ[1], XYZ[2]]
return self._gridFz
return locals()
@@ -118,10 +123,10 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
if self._gridEx is None:
N = self.r(self.gridN, 'N', 'N', 'M')
if self.dim == 2:
- XY = [utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N]
+ XY = [Utils.mkvc(0.5 * (n[:-1, :] + n[1:, :])) for n in N]
self._gridEx = np.c_[XY[0], XY[1]]
elif self.dim == 3:
- XYZ = [utils.mkvc(0.5 * (n[:-1, :, :] + n[1:, :, :])) for n in N]
+ XYZ = [Utils.mkvc(0.5 * (n[:-1, :, :] + n[1:, :, :])) for n in N]
self._gridEx = np.c_[XYZ[0], XYZ[1], XYZ[2]]
return self._gridEx
return locals()
@@ -135,10 +140,10 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
if self._gridEy is None:
N = self.r(self.gridN, 'N', 'N', 'M')
if self.dim == 2:
- XY = [utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
+ XY = [Utils.mkvc(0.5 * (n[:, :-1] + n[:, 1:])) for n in N]
self._gridEy = np.c_[XY[0], XY[1]]
elif self.dim == 3:
- XYZ = [utils.mkvc(0.5 * (n[:, :-1, :] + n[:, 1:, :])) for n in N]
+ XYZ = [Utils.mkvc(0.5 * (n[:, :-1, :] + n[:, 1:, :])) for n in N]
self._gridEy = np.c_[XYZ[0], XYZ[1], XYZ[2]]
return self._gridEy
return locals()
@@ -151,7 +156,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
def fget(self):
if self._gridEz is None and self.dim == 3:
N = self.r(self.gridN, 'N', 'N', 'M')
- XYZ = [utils.mkvc(0.5 * (n[:, :, :-1] + n[:, :, 1:])) for n in N]
+ XYZ = [Utils.mkvc(0.5 * (n[:, :, :-1] + n[:, :, 1:])) for n in N]
self._gridEz = np.c_[XYZ[0], XYZ[1], XYZ[2]]
return self._gridEz
return locals()
@@ -194,25 +199,25 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
def fget(self):
if(self._vol is None):
if self.dim == 2:
- A, B, C, D = utils.indexCube('ABCD', self.n+1)
- normal, area = utils.faceInfo(np.c_[self.gridN, np.zeros((self.nN, 1))], A, B, C, D)
+ A, B, C, D = Utils.indexCube('ABCD', self.nCv+1)
+ normal, area = Utils.faceInfo(np.c_[self.gridN, np.zeros((self.nN, 1))], A, B, C, D)
self._vol = area
elif self.dim == 3:
# Each polyhedron can be decomposed into 5 tetrahedrons
# However, this presents a choice so we may as well divide in two ways and average.
- A, B, C, D, E, F, G, H = utils.indexCube('ABCDEFGH', self.n+1)
+ A, B, C, D, E, F, G, H = Utils.indexCube('ABCDEFGH', self.nCv+1)
- vol1 = (utils.volTetra(self.gridN, A, B, D, E) + # cutted edge top
- utils.volTetra(self.gridN, B, E, F, G) + # cutted edge top
- utils.volTetra(self.gridN, B, D, E, G) + # middle
- utils.volTetra(self.gridN, B, C, D, G) + # cutted edge bottom
- utils.volTetra(self.gridN, D, E, G, H)) # cutted edge bottom
+ vol1 = (Utils.volTetra(self.gridN, A, B, D, E) + # cutted edge top
+ Utils.volTetra(self.gridN, B, E, F, G) + # cutted edge top
+ Utils.volTetra(self.gridN, B, D, E, G) + # middle
+ Utils.volTetra(self.gridN, B, C, D, G) + # cutted edge bottom
+ Utils.volTetra(self.gridN, D, E, G, H)) # cutted edge bottom
- vol2 = (utils.volTetra(self.gridN, A, F, B, C) + # cutted edge top
- utils.volTetra(self.gridN, A, E, F, H) + # cutted edge top
- utils.volTetra(self.gridN, A, H, F, C) + # middle
- utils.volTetra(self.gridN, C, H, D, A) + # cutted edge bottom
- utils.volTetra(self.gridN, C, G, H, F)) # cutted edge bottom
+ vol2 = (Utils.volTetra(self.gridN, A, F, B, C) + # cutted edge top
+ Utils.volTetra(self.gridN, A, E, F, H) + # cutted edge top
+ Utils.volTetra(self.gridN, A, H, F, C) + # middle
+ Utils.volTetra(self.gridN, C, H, D, A) + # cutted edge bottom
+ Utils.volTetra(self.gridN, C, G, H, F)) # cutted edge bottom
self._vol = (vol1 + vol2)/2
return self._vol
@@ -228,30 +233,30 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
# Compute areas of cell faces
if(self.dim == 2):
xy = self.gridN
- A, B = utils.indexCube('AB', self.n+1, np.array([self.nNx, self.nCy]))
+ A, B = Utils.indexCube('AB', self.nCv+1, np.array([self.nNx, self.nCy]))
edge1 = xy[B, :] - xy[A, :]
normal1 = np.c_[edge1[:, 1], -edge1[:, 0]]
area1 = length2D(edge1)
- A, D = utils.indexCube('AD', self.n+1, np.array([self.nCx, self.nNy]))
+ A, D = Utils.indexCube('AD', self.nCv+1, np.array([self.nCx, self.nNy]))
# Note that we are doing A-D to make sure the normal points the right way.
# Think about it. Look at the picture. Normal points towards C iff you do this.
edge2 = xy[A, :] - xy[D, :]
normal2 = np.c_[edge2[:, 1], -edge2[:, 0]]
area2 = length2D(edge2)
- self._area = np.r_[utils.mkvc(area1), utils.mkvc(area2)]
+ self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2)]
self._normals = [normalize2D(normal1), normalize2D(normal2)]
elif(self.dim == 3):
- A, E, F, B = utils.indexCube('AEFB', self.n+1, np.array([self.nNx, self.nCy, self.nCz]))
- normal1, area1 = utils.faceInfo(self.gridN, A, E, F, B, average=False, normalizeNormals=False)
+ A, E, F, B = Utils.indexCube('AEFB', self.nCv+1, np.array([self.nNx, self.nCy, self.nCz]))
+ normal1, area1 = Utils.faceInfo(self.gridN, A, E, F, B, average=False, normalizeNormals=False)
- A, D, H, E = utils.indexCube('ADHE', self.n+1, np.array([self.nCx, self.nNy, self.nCz]))
- normal2, area2 = utils.faceInfo(self.gridN, A, D, H, E, average=False, normalizeNormals=False)
+ A, D, H, E = Utils.indexCube('ADHE', self.nCv+1, np.array([self.nCx, self.nNy, self.nCz]))
+ normal2, area2 = Utils.faceInfo(self.gridN, A, D, H, E, average=False, normalizeNormals=False)
- A, B, C, D = utils.indexCube('ABCD', self.n+1, np.array([self.nCx, self.nCy, self.nNz]))
- normal3, area3 = utils.faceInfo(self.gridN, A, B, C, D, average=False, normalizeNormals=False)
+ A, B, C, D = Utils.indexCube('ABCD', self.nCv+1, np.array([self.nCx, self.nCy, self.nNz]))
+ normal3, area3 = Utils.faceInfo(self.gridN, A, B, C, D, average=False, normalizeNormals=False)
- self._area = np.r_[utils.mkvc(area1), utils.mkvc(area2), utils.mkvc(area3)]
+ self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2), Utils.mkvc(area3)]
self._normals = [normal1, normal2, normal3]
return self._area
return locals()
@@ -291,21 +296,21 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
if(self._edge is None or self._tangents is None):
if(self.dim == 2):
xy = self.gridN
- A, D = utils.indexCube('AD', self.n+1, np.array([self.nCx, self.nNy]))
+ A, D = Utils.indexCube('AD', self.nCv+1, np.array([self.nCx, self.nNy]))
edge1 = xy[D, :] - xy[A, :]
- A, B = utils.indexCube('AB', self.n+1, np.array([self.nNx, self.nCy]))
+ A, B = Utils.indexCube('AB', self.nCv+1, np.array([self.nNx, self.nCy]))
edge2 = xy[B, :] - xy[A, :]
- self._edge = np.r_[utils.mkvc(length2D(edge1)), utils.mkvc(length2D(edge2))]
+ self._edge = np.r_[Utils.mkvc(length2D(edge1)), Utils.mkvc(length2D(edge2))]
self._tangents = np.r_[edge1, edge2]/np.c_[self._edge, self._edge]
elif(self.dim == 3):
xyz = self.gridN
- A, D = utils.indexCube('AD', self.n+1, np.array([self.nCx, self.nNy, self.nNz]))
+ A, D = Utils.indexCube('AD', self.nCv+1, np.array([self.nCx, self.nNy, self.nNz]))
edge1 = xyz[D, :] - xyz[A, :]
- A, B = utils.indexCube('AB', self.n+1, np.array([self.nNx, self.nCy, self.nNz]))
+ A, B = Utils.indexCube('AB', self.nCv+1, np.array([self.nNx, self.nCy, self.nNz]))
edge2 = xyz[B, :] - xyz[A, :]
- A, E = utils.indexCube('AE', self.n+1, np.array([self.nNx, self.nNy, self.nCz]))
+ A, E = Utils.indexCube('AE', self.nCv+1, np.array([self.nNx, self.nNy, self.nCz]))
edge3 = xyz[E, :] - xyz[A, :]
- self._edge = np.r_[utils.mkvc(length3D(edge1)), utils.mkvc(length3D(edge2)), utils.mkvc(length3D(edge3))]
+ self._edge = np.r_[Utils.mkvc(length3D(edge1)), Utils.mkvc(length3D(edge2)), Utils.mkvc(length3D(edge3))]
self._tangents = np.r_[edge1, edge2, edge3]/np.c_[self._edge, self._edge, self._edge]
return self._edge
return locals()
@@ -331,10 +336,10 @@ if __name__ == '__main__':
h3 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
dee3 = True
if dee3:
- X, Y, Z = utils.ndgrid(h1, h2, h3, vector=False)
+ X, Y, Z = Utils.ndgrid(h1, h2, h3, vector=False)
M = LogicallyOrthogonalMesh([X, Y, Z])
else:
- X, Y = utils.ndgrid(h1, h2, vector=False)
+ X, Y = Utils.ndgrid(h1, h2, vector=False)
M = LogicallyOrthogonalMesh([X, Y])
print M.r(M.normals, 'F', 'Fx', 'V')
diff --git a/SimPEG/mesh/LomView.py b/SimPEG/Mesh/LomView.py
similarity index 92%
rename from SimPEG/mesh/LomView.py
rename to SimPEG/Mesh/LomView.py
index 2a9b242b..4eceb918 100644
--- a/SimPEG/mesh/LomView.py
+++ b/SimPEG/Mesh/LomView.py
@@ -2,7 +2,7 @@ import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from mpl_toolkits.mplot3d import Axes3D
-from SimPEG.utils import mkvc
+from SimPEG.Utils import mkvc
class LomView(object):
@@ -15,10 +15,18 @@ class LomView(object):
def __init__(self):
pass
- def plotGrid(self, length=0.05):
+ def plotGrid(self, length=0.05, showIt=False):
"""Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions.
- .. plot:: examples/mesh/plot_LogicallyOrthogonalMesh.py
+
+ .. plot::
+ :include-source:
+
+ from SimPEG import Mesh, Utils
+ X, Y = Utils.exampleLomGird([3,3],'rotate')
+ M = Mesh.LogicallyOrthogonalMesh([X, Y])
+ M.plotGrid(showIt=True)
+
"""
NN = self.r(self.gridN, 'N', 'N', 'M')
if self.dim == 2:
@@ -92,4 +100,5 @@ class LomView(object):
ax.hold(False)
ax.set_xlabel('x1')
ax.set_ylabel('x2')
- fig.show()
+
+ if showIt: plt.show()
diff --git a/SimPEG/mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py
similarity index 85%
rename from SimPEG/mesh/TensorMesh.py
rename to SimPEG/Mesh/TensorMesh.py
index d4ab86b0..341cfba7 100644
--- a/SimPEG/mesh/TensorMesh.py
+++ b/SimPEG/Mesh/TensorMesh.py
@@ -1,4 +1,4 @@
-from SimPEG import utils, np, sp
+from SimPEG import Utils, np, sp
from BaseMesh import BaseMesh
from TensorView import TensorView
from DiffOperators import DiffOperators
@@ -17,23 +17,23 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
hy = np.array([1,2])
hz = np.array([1,1,1,1])
- mesh = TensorMesh([hx, hy, hz])
+ mesh = Mesh.TensorMesh([hx, hy, hz])
Example of a padded tensor mesh:
.. plot::
- from SimPEG import mesh, utils
- M = mesh.TensorMesh(utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0))))
+ from SimPEG import Mesh, Utils
+ M = Mesh.TensorMesh(Utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0))))
M.plotGrid()
For a quick tensor mesh on a (10x12x15) unit cube::
- mesh = TensorMesh([10, 12, 15])
+ mesh = Mesh.TensorMesh([10, 12, 15])
"""
- __metaclass__ = utils.Save.Savable
+ __metaclass__ = Utils.Save.Savable
_meshType = 'TENSOR'
@@ -52,7 +52,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
assert len(h) == len(self.x0), "Dimension mismatch. x0 != len(h)"
# Ensure h contains 1D vectors
- self._h = [utils.mkvc(x.astype(float)) for x in h]
+ self._h = [Utils.mkvc(x.astype(float)) for x in h]
def __str__(self):
outStr = ' ---- {0:d}-D TensorMesh ---- '.format(self.dim)
@@ -170,7 +170,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridCC is None:
- self._gridCC = utils.ndgrid(self.getTensor('CC'))
+ self._gridCC = Utils.ndgrid(self.getTensor('CC'))
return self._gridCC
return locals()
_gridCC = None # Store grid by default
@@ -181,7 +181,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridN is None:
- self._gridN = utils.ndgrid(self.getTensor('N'))
+ self._gridN = Utils.ndgrid(self.getTensor('N'))
return self._gridN
return locals()
_gridN = None # Store grid by default
@@ -192,7 +192,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridFx is None:
- self._gridFx = utils.ndgrid(self.getTensor('Fx'))
+ self._gridFx = Utils.ndgrid(self.getTensor('Fx'))
return self._gridFx
return locals()
_gridFx = None # Store grid by default
@@ -203,7 +203,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridFy is None and self.dim > 1:
- self._gridFy = utils.ndgrid(self.getTensor('Fy'))
+ self._gridFy = Utils.ndgrid(self.getTensor('Fy'))
return self._gridFy
return locals()
_gridFy = None # Store grid by default
@@ -214,7 +214,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridFz is None and self.dim > 2:
- self._gridFz = utils.ndgrid(self.getTensor('Fz'))
+ self._gridFz = Utils.ndgrid(self.getTensor('Fz'))
return self._gridFz
return locals()
_gridFz = None # Store grid by default
@@ -225,7 +225,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridEx is None:
- self._gridEx = utils.ndgrid(self.getTensor('Ex'))
+ self._gridEx = Utils.ndgrid(self.getTensor('Ex'))
return self._gridEx
return locals()
_gridEx = None # Store grid by default
@@ -236,7 +236,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridEy is None and self.dim > 1:
- self._gridEy = utils.ndgrid(self.getTensor('Ey'))
+ self._gridEy = Utils.ndgrid(self.getTensor('Ey'))
return self._gridEy
return locals()
_gridEy = None # Store grid by default
@@ -247,7 +247,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
def fget(self):
if self._gridEz is None and self.dim > 2:
- self._gridEz = utils.ndgrid(self.getTensor('Ez'))
+ self._gridEz = Utils.ndgrid(self.getTensor('Ez'))
return self._gridEz
return locals()
_gridEz = None # Store grid by default
@@ -262,13 +262,13 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
vh = self.h
# Compute cell volumes
if(self.dim == 1):
- self._vol = utils.mkvc(vh[0])
+ self._vol = Utils.mkvc(vh[0])
elif(self.dim == 2):
# Cell sizes in each direction
- self._vol = utils.mkvc(np.outer(vh[0], vh[1]))
+ self._vol = Utils.mkvc(np.outer(vh[0], vh[1]))
elif(self.dim == 3):
# Cell sizes in each direction
- self._vol = utils.mkvc(np.outer(utils.mkvc(np.outer(vh[0], vh[1])), vh[2]))
+ self._vol = Utils.mkvc(np.outer(Utils.mkvc(np.outer(vh[0], vh[1])), vh[2]))
return self._vol
return locals()
_vol = None
@@ -282,19 +282,19 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
# Ensure that we are working with column vectors
vh = self.h
# The number of cell centers in each direction
- n = self.n
+ n = self.nCv
# Compute areas of cell faces
if(self.dim == 1):
self._area = np.ones(n[0]+1)
elif(self.dim == 2):
area1 = np.outer(np.ones(n[0]+1), vh[1])
area2 = np.outer(vh[0], np.ones(n[1]+1))
- self._area = np.r_[utils.mkvc(area1), utils.mkvc(area2)]
+ self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2)]
elif(self.dim == 3):
- area1 = np.outer(np.ones(n[0]+1), utils.mkvc(np.outer(vh[1], vh[2])))
- area2 = np.outer(vh[0], utils.mkvc(np.outer(np.ones(n[1]+1), vh[2])))
- area3 = np.outer(vh[0], utils.mkvc(np.outer(vh[1], np.ones(n[2]+1))))
- self._area = np.r_[utils.mkvc(area1), utils.mkvc(area2), utils.mkvc(area3)]
+ area1 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(vh[1], vh[2])))
+ area2 = np.outer(vh[0], Utils.mkvc(np.outer(np.ones(n[1]+1), vh[2])))
+ area3 = np.outer(vh[0], Utils.mkvc(np.outer(vh[1], np.ones(n[2]+1))))
+ self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2), Utils.mkvc(area3)]
return self._area
return locals()
_area = None
@@ -308,19 +308,19 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
# Ensure that we are working with column vectors
vh = self.h
# The number of cell centers in each direction
- n = self.n
+ n = self.nCv
# Compute edge lengths
if(self.dim == 1):
- self._edge = utils.mkvc(vh[0])
+ self._edge = Utils.mkvc(vh[0])
elif(self.dim == 2):
l1 = np.outer(vh[0], np.ones(n[1]+1))
l2 = np.outer(np.ones(n[0]+1), vh[1])
- self._edge = np.r_[utils.mkvc(l1), utils.mkvc(l2)]
+ self._edge = np.r_[Utils.mkvc(l1), Utils.mkvc(l2)]
elif(self.dim == 3):
- l1 = np.outer(vh[0], utils.mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1))))
- l2 = np.outer(np.ones(n[0]+1), utils.mkvc(np.outer(vh[1], np.ones(n[2]+1))))
- l3 = np.outer(np.ones(n[0]+1), utils.mkvc(np.outer(np.ones(n[1]+1), vh[2])))
- self._edge = np.r_[utils.mkvc(l1), utils.mkvc(l2), utils.mkvc(l3)]
+ l1 = np.outer(vh[0], Utils.mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1))))
+ l2 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(vh[1], np.ones(n[2]+1))))
+ l3 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(np.ones(n[1]+1), vh[2])))
+ self._edge = np.r_[Utils.mkvc(l1), Utils.mkvc(l2), Utils.mkvc(l3)]
return self._edge
return locals()
_edge = None
@@ -410,11 +410,11 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
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.nFv if 'F' in locType else self.nEv
- components = [utils.spzeros(loc.shape[0], n) for n in nF_nE]
- components[ind] = utils.interpmat(loc, *self.getTensor(locType))
+ components = [Utils.spzeros(loc.shape[0], n) for n in nF_nE]
+ components[ind] = Utils.interpmat(loc, *self.getTensor(locType))
Q = sp.hstack(components)
elif locType in ['CC', 'N']:
- Q = utils.interpmat(loc, *self.getTensor(locType))
+ Q = Utils.interpmat(loc, *self.getTensor(locType))
else:
raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim))
return Q
diff --git a/SimPEG/mesh/TensorView.py b/SimPEG/Mesh/TensorView.py
similarity index 92%
rename from SimPEG/mesh/TensorView.py
rename to SimPEG/Mesh/TensorView.py
index a745c1cd..7446aecd 100644
--- a/SimPEG/mesh/TensorView.py
+++ b/SimPEG/Mesh/TensorView.py
@@ -2,7 +2,7 @@ import numpy as np
import matplotlib.pyplot as plt
import matplotlib
from mpl_toolkits.mplot3d import Axes3D
-from SimPEG.utils import mkvc, animate
+from SimPEG.Utils import mkvc, animate
class TensorView(object):
@@ -34,11 +34,22 @@ class TensorView(object):
:param str annotationColor: color of annotation, e.g. 'w', 'k', 'b'
:param bool showIt: call plt.show()
- .. plot:: examples/mesh/plot_image_2D.py
- :include-source:
+ .. plot::
+ :include-source:
+
+ from SimPEG import Mesh, np
+ M = Mesh.TensorMesh([20, 20])
+ I = np.sin(M.gridCC[:,0]*2*np.pi)*np.sin(M.gridCC[:,1]*2*np.pi)
+ M.plotImage(I, showIt=True)
+
+ .. plot::
+ :include-source:
+
+ from SimPEG import Mesh, np
+ M = Mesh.TensorMesh([20,20,20])
+ I = np.sin(M.gridCC[:,0]*2*np.pi)*np.sin(M.gridCC[:,1]*2*np.pi)*np.sin(M.gridCC[:,2]*2*np.pi)
+ M.plotImage(I, annotationColor='k', showIt=True)
- .. plot:: examples/mesh/plot_image_3D.py
- :include-source:
"""
assert type(I) == np.ndarray, "I must be a numpy array"
assert type(numbering) == bool, "numbering must be a bool"
@@ -124,9 +135,9 @@ class TensorView(object):
ax.axis('tight')
elif self.dim == 2:
if imageType == 'CC':
- C = I[:].reshape(self.n, order='F')
+ C = I[:].reshape(self.nCv, order='F')
elif imageType == 'N':
- C = I[:].reshape(self.n+1, order='F')
+ C = I[:].reshape(self.nNv, order='F')
C = 0.25*(C[:-1, :-1] + C[1:, :-1] + C[:-1, 1:] + C[1:, 1:])
elif imageType == 'Fx':
C = I[:].reshape(self.nFx, order='F')
@@ -153,9 +164,9 @@ class TensorView(object):
# get copy of image and average to cell-centres is necessary
if imageType == 'CC':
- Ic = I[:].reshape(self.n, order='F')
+ Ic = I[:].reshape(self.nCv, order='F')
elif imageType == 'N':
- Ic = I[:].reshape(self.n+1, order='F')
+ Ic = I[:].reshape(self.nNv, order='F')
Ic = .125*(Ic[:-1,:-1,:-1]+Ic[1:,:-1,:-1] + Ic[:-1,1:,:-1]+ Ic[1:,1:,:-1]+ Ic[:-1,:-1,1:]+Ic[1:,:-1,1:] + Ic[:-1,1:,1:]+ Ic[1:,1:,1:] )
elif imageType == 'Fx':
Ic = I[:].reshape(self.nFx, order='F')
@@ -237,11 +248,25 @@ class TensorView(object):
:param bool lines: plot lines connecting nodes
:param bool showIt: call plt.show()
- .. plot:: examples/mesh/plot_grid_2D.py
+ .. plot::
:include-source:
- .. plot:: examples/mesh/plot_grid_3D.py
+ from SimPEG import Mesh, np
+ h1 = np.linspace(.1,.5,3)
+ h2 = np.linspace(.1,.5,5)
+ mesh = Mesh.TensorMesh([h1, h2])
+ mesh.plotGrid(nodes=True, faces=True, centers=True, lines=True, showIt=True)
+
+ .. plot::
:include-source:
+
+ from SimPEG import Mesh, np
+ h1 = np.linspace(.1,.5,3)
+ h2 = np.linspace(.1,.5,5)
+ h3 = np.linspace(.1,.5,3)
+ mesh = Mesh.TensorMesh([h1,h2,h3])
+ mesh.plotGrid(nodes=True, faces=True, centers=True, lines=True, showIt=True)
+
"""
if self.dim == 1:
fig = plt.figure(1)
diff --git a/SimPEG/mesh/__init__.py b/SimPEG/Mesh/__init__.py
similarity index 100%
rename from SimPEG/mesh/__init__.py
rename to SimPEG/Mesh/__init__.py
diff --git a/SimPEG/Model.py b/SimPEG/Model.py
new file mode 100644
index 00000000..33cfeef8
--- /dev/null
+++ b/SimPEG/Model.py
@@ -0,0 +1,129 @@
+import Utils, Parameters, numpy as np, scipy.sparse as sp
+
+
+class BaseModel(object):
+ """
+ SimPEG Model
+
+ """
+
+ __metaclass__ = Utils.Save.Savable
+
+ counter = None #: A SimPEG.Utils.Counter object
+ mesh = None #: A SimPEG Mesh
+
+ def __init__(self, mesh):
+ self.mesh = mesh
+
+ def transform(self, m):
+ """
+ :param numpy.array m: model
+ :rtype: numpy.array
+ :return: transformed model
+
+ The *transform* changes the model into the physical property.
+
+ """
+ return m
+
+ def transformInverse(self, D):
+ """
+ :param numpy.array D: physical property
+ :rtype: numpy.array
+ :return: model
+
+ The *transformInverse* changes the physical property into the model.
+
+ .. note:: The *transformInverse* may not be easy to create in general.
+
+ """
+ raise NotImplementedError('The transformInverse is not implemented.')
+
+ def transformDeriv(self, m):
+ """
+ :param numpy.array m: model
+ :rtype: scipy.csr_matrix
+ :return: derivative of transformed model
+
+ The *transform* changes the model into the physical property.
+ The *transformDeriv* provides the derivative of the *transform*.
+ """
+ return sp.identity(m.size)
+
+ @property
+ def nP(self):
+ """Number of parameters in the model."""
+ return self.mesh.nC
+
+ def example(self, modelType=None):
+ return np.random.rand(self.mesh.nC)
+
+
+
+class LogModel(BaseModel):
+ """SimPEG LogModel"""
+
+ def __init__(self, mesh, **kwargs):
+ BaseModel.__init__(self, mesh, **kwargs)
+
+ def transform(self, m):
+ """
+ :param numpy.array m: model
+ :rtype: numpy.array
+ :return: transformed model
+
+ The *transform* 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(Utils.mkvc(m))
+
+
+ def transformInverse(self, D):
+ """
+ :param numpy.array D: physical property
+ :rtype: numpy.array
+ :return: model
+
+ The *transformInverse* changes the physical property into the model.
+
+ .. math::
+
+ m = \log{\sigma}
+
+ """
+ return np.log(Utils.mkvc(D))
+
+
+ def transformDeriv(self, m):
+ """
+ :param numpy.array m: model
+ :rtype: scipy.csr_matrix
+ :return: derivative of transformed model
+
+ The *transform* changes the model into the physical property.
+ The *transformDeriv* provides the derivative of the *transform*.
+
+ 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 Utils.sdiag(np.exp(Utils.mkvc(m)))
diff --git a/SimPEG/ObjFunction.py b/SimPEG/ObjFunction.py
new file mode 100644
index 00000000..ea2fe831
--- /dev/null
+++ b/SimPEG/ObjFunction.py
@@ -0,0 +1,215 @@
+import Utils, Parameters, numpy as np, scipy.sparse as sp
+
+class BaseObjFunction(object):
+ """BaseObjFunction(data, reg, **kwargs)"""
+
+ __metaclass__ = Utils.Save.Savable
+
+ beta = Parameters.ParameterProperty('beta', default=1, doc='Regularization trade-off parameter')
+
+ debug = False #: Print debugging information
+ counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things
+
+ name = 'Base Objective Function' #: Name of the objective function
+
+ u_current = None #: The most current evaluated field
+ m_current = None #: The most current model
+
+ @property
+ def parent(self):
+ """This is the parent of the objective function."""
+ return getattr(self,'_parent',None)
+ @parent.setter
+ def parent(self, p):
+ if getattr(self,'_parent',None) is not None:
+ print 'Objective function has switched to a new parent!'
+ self._parent = p
+
+ @property
+ def inv(self): return self.parent
+ @property
+ def objFunc(self): return self
+ @property
+ def opt(self): return getattr(self.parent,'opt',None)
+ @property
+ def prob(self): return self.data.prob
+ @property
+ def mesh(self): return self.data.prob.mesh
+ @property
+ def model(self): return self.data.prob.model
+
+
+ def __init__(self, data, reg, **kwargs):
+ Utils.setKwargs(self, **kwargs)
+
+ self.data = data
+
+ self.reg = reg
+ self.reg.parent = self
+
+
+ @Utils.callHooks('startup')
+ def startup(self, m0):
+ """startup(m0)
+
+ Called when inversion is first starting.
+ """
+ if self.debug: print 'Calling ObjFunction.startup'
+
+ if self.reg.mref is None:
+ print 'Regularization has not set mref. SimPEG will set it to m0.'
+ self.reg.mref = m0
+
+ self.phi_d = np.nan
+ self.phi_m = np.nan
+
+ self.m_current = m0
+
+ @Utils.timeIt
+ def evalFunction(self, m, return_g=True, return_H=True):
+ """evalFunction(m, return_g=True, return_H=True)
+ """
+
+ self.u_current = None
+ self.m_current = m
+
+ u = self.data.prob.field(m)
+ self.u_current = u
+
+ phi_d = self.dataObj(m, u=u)
+ phi_m = self.reg.modelObj(m)
+
+ self.dpred = self.data.dpred(m, u=u) # This is a cheap matrix vector calculation.
+
+ self.phi_d, self.phi_d_last = phi_d, self.phi_d
+ self.phi_m, self.phi_m_last = phi_m, self.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()*v
+
+ return phi_d2Deriv + self.beta * phi_m2Deriv
+
+ operator = sp.linalg.LinearOperator( (m.size, m.size), H_fun, dtype=m.dtype )
+ out += (operator,)
+ return out if len(out) > 1 else out[0]
+
+ @Utils.timeIt
+ def dataObj(self, m, u=None):
+ """dataObj(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.data.residualWeighted(m, u=u)
+ return 0.5*np.vdot(R, R)
+
+ @Utils.timeIt
+ def dataObjDeriv(self, m, u=None):
+ """dataObjDeriv(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.data.prob.field(m)
+
+ R = self.data.residualWeighted(m, u=u)
+
+ dmisfit = self.data.prob.Jt(m, self.data.Wd * R, u=u)
+
+ return dmisfit
+
+ @Utils.timeIt
+ def dataObj2Deriv(self, m, v, u=None):
+ """dataObj2Deriv(m, v, u=None)
+
+ :param numpy.array m: geophysical model
+ :param numpy.array v: vector to multiply
+ :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.data.prob.field(m)
+
+ R = self.data.residualWeighted(m, u=u)
+
+ # TODO: abstract to different norms a little cleaner.
+ # \/ it goes here. in l2 it is the identity.
+ dmisfit = self.data.prob.Jt_approx(m, self.data.Wd * self.data.Wd * self.data.prob.J_approx(m, v, u=u), u=u)
+
+ return dmisfit
diff --git a/SimPEG/inverse/Optimize.py b/SimPEG/Optimization.py
similarity index 90%
rename from SimPEG/inverse/Optimize.py
rename to SimPEG/Optimization.py
index 9d5edcc2..4e5963f6 100644
--- a/SimPEG/inverse/Optimize.py
+++ b/SimPEG/Optimization.py
@@ -1,5 +1,5 @@
-from SimPEG import Solver, utils, sp, np
-import matplotlib.pyplot as plt
+import Utils, numpy as np, scipy.sparse as sp
+from Solver import Solver
norm = np.linalg.norm
@@ -10,11 +10,11 @@ class StoppingCriteria(object):
"""docstring for StoppingCriteria"""
iteration = { "str": "%d : maxIter = %3d <= iter = %3d",
- "left": lambda M: M.maxIter, "right": lambda M: M._iter,
+ "left": lambda M: M.maxIter, "right": lambda M: M.iter,
"stopType": "critical"}
iterationLS = { "str": "%d : maxIterLS = %3d <= iterLS = %3d",
- "left": lambda M: M.maxIterLS, "right": lambda M: M._iterLS,
+ "left": lambda M: M.maxIterLS, "right": lambda M: M.iterLS,
"stopType": "critical"}
armijoGoldstein = { "str": "%d : ft = %1.4e <= alp*descent = %1.4e",
@@ -22,11 +22,11 @@ class StoppingCriteria(object):
"stopType": "optimal"}
tolerance_f = { "str": "%d : |fc-fOld| = %1.4e <= tolF*(1+|f0|) = %1.4e",
- "left": lambda M: 1 if M._iter==0 else abs(M.f-M.f_last), "right": lambda M: 0 if M._iter==0 else M.tolF*(1+abs(M.f0)),
+ "left": lambda M: 1 if M.iter==0 else abs(M.f-M.f_last), "right": lambda M: 0 if M.iter==0 else M.tolF*(1+abs(M.f0)),
"stopType": "optimal"}
moving_x = { "str": "%d : |xc-x_last| = %1.4e <= tolX*(1+|x0|) = %1.4e",
- "left": lambda M: 1 if M._iter==0 else norm(M.xc-M.x_last), "right": lambda M: 0 if M._iter==0 else M.tolX*(1+norm(M.x0)),
+ "left": lambda M: 1 if M.iter==0 else norm(M.xc-M.x_last), "right": lambda M: 0 if M.iter==0 else M.tolX*(1+norm(M.x0)),
"stopType": "optimal"}
tolerance_g = { "str": "%d : |proj(x-g)-x| = %1.4e <= tolG = %1.4e",
@@ -57,12 +57,12 @@ class StoppingCriteria(object):
class IterationPrinters(object):
"""docstring for IterationPrinters"""
- iteration = {"title": "#", "value": lambda M: M._iter, "width": 5, "format": "%3d"}
+ iteration = {"title": "#", "value": lambda M: M.iter, "width": 5, "format": "%3d"}
f = {"title": "f", "value": lambda M: M.f, "width": 10, "format": "%1.2e"}
norm_g = {"title": "|proj(x-g)-x|", "value": lambda M: norm(M.projection(M.xc - M.g) - M.xc), "width": 15, "format": "%1.2e"}
- totalLS = {"title": "LS", "value": lambda M: M._iterLS, "width": 5, "format": "%d"}
+ totalLS = {"title": "LS", "value": lambda M: M.iterLS, "width": 5, "format": "%d"}
- iterationLS = {"title": "#", "value": lambda M: (M._iter, M._iterLS), "width": 5, "format": "%3d.%d"}
+ iterationLS = {"title": "#", "value": lambda M: (M.iter, M.iterLS), "width": 5, "format": "%3d.%d"}
LS_ft = {"title": "ft", "value": lambda M: M._LS_ft, "width": 10, "format": "%1.2e"}
LS_t = {"title": "t", "value": lambda M: M._LS_t, "width": 10, "format": "%0.5f"}
LS_armijoGoldstein = {"title": "f + alp*g.T*p", "value": lambda M: M.f + M.LSreduction*M._LS_descent, "width": 16, "format": "%1.2e"}
@@ -72,9 +72,9 @@ class IterationPrinters(object):
bSet = {"title": "bSet", "value": lambda M: np.sum(M.bindingSet(M.xc)), "width": 8, "format": "%d"}
comment = {"title": "Comment", "value": lambda M: M.comment, "width": 12, "format": "%s"}
- beta = {"title": "beta", "value": lambda M: M.parent._beta, "width": 10, "format": "%1.2e"}
- phi_d = {"title": "phi_d", "value": lambda M: M.parent.phi_d, "width": 10, "format": "%1.2e"}
- phi_m = {"title": "phi_m", "value": lambda M: M.parent.phi_m, "width": 10, "format": "%1.2e"}
+ beta = {"title": "beta", "value": lambda M: M.parent.objFunc.beta, "width": 10, "format": "%1.2e"}
+ phi_d = {"title": "phi_d", "value": lambda M: M.parent.objFunc.phi_d, "width": 10, "format": "%1.2e"}
+ phi_m = {"title": "phi_m", "value": lambda M: M.parent.objFunc.phi_m, "width": 10, "format": "%1.2e"}
class Minimize(object):
@@ -82,7 +82,7 @@ class Minimize(object):
Minimize is a general class for derivative based optimization.
"""
- __metaclass__ = utils.Save.Savable
+ __metaclass__ = Utils.Save.Savable
name = "General Optimization Algorithm" #: The name of the optimization algorithm
@@ -100,7 +100,8 @@ class Minimize(object):
debugLS = False #: Print debugging information for the line-search
comment = '' #: Used by some functions to indicate what is going on in the algorithm
- counter = None #: Set this to a SimPEG.utils.Counter() if you want to count things
+ counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things
+ parent = None #: This is the parent of the optimization routine.
def __init__(self, **kwargs):
self.stoppers = [StoppingCriteria.tolerance_f, StoppingCriteria.moving_x, StoppingCriteria.tolerance_g, StoppingCriteria.norm_g, StoppingCriteria.iteration]
@@ -109,9 +110,9 @@ class Minimize(object):
self.printers = [IterationPrinters.iteration, IterationPrinters.f, IterationPrinters.norm_g, IterationPrinters.totalLS]
self.printersLS = [IterationPrinters.iterationLS, IterationPrinters.LS_ft, IterationPrinters.LS_t, IterationPrinters.LS_armijoGoldstein]
- utils.setKwargs(self, **kwargs)
+ Utils.setKwargs(self, **kwargs)
- @utils.timeIt
+ @Utils.timeIt
def minimize(self, evalFunction, x0):
"""minimize(evalFunction, x0)
@@ -179,17 +180,7 @@ class Minimize(object):
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
-
- @utils.callHooks('startup')
+ @Utils.callHooks('startup')
def startup(self, x0):
"""
**startup** is called at the start of any new minimize call.
@@ -198,15 +189,15 @@ class Minimize(object):
x0 = x0
xc = x0
- _iter = _iterLS = 0
+ iter = iterLS = 0
:param numpy.ndarray x0: initial x
:rtype: None
:return: None
"""
- self._iter = 0
- self._iterLS = 0
+ self.iter = 0
+ self.iterLS = 0
x0 = self.projection(x0) # ensure that we start of feasible.
self.x0 = x0
@@ -214,8 +205,8 @@ class Minimize(object):
self.f_last = np.nan
self.x_last = x0
- @utils.count
- @utils.callHooks('doStartIteration')
+ @Utils.count
+ @Utils.callHooks('doStartIteration')
def doStartIteration(self):
"""doStartIteration()
@@ -237,9 +228,9 @@ class Minimize(object):
"""
pad = ' '*10 if inLS else ''
name = self.name if not inLS else self.nameLS
- utils.printTitles(self, self.printers if not inLS else self.printersLS, name, pad)
+ Utils.printTitles(self, self.printers if not inLS else self.printersLS, name, pad)
- @utils.callHooks('printIter')
+ @Utils.callHooks('printIter')
def printIter(self, inLS=False):
"""
**printIter** is called directly after function evaluations.
@@ -249,7 +240,7 @@ class Minimize(object):
"""
pad = ' '*10 if inLS else ''
- utils.printLine(self, self.printers if not inLS else self.printersLS, pad=pad)
+ Utils.printLine(self, self.printers if not inLS else self.printersLS, pad=pad)
def printDone(self, inLS=False):
"""
@@ -262,10 +253,10 @@ class Minimize(object):
pad = ' '*10 if inLS else ''
stop, done = (' STOP! ', ' DONE! ') if not inLS else ('----------------', ' End Linesearch ')
stoppers = self.stoppers if not inLS else self.stoppersLS
- utils.printStoppers(self, stoppers, pad='', stop=stop, done=done)
+ Utils.printStoppers(self, stoppers, pad='', stop=stop, done=done)
- @utils.callHooks('finish')
+ @Utils.callHooks('finish')
def finish(self):
"""finish()
@@ -278,13 +269,13 @@ class Minimize(object):
pass
def stoppingCriteria(self, inLS=False):
- if self._iter == 0:
+ if self.iter == 0:
self.f0 = self.f
self.g0 = self.g
- return utils.checkStoppers(self, self.stoppers if not inLS else self.stoppersLS)
+ return Utils.checkStoppers(self, self.stoppers if not inLS else self.stoppersLS)
- @utils.timeIt
- @utils.callHooks('projection')
+ @Utils.timeIt
+ @Utils.callHooks('projection')
def projection(self, p):
"""projection(p)
@@ -298,7 +289,7 @@ class Minimize(object):
"""
return p
- @utils.timeIt
+ @Utils.timeIt
def findSearchDirection(self):
"""findSearchDirection()
@@ -329,7 +320,7 @@ class Minimize(object):
"""
return -self.g
- @utils.count
+ @Utils.count
def scaleSearchDirection(self, p):
"""scaleSearchDirection(p)
@@ -348,7 +339,7 @@ class Minimize(object):
nameLS = "Armijo linesearch" #: The line-search name
- @utils.timeIt
+ @Utils.timeIt
def modifySearchDirection(self, p):
"""modifySearchDirection(p)
@@ -370,23 +361,23 @@ class Minimize(object):
"""
# Projected Armijo linesearch
self._LS_t = 1
- self._iterLS = 0
- while self._iterLS < self.maxIterLS:
+ self.iterLS = 0
+ while self.iterLS < self.maxIterLS:
self._LS_xt = self.projection(self.xc + self._LS_t*p)
self._LS_ft = self.evalFunction(self._LS_xt, return_g=False, return_H=False)
self._LS_descent = np.inner(self.g, self._LS_xt - self.xc) # this takes into account multiplying by t, but is important for projection.
if self.stoppingCriteria(inLS=True): break
- self._iterLS += 1
+ self.iterLS += 1
self._LS_t = self.LSshorten*self._LS_t
if self.debugLS:
- if self._iterLS == 1: self.printInit(inLS=True)
+ if self.iterLS == 1: self.printInit(inLS=True)
self.printIter(inLS=True)
- if self.debugLS and self._iterLS > 0: self.printDone(inLS=True)
+ if self.debugLS and self.iterLS > 0: self.printDone(inLS=True)
- return self._LS_xt, self._iterLS < self.maxIterLS
+ return self._LS_xt, self.iterLS < self.maxIterLS
- @utils.count
+ @Utils.count
def modifySearchDirectionBreak(self, p):
"""modifySearchDirectionBreak(p)
@@ -408,8 +399,8 @@ class Minimize(object):
print 'The linesearch got broken. Boo.'
return p, False
- @utils.count
- @utils.callHooks('doEndIteration')
+ @Utils.count
+ @Utils.callHooks('doEndIteration')
def doEndIteration(self, xt):
"""doEndIteration(xt)
@@ -426,7 +417,7 @@ class Minimize(object):
# store old values
self.f_last = self.f
self.x_last, self.xc = self.xc, xt
- self._iter += 1
+ self.iter += 1
if self.debug: self.printDone()
@@ -527,7 +518,7 @@ class ProjectedGradient(Minimize, Remember):
self.aSet_prev = self.activeSet(x0)
- @utils.count
+ @Utils.count
def projection(self, x):
"""projection(x)
@@ -536,7 +527,7 @@ class ProjectedGradient(Minimize, Remember):
"""
return np.median(np.c_[self.lower,x,self.upper],axis=1)
- @utils.count
+ @Utils.count
def activeSet(self, x):
"""activeSet(x)
@@ -545,7 +536,7 @@ class ProjectedGradient(Minimize, Remember):
"""
return np.logical_or(x == self.lower, x == self.upper)
- @utils.count
+ @Utils.count
def inactiveSet(self, x):
"""inactiveSet(x)
@@ -554,7 +545,7 @@ class ProjectedGradient(Minimize, Remember):
"""
return np.logical_not(self.activeSet(x))
- @utils.count
+ @Utils.count
def bindingSet(self, x):
"""bindingSet(x)
@@ -567,7 +558,7 @@ class ProjectedGradient(Minimize, Remember):
bind_low = np.logical_and(x == self.upper, self.g <= 0)
return np.logical_or(bind_up, bind_low)
- @utils.timeIt
+ @Utils.timeIt
def findSearchDirection(self):
"""findSearchDirection()
@@ -612,7 +603,7 @@ class ProjectedGradient(Minimize, Remember):
# aSet_after = self.activeSet(self.xc+p)
return p
- @utils.timeIt
+ @Utils.timeIt
def _doEndIteration_ProjectedGradient(self, xt):
"""_doEndIteration_ProjectedGradient(xt)"""
aSet = self.activeSet(xt)
@@ -623,7 +614,7 @@ class ProjectedGradient(Minimize, Remember):
f_current_decrease = self.f_last - self.f
self.comment = ''
- if self._iter < 1:
+ if self.iter < 1:
# Note that this is reset on every CG iteration.
self.f_decrease_max = -np.inf
else:
@@ -694,7 +685,7 @@ class BFGS(Minimize, Remember):
return self.bfgs(-self.g)
def _doEndIteration_BFGS(self, xt):
- if self._iter is 0:
+ if self.iter is 0:
self.g_last = self.g
return
@@ -718,7 +709,7 @@ class GaussNewton(Minimize, Remember):
def __init__(self, **kwargs):
Minimize.__init__(self, **kwargs)
- @utils.timeIt
+ @Utils.timeIt
def findSearchDirection(self):
return Solver(self.H).solve(-self.g)
@@ -765,7 +756,7 @@ class InexactGaussNewton(BFGS, Minimize, Remember):
def approxHinv(self, value):
self._approxHinv = value
- @utils.timeIt
+ @Utils.timeIt
def findSearchDirection(self):
Hinv = Solver(self.H, doDirect=False, options={'iterSolver': 'CG', 'M': self.approxHinv, 'tol': self.tolCG, 'maxIter': self.maxIterCG})
p = Hinv.solve(-self.g)
@@ -778,7 +769,7 @@ class SteepestDescent(Minimize, Remember):
def __init__(self, **kwargs):
Minimize.__init__(self, **kwargs)
- @utils.timeIt
+ @Utils.timeIt
def findSearchDirection(self):
return -self.g
@@ -811,7 +802,7 @@ class NewtonRoot(object):
doLS = True
def __init__(self, **kwargs):
- utils.setKwargs(self, **kwargs)
+ Utils.setKwargs(self, **kwargs)
def root(self, fun, x):
"""root(fun, x)
@@ -827,7 +818,7 @@ class NewtonRoot(object):
"""
if self.comments: print 'Newton Method:\n'
- self._iter = 0
+ self.iter = 0
while True:
r, J = fun(x, return_g=True)
@@ -861,10 +852,10 @@ class NewtonRoot(object):
rt = fun(xt, return_g=False)
x = xt
- self._iter += 1
+ self.iter += 1
if norm(rt) < self.tol:
break
- if self._iter > self.maxIter:
+ if self.iter > self.maxIter:
print 'NewtonRoot stopped by maxIters. norm: %4.4e' % norm(rt)
break
@@ -873,7 +864,7 @@ class NewtonRoot(object):
if __name__ == '__main__':
- from SimPEG.tests import Rosenbrock, checkDerivative
+ from SimPEG.Tests import Rosenbrock, checkDerivative
import matplotlib.pyplot as plt
x0 = np.array([2.6, 3.7])
checkDerivative(Rosenbrock, x0, plotIt=False)
@@ -885,7 +876,7 @@ if __name__ == '__main__':
print 'test the newtonRoot finding.'
- fun = lambda x, return_g=True: np.sin(x) if not return_g else ( np.sin(x), utils.sdiag( np.cos(x) ) )
+ fun = lambda x, return_g=True: np.sin(x) if not return_g else ( np.sin(x), Utils.sdiag( np.cos(x) ) )
x = np.array([np.pi-0.3, np.pi+0.1, 0])
pnt = NewtonRoot(comments=True).root(fun,x)
print pnt
diff --git a/SimPEG/Parameters.py b/SimPEG/Parameters.py
new file mode 100644
index 00000000..2a07aa69
--- /dev/null
+++ b/SimPEG/Parameters.py
@@ -0,0 +1,164 @@
+import Utils, numpy as np
+
+
+class Parameter(object):
+ """Parameter"""
+
+ debug = False #: Print debugging information
+
+ current = None #: This hold
+ currentIter = 0
+
+ def __init__(self, **kwargs):
+ Utils.setKwargs(self, **kwargs)
+
+ @property
+ def parent(self):
+ """This is the parent of the Parameter instance."""
+ return getattr(self,'_parent',None)
+ @parent.setter
+ def parent(self, p):
+ startupName = '_startup_paramProperty_'+self._propertyName
+ if getattr(self,'_parent',None) is not None:
+ delattr(self._parent,startupName)
+ print 'Warning: Parameter %s has switched to a new parent.' % self._propertyName
+ if self.debug: print '%s function has been deleted' % startupName
+ self._parent = p
+
+ prop = self
+ def _startup_paramProperty(self, *args):
+ if prop.debug: print 'initializing %s' % prop._propertyName
+ prop.initialize()
+
+ Utils.hook(self._parent, _startup_paramProperty, name=startupName, overwrite=True)
+
+ @property
+ def inv(self): return self.parent.inv
+ @property
+ def objFunc(self): return self.parent.objFunc
+ @property
+ def opt(self): return self.parent.opt
+ @property
+ def reg(self): return self.parent.reg
+ @property
+ def data(self): return self.parent.data
+ @property
+ def prob(self): return self.parent.prob
+ @property
+ def model(self): return self.parent.model
+ @property
+ def mesh(self): return self.parent.mesh
+
+ def initialize(self):
+ pass
+
+ def get(self):
+ if (self.current is None or
+ not self.opt.iter == self.currentIter):
+ self.current = self.nextIter()
+ self.currentIter = self.opt.iter
+ return self.current
+
+ def nextIter(self):
+ raise NotImplementedError('Getting the Parameter is not yet implemented.')
+
+
+def ParameterProperty(name, default=None, doc=""):
+ def getter(self):
+ out = getattr(self,'_'+name,default)
+ if isinstance(out, Parameter):
+ out = out.get()
+ return out
+ def setter(self, value):
+ if isinstance(value, Parameter):
+ value._propertyName = name
+ value.parent = self
+ setattr(self, '_'+name, value)
+
+ return property(fget=getter, fset=setter, doc=doc)
+
+
+class BetaEstimate(Parameter):
+ """BetaEstimate"""
+
+ beta0 = 'guess' #: The initial Beta (regularization parameter)
+ beta0_ratio = 0.1 #: When beta0 is set to 'guess', estimateBeta0 is used with this ratio
+
+ beta = None #: Beta parameter
+
+ def __init__(self, **kwargs):
+ Parameter.__init__(self, **kwargs)
+
+ def initialize(self):
+ self.beta = self.beta0
+
+ @Utils.requires('parent')
+ def nextIter(self):
+ if self.beta is 'guess':
+ if self.debug: print 'BetaSchedule is estimating Beta0.'
+ self.beta = self.estimateBeta0()
+ return self.beta
+
+ @Utils.requires('parent')
+ def estimateBeta0(self):
+ """estimateBeta0(u=None)
+
+ The initial beta is calculated by comparing the estimated
+ eigenvalues of JtJ and WtW.
+
+ To estimate the eigenvector of **A**, we will use one iteration
+ of the *Power Method*:
+
+ .. math::
+
+ \mathbf{x_1 = A x_0}
+
+ Given this (very course) approximation of the eigenvector,
+ we can use the *Rayleigh quotient* to approximate the largest eigenvalue.
+
+ .. math::
+
+ \lambda_0 = \\frac{\mathbf{x^\\top A x}}{\mathbf{x^\\top x}}
+
+ We will approximate the largest eigenvalue for both JtJ and WtW, and
+ use some ratio of the quotient to estimate beta0.
+
+ .. math::
+
+ \\beta_0 = \gamma \\frac{\mathbf{x^\\top J^\\top J x}}{\mathbf{x^\\top W^\\top W x}}
+
+ :rtype: float
+ :return: beta0
+ """
+ objFunc = self.parent
+ data = objFunc.data
+
+ m = objFunc.m_current
+ u = objFunc.u_current
+
+ if u is None:
+ u = data.prob.field(m)
+
+ x0 = np.random.rand(*m.shape)
+ t = x0.dot(objFunc.dataObj2Deriv(m,x0,u=u))
+ b = x0.dot(objFunc.reg.modelObj2Deriv()*x0)
+ return self.beta0_ratio*(t/b)
+
+
+class BetaSchedule(BetaEstimate):
+ """BetaSchedule"""
+
+ coolingFactor = 2.
+ coolingRate = 3
+
+ @Utils.requires('parent')
+ def nextIter(self):
+ if self.beta is 'guess':
+ if self.debug: print 'BetaSchedule is estimating Beta0.'
+ self.beta = self.estimateBeta0()
+
+ if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0:
+ if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
+ self.beta /= self.coolingFactor
+
+ return self.beta
diff --git a/SimPEG/forward/Problem.py b/SimPEG/Problem.py
similarity index 58%
rename from SimPEG/forward/Problem.py
rename to SimPEG/Problem.py
index 1412ff5b..386fac9d 100644
--- a/SimPEG/forward/Problem.py
+++ b/SimPEG/Problem.py
@@ -1,8 +1,7 @@
-from SimPEG import utils, data, np, sp
-norm = np.linalg.norm
+import Utils, Data, numpy as np, scipy.sparse as sp
-class Problem(object):
+class BaseProblem(object):
"""
Problem is the base class for all geophysical forward problems in SimPEG.
@@ -35,71 +34,42 @@ class Problem(object):
to (locally) find how model parameters change the data, and optimize!
"""
- __metaclass__ = utils.Save.Savable
+ __metaclass__ = Utils.Save.Savable
- counter = None #: A SimPEG.utils.Counter object
+ counter = None #: A SimPEG.Utils.Counter object
+ dataPair = Data.BaseData
- def __init__(self, mesh, *args, **kwargs):
- utils.setKwargs(self, **kwargs)
+ def __init__(self, mesh, model, *args, **kwargs):
+ Utils.setKwargs(self, **kwargs)
self.mesh = mesh
+ self.model = model
@property
- def RHS(self):
+ def data(self):
"""
- Source matrix.
+ The data object for this problem.
"""
- return self._RHS
- @RHS.setter
- def RHS(self, value):
- self._RHS = value
+ return getattr(self, '_data', None)
+
+ def pair(self, d):
+ """Bind a data to this problem instance using pointers."""
+ assert isinstance(d, self.dataPair), "Data object must be an instance of a %s class."%(self.dataPair.__name__)
+ if d.ispaired:
+ raise Exception("The data object is already paired to a problem. Use data.unpair()")
+ self._data = d
+ d._prob = self
+
+ def unpair(self):
+ """Unbind a data from this problem instance."""
+ if not self.ispaired: return
+ self.data._prob = None
+ self._data = None
@property
- def P(self):
- """
- Projection matrix.
+ def ispaired(self): return self.data is not None
- .. math::
- d_\\text{pred} = Pu(m)
- """
- return self._P
- @P.setter
- def P(self, value):
- self._P = value
-
- @utils.count
- 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
-
- @utils.count
- def dataResidual(self, m, data, 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) - data.dobs
-
- @utils.timeIt
+ @Utils.timeIt
def J(self, m, v, u=None):
"""
:param numpy.array m: model
@@ -129,7 +99,7 @@ class Problem(object):
"""
raise NotImplementedError('J is not yet implemented.')
- @utils.timeIt
+ @Utils.timeIt
def Jt(self, m, v, u=None):
"""
:param numpy.array m: model
@@ -143,7 +113,7 @@ class Problem(object):
raise NotImplementedError('Jt is not yet implemented.')
- @utils.timeIt
+ @Utils.timeIt
def J_approx(self, m, v, u=None):
"""
@@ -158,7 +128,7 @@ class Problem(object):
"""
return self.J(m, v, u)
- @utils.timeIt
+ @Utils.timeIt
def Jt_approx(self, m, v, u=None):
"""
:param numpy.array m: model
@@ -182,33 +152,7 @@ class Problem(object):
"""
pass
- 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:
-
- """
- return 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.
- """
- return sp.identity(m.size)
-
- def createSyntheticData(self, m, std=0.05, u=None):
+ def createSyntheticData(self, m, std=0.05, u=None, **geometry_kwargs):
"""
Create synthetic data given a model, and a standard deviation.
@@ -220,11 +164,13 @@ class Problem(object):
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.
"""
- dtrue = self.dpred(m,u=u)
- noise = std*abs(dtrue)*np.random.randn(*dtrue.shape)
- dobs = dtrue+noise
- stdev = dobs*0 + std
- return data.SimPEGData(self, dobs=dobs, std=stdev, dtrue=dtrue, mtrue=m)
+ data = self.dataPair(mtrue=m, **geometry_kwargs)
+ data.pair(self)
+ data.dtrue = data.dpred(m, u=u)
+ noise = std*abs(data.dtrue)*np.random.randn(*data.dtrue.shape)
+ data.dobs = data.dtrue+noise
+ data.std = data.dobs*0 + std
+ return data
diff --git a/SimPEG/inverse/Regularization.py b/SimPEG/Regularization.py
similarity index 70%
rename from SimPEG/inverse/Regularization.py
rename to SimPEG/Regularization.py
index 8c6f19ef..aa62449a 100644
--- a/SimPEG/inverse/Regularization.py
+++ b/SimPEG/Regularization.py
@@ -1,7 +1,109 @@
-from SimPEG import utils, np, sp
+import Utils, Model, Parameters, numpy as np, scipy.sparse as sp
-class Regularization(object):
- """**Regularization**
+class BaseRegularization(object):
+ """
+ **Base Regularization Class**
+
+ This is used to regularize the model space::
+
+ reg = Regularization(mesh, model)
+
+ """
+
+ __metaclass__ = Utils.Save.Savable
+
+ modelPair = Model.BaseModel #: Some regularizations only work on specific models
+
+ model = None #: A SimPEG.Model instance.
+
+ counter = None
+
+ def __init__(self, model, **kwargs):
+ Utils.setKwargs(self, **kwargs)
+ assert isinstance(model, self.modelPair), "Incorrect model for this regularization"
+ self.model = model
+
+ mref = Parameters.ParameterProperty('mref', default=None, doc='Reference model.')
+
+ @property
+ def parent(self):
+ """This is the parent of the regularization."""
+ return getattr(self,'_parent',None)
+ @parent.setter
+ def parent(self, p):
+ if getattr(self,'_parent',None) is not None:
+ print 'Regularization has switched to a new parent!'
+ self._parent = p
+
+ @property
+ def inv(self): return self.parent.inv
+ @property
+ def objFunc(self): return self.parent
+ @property
+ def reg(self): return self
+ @property
+ def opt(self): return self.parent.opt
+ @property
+ def prob(self): return self.parent.prob
+ @property
+ def data(self): return self.parent.data
+ @property
+ def mesh(self): return self.model.mesh
+
+
+ @property
+ def W(self):
+ """Full regularization weighting matrix W."""
+ return sp.identity(self.model.nP)
+
+
+ @Utils.timeIt
+ def modelObj(self, m):
+ r = self.W * (m - self.mref)
+ return 0.5*r.dot(r)
+
+ @Utils.timeIt
+ def modelObjDeriv(self, m):
+ """
+
+ The regularization is:
+
+ .. math::
+
+ R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
+
+ So the derivative is straight forward:
+
+ .. math::
+
+ R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
+
+ """
+ return self.W.T * ( self.W * (m - self.mref) )
+
+ @Utils.timeIt
+ def modelObj2Deriv(self):
+ """
+
+ The regularization is:
+
+ .. math::
+
+ R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
+
+ So the second derivative is straight forward:
+
+ .. math::
+
+ R(m) = \mathbf{W^\\top W}
+
+ """
+ return self.W.T * self.W
+
+
+
+class Tikhonov(BaseRegularization):
+ """**Tikhonov Regularization**
Here we will define regularization of a model, m, in general however, this should be thought of as (m-m_ref) but otherwise it is exactly the same:
@@ -83,36 +185,22 @@ class Regularization(object):
"""
- __metaclass__ = utils.Save.Savable
+ alpha_s = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Ws'], "Smallness weight")
+ alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
+ alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
+ alpha_z = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
+ alpha_xx = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction")
+ alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction")
+ alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction")
- alpha_s = utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Ws'], "Smallness weight")
- alpha_x = utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
- alpha_y = utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
- alpha_z = utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
- alpha_xx = utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction")
- alpha_yy = utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction")
- alpha_zz = utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction")
-
- counter = None
-
- def __init__(self, mesh, **kwargs):
- utils.setKwargs(self, **kwargs)
- self.mesh = mesh
-
- @property
- def mref(self):
- if getattr(self, '_mref', None) is None:
- return np.zeros(self.mesh.nC);
- return self._mref
- @mref.setter
- def mref(self, value):
- self._mref = value
+ def __init__(self, model, **kwargs):
+ BaseRegularization.__init__(self, model, **kwargs)
@property
def Ws(self):
"""Regularization matrix Ws"""
if getattr(self,'_Ws', None) is None:
- self._Ws = utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
+ self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
return self._Ws
@property
@@ -120,7 +208,7 @@ class Regularization(object):
"""Regularization matrix Wx"""
if getattr(self, '_Wx', None) is None:
Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFv[0]].T*self.mesh.vol
- self._Wx = utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx
+ self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx
return self._Wx
@property
@@ -128,7 +216,7 @@ class Regularization(object):
"""Regularization matrix Wy"""
if getattr(self, '_Wy', None) is None:
Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFv[0]:np.sum(self.mesh.nFv[:2])].T*self.mesh.vol
- self._Wy = utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady
+ self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady
return self._Wy
@property
@@ -136,31 +224,30 @@ class Regularization(object):
"""Regularization matrix Wz"""
if getattr(self, '_Wz', None) is None:
Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.nFv[:2]):].T*self.mesh.vol
- self._Wz = utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz
+ self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz
return self._Wz
@property
def Wxx(self):
"""Regularization matrix Wxx"""
if getattr(self, '_Wxx', None) is None:
- self._Wxx = utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx
+ self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx
return self._Wxx
@property
def Wyy(self):
"""Regularization matrix Wyy"""
if getattr(self, '_Wyy', None) is None:
- self._Wyy = utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady
+ self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady
return self._Wyy
@property
def Wzz(self):
"""Regularization matrix Wzz"""
if getattr(self, '_Wzz', None) is None:
- self._Wzz = utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz
+ self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz
return self._Wzz
-
@property
def W(self):
"""Full regularization matrix W"""
@@ -173,47 +260,3 @@ class Regularization(object):
self._W = sp.vstack(wlist)
return self._W
-
- @utils.timeIt
- def modelObj(self, m):
- r = self.W * (m - self.mref)
- return 0.5*r.dot(r)
-
- @utils.timeIt
- def modelObjDeriv(self, m):
- """
-
- The regularization is:
-
- .. math::
-
- R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
-
- So the derivative is straight forward:
-
- .. math::
-
- R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
-
- """
- return self.W.T * ( self.W * (m - self.mref) )
-
- @utils.timeIt
- def modelObj2Deriv(self):
- """
-
- The regularization is:
-
- .. math::
-
- R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
-
- So the second derivative is straight forward:
-
- .. math::
-
- R(m) = \mathbf{W^\\top W}
-
- """
- return self.W.T * self.W
-
diff --git a/SimPEG/utils/Solver.py b/SimPEG/Solver.py
similarity index 98%
rename from SimPEG/utils/Solver.py
rename to SimPEG/Solver.py
index e0f1e017..3e9c2d02 100644
--- a/SimPEG/utils/Solver.py
+++ b/SimPEG/Solver.py
@@ -1,14 +1,15 @@
import numpy as np
import scipy.sparse as sp
import scipy.sparse.linalg as linalg
-from SimPEG.utils import mkvc, sdiag
+from Utils.matutils import mkvc
+from Utils.sputils import sdiag
import warnings
DEFAULTS = {'direct':'scipy', 'iter':'scipy', 'triangular':'fortran', 'diagonal':'python'}
OPTIONS = {'direct':['scipy'], 'iter':['scipy'], 'triangular':['python'], 'diagonal':['python']}
try:
- import TriSolve
+ import Utils.TriSolve as TriSolve
OPTIONS['triangular'].append('fortran')
except Exception, e:
print 'Warning: Python backend is being used for solver. Run setup.py from the command line.'
@@ -319,7 +320,7 @@ class Solver(object):
if __name__ == '__main__':
- from SimPEG.mesh import TensorMesh
+ from SimPEG.Mesh import TensorMesh
from time import time
h1 = np.ones(20)*100.
h2 = np.ones(20)*100.
diff --git a/SimPEG/tests/TestUtils.py b/SimPEG/Tests/TestUtils.py
similarity index 95%
rename from SimPEG/tests/TestUtils.py
rename to SimPEG/Tests/TestUtils.py
index 2723b0b7..d84408d4 100644
--- a/SimPEG/tests/TestUtils.py
+++ b/SimPEG/Tests/TestUtils.py
@@ -1,9 +1,9 @@
import numpy as np
import matplotlib.pyplot as plt
-from pylab import norm
-from SimPEG.utils import mkvc, sdiag
-from SimPEG import utils
-from SimPEG.mesh import TensorMesh, LogicallyOrthogonalMesh
+from numpy.linalg import norm
+from SimPEG.Utils import mkvc, sdiag
+from SimPEG import Utils
+from SimPEG.Mesh import TensorMesh, LogicallyOrthogonalMesh
import numpy as np
import scipy.sparse as sp
import unittest
@@ -112,10 +112,10 @@ class OrderTest(unittest.TestCase):
else:
raise Exception('Unexpected meshType')
if self.meshDimension == 2:
- X, Y = utils.exampleLomGird([nc, nc], kwrd)
+ X, Y = Utils.exampleLomGird([nc, nc], kwrd)
self.M = LogicallyOrthogonalMesh([X, Y])
if self.meshDimension == 3:
- X, Y, Z = utils.exampleLomGird([nc, nc, nc], kwrd)
+ X, Y, Z = Utils.exampleLomGird([nc, nc, nc], kwrd)
self.M = LogicallyOrthogonalMesh([X, Y, Z])
return 1./nc
@@ -211,12 +211,10 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tole
.. plot::
:include-source:
- from SimPEG.tests import checkDerivative
- from SimPEG.utils import sdiag
- import numpy as np
+ from SimPEG import Tests, Utils, np
def simplePass(x):
- return np.sin(x), sdiag(np.cos(x))
- checkDerivative(simplePass, np.random.randn(5))
+ return np.sin(x), Utils.sdiag(np.cos(x))
+ Tests.checkDerivative(simplePass, np.random.randn(5))
"""
print "%s checkDerivative %s" % ('='*20, '='*20)
diff --git a/SimPEG/Tests/__init__.py b/SimPEG/Tests/__init__.py
new file mode 100644
index 00000000..32112042
--- /dev/null
+++ b/SimPEG/Tests/__init__.py
@@ -0,0 +1,14 @@
+from TestUtils import checkDerivative, Rosenbrock, OrderTest, getQuadratic
+
+
+if __name__ == '__main__':
+ import os
+ import glob
+ import unittest
+ test_file_strings = glob.glob('test_*.py')
+ module_strings = [str[0:len(str)-3] for str in test_file_strings]
+ suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
+ in module_strings]
+ testSuite = unittest.TestSuite(suites)
+
+ unittest.TextTestRunner(verbosity=2).run(testSuite)
diff --git a/SimPEG/tests/test_LogicallyOrthogonalMesh.py b/SimPEG/Tests/test_LogicallyOrthogonalMesh.py
similarity index 98%
rename from SimPEG/tests/test_LogicallyOrthogonalMesh.py
rename to SimPEG/Tests/test_LogicallyOrthogonalMesh.py
index d760e230..241b2da8 100644
--- a/SimPEG/tests/test_LogicallyOrthogonalMesh.py
+++ b/SimPEG/Tests/test_LogicallyOrthogonalMesh.py
@@ -1,7 +1,7 @@
import numpy as np
import unittest
-from SimPEG.mesh import TensorMesh, LogicallyOrthogonalMesh
-from SimPEG.utils import ndgrid
+from SimPEG.Mesh import TensorMesh, LogicallyOrthogonalMesh
+from SimPEG.Utils import ndgrid
class BasicLOMTests(unittest.TestCase):
diff --git a/SimPEG/tests/test_Solver.py b/SimPEG/Tests/test_Solver.py
similarity index 98%
rename from SimPEG/tests/test_Solver.py
rename to SimPEG/Tests/test_Solver.py
index f8c2964d..f8a8983f 100644
--- a/SimPEG/tests/test_Solver.py
+++ b/SimPEG/Tests/test_Solver.py
@@ -1,7 +1,7 @@
import unittest
from SimPEG import Solver
-from SimPEG.mesh import TensorMesh
-from SimPEG.utils import sdiag
+from SimPEG.Mesh import TensorMesh
+from SimPEG.Utils import sdiag
import numpy as np
import scipy.sparse as sparse
diff --git a/SimPEG/tests/test_basemesh.py b/SimPEG/Tests/test_basemesh.py
similarity index 93%
rename from SimPEG/tests/test_basemesh.py
rename to SimPEG/Tests/test_basemesh.py
index d81a8553..75b67347 100644
--- a/SimPEG/tests/test_basemesh.py
+++ b/SimPEG/Tests/test_basemesh.py
@@ -1,6 +1,6 @@
import unittest
import sys
-from SimPEG.mesh import BaseMesh
+from SimPEG.Mesh import BaseMesh
import numpy as np
@@ -13,7 +13,7 @@ class TestBaseMesh(unittest.TestCase):
self.assertTrue(self.mesh.dim, 3)
def test_mesh_nc(self):
- self.assertTrue(np.all(self.mesh.n == [6, 2, 3]))
+ self.assertTrue(np.all(self.mesh.nCv == [6, 2, 3]))
def test_mesh_nc_xyz(self):
x = np.all(self.mesh.nCx == 6)
@@ -106,9 +106,9 @@ class TestBaseMesh(unittest.TestCase):
g[:, 1] = 2
g[:, 2] = 3
Xc, Yc, Zc = self.mesh.r(g, 'CC', 'CC', 'M')
- self.assertTrue(np.all(Xc.shape == self.mesh.n))
- self.assertTrue(np.all(Yc.shape == self.mesh.n))
- self.assertTrue(np.all(Zc.shape == self.mesh.n))
+ self.assertTrue(np.all(Xc.shape == self.mesh.nCv))
+ self.assertTrue(np.all(Yc.shape == self.mesh.nCv))
+ self.assertTrue(np.all(Zc.shape == self.mesh.nCv))
self.assertTrue(np.all(Xc == 1))
self.assertTrue(np.all(Yc == 2))
self.assertTrue(np.all(Zc == 3))
@@ -123,7 +123,7 @@ class TestMeshNumbers2D(unittest.TestCase):
self.assertTrue(self.mesh.dim, 2)
def test_mesh_nc(self):
- self.assertTrue(np.all(self.mesh.n == [6, 2]))
+ self.assertTrue(np.all(self.mesh.nCv == [6, 2]))
def test_mesh_nc_xyz(self):
x = np.all(self.mesh.nCx == 6)
@@ -203,8 +203,8 @@ class TestMeshNumbers2D(unittest.TestCase):
g = np.ones((self.mesh.nC, 2))
g[:, 1] = 2
Xc, Yc = self.mesh.r(g, 'CC', 'CC', 'M')
- self.assertTrue(np.all(Xc.shape == self.mesh.n))
- self.assertTrue(np.all(Yc.shape == self.mesh.n))
+ self.assertTrue(np.all(Xc.shape == self.mesh.nCv))
+ self.assertTrue(np.all(Yc.shape == self.mesh.nCv))
self.assertTrue(np.all(Xc == 1))
self.assertTrue(np.all(Yc == 2))
diff --git a/SimPEG/Tests/test_forward_DCproblem.py b/SimPEG/Tests/test_forward_DCproblem.py
new file mode 100644
index 00000000..1a44b350
--- /dev/null
+++ b/SimPEG/Tests/test_forward_DCproblem.py
@@ -0,0 +1,85 @@
+# import numpy as np
+# import unittest
+# from SimPEG.mesh import TensorMesh
+# from SimPEG.Utils import ModelBuilder, sdiag
+# from SimPEG.forward import Problem
+# from SimPEG.examples.DC import *
+# from TestUtils import checkDerivative
+# from scipy.sparse.linalg import dsolve
+# from SimPEG import inverse
+
+
+# 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 = genTxRxmat(nelec, spacelec, surfloc, elecini, mesh)
+# P = Q.T
+
+# # Create some data
+
+# problem = DCProblem(mesh)
+# problem.P = P
+# problem.RHS = q
+# data = problem.createSyntheticData(mSynth, std=0.05)
+
+# # Now set up the problem to do some minimization
+# opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
+# reg = inverse.Regularization(mesh)
+# inv = inverse.Inversion(problem, reg, opt, data, beta0=1e4)
+
+# self.inv = inv
+# self.reg = reg
+# self.p = problem
+# self.mesh = mesh
+# self.m0 = mSynth
+# self.data = data
+
+# def test_misfit(self):
+# derChk = lambda m: [self.p.dpred(m), lambda mx: self.p.J(self.m0, mx)]
+# passed = checkDerivative(derChk, self.m0, plotIt=False)
+# self.assertTrue(passed)
+
+# def test_adjoint(self):
+# # Adjoint Test
+# u = np.random.rand(self.mesh.nC*self.p.RHS.shape[1])
+# v = np.random.rand(self.mesh.nC)
+# w = np.random.rand(self.data.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)
+
+# def test_dataObj(self):
+# derChk = lambda m: [self.inv.dataObj(m), self.inv.dataObjDeriv(m)]
+# checkDerivative(derChk, self.m0, plotIt=False)
+
+# def test_modelObj(self):
+# derChk = lambda m: [self.reg.modelObj(m), self.reg.modelObjDeriv(m)]
+# checkDerivative(derChk, self.m0, plotIt=False)
+
+
+# if __name__ == '__main__':
+# unittest.main()
diff --git a/SimPEG/tests/test_interpolation.py b/SimPEG/Tests/test_interpolation.py
similarity index 99%
rename from SimPEG/tests/test_interpolation.py
rename to SimPEG/Tests/test_interpolation.py
index 3d751dde..d01dcbc3 100644
--- a/SimPEG/tests/test_interpolation.py
+++ b/SimPEG/Tests/test_interpolation.py
@@ -1,7 +1,7 @@
import numpy as np
import unittest
from TestUtils import OrderTest
-from SimPEG.utils import mkvc
+from SimPEG.Utils import mkvc
MESHTYPES = ['uniformTensorMesh', 'randomTensorMesh']
TOLERANCES = [0.9, 0.55]
diff --git a/SimPEG/tests/test_massMatrices.py b/SimPEG/Tests/test_massMatrices.py
similarity index 100%
rename from SimPEG/tests/test_massMatrices.py
rename to SimPEG/Tests/test_massMatrices.py
diff --git a/SimPEG/Tests/test_model.py b/SimPEG/Tests/test_model.py
new file mode 100644
index 00000000..0c3b2c54
--- /dev/null
+++ b/SimPEG/Tests/test_model.py
@@ -0,0 +1,27 @@
+import numpy as np
+import unittest
+from SimPEG import *
+from TestUtils import checkDerivative
+from scipy.sparse.linalg import dsolve
+
+
+class ModelTests(unittest.TestCase):
+
+ def setUp(self):
+
+ a = np.array([1, 1, 1])
+ b = np.array([1, 2])
+ c = np.array([1, 4])
+ self.mesh2 = Mesh.TensorMesh([a, b], np.array([3, 5]))
+
+ def test_modelTransforms(self):
+ print 'SimPEG.Model.BaseModel: Testing Model Transform'
+ for M in dir(Model):
+ if 'Model' not in M: continue
+ model = getattr(Model, M)(self.mesh2)
+ m = model.example()
+ passed = checkDerivative(lambda m : [model.transform(m), model.transformDeriv(m)], m, plotIt=False)
+ self.assertTrue(passed)
+
+if __name__ == '__main__':
+ unittest.main()
diff --git a/SimPEG/tests/test_operators.py b/SimPEG/Tests/test_operators.py
similarity index 100%
rename from SimPEG/tests/test_operators.py
rename to SimPEG/Tests/test_operators.py
diff --git a/SimPEG/tests/test_optimizers.py b/SimPEG/Tests/test_optimizers.py
similarity index 82%
rename from SimPEG/tests/test_optimizers.py
rename to SimPEG/Tests/test_optimizers.py
index 0253852c..3132beb9 100644
--- a/SimPEG/tests/test_optimizers.py
+++ b/SimPEG/Tests/test_optimizers.py
@@ -1,11 +1,11 @@
import unittest
from SimPEG import Solver
-from SimPEG.mesh import TensorMesh
-from SimPEG.utils import sdiag
+from SimPEG.Mesh import TensorMesh
+from SimPEG.Utils import sdiag
import numpy as np
import scipy.sparse as sp
-from SimPEG import inverse
-from SimPEG.tests import getQuadratic, Rosenbrock
+from SimPEG import Optimization
+from SimPEG.Tests import getQuadratic, Rosenbrock
TOL = 1e-2
@@ -16,7 +16,7 @@ class TestOptimizers(unittest.TestCase):
self.b = np.array([-5,-5])
def test_GN_Rosenbrock(self):
- GN = inverse.GaussNewton()
+ GN = Optimization.GaussNewton()
xopt = GN.minimize(Rosenbrock,np.array([0,0]))
x_true = np.array([1.,1.])
print 'xopt: ', xopt
@@ -24,7 +24,7 @@ class TestOptimizers(unittest.TestCase):
self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True)
def test_GN_quadratic(self):
- GN = inverse.GaussNewton()
+ GN = Optimization.GaussNewton()
xopt = GN.minimize(getQuadratic(self.A,self.b),np.array([0,0]))
x_true = np.array([5.,5.])
print 'xopt: ', xopt
@@ -32,7 +32,7 @@ class TestOptimizers(unittest.TestCase):
self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True)
def test_ProjGradient_quadraticBounded(self):
- PG = inverse.ProjectedGradient(debug=True)
+ PG = Optimization.ProjectedGradient(debug=True)
PG.lower, PG.upper = -2, 2
xopt = PG.minimize(getQuadratic(self.A,self.b),np.array([0,0]))
x_true = np.array([2.,2.])
@@ -42,7 +42,7 @@ class TestOptimizers(unittest.TestCase):
def test_ProjGradient_quadratic1Bound(self):
myB = np.array([-5,1])
- PG = inverse.ProjectedGradient()
+ PG = Optimization.ProjectedGradient()
PG.lower, PG.upper = -2, 2
xopt = PG.minimize(getQuadratic(self.A,myB),np.array([0,0]))
x_true = np.array([2.,-1.])
@@ -53,7 +53,7 @@ class TestOptimizers(unittest.TestCase):
def test_NewtonRoot(self):
fun = lambda x, return_g=True: np.sin(x) if not return_g else ( np.sin(x), sdiag( np.cos(x) ) )
x = np.array([np.pi-0.3, np.pi+0.1, 0])
- xopt = inverse.NewtonRoot(comments=False).root(fun,x)
+ xopt = Optimization.NewtonRoot(comments=False).root(fun,x)
x_true = np.array([np.pi,np.pi,0])
print 'Newton Root Finding'
print 'xopt: ', xopt
diff --git a/SimPEG/Tests/test_regularization.py b/SimPEG/Tests/test_regularization.py
new file mode 100644
index 00000000..3df9180c
--- /dev/null
+++ b/SimPEG/Tests/test_regularization.py
@@ -0,0 +1,31 @@
+import numpy as np
+import unittest
+from SimPEG import *
+from TestUtils import checkDerivative
+from scipy.sparse.linalg import dsolve
+import inspect
+
+
+class RegularizationTests(unittest.TestCase):
+
+ def setUp(self):
+ self.mesh2 = Mesh.TensorMesh([3, 2])
+
+ def test_regularization(self):
+ for R in dir(Regularization):
+ r = getattr(Regularization, R)
+ if not inspect.isclass(r): continue
+ if not issubclass(r, Regularization.BaseRegularization):
+ continue
+ # if 'Regularization' not in R: continue
+ print 'Check:', R
+ model = r.modelPair(self.mesh2)
+ reg = r(model)
+ m = model.example()
+ reg.mref = model.example()*0
+ passed = checkDerivative(lambda m : [reg.modelObj(m), reg.modelObjDeriv(m)], m, plotIt=False)
+ self.assertTrue(passed)
+
+
+if __name__ == '__main__':
+ unittest.main()
diff --git a/SimPEG/tests/test_tensorMesh.py b/SimPEG/Tests/test_tensorMesh.py
similarity index 98%
rename from SimPEG/tests/test_tensorMesh.py
rename to SimPEG/Tests/test_tensorMesh.py
index 9544dab6..3e01181b 100644
--- a/SimPEG/tests/test_tensorMesh.py
+++ b/SimPEG/Tests/test_tensorMesh.py
@@ -1,6 +1,6 @@
import numpy as np
import unittest
-from SimPEG.mesh import TensorMesh
+from SimPEG.Mesh import TensorMesh
from TestUtils import OrderTest
from scipy.sparse.linalg import dsolve
diff --git a/SimPEG/tests/test_utils.py b/SimPEG/Tests/test_utils.py
similarity index 97%
rename from SimPEG/tests/test_utils.py
rename to SimPEG/Tests/test_utils.py
index 058dba56..fea231f2 100644
--- a/SimPEG/tests/test_utils.py
+++ b/SimPEG/Tests/test_utils.py
@@ -1,7 +1,7 @@
import numpy as np
import unittest
-from SimPEG.utils import mkvc, ndgrid, indexCube, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
-from SimPEG.tests import checkDerivative
+from SimPEG.Utils import mkvc, ndgrid, indexCube, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
+from SimPEG.Tests import checkDerivative
class TestCheckDerivative(unittest.TestCase):
diff --git a/SimPEG/utils/Geophysics/__init__.py b/SimPEG/Utils/Geophysics/__init__.py
similarity index 100%
rename from SimPEG/utils/Geophysics/__init__.py
rename to SimPEG/Utils/Geophysics/__init__.py
diff --git a/SimPEG/utils/Geophysics/emSources/__init__.py b/SimPEG/Utils/Geophysics/emSources/__init__.py
similarity index 100%
rename from SimPEG/utils/Geophysics/emSources/__init__.py
rename to SimPEG/Utils/Geophysics/emSources/__init__.py
diff --git a/SimPEG/utils/Geophysics/emSources/emSources.py b/SimPEG/Utils/Geophysics/emSources/emSources.py
similarity index 100%
rename from SimPEG/utils/Geophysics/emSources/emSources.py
rename to SimPEG/Utils/Geophysics/emSources/emSources.py
diff --git a/SimPEG/utils/ModelBuilder.py b/SimPEG/Utils/ModelBuilder.py
similarity index 94%
rename from SimPEG/utils/ModelBuilder.py
rename to SimPEG/Utils/ModelBuilder.py
index 969551ed..c70c7019 100644
--- a/SimPEG/utils/ModelBuilder.py
+++ b/SimPEG/Utils/ModelBuilder.py
@@ -27,6 +27,9 @@ def getIndecesBlock(p0,p1,ccMesh):
dimMesh = np.size(ccMesh[0,:])
assert len(p0) == dimMesh, "Dimension mismatch. len(p0) != dimMesh"
+ for ii in range(len(p0)):
+ p0[ii], p1[ii] = np.min([p0[ii], p1[ii]]), np.max([p0[ii], p1[ii]])
+
if dimMesh == 1:
# Define the reference points
x1 = p0[0]
@@ -67,7 +70,7 @@ def getIndecesBlock(p0,p1,ccMesh):
# Return a tuple
return ind
-def defineBlockConductivity(p0,p1,ccMesh,condVals):
+def defineBlockConductivity(ccMesh,p0,p1,condVals):
"""
Build a block with the conductivity specified by condVal. Returns an array.
condVals[0] conductivity of the block
@@ -80,7 +83,7 @@ def defineBlockConductivity(p0,p1,ccMesh,condVals):
return sigma
-def defineTwoLayeredConductivity(depth,ccMesh,condVals):
+def defineTwoLayeredConductivity(ccMesh,depth,condVals):
"""
Define a two layered model. Depth of the first layer must be specified.
CondVals vector with the conductivity values of the layers. Eg:
@@ -151,7 +154,7 @@ def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]):
.. plot::
import matplotlib.pyplot as plt
- import SimPEG.utils.ModelBuilder as MB
+ import SimPEG.Utils.ModelBuilder as MB
plt.colorbar(plt.imshow(MB.randomModel((50,50),bounds=[-4,0])))
plt.title('A very cool, yet completely random model.')
plt.show()
@@ -196,7 +199,7 @@ def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]):
if __name__ == '__main__':
- from SimPEG.mesh import TensorMesh
+ from SimPEG.Mesh import TensorMesh
from matplotlib import pyplot as plt
# Define the mesh
@@ -229,7 +232,7 @@ if __name__ == '__main__':
p1 = np.array([1.0,1.0,1.0])[:testDim]
condVals = np.array([100,1e-6])
- sigma = defineBlockConductivity(p0,p1,ccMesh,condVals)
+ sigma = defineBlockConductivity(ccMesh,p0,p1,condVals)
# Plot sigma model
print sigma.shape
@@ -242,7 +245,7 @@ if __name__ == '__main__':
condVals = np.array([100,1e-5]);
depth = 1.0;
- sigma = defineTwoLayeredConductivity(depth,ccMesh,condVals)
+ sigma = defineTwoLayeredConductivity(ccMesh,depth,condVals)
M.plotImage(sigma)
print sigma
diff --git a/SimPEG/utils/Save.py b/SimPEG/Utils/Save.py
similarity index 98%
rename from SimPEG/utils/Save.py
rename to SimPEG/Utils/Save.py
index 25fd9b9e..a9c77191 100644
--- a/SimPEG/utils/Save.py
+++ b/SimPEG/Utils/Save.py
@@ -5,7 +5,7 @@ import re
try:
import h5py
except Exception, e:
- print 'Warning: SimPEG.utils.Save needs h5py to be installed.'
+ print 'Warning: SimPEG.Utils.Save needs h5py to be installed.'
SAVEABLES = {}
@@ -57,7 +57,7 @@ class SimPEGTable:
# At the start of every iteration we will create a inversion iteration node.
def _doStartIteration_hdf5_inv(invObj):
- invObj._invNodeIt = invObj._invNode.addGroup('%d'%(invObj._iter+1))
+ invObj._invNodeIt = invObj._invNode.addGroup('%d'%(invObj.iter+1))
preIteration(invObj._invNodeIt)
invObj.hook(_doStartIteration_hdf5_inv, overwrite=True)
@@ -78,7 +78,7 @@ class SimPEGTable:
invObj.hook(_finish_hdf5_inv, overwrite=True)
def _doStartIteration_hdf5_opt(optObj):
- optObj._optNodeIt = optObj.parent._invNode.addGroup('%d.%d'%(optObj.parent._iter, optObj._iter))
+ optObj._optNodeIt = optObj.parent._invNode.addGroup('%d.%d'%(optObj.parent.iter, optObj.iter))
preIteration(optObj._optNodeIt)
invObj.opt.hook(_doStartIteration_hdf5_opt, overwrite=True)
@@ -347,6 +347,6 @@ def loadSavable(node, pointers=None):
print 'KWARGS: ', KWARGS
return (cls, ARGS, KWARGS, node)
else:
- print 'Warning: %s Class not found in SimPEG.utils.Save.SAVABLES' % cls
+ print 'Warning: %s Class not found in SimPEG.Utils.Save.SAVABLES' % cls
return (cls, ARGS, KWARGS, node)
diff --git a/SimPEG/utils/TriSolve.f b/SimPEG/Utils/TriSolve.f
similarity index 100%
rename from SimPEG/utils/TriSolve.f
rename to SimPEG/Utils/TriSolve.f
diff --git a/SimPEG/utils/__init__.py b/SimPEG/Utils/__init__.py
similarity index 85%
rename from SimPEG/utils/__init__.py
rename to SimPEG/Utils/__init__.py
index d20c2cec..1daa2e29 100644
--- a/SimPEG/utils/__init__.py
+++ b/SimPEG/Utils/__init__.py
@@ -1,19 +1,12 @@
-import matutils
-import sputils
-import lomutils
-import interputils
-import ModelBuilder
-import meshutils
from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind
from sputils import spzeros, kron3, speye, sdiag, ddx, av, avExtrap
from meshutils import exampleLomGird, meshTensors
from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube
from interputils import interpmat
-from ipythonUtils import easyAnimate as animate
-import Solver
-from Solver import Solver
+from ipythonutils import easyAnimate as animate
import Save
import Geophysics
+import ModelBuilder
import types
import time
@@ -98,7 +91,7 @@ def callHooks(match, mainFirst=False):
def doEndIteration(self):
pass
- This will call everything named _doEndIteration* at the beginning of the function call.
+ This will call everything named _doEndIteration* at the beginning of the function call.
By default the master method (doEndIteration) is run after all of the sub methods (_doEndIteration*).
This can be reversed by adding the mainFirst=True kwarg.
"""
@@ -131,10 +124,8 @@ def callHooks(match, mainFirst=False):
Where the * can be any string. If present, _%s* will be called at the start of the default %s call.
You may also completely overwrite this function.
""" % (match, match, match, match)
- try:
- wrapper.__doc__ += extra
- except Exception, e:
- pass
+ doc = wrapper.__doc__
+ wrapper.__doc__ = ('' if doc is None else doc) + extra
return wrapper
return callHooksWrap
@@ -147,6 +138,45 @@ def dependentProperty(name, value, children, doc):
setattr(self, name, val)
return property(fget=fget, fset=fset, doc=doc)
+def requires(var):
+ """
+ Use this to wrap a funciton::
+
+ @requires('prob')
+ def dpred(self):
+ pass
+
+ This wrapper will ensure that a problem has been bound to the data.
+ If a problem is not bound an Exception will be raised, and an nice error message printed.
+ """
+ def requiresVar(f):
+ if var is 'prob':
+ extra = """
+
+ .. note::
+
+ To use data.%s(), SimPEG requires that a problem be bound to the data.
+ If a problem has not been bound, an Exception will be raised.
+ To bind a problem to the Data object::
+
+ data.pair(myProblem)
+
+ """ % f.__name__
+ else:
+ extra = """
+ To use *%s* method, SimPEG requires that the %s be specified.
+ """ % (f.__name__, var)
+ @wraps(f)
+ def requiresVarWrapper(self,*args,**kwargs):
+ if getattr(self, var, None) is None:
+ raise Exception(extra)
+ return f(self,*args,**kwargs)
+
+ doc = requiresVarWrapper.__doc__
+ requiresVarWrapper.__doc__ = ('' if doc is None else doc) + extra
+
+ return requiresVarWrapper
+ return requiresVar
class Counter(object):
"""
diff --git a/SimPEG/utils/interputils.py b/SimPEG/Utils/interputils.py
similarity index 96%
rename from SimPEG/utils/interputils.py
rename to SimPEG/Utils/interputils.py
index c8bcb4ec..38308492 100644
--- a/SimPEG/utils/interputils.py
+++ b/SimPEG/Utils/interputils.py
@@ -14,10 +14,10 @@ def _interp_point_1D(x, xr_i):
"""
# TODO: This fails if the point is on the outside of the mesh. We may want to replace this by extrapolation?
im = np.argmin(abs(x-xr_i))
- if xr_i - x[im] >= 0: # Point on the left
+ if xr_i - x[im] >= 0: # Point on the left
ind_x1 = im
ind_x2 = im+1
- elif xr_i - x[im] < 0: # Point on the right
+ elif xr_i - x[im] < 0: # Point on the right
ind_x1 = im-1
ind_x2 = im
dx1 = xr_i - x[ind_x1]
@@ -45,7 +45,7 @@ def interpmat(locs, x, y=None, z=None):
x = np.linspace(0,1,7)
dense = np.linspace(0,1,200)
fun = lambda x: np.cos(2*np.pi*x)
- Q = SimPEG.utils.interpmat(locs, x)
+ Q = SimPEG.Utils.interpmat(locs, x)
plt.plot(x, fun(x), 'bs-')
plt.plot(dense, fun(dense), 'y:')
plt.plot(locs, Q*fun(x), 'mo')
@@ -173,7 +173,7 @@ if __name__ == '__main__':
x = np.linspace(0,1,7)
dense = np.linspace(0,1,200)
fun = lambda x: np.cos(2*np.pi*x)
- Q = SimPEG.utils.interpmat(locs, x)
+ Q = SimPEG.Utils.interpmat(locs, x)
plt.plot(x, fun(x), 'bs-')
plt.plot(dense, fun(dense), 'y:')
plt.plot(locs, Q*fun(x), 'mo')
diff --git a/SimPEG/utils/ipythonUtils.py b/SimPEG/Utils/ipythonutils.py
similarity index 100%
rename from SimPEG/utils/ipythonUtils.py
rename to SimPEG/Utils/ipythonutils.py
diff --git a/SimPEG/utils/lomutils.py b/SimPEG/Utils/lomutils.py
similarity index 100%
rename from SimPEG/utils/lomutils.py
rename to SimPEG/Utils/lomutils.py
diff --git a/SimPEG/utils/matutils.py b/SimPEG/Utils/matutils.py
similarity index 100%
rename from SimPEG/utils/matutils.py
rename to SimPEG/Utils/matutils.py
diff --git a/SimPEG/utils/meshutils.py b/SimPEG/Utils/meshutils.py
similarity index 95%
rename from SimPEG/utils/meshutils.py
rename to SimPEG/Utils/meshutils.py
index 429db04b..e33c88f7 100644
--- a/SimPEG/utils/meshutils.py
+++ b/SimPEG/Utils/meshutils.py
@@ -34,8 +34,8 @@ def meshTensors(*args):
.. plot::
- from SimPEG import mesh, utils
- M = mesh.TensorMesh(utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0))))
+ from SimPEG import Mesh, Utils
+ M = Mesh.TensorMesh(Utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0))))
M.plotGrid()
"""
diff --git a/SimPEG/utils/sputils.py b/SimPEG/Utils/sputils.py
similarity index 100%
rename from SimPEG/utils/sputils.py
rename to SimPEG/Utils/sputils.py
diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py
index f421bc82..4967608e 100644
--- a/SimPEG/__init__.py
+++ b/SimPEG/__init__.py
@@ -1,13 +1,19 @@
import numpy as np
import scipy.sparse as sp
-import utils
-from utils import Solver
-import mesh
-import data
-import forward
-import inverse
-import visualize
-import examples
+import Utils
+from Solver import Solver
+import Mesh
+import Model
+import Problem
+import Data
+import Regularization
+import ObjFunction
+import Optimization
+import Inversion
+import Parameters
+import Examples
+import Tests
+
import scipy.version as _v
if _v.version < '0.13.0':
diff --git a/SimPEG/data/__init__.py b/SimPEG/data/__init__.py
deleted file mode 100644
index 86123003..00000000
--- a/SimPEG/data/__init__.py
+++ /dev/null
@@ -1,20 +0,0 @@
-from SimPEG import utils
-
-
-class SimPEGData(object):
- """Data holds the observed data, and the standard deviations."""
-
- __metaclass__ = utils.Save.Savable
-
- std = None #: Estimated Standard Deviations
- dobs = None #: Observed data
- dtrue = None #: True data, if data is synthetic
- mtrue = None #: True model, if data is synthetic
-
- def __init__(self, prob, **kwargs):
- utils.setKwargs(self, **kwargs)
- self.prob = prob
-
- def isSynthetic(self):
- "Check if the data is synthetic."
- return self.mtrue is not None
diff --git a/SimPEG/forward/ModelTransforms.py b/SimPEG/forward/ModelTransforms.py
deleted file mode 100644
index ea89b974..00000000
--- a/SimPEG/forward/ModelTransforms.py
+++ /dev/null
@@ -1,49 +0,0 @@
-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/__init__.py b/SimPEG/forward/__init__.py
deleted file mode 100644
index ce6f7e52..00000000
--- a/SimPEG/forward/__init__.py
+++ /dev/null
@@ -1,2 +0,0 @@
-from Problem import *
-import ModelTransforms
diff --git a/SimPEG/inverse/BetaSchedule.py b/SimPEG/inverse/BetaSchedule.py
deleted file mode 100644
index c225c0c0..00000000
--- a/SimPEG/inverse/BetaSchedule.py
+++ /dev/null
@@ -1,12 +0,0 @@
-
-
-class Cooling(object):
- """Simple Beta Schedule"""
-
- beta0 = None #: The initial beta value, set to none means that it will be approximated in the first iteration.
- beta_coolingFactor = 2.
-
- def getBeta(self):
- if self._beta is None:
- return self.beta0
- return self._beta / self.beta_coolingFactor
diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py
deleted file mode 100644
index 926f719a..00000000
--- a/SimPEG/inverse/Inversion.py
+++ /dev/null
@@ -1,396 +0,0 @@
-import SimPEG
-from SimPEG import utils, sp, np
-from Optimize import Remember
-from BetaSchedule import Cooling
-from SimPEG.inverse import IterationPrinters, StoppingCriteria
-
-class BaseInversion(object):
- """BaseInversion(prob, reg, opt, data, **kwargs)
- """
-
- __metaclass__ = utils.Save.Savable
-
- maxIter = 1 #: Maximum number of iterations
- name = 'BaseInversion'
-
- debug = False #: Print debugging information
-
- comment = '' #: Used by some functions to indicate what is going on in the algorithm
- counter = None #: Set this to a SimPEG.utils.Counter() if you want to count things
-
- beta0 = None #: The initial Beta (regularization parameter)
- beta0_ratio = 0.1 #: When beta0 is set to None, estimateBeta0 is used with this ratio
-
- def __init__(self, prob, reg, opt, data, **kwargs):
- utils.setKwargs(self, **kwargs)
- self.prob = prob
- self.reg = reg
- self.opt = opt
- self.data = data
- self.opt.parent = self
-
- self.stoppers = [StoppingCriteria.iteration]
-
- # Check if we have inserted printers into the optimization
- if IterationPrinters.phi_d not in self.opt.printers:
- self.opt.printers.insert(1,IterationPrinters.beta)
- self.opt.printers.insert(2,IterationPrinters.phi_d)
- self.opt.printers.insert(3,IterationPrinters.phi_m)
-
- if not hasattr(opt, '_bfgsH0') and hasattr(opt, 'bfgsH0'): # Check if it has been set by the user and the default is not being used.
- print 'Setting bfgsH0 to the inverse of the modelObj2Deriv. Done using direct methods.'
- opt.bfgsH0 = SimPEG.Solver(reg.modelObj2Deriv())
-
-
- @property
- def Wd(self):
- """
- Standard deviation weighting matrix.
- """
- if getattr(self,'_Wd',None) is None:
- eps = np.linalg.norm(utils.mkvc(self.data.dobs),2)*1e-5
- self._Wd = 1/(abs(self.data.dobs)*self.data.std+eps)
- return self._Wd
- @Wd.setter
- def Wd(self, value):
- self._Wd = value
-
- @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.data.dobs.size #
- return self._phi_d_target
-
- @phi_d_target.setter
- def phi_d_target(self, value):
- self._phi_d_target = value
-
- @utils.timeIt
- def run(self, m0):
- """run(m0)
-
- Runs the inversion!
-
- """
- self.startup(m0)
- while True:
- self.doStartIteration()
- self.m = self.opt.minimize(self.evalFunction, self.m)
- self.doEndIteration()
- if self.stoppingCriteria(): break
-
- self.printDone()
- self.finish()
-
- return self.m
-
- @utils.callHooks('startup')
- def startup(self, m0):
- """
- **startup** is called at the start of any new run call.
-
- :param numpy.ndarray x0: initial x
- :rtype: None
- :return: None
- """
-
- if not hasattr(self.reg, '_mref'):
- print 'Regularization has not set mref. SimPEG will set it to m0.'
- self.reg.mref = m0
-
- self.m = m0
- self._iter = 0
- self._beta = None
- self.phi_d_last = np.nan
- self.phi_m_last = np.nan
-
- @utils.callHooks('doStartIteration')
- def doStartIteration(self):
- """
- **doStartIteration** is called at the end of each run iteration.
-
- :rtype: None
- :return: None
- """
- self._beta = self.getBeta()
-
-
- @utils.callHooks('doEndIteration')
- def doEndIteration(self):
- """
- **doEndIteration** is called at the end of each run iteration.
-
- :rtype: None
- :return: None
- """
- # store old values
- self.phi_d_last = self.phi_d
- self.phi_m_last = self.phi_m
- self._iter += 1
-
- def getBeta(self):
- return self.beta0
-
- def estimateBeta0(self, u=None, ratio=0.1):
- """estimateBeta0(u=None, ratio=0.1)
-
- The initial beta is calculated by comparing the estimated
- eigenvalues of JtJ and WtW.
-
- To estimate the eigenvector of **A**, we will use one iteration
- of the *Power Method*:
-
- .. math::
-
- \mathbf{x_1 = A x_0}
-
- Given this (very course) approximation of the eigenvector,
- we can use the *Rayleigh quotient* to approximate the largest eigenvalue.
-
- .. math::
-
- \lambda_0 = \\frac{\mathbf{x^\\top A x}}{\mathbf{x^\\top x}}
-
- We will approximate the largest eigenvalue for both JtJ and WtW, and
- use some ratio of the quotient to estimate beta0.
-
- .. math::
-
- \\beta_0 = \gamma \\frac{\mathbf{x^\\top J^\\top J x}}{\mathbf{x^\\top W^\\top W x}}
-
-
- :param numpy.array u: fields
- :param float ratio: desired ratio of the eigenvalues, default is 0.1
- :rtype: float
- :return: beta0
- """
- if u is None:
- u = self.prob.field(self.m)
-
- x0 = np.random.rand(*self.m.shape)
- t = x0.dot(self.dataObj2Deriv(self.m,x0,u=u))
- b = x0.dot(self.reg.modelObj2Deriv()*x0)
- return ratio*(t/b)
-
- def stoppingCriteria(self):
- if self.debug: print 'checking stoppingCriteria'
- return utils.checkStoppers(self, self.stoppers)
-
-
- def printDone(self):
- """
- **printDone** is called at the end of the inversion routine.
-
- """
- utils.printStoppers(self, self.stoppers)
-
- @utils.callHooks('finish')
- def finish(self):
- """finish()
-
- **finish** is called at the end of the optimization.
- """
- pass
-
- @utils.timeIt
- def evalFunction(self, m, return_g=True, return_H=True):
- """evalFunction(m, return_g=True, return_H=True)
-
-
- """
-
- u = self.prob.field(m)
-
- if self._iter is 0 and self._beta is None:
- self._beta = self.beta0 = self.estimateBeta0(u=u,ratio=self.beta0_ratio)
-
- phi_d = self.dataObj(m, u)
- phi_m = self.reg.modelObj(m)
-
- self.dpred = self.prob.dpred(m, u=u) # This is a cheap matrix vector calculation.
- self.phi_d = phi_d
- self.phi_m = 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()*v
-
- return phi_d2Deriv + self._beta * phi_m2Deriv
-
- operator = sp.linalg.LinearOperator( (m.size, m.size), H_fun, dtype=m.dtype )
- out += (operator,)
- return out if len(out) > 1 else out[0]
-
- @utils.timeIt
- def dataObj(self, m, u=None):
- """dataObj(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, self.data, u=u)
- R = utils.mkvc(R)
- return 0.5*np.vdot(R, R)
-
- @utils.timeIt
- def dataObjDeriv(self, m, u=None):
- """dataObjDeriv(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, self.data, u=u)
-
- dmisfit = self.prob.Jt(m, self.Wd * R, u=u)
-
- return dmisfit
-
- @utils.timeIt
- def dataObj2Deriv(self, m, v, u=None):
- """dataObj2Deriv(m, v, u=None)
-
- :param numpy.array m: geophysical model
- :param numpy.array v: vector to multiply
- :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, self.data, 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
-
- def save(self, group):
- group.attrs['phi_d'] = self.phi_d
- group.attrs['phi_m'] = self.phi_m
- group.setArray('m', self.m)
- group.setArray('dpred', self.dpred)
-
-class Inversion(Cooling, Remember, BaseInversion):
-
- maxIter = 10
- name = "SimPEG Inversion"
-
- def __init__(self, prob, reg, opt, data, **kwargs):
- BaseInversion.__init__(self, prob, reg, opt, data, **kwargs)
-
- self.stoppers.append(StoppingCriteria.phi_d_target_Inversion)
-
- if StoppingCriteria.phi_d_target_Minimize not in self.opt.stoppers:
- self.opt.stoppers.append(StoppingCriteria.phi_d_target_Minimize)
-
-class TimeSteppingInversion(Remember, BaseInversion):
- """
- A slightly different view on regularization parameters,
- let Beta be viewed as 1/dt, and timestep by updating the
- reference model every optimization iteration.
- """
- maxIter = 1
- name = "Time-Stepping SimPEG Inversion"
-
- def __init__(self, prob, reg, opt, data, **kwargs):
- BaseInversion.__init__(self, prob, reg, opt, data, **kwargs)
-
- self.stoppers.append(StoppingCriteria.phi_d_target_Inversion)
-
- if StoppingCriteria.phi_d_target_Minimize not in self.opt.stoppers:
- self.opt.stoppers.append(StoppingCriteria.phi_d_target_Minimize)
-
- def _startup_TimeSteppingInversion(self, m0):
-
- def _doEndIteration_updateMref(self, xt):
- if self.debug: 'Updating the reference model.'
- self.parent.reg.mref = self.xc
-
- self.opt.hook(_doEndIteration_updateMref, overwrite=True)
diff --git a/SimPEG/inverse/__init__.py b/SimPEG/inverse/__init__.py
deleted file mode 100644
index d83a8269..00000000
--- a/SimPEG/inverse/__init__.py
+++ /dev/null
@@ -1,4 +0,0 @@
-from Optimize import *
-from Inversion import *
-from Regularization import Regularization
-import BetaSchedule
diff --git a/SimPEG/license.txt b/SimPEG/license.txt
deleted file mode 100644
index 3110a4fb..00000000
--- a/SimPEG/license.txt
+++ /dev/null
@@ -1,21 +0,0 @@
-The MIT License (MIT)
-
-Copyright (c) 2013 SimPEG Developers
-
-Permission is hereby granted, free of charge, to any person obtaining a copy
-of this software and associated documentation files (the "Software"), to deal
-in the Software without restriction, including without limitation the rights
-to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-copies of the Software, and to permit persons to whom the Software is
-furnished to do so, subject to the following conditions:
-
-The above copyright notice and this permission notice shall be included in
-all copies or substantial portions of the Software.
-
-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
-THE SOFTWARE.
diff --git a/SimPEG/setup.py b/SimPEG/setup.py
index c421cb4b..5afc57a4 100644
--- a/SimPEG/setup.py
+++ b/SimPEG/setup.py
@@ -1,8 +1,8 @@
import os
print 'Compiling TriSolve.'
-os.system('f2py -c utils/TriSolve.f -m TriSolve')
+os.system('f2py -c Utils/TriSolve.f -m TriSolve')
print 'TriSolve Compiled! yay.'
print 'Moving TriSolve into Utils.'
-os.system('mv TriSolve.so utils/TriSolve.so')
+os.system('mv TriSolve.so Utils/TriSolve.so')
print 'Thats it. Well Done Computer.'
diff --git a/SimPEG/tests/HTMLTestRunner.py b/SimPEG/tests/HTMLTestRunner.py
deleted file mode 100644
index 05ae09df..00000000
--- a/SimPEG/tests/HTMLTestRunner.py
+++ /dev/null
@@ -1,825 +0,0 @@
-"""
-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
-
-