Merge pull request #44 from simpeg/DCproblem

SimPEG Framework

.gitignore

@@ -3,4 +3,3 @@
 *.sublime-project
 *.sublime-workspace
 docs/_build/
-myNotebooks/*

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.

README.md

@@ -1,6 +1,6 @@
 ![SimPEG](https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png)
 
-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:
 [![Build Status](https://travis-ci.org/simpeg/simpeg.png)](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)

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()

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()
 

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)

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)

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()

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"

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):

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
 
 

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')

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()

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

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)

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)))

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

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

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

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
 
 
 

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
-

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.

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)

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)

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):

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
 

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))
 

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()

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]

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()

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

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()

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
 

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):

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

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)
 

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):
     """

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')

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()
 
     """

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':

SimPEG/data/__init__.py

@@ -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

SimPEG/forward/ModelTransforms.py

@@ -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)))

SimPEG/forward/__init__.py

@@ -1,2 +0,0 @@
-from Problem import *
-import ModelTransforms

SimPEG/inverse/BetaSchedule.py

@@ -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

SimPEG/inverse/Inversion.py

@@ -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)

SimPEG/inverse/__init__.py

@@ -1,4 +0,0 @@
-from Optimize import *
-from Inversion import *
-from Regularization import Regularization
-import BetaSchedule

SimPEG/license.txt

@@ -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.

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.'
 

SimPEG/tests/HTMLTestRunner.py

@@ -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 = '<link rel="stylesheet" href="my_stylesheet.css" type="text/css">'
-
-    # 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 <script> block as CDATA.
-
-Version in 0.7.1
-* Back port to Python 2.3 (Frank Horowitz).
-* Fix missing scroll bars in detail log (Podi).
-"""
-
-# TODO: color stderr
-# TODO: simplify javascript using ,ore than 1 class in the class attribute?
-
-import datetime
-import StringIO
-import sys
-import time
-import unittest
-from xml.sax import saxutils
-
-
-# ------------------------------------------------------------------------
-# The redirectors below are used to capture output during testing. Output
-# sent to sys.stdout and sys.stderr are automatically captured. However
-# in some cases sys.stdout is already cached before HTMLTestRunner is
-# invoked (e.g. calling logging.basicConfig). In order to capture those
-# output, use the redirectors for the cached stream.
-#
-# e.g.
-#   >>> logging.basicConfig(stream=HTMLTestRunner.stdout_redirector)
-#   >>>
-
-class OutputRedirector(object):
-    """ Wrapper to redirect stdout or stderr """
-    def __init__(self, fp):
-        self.fp = fp
-
-    def write(self, s):
-        self.fp.write(s)
-
-    def writelines(self, lines):
-        self.fp.writelines(lines)
-
-    def flush(self):
-        self.fp.flush()
-
-stdout_redirector = OutputRedirector(sys.stdout)
-stderr_redirector = OutputRedirector(sys.stderr)
-
-
-
-# ----------------------------------------------------------------------
-# Template
-
-class Template_mixin(object):
-    """
-    Define a HTML template for report customerization and generation.
-
-    Overall structure of an HTML report
-
-    HTML
-    +------------------------+
-    |<html>                  |
-    |  <head>                |
-    |                        |
-    |   STYLESHEET           |
-    |   +----------------+   |
-    |   |                |   |
-    |   +----------------+   |
-    |                        |
-    |  </head>               |
-    |                        |
-    |  <body>                |
-    |                        |
-    |   HEADING              |
-    |   +----------------+   |
-    |   |                |   |
-    |   +----------------+   |
-    |                        |
-    |   REPORT               |
-    |   +----------------+   |
-    |   |                |   |
-    |   +----------------+   |
-    |                        |
-    |   ENDING               |
-    |   +----------------+   |
-    |   |                |   |
-    |   +----------------+   |
-    |                        |
-    |  </body>               |
-    |</html>                 |
-    +------------------------+
-    """
-
-    STATUS = {
-    0: 'pass',
-    1: 'fail',
-    2: 'error',
-    }
-
-    DEFAULT_TITLE = 'Unit Test Report'
-    DEFAULT_DESCRIPTION = ''
-
-    # ------------------------------------------------------------------------
-    # HTML Template
-
-    HTML_TMPL = r"""<?xml version="1.0" encoding="UTF-8"?>
-<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Strict//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-strict.dtd">
-<html xmlns="http://www.w3.org/1999/xhtml">
-<head>
-    <title>%(title)s</title>
-    <meta name="generator" content="%(generator)s"/>
-    <meta http-equiv="Content-Type" content="text/html; charset=UTF-8"/>
-    %(stylesheet)s
-</head>
-<body>
-<script language="javascript" type="text/javascript"><!--
-output_list = Array();
-
-/* level - 0:Summary; 1:Failed; 2:All */
-function showCase(level) {
-    trs = document.getElementsByTagName("tr");
-    for (var i = 0; i < trs.length; i++) {
-        tr = trs[i];
-        id = tr.id;
-        if (id.substr(0,2) == 'ft') {
-            if (level < 1) {
-                tr.className = 'hiddenRow';
-            }
-            else {
-                tr.className = '';
-            }
-        }
-        if (id.substr(0,2) == 'pt') {
-            if (level > 1) {
-                tr.className = '';
-            }
-            else {
-                tr.className = 'hiddenRow';
-            }
-        }
-    }
-}
-
-
-function showClassDetail(cid, count) {
-    var id_list = Array(count);
-    var toHide = 1;
-    for (var i = 0; i < count; i++) {
-        tid0 = 't' + cid.substr(1) + '.' + (i+1);
-        tid = 'f' + tid0;
-        tr = document.getElementById(tid);
-        if (!tr) {
-            tid = 'p' + tid0;
-            tr = document.getElementById(tid);
-        }
-        id_list[i] = tid;
-        if (tr.className) {
-            toHide = 0;
-        }
-    }
-    for (var i = 0; i < count; i++) {
-        tid = id_list[i];
-        if (toHide) {
-            var divTid = document.getElementById('div_'+tid);
-            if(divTid !== null){divTid.style.display = 'none';}
-            document.getElementById(tid).className = 'hiddenRow';
-        }
-        else {
-            document.getElementById(tid).className = '';
-        }
-    }
-}
-
-
-function showTestDetail(div_id){
-    var details_div = document.getElementById(div_id)
-    var displayState = details_div.style.display
-    // alert(displayState)
-    if (displayState != 'block' ) {
-        displayState = 'block'
-        details_div.style.display = 'block'
-    }
-    else {
-        details_div.style.display = 'none'
-    }
-}
-
-
-function html_escape(s) {
-    s = s.replace(/&/g,'&amp;');
-    s = s.replace(/</g,'&lt;');
-    s = s.replace(/>/g,'&gt;');
-    return s;
-}
-
-/* obsoleted by detail in <div>
-function showOutput(id, name) {
-    var w = window.open("", //url
-                    name,
-                    "resizable,scrollbars,status,width=800,height=450");
-    d = w.document;
-    d.write("<pre>");
-    d.write(html_escape(output_list[id]));
-    d.write("\n");
-    d.write("<a href='javascript:window.close()'>close</a>\n");
-    d.write("</pre>\n");
-    d.close();
-}
-*/
---></script>
-
-%(heading)s
-%(report)s
-%(ending)s
-
-</body>
-</html>
-"""
-    # variables: (title, generator, stylesheet, heading, report, ending)
-
-
-    # ------------------------------------------------------------------------
-    # Stylesheet
-    #
-    # alternatively use a <link> for external style sheet, e.g.
-    #   <link rel="stylesheet" href="$url" type="text/css">
-
-    STYLESHEET_TMPL = """
-<style type="text/css" media="screen">
-body        { font-family: verdana, arial, helvetica, sans-serif; font-size: 80%; }
-table       { font-size: 100%; }
-pre         { }
-
-/* -- heading ---------------------------------------------------------------------- */
-h1 {
-    font-size: 16pt;
-    color: gray;
-}
-.heading {
-    margin-top: 0ex;
-    margin-bottom: 1ex;
-}
-
-.heading .attribute {
-    margin-top: 1ex;
-    margin-bottom: 0;
-}
-
-.heading .description {
-    margin-top: 4ex;
-    margin-bottom: 6ex;
-}
-
-/* -- css div popup ------------------------------------------------------------------------ */
-a.popup_link {
-}
-
-a.popup_link:hover {
-    color: red;
-}
-
-.popup_window {
-    display: none;
-    position: relative;
-    left: 0px;
-    top: 0px;
-    /*border: solid #627173 1px; */
-    padding: 10px;
-    background-color: #E6E6D6;
-    font-family: "Lucida Console", "Courier New", Courier, monospace;
-    text-align: left;
-    font-size: 8pt;
-    width: 500px;
-}
-
-}
-/* -- report ------------------------------------------------------------------------ */
-#show_detail_line {
-    margin-top: 3ex;
-    margin-bottom: 1ex;
-}
-#result_table {
-    width: 80%;
-    border-collapse: collapse;
-    border: 1px solid #777;
-}
-#header_row {
-    font-weight: bold;
-    color: white;
-    background-color: #777;
-}
-#result_table td {
-    border: 1px solid #777;
-    padding: 2px;
-}
-#total_row  { font-weight: bold; }
-.passClass  { background-color: #6c6; }
-.failClass  { background-color: #c60; }
-.errorClass { background-color: #c00; }
-.passCase   { color: #6c6; }
-.failCase   { color: #c60; font-weight: bold; }
-.errorCase  { color: #c00; font-weight: bold; }
-.hiddenRow  { display: none; }
-.testcase   { margin-left: 2em; }
-
-
-/* -- ending ---------------------------------------------------------------------- */
-#ending {
-}
-
-</style>
-"""
-
-
-
-    # ------------------------------------------------------------------------
-    # Heading
-    #
-
-    HEADING_TMPL = """<div class='heading'>
-<h1>%(title)s</h1>
-%(parameters)s
-<p class='description'>%(description)s</p>
-</div>
-
-""" # variables: (title, parameters, description)
-
-    HEADING_ATTRIBUTE_TMPL = """<p class='attribute'><strong>%(name)s:</strong> %(value)s</p>
-""" # variables: (name, value)
-
-
-
-    # ------------------------------------------------------------------------
-    # Report
-    #
-
-    REPORT_TMPL = """
-<p id='show_detail_line'>Show
-<a href='javascript:showCase(0)'>Summary</a>
-<a href='javascript:showCase(1)'>Failed</a>
-<a href='javascript:showCase(2)'>All</a>
-</p>
-<table id='result_table'>
-<colgroup>
-<col align='left' />
-<col align='right' />
-<col align='right' />
-<col align='right' />
-<col align='right' />
-<col align='right' />
-</colgroup>
-<tr id='header_row'>
-    <td>Test Group/Test case</td>
-    <td>Count</td>
-    <td>Pass</td>
-    <td>Fail</td>
-    <td>Error</td>
-    <td>View</td>
-</tr>
-%(test_list)s
-<tr id='total_row'>
-    <td>Total</td>
-    <td>%(count)s</td>
-    <td>%(Pass)s</td>
-    <td>%(fail)s</td>
-    <td>%(error)s</td>
-    <td>&nbsp;</td>
-</tr>
-</table>
-""" # variables: (test_list, count, Pass, fail, error)
-
-    REPORT_CLASS_TMPL = r"""
-<tr class='%(style)s'>
-    <td>%(desc)s</td>
-    <td>%(count)s</td>
-    <td>%(Pass)s</td>
-    <td>%(fail)s</td>
-    <td>%(error)s</td>
-    <td><a href="javascript:showClassDetail('%(cid)s',%(count)s)">Detail</a></td>
-</tr>
-""" # variables: (style, desc, count, Pass, fail, error, cid)
-
-
-    REPORT_TEST_WITH_OUTPUT_TMPL = r"""
-<tr id='%(tid)s' class='%(Class)s'>
-    <td class='%(style)s'><div class='testcase'>%(desc)s</div></td>
-    <td colspan='5' align='center'>
-
-    <!--css div popup start-->
-    <a class="popup_link" onfocus='this.blur();' href="javascript:showTestDetail('div_%(tid)s')" >
-        %(status)s</a>
-
-    <div id='div_%(tid)s' class="popup_window">
-        <div style='text-align: right; color:red;cursor:pointer'>
-        <a onfocus='this.blur();' onclick="document.getElementById('div_%(tid)s').style.display = 'none' " >
-           [x]</a>
-        </div>
-        <pre>
-        %(script)s
-        </pre>
-    </div>
-    <!--css div popup end-->
-
-    </td>
-</tr>
-""" # variables: (tid, Class, style, desc, status)
-
-
-    REPORT_TEST_NO_OUTPUT_TMPL = r"""
-<tr id='%(tid)s' class='%(Class)s'>
-    <td class='%(style)s'><div class='testcase'>%(desc)s</div></td>
-    <td colspan='5' align='center'>%(status)s</td>
-</tr>
-""" # variables: (tid, Class, style, desc, status)
-
-
-    REPORT_TEST_OUTPUT_TMPL = r"""
-%(id)s: %(output)s
-""" # variables: (id, output)
-
-
-
-    # ------------------------------------------------------------------------
-    # ENDING
-    #
-
-    ENDING_TMPL = """<div id='ending'>&nbsp;</div>"""
-
-# -------------------- The end of the Template class -------------------
-
-
-TestResult = unittest.TestResult
-
-class _TestResult(TestResult):
-    # note: _TestResult is a pure representation of results.
-    # It lacks the output and reporting ability compares to unittest._TextTestResult.
-
-    def __init__(self, verbosity=1):
-        TestResult.__init__(self)
-        self.stdout0 = None
-        self.stderr0 = None
-        self.success_count = 0
-        self.failure_count = 0
-        self.error_count = 0
-        self.verbosity = verbosity
-
-        # result is a list of result in 4 tuple
-        # (
-        #   result code (0: success; 1: fail; 2: error),
-        #   TestCase object,
-        #   Test output (byte string),
-        #   stack trace,
-        # )
-        self.result = []
-
-
-    def startTest(self, test):
-        TestResult.startTest(self, test)
-        # just one buffer for both stdout and stderr
-        self.outputBuffer = StringIO.StringIO()
-        stdout_redirector.fp = self.outputBuffer
-        stderr_redirector.fp = self.outputBuffer
-        self.stdout0 = sys.stdout
-        self.stderr0 = sys.stderr
-        sys.stdout = stdout_redirector
-        sys.stderr = stderr_redirector
-
-
-    def complete_output(self):
-        """
-        Disconnect output redirection and return buffer.
-        Safe to call multiple times.
-        """
-        if self.stdout0:
-            sys.stdout = self.stdout0
-            sys.stderr = self.stderr0
-            self.stdout0 = None
-            self.stderr0 = None
-        return self.outputBuffer.getvalue()
-
-
-    def stopTest(self, test):
-        # Usually one of addSuccess, addError or addFailure would have been called.
-        # But there are some path in unittest that would bypass this.
-        # We must disconnect stdout in stopTest(), which is guaranteed to be called.
-        self.complete_output()
-
-
-    def addSuccess(self, test):
-        self.success_count += 1
-        TestResult.addSuccess(self, test)
-        output = self.complete_output()
-        self.result.append((0, test, output, ''))
-        if self.verbosity > 1:
-            sys.stderr.write('ok ')
-            sys.stderr.write(str(test))
-            sys.stderr.write('\n')
-        else:
-            sys.stderr.write('.')
-
-    def addError(self, test, err):
-        self.error_count += 1
-        TestResult.addError(self, test, err)
-        _, _exc_str = self.errors[-1]
-        output = self.complete_output()
-        self.result.append((2, test, output, _exc_str))
-        if self.verbosity > 1:
-            sys.stderr.write('E  ')
-            sys.stderr.write(str(test))
-            sys.stderr.write('\n')
-        else:
-            sys.stderr.write('E')
-
-    def addFailure(self, test, err):
-        self.failure_count += 1
-        TestResult.addFailure(self, test, err)
-        _, _exc_str = self.failures[-1]
-        output = self.complete_output()
-        self.result.append((1, test, output, _exc_str))
-        if self.verbosity > 1:
-            sys.stderr.write('F  ')
-            sys.stderr.write(str(test))
-            sys.stderr.write('\n')
-        else:
-            sys.stderr.write('F')
-
-
-class HTMLTestRunner(Template_mixin):
-    """
-    """
-    def __init__(self, stream=sys.stdout, verbosity=1, title=None, description=None):
-        self.stream = stream
-        self.verbosity = verbosity
-        if title is None:
-            self.title = self.DEFAULT_TITLE
-        else:
-            self.title = title
-        if description is None:
-            self.description = self.DEFAULT_DESCRIPTION
-        else:
-            self.description = description
-
-        self.startTime = datetime.datetime.now()
-
-
-    def run(self, test):
-        "Run the given test case or test suite."
-        result = _TestResult(self.verbosity)
-        test(result)
-        self.stopTime = datetime.datetime.now()
-        self.generateReport(test, result)
-        print >>sys.stderr, '\nTime Elapsed: %s' % (self.stopTime-self.startTime)
-        return result
-
-
-    def sortResult(self, result_list):
-        # unittest does not seems to run in any particular order.
-        # Here at least we want to group them together by class.
-        rmap = {}
-        classes = []
-        for n,t,o,e in result_list:
-            cls = t.__class__
-            if not rmap.has_key(cls):
-                rmap[cls] = []
-                classes.append(cls)
-            rmap[cls].append((n,t,o,e))
-        r = [(cls, rmap[cls]) for cls in classes]
-        return r
-
-
-    def getReportAttributes(self, result):
-        """
-        Return report attributes as a list of (name, value).
-        Override this to add custom attributes.
-        """
-        startTime = str(self.startTime)[:19]
-        duration = str(self.stopTime - self.startTime)
-        status = []
-        if result.success_count: status.append('Pass %s'    % result.success_count)
-        if result.failure_count: status.append('Failure %s' % result.failure_count)
-        if result.error_count:   status.append('Error %s'   % result.error_count  )
-        if status:
-            status = ' '.join(status)
-        else:
-            status = 'none'
-        return [
-            ('Start Time', startTime),
-            ('Duration', duration),
-            ('Status', status),
-        ]
-
-
-    def generateReport(self, test, result):
-        report_attrs = self.getReportAttributes(result)
-        generator = 'HTMLTestRunner %s' % __version__
-        stylesheet = self._generate_stylesheet()
-        heading = self._generate_heading(report_attrs)
-        report = self._generate_report(result)
-        ending = self._generate_ending()
-        output = self.HTML_TMPL % dict(
-            title = saxutils.escape(self.title),
-            generator = generator,
-            stylesheet = stylesheet,
-            heading = heading,
-            report = report,
-            ending = ending,
-        )
-        self.stream.write(output.encode('utf8'))
-
-
-    def _generate_stylesheet(self):
-        return self.STYLESHEET_TMPL
-
-
-    def _generate_heading(self, report_attrs):
-        a_lines = []
-        for name, value in report_attrs:
-            line = self.HEADING_ATTRIBUTE_TMPL % dict(
-                    name = saxutils.escape(name),
-                    value = saxutils.escape(value),
-                )
-            a_lines.append(line)
-        heading = self.HEADING_TMPL % dict(
-            title = saxutils.escape(self.title),
-            parameters = ''.join(a_lines),
-            description = saxutils.escape(self.description),
-        )
-        return heading
-
-
-    def _generate_report(self, result):
-        rows = []
-        sortedResult = self.sortResult(result.result)
-        for cid, (cls, cls_results) in enumerate(sortedResult):
-            # subtotal for a class
-            np = nf = ne = 0
-            for n,t,o,e in cls_results:
-                if n == 0: np += 1
-                elif n == 1: nf += 1
-                else: ne += 1
-
-            # format class description
-            if cls.__module__ == "__main__":
-                name = cls.__name__
-            else:
-                name = "%s.%s" % (cls.__module__, cls.__name__)
-            doc = cls.__doc__ and cls.__doc__.split("\n")[0] or ""
-            desc = doc and '%s: %s' % (name, doc) or name
-
-            row = self.REPORT_CLASS_TMPL % dict(
-                style = ne > 0 and 'errorClass' or nf > 0 and 'failClass' or 'passClass',
-                desc = desc,
-                count = np+nf+ne,
-                Pass = np,
-                fail = nf,
-                error = ne,
-                cid = 'c%s' % (cid+1),
-            )
-            rows.append(row)
-
-            for tid, (n,t,o,e) in enumerate(cls_results):
-                self._generate_report_test(rows, cid, tid, n, t, o, e)
-
-        report = self.REPORT_TMPL % dict(
-            test_list = ''.join(rows),
-            count = str(result.success_count+result.failure_count+result.error_count),
-            Pass = str(result.success_count),
-            fail = str(result.failure_count),
-            error = str(result.error_count),
-        )
-        return report
-
-
-    def _generate_report_test(self, rows, cid, tid, n, t, o, e):
-        # e.g. 'pt1.1', 'ft1.1', etc
-        has_output = bool(o or e)
-        tid = (n == 0 and 'p' or 'f') + 't%s.%s' % (cid+1,tid+1)
-        name = t.id().split('.')[-1]
-        doc = t.shortDescription() or ""
-        desc = doc and ('%s: %s' % (name, doc)) or name
-        tmpl = has_output and self.REPORT_TEST_WITH_OUTPUT_TMPL or self.REPORT_TEST_NO_OUTPUT_TMPL
-
-        # o and e should be byte string because they are collected from stdout and stderr?
-        if isinstance(o,str):
-            # TODO: some problem with 'string_escape': it escape \n and mess up formating
-            # uo = unicode(o.encode('string_escape'))
-            uo = o.decode('latin-1')
-        else:
-            uo = o
-        if isinstance(e,str):
-            # TODO: some problem with 'string_escape': it escape \n and mess up formating
-            # ue = unicode(e.encode('string_escape'))
-            ue = e.decode('latin-1')
-        else:
-            ue = e
-
-        script = self.REPORT_TEST_OUTPUT_TMPL % dict(
-            id = tid,
-            output = saxutils.escape(uo+ue),
-        )
-
-        row = tmpl % dict(
-            tid = tid,
-            Class = (n == 0 and 'hiddenRow' or 'none'),
-            style = n == 2 and 'errorCase' or (n == 1 and 'failCase' or 'none'),
-            desc = desc,
-            script = script,
-            status = self.STATUS[n],
-        )
-        rows.append(row)
-        if not has_output:
-            return
-
-    def _generate_ending(self):
-        return self.ENDING_TMPL
-
-
-##############################################################################
-# Facilities for running tests from the command line
-##############################################################################
-
-# Note: Reuse unittest.TestProgram to launch test. In the future we may
-# build our own launcher to support more specific command line
-# parameters like test title, CSS, etc.
-class TestProgram(unittest.TestProgram):
-    """
-    A variation of the unittest.TestProgram. Please refer to the base
-    class for command line parameters.
-    """
-    def runTests(self):
-        # Pick HTMLTestRunner as the default test runner.
-        # base class's testRunner parameter is not useful because it means
-        # we have to instantiate HTMLTestRunner before we know self.verbosity.
-        if self.testRunner is None:
-            self.testRunner = HTMLTestRunner(verbosity=self.verbosity)
-        unittest.TestProgram.runTests(self)
-
-main = TestProgram
-
-##############################################################################
-# Executing this module from the command line
-##############################################################################
-
-if __name__ == "__main__":
-    main(module=None)

SimPEG/tests/__init__.py

@@ -1,2 +0,0 @@
-import TestUtils
-from TestUtils import checkDerivative, Rosenbrock, OrderTest, getQuadratic

SimPEG/tests/runTests.py

@@ -1,57 +0,0 @@
-import os
-import glob
-import unittest
-import HTMLTestRunner
-
-# This code will run all tests in directory named test_*.py
-def main(html=False):
-    TITLE = 'Test Results'
-    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)
-
-    if not html:
-        unittest.TextTestRunner(verbosity=2).run(testSuite)
-        return
-
-
-    outfile = open("report.html", "w")
-    runner = HTMLTestRunner.HTMLTestRunner(
-                    stream=outfile,
-                    title=TITLE,
-                    description='SimPEG Test Report was automatically generated.',
-                    verbosity=2
-                    )
-
-    runner.run(testSuite)
-    outfile.close()
-
-    reader = open("report.html", "r")
-    writer = open("../../docs/api_TestResults.rst", "w")
-
-    writer.write('.. _api_TestResults:\n\nTest Results\n============\n\n.. raw:: html\n\n')
-
-    go = False
-    for line in reader:
-        skip = False
-        if line == '<style type="text/css" media="screen">\n':
-            go = True
-        elif line == "<div id='ending'>&nbsp;</div>\n":
-            go = False
-        elif line == '</head>\n':
-            skip = True
-        elif line == '<h1>'+TITLE+'</h1>\n':
-            skip = True
-        elif line == '<body>\n':
-            skip = True
-        if go and not skip:
-            writer.write('    '+line)
-
-    writer.close()
-    reader.close()
-    os.remove("report.html")
-
-if __name__ == '__main__':
-    main(True)

SimPEG/tests/runTests.sh

@@ -1,3 +0,0 @@
-#!/bin/sh
-
-python -m unittest discover

SimPEG/tests/test_forward_DCproblem.py

@@ -1,85 +0,0 @@
-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()

SimPEG/tests/test_forward_problem.py

@@ -1,34 +0,0 @@
-import numpy as np
-import unittest
-from SimPEG import mesh, forward, inverse
-from TestUtils import checkDerivative
-from scipy.sparse.linalg import dsolve
-
-
-class ProblemTests(unittest.TestCase):
-
-    def setUp(self):
-
-        a = np.array([1, 1, 1])
-        b = np.array([1, 2])
-        c = np.array([1, 4])
-        self.mesh2 = mesh.TensorMesh([a, b], np.array([3, 5]))
-        self.p2 = forward.Problem(self.mesh2)
-        self.reg = inverse.Regularization(self.mesh2)
-
-    def test_modelTransform(self):
-        print 'SimPEG.forward.Problem: Testing Model Transform'
-        m = np.random.rand(self.mesh2.nC)
-        passed = checkDerivative(lambda m : [self.p2.modelTransform(m), self.p2.modelTransformDeriv(m)], m, plotIt=False)
-        self.assertTrue(passed)
-
-    def test_regularization(self):
-        derChk = lambda m: [self.reg.modelObj(m), self.reg.modelObjDeriv(m)]
-        mSynth = np.random.randn(self.mesh2.nC)
-        checkDerivative(derChk, mSynth, plotIt=False)
-
-
-
-
-if __name__ == '__main__':
-    unittest.main()

SimPEG/visualize/vtk/vtkTools.py

@@ -3,7 +3,7 @@ try:
 	import vtk, vtk.util.numpy_support as npsup, pdb
 except Exception, e:
 	print 'VTK import error. Please ensure you have VTK installed to use this visualization package.'
-from SimPEG.utils import mkvc
+from SimPEG.Utils import mkvc
 
 
 class vtkTools(object):

docs/api_BaseMesh.rst

@@ -3,6 +3,6 @@
 Base Mesh
 *********
 
-.. automodule:: SimPEG.mesh.BaseMesh
+.. automodule:: SimPEG.Mesh.BaseMesh
     :members:
     :undoc-members:

docs/api_Cyl1DMesh.rst

@@ -3,6 +3,6 @@
 Cylindrical 1D Mesh
 *******************
 
-.. automodule:: SimPEG.mesh.Cyl1DMesh
+.. automodule:: SimPEG.Mesh.Cyl1DMesh
     :members:
-    :undoc-members:
+    :undoc-members:

docs/api_DiffOperators.rst

@@ -3,6 +3,6 @@
 Differential Operators
 **********************
 
-.. automodule:: SimPEG.mesh.DiffOperators
+.. automodule:: SimPEG.Mesh.DiffOperators
     :members:
     :undoc-members:

docs/api_Forward.rst

@@ -1,32 +1,30 @@
-.. _api_Problem:
+.. _api_Forward:
 
 
+Model
+*****
+
+.. automodule:: SimPEG.Model
+    :show-inheritance:
+    :members:
+    :undoc-members:
+    :inherited-members:
+
+Data
+****
+
+.. automodule:: SimPEG.Data
+    :show-inheritance:
+    :members:
+    :undoc-members:
+    :inherited-members:
+
 Problem
 *******
 
-.. automodule:: SimPEG.forward.Problem
+.. automodule:: SimPEG.Problem
     :show-inheritance:
     :members:
     :undoc-members:
     :inherited-members:
 
-
-DCProblem
-*********
-
-.. automodule:: SimPEG.examples.DC
-    :show-inheritance:
-    :members:
-    :undoc-members:
-    :inherited-members:
-
-
-
-Linear Problem
-**************
-
-.. automodule:: SimPEG.examples.Linear
-    :show-inheritance:
-    :members:
-    :undoc-members:
-    :inherited-members:

docs/api_InnerProducts.rst

@@ -3,6 +3,6 @@
 Inner Products
 **************
 
-.. automodule:: SimPEG.mesh.InnerProducts
+.. automodule:: SimPEG.Mesh.InnerProducts
     :members:
     :undoc-members:

docs/api_Inverse.rst

@@ -1,35 +1,38 @@
 .. _api_Inverse:
 
-Optimize
-********
-
-.. automodule:: SimPEG.inverse.Optimize
-    :show-inheritance:
-    :private-members:
-    :members:
-    :undoc-members:
-
-
-Inversion
-*********
-
-.. automodule:: SimPEG.inverse.Inversion
-    :show-inheritance:
-    :members:
-    :undoc-members:
-
-Beta Schedule
-*************
-
-.. automodule:: SimPEG.inverse.BetaSchedule
-    :members:
-    :undoc-members:
-
 
 Regularization
 **************
 
-.. automodule:: SimPEG.inverse.Regularization
+.. automodule:: SimPEG.Regularization
     :show-inheritance:
     :members:
     :undoc-members:
+
+
+Objective Function
+******************
+
+.. automodule:: SimPEG.ObjFunction
+    :members:
+    :undoc-members:
+
+
+Optimize
+********
+
+.. automodule:: SimPEG.Optimization
+    :show-inheritance:
+    :private-members:
+    :members:
+    :undoc-members:
+
+Inversion
+*********
+
+.. automodule:: SimPEG.Inversion
+    :show-inheritance:
+    :members:
+    :undoc-members:
+
+

docs/api_LogicallyOrthogonalMesh.rst

@@ -3,16 +3,8 @@
 Logically Orthogonal Mesh
 *************************
 
-.. automodule:: SimPEG.mesh.LogicallyOrthogonalMesh
+.. automodule:: SimPEG.Mesh.LogicallyOrthogonalMesh
     :show-inheritance:
     :members:
     :undoc-members:
     :inherited-members:
-
-
-LOM View
-********
-
-.. automodule:: SimPEG.mesh.LomView
-    :members:
-    :undoc-members:

docs/api_Parameters.rst

@@ -0,0 +1,11 @@
+.. _api_Parameters:
+
+
+Parameters
+**********
+
+.. automodule:: SimPEG.Parameters
+    :show-inheritance:
+    :members:
+    :undoc-members:
+    :inherited-members:

docs/api_TensorMesh.rst

@@ -3,15 +3,8 @@
 Tensor Mesh
 ***********
 
-.. automodule:: SimPEG.mesh.TensorMesh
+.. automodule:: SimPEG.Mesh.TensorMesh
     :show-inheritance:
     :members:
     :undoc-members:
     :inherited-members:
-
-Tensor View
-***********
-
-.. automodule:: SimPEG.mesh.TensorView
-    :members:
-    :undoc-members:

docs/api_Tests.rst

@@ -3,6 +3,6 @@
 Testing SimPEG
 **************
 
-.. automodule:: SimPEG.tests.TestUtils
+.. automodule:: SimPEG.Tests.TestUtils
     :members:
     :undoc-members:

docs/api_Utils.rst

@@ -4,7 +4,7 @@
 Solver
 ******
 
-.. automodule:: SimPEG.utils.Solver
+.. automodule:: SimPEG.Solver
     :members:
     :undoc-members:
 
@@ -12,49 +12,49 @@ Solver
 Utilities
 *********
 
-.. automodule:: SimPEG.utils
+.. automodule:: SimPEG.Utils
     :members:
     :undoc-members:
 
 Matrix Utilities
 ****************
 
-.. automodule:: SimPEG.utils.matutils
+.. automodule:: SimPEG.Utils.matutils
     :members:
     :undoc-members:
 
 Sparse Utilities
 ****************
 
-.. automodule:: SimPEG.utils.sputils
+.. automodule:: SimPEG.Utils.sputils
     :members:
     :undoc-members:
 
 LOM Utilities
 *************
 
-.. automodule:: SimPEG.utils.lomutils
+.. automodule:: SimPEG.Utils.lomutils
     :members:
     :undoc-members:
 
 Mesh Utilities
 **************
 
-.. automodule:: SimPEG.utils.meshutils
+.. automodule:: SimPEG.Utils.meshutils
     :members:
     :undoc-members:
 
 Model Builder Utilities
 ***********************
 
-.. automodule:: SimPEG.utils.ModelBuilder
+.. automodule:: SimPEG.Utils.ModelBuilder
     :members:
     :undoc-members:
 
 Interpolation Utilities
 ***********************
 
-.. automodule:: SimPEG.utils.interputils
+.. automodule:: SimPEG.Utils.interputils
     :members:
     :undoc-members:
 

docs/examples/mesh/plot_LogicallyOrthogonalMesh.py

@@ -1,7 +0,0 @@
-from SimPEG.mesh import LogicallyOrthogonalMesh
-from SimPEG import utils
-import matplotlib.pyplot as plt
-X, Y = utils.exampleLomGird([3,3],'rotate')
-M = LogicallyOrthogonalMesh([X, Y])
-M.plotGrid()
-plt.show()

docs/examples/mesh/plot_grid_2D.py

@@ -1,11 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-from SimPEG.mesh import TensorMesh
-
-h1 = np.linspace(.1,.5,3)
-h2 = np.linspace(.1,.5,5)
-mesh = TensorMesh([h1, h2])
-mesh.plotGrid(nodes=True, faces=True, centers=True, lines=True)
-
-plt.show()
-

docs/examples/mesh/plot_grid_3D.py

@@ -1,12 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-from SimPEG.mesh import TensorMesh
-
-h1 = np.linspace(.1,.5,3)
-h2 = np.linspace(.1,.5,5)
-h3 = np.linspace(.1,.5,3)
-mesh = TensorMesh([h1,h2,h3])
-mesh.plotGrid(nodes=True, faces=True, centers=True, lines=True)
-
-plt.show()
-

docs/examples/mesh/plot_image_2D.py

@@ -1,11 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-from SimPEG.mesh import TensorMesh
-
-n = 20
-h = np.ones(n)/n
-M = TensorMesh([h, h])
-I = np.sin(M.gridCC[:,0]*2*np.pi)*np.sin(M.gridCC[:,1]*2*np.pi)
-M.plotImage(I)
-
-plt.show()

docs/examples/mesh/plot_image_3D.py

@@ -1,12 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-from SimPEG.mesh import TensorMesh
-
-n = 20
-h = np.ones(n)/n
-M = TensorMesh([h,h,h])
-
-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')
-
-plt.show()

docs/examples/mesh/plot_nodes.py

@@ -1,13 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-from SimPEG.mesh import TensorMesh
-
-x0 = np.zeros(2)
-h1 = np.linspace(.1,.5,3)
-h2 = np.linspace(.1,.5,5)
-M = TensorMesh([h1,h2],x0)
-M.plotGrid()
-plt.hold()
-plt.plot(M.gridN[:,0], M.gridN[:,1], 'ks', markersize=10)
-plt.show()
-

docs/index.rst

@@ -3,7 +3,7 @@
    :alt: SimPEG
    :align: center
 
-SimPEG (Simulation and Parameter Estimation in Geoscience) is a python
+SimPEG (Simulation and Parameter Estimation in Geophysics) is a python
 package for simulation and gradient based parameter estimation in the
 context of geoscience applications.
 
@@ -16,11 +16,10 @@ these goals, this package has the following features:
 * provides a framework for geophysical and hydrogeologic problems
 * supports 1D, 2D and 3D problems
 
-
-.. raw:: html
-
-   <iframe src="http://row1.ca/labs/simpegvis" width="100%" height="500px" style="border:none;background:#eee;margin:50px 0;"></iframe>
-
+.. image:: simpeg-framework.png
+   :width: 400 px
+   :alt: Framework
+   :align: center
 
 Meshing & Operators
 ===================
@@ -41,7 +40,7 @@ Forward Problems
 .. toctree::
    :maxdepth: 2
 
-   api_Problem
+   api_Forward
 
 Inversion
 =========
@@ -49,7 +48,8 @@ Inversion
 .. toctree::
    :maxdepth: 2
 
-   api_Optimize
+   api_Inverse
+   api_Parameters
 
 Testing SimPEG
 ==============
@@ -58,22 +58,18 @@ Testing SimPEG
    :maxdepth: 2
 
    api_Tests
-   api_TestResults
-
-Build Results
-=============
 
 * Master Branch
-.. image:: https://travis-ci.org/simpeg/simpeg.png?branch=master
-   :target: https://travis-ci.org/simpeg/simpeg
-   :alt: Master Branch
-   :align: center
+   .. image:: https://travis-ci.org/simpeg/simpeg.png?branch=master
+      :target: https://travis-ci.org/simpeg/simpeg
+      :alt: Master Branch
+      :align: center
 
 * Develop Branch
-.. image:: https://travis-ci.org/simpeg/simpeg.png?branch=develop
-   :target: https://travis-ci.org/simpeg/simpeg
-   :alt: Develop Branch
-   :align: center
+   .. image:: https://travis-ci.org/simpeg/simpeg.png?branch=develop
+      :target: https://travis-ci.org/simpeg/simpeg
+      :alt: Develop Branch
+      :align: center
 
 
 Utility Codes

notebooks/BFGS.ipynb

@@ -1,158 +0,0 @@
-{
- "metadata": {
-  "name": ""
- },
- "nbformat": 3,
- "nbformat_minor": 0,
- "worksheets": [
-  {
-   "cells": [
-    {
-     "cell_type": "code",
-     "collapsed": false,
-     "input": [
-      "import SimPEG\n",
-      "from SimPEG import Solver\n",
-      "from SimPEG.mesh import TensorMesh\n",
-      "from SimPEG.regularization import Regularization\n",
-      "import SimPEG.inverse as inverse\n",
-      "from SimPEG.inverse import Minimize, Remember, IterationPrinters\n",
-      "import numpy as np\n",
-      "import scipy.sparse as sp"
-     ],
-     "language": "python",
-     "metadata": {},
-     "outputs": [],
-     "prompt_number": 2
-    },
-    {
-     "cell_type": "code",
-     "collapsed": false,
-     "input": [
-      "FUN = SimPEG.tests.Rosenbrock\n",
-      "FUN = SimPEG.tests.getQuadratic(sp.csr_matrix(([100,1],([0,1],[0,1])),shape=(2,2)),np.array([-5,-5]),100)\n",
-      "\n",
-      "x0 = np.array([1,0])\n",
-      "opt = inverse.BFGS()\n",
-      "xopt = opt.minimize(FUN,x0)\n",
-      "print xopt\n",
-      "opt = inverse.GaussNewton()\n",
-      "xopt = opt.minimize(FUN,x0)\n",
-      "print xopt\n",
-      "opt = inverse.SteepestDescent()\n",
-      "xopt = opt.minimize(FUN,x0)\n",
-      "print xopt"
-     ],
-     "language": "python",
-     "metadata": {},
-     "outputs": [
-      {
-       "output_type": "stream",
-       "stream": "stdout",
-       "text": [
-        "===================== BFGS =====================\n",
-        "  #      f      |proj(x-g)-x|  LS    Comment   \n",
-        "-----------------------------------------------\n",
-        "   0  1.45e+02    9.51e+01      0              \n",
-        "   1  1.14e+02    5.37e+01      6              \n",
-        "   2  1.04e+02    3.04e+01      6              \n",
-        "   3  8.83e+01    1.37e+01      0              \n",
-        "   4  8.76e+01    5.97e+00      0   Skip BFGS  \n",
-        "   5  8.74e+01    2.61e+00      0   Skip BFGS  \n",
-        "   6  8.74e+01    1.14e+00      0   Skip BFGS  \n",
-        "   7  8.74e+01    5.01e-01      0   Skip BFGS  \n",
-        "   8  8.74e+01    2.19e-01      0   Skip BFGS  \n",
-        "   9  8.74e+01    9.60e-02      0   Skip BFGS  \n",
-        "------------------------- STOP! -------------------------\n",
-        "1 : |fc-fOld| = 1.9437e-04 <= tolF*(1+|f0|) = 1.4600e+01\n",
-        "1 : |xc-x_last| = 1.2663e-03 <= tolX*(1+|x0|) = 2.0000e-01\n",
-        "1 : |proj(x-g)-x|    = 9.5952e-02 <= tolG          = 1.0000e-01\n",
-        "0 : |proj(x-g)-x|    = 9.5952e-02 <= 1e3*eps       = 1.0000e-02\n",
-        "0 : maxIter   =      20    <= iter          =      9\n",
-        "------------------------- DONE! -------------------------\n",
-        "[ 0.05095952  4.99977449]\n",
-        "=========== Gauss Newton ===========\n",
-        "  #      f      |proj(x-g)-x|  LS  \n",
-        "-----------------------------------\n",
-        "   0  1.45e+02    9.51e+01      0  \n",
-        "   1  8.74e+01    4.44e-15      0  \n",
-        "------------------------- STOP! -------------------------\n",
-        "0 : |fc-fOld| = 5.7625e+01 <= tolF*(1+|f0|) = 1.4600e+01\n",
-        "0 : |xc-x_last| = 5.0894e+00 <= tolX*(1+|x0|) = 2.0000e-01\n",
-        "1 : |proj(x-g)-x|    = 4.4409e-15 <= tolG          = 1.0000e-01\n",
-        "1 : |proj(x-g)-x|    = 4.4409e-15 <= 1e3*eps       = 1.0000e-02\n",
-        "0 : maxIter   =      20    <= iter          =      1\n",
-        "------------------------- DONE! -------------------------\n",
-        "[ 0.05  5.  ]\n",
-        "========= Steepest Descent =========\n",
-        "  #      f      |proj(x-g)-x|  LS  \n",
-        "-----------------------------------\n",
-        "   0  1.45e+02    9.51e+01      0  \n",
-        "   1  1.14e+02    5.37e+01      6  \n",
-        "   2  1.04e+02    3.04e+01      6  \n",
-        "   3  1.00e+02    1.76e+01      6  \n",
-        "   4  9.88e+01    1.06e+01      6  \n",
-        "   5  9.82e+01    7.07e+00      6  \n",
-        "   6  9.80e+01    1.22e+01      5  \n",
-        "   7  9.73e+01    7.77e+00      6  \n",
-        "   8  9.68e+01    5.64e+00      6  \n",
-        "   9  9.65e+01    8.72e+00      5  \n",
-        "  10  9.60e+01    5.97e+00      6  \n",
-        "  11  9.58e+01    9.98e+00      5  \n",
-        "  12  9.53e+01    6.48e+00      6  \n",
-        "  13  9.53e+01    1.16e+01      5  \n",
-        "  14  9.46e+01    7.20e+00      6  \n",
-        "  15  9.43e+01    5.07e+00      6  \n",
-        "  16  9.41e+01    8.17e+00      5  \n",
-        "  17  9.37e+01    5.43e+00      6  \n",
-        "  18  9.36e+01    9.42e+00      5  \n",
-        "  19  9.32e+01    5.98e+00      6  \n",
-        "  20  9.29e+01    4.32e+00      6  \n",
-        "------------------------- STOP! -------------------------\n",
-        "1 : |fc-fOld| = 2.5913e-01 <= tolF*(1+|f0|) = 1.4600e+01\n",
-        "1 : |xc-x_last| = 9.3379e-02 <= tolX*(1+|x0|) = 2.0000e-01\n",
-        "0 : |proj(x-g)-x|    = 4.3246e+00 <= tolG          = 1.0000e-01\n",
-        "0 : |proj(x-g)-x|    = 4.3246e+00 <= 1e3*eps       = 1.0000e-02\n",
-        "1 : maxIter   =      20    <= iter          =     20\n",
-        "------------------------- DONE! -------------------------\n",
-        "[ 0.07777107  1.6849632 ]\n"
-       ]
-      },
-      {
-       "output_type": "stream",
-       "stream": "stderr",
-       "text": [
-        "/Users/rowan/git/simpeg/SimPEG/inverse/Optimize.py:664: RuntimeWarning: divide by zero encountered in remainder\n",
-        "  khat    = np.mod(n-nn+k,nn)\n"
-       ]
-      }
-     ],
-     "prompt_number": 3
-    },
-    {
-     "cell_type": "code",
-     "collapsed": false,
-     "input": [
-      "A = sp.identity(2)\n",
-      "S = Solver(A)\n",
-      "\n",
-      "assert type(S) is Solver"
-     ],
-     "language": "python",
-     "metadata": {},
-     "outputs": [],
-     "prompt_number": 6
-    },
-    {
-     "cell_type": "code",
-     "collapsed": false,
-     "input": [],
-     "language": "python",
-     "metadata": {},
-     "outputs": []
-    }
-   ],
-   "metadata": {}
-  }
- ]
-}

notebooks/VisualizeWithvtkView-updated.ipynb

@@ -1,142 +0,0 @@
-{
- "metadata": {
-  "name": "VisualizeWithvtkView-updated"
- },
- "nbformat": 3,
- "nbformat_minor": 0,
- "worksheets": [
-  {
-   "cells": [
-    {
-     "cell_type": "code",
-     "collapsed": false,
-     "input": [
-      "import SimPEG as simpeg, matplotlib as mpl"
-     ],
-     "language": "python",
-     "metadata": {},
-     "outputs": [
-      {
-       "output_type": "stream",
-       "stream": "stdout",
-       "text": [
-        "The history saving thread hit an unexpected error (OperationalError('disk I/O error',)).History will not be written to the database.\n",
-        "Warning: mumps solver not available."
-       ]
-      },
-      {
-       "output_type": "stream",
-       "stream": "stdout",
-       "text": [
-        "\n"
-       ]
-      }
-     ],
-     "prompt_number": 1
-    },
-    {
-     "cell_type": "markdown",
-     "metadata": {},
-     "source": [
-      "Simple notebook of how to use vtkView to visualize SimPEG models. It will pop-up external vtk windows."
-     ]
-    },
-    {
-     "cell_type": "code",
-     "collapsed": false,
-     "input": [
-      "# Make a mesh and model\n",
-      "x0 = np.zeros(3)\n",
-      "h1 = np.ones(60)*50\n",
-      "h2 = np.ones(60)*100\n",
-      "h3 = np.ones(50)*200\n",
-      "\n",
-      "mesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n",
-      "\n",
-      "# Make a models that correspond to the cells, faces and edges.\n",
-      "t = np.ones(mesh.nC)\n",
-      "t[10000:50000] = 100\n",
-      "t[100000:120000] = 100\n",
-      "t[100000:120000] = 50\n",
-      "# Make models called 'Test' for all with a range. \n",
-      "models = {'C':{'Test':np.arange(0,mesh.nC),'Model':t},'F':{'Test':np.arange(0,mesh.nF)},'E':{'Test':np.arange(0,mesh.nE)}}\n"
-     ],
-     "language": "python",
-     "metadata": {},
-     "outputs": [],
-     "prompt_number": 2
-    },
-    {
-     "cell_type": "code",
-     "collapsed": false,
-     "input": [
-      "# Make the vtk viewer object.\n",
-      "vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models)\n",
-      "# Set the .viewprop for which model to view\n",
-      "vtkViewer.viewprop = {'F':'Test'}\n",
-      "# Show the image\n",
-      "vtkViewer.Show()\n",
-      "\n"
-     ],
-     "language": "python",
-     "metadata": {},
-     "outputs": [],
-     "prompt_number": 3
-    },
-    {
-     "cell_type": "code",
-     "collapsed": false,
-     "input": [
-      "# Set subset of the mesh to view (remove padding)\n",
-      "vtkViewer.extent = [4,14,0,7,0,3]\n",
-      "vtkViewer.Show()\n",
-      "\n",
-      "# Change viewing property \n",
-      "vtkViewer.viewprop = {'C':'Model'}\n",
-      "# Set the color range\n",
-      "# Reset extent. Error check will reset the limits correctly.\n",
-      "vtkViewer.extent = [-1,1000,-1,1000,-1,1000]\n",
-      "# Set the range\n",
-      "vtkViewer.range = [0.,100.]\n",
-      "# Show\n",
-      "vtkViewer.Show()\n"
-     ],
-     "language": "python",
-     "metadata": {},
-     "outputs": [
-      {
-       "output_type": "stream",
-       "stream": "stderr",
-       "text": [
-        "/home/Gudni/Codes/python/simpeg/SimPEG/visualize/vtk/vtkView.py:116: UserWarning: Lower bounds smaller then 0\n",
-        "  warnings.warn('Lower bounds smaller then 0')\n",
-        "/home/Gudni/Codes/python/simpeg/SimPEG/visualize/vtk/vtkView.py:128: UserWarning: Upper bounds greater then number of cells\n",
-        "  warnings.warn('Upper bounds greater then number of cells')\n",
-        "/home/Gudni/Codes/python/simpeg/SimPEG/visualize/vtk/vtkView.py:137: UserWarning: Changed given extent from [-1, 1000, -1, 1000, -1, 1000] to [0, 59, 0, 59, 0, 49]\n",
-        "  warnings.warn('Changed given extent from {:s} to {:s}'.format(value,valnp.tolist()))\n"
-       ]
-      }
-     ],
-     "prompt_number": 4
-    },
-    {
-     "cell_type": "code",
-     "collapsed": false,
-     "input": [
-      "# Change color scale, has to be set to bytes=True.\n",
-      "vtkViewer.cmap = mpl.cm.copper(np.arange(0.,1.,0.01),bytes=True)\n",
-      "vtkViewer.Show()\n",
-      "# Set limits of values to view \n",
-      "vtkViewer.limits = [5.0,100.0]\n",
-      "vtkViewer.Show()"
-     ],
-     "language": "python",
-     "metadata": {},
-     "outputs": [],
-     "prompt_number": 5
-    }
-   ],
-   "metadata": {}
-  }
- ]
-}