Merge pull request #42 from simpeg/hdf5

Hdf5 saving ability.

README.md

@@ -1,6 +1,6 @@
 ![SimPEG](https://raw.github.com/simpeg/simpeg/master/docs/simpeg-logo.png)
 
-Simulation and Parameter Estimation in Geophysics  -  A python package for simulation and gradient based parameter estimation in the context of geophysical applications.
+Simulation and Parameter Estimation in Geoscience  -  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:
 

SimPEG/__init__.py

@@ -1,10 +1,12 @@
+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 forward
-import regularization
 import examples
 
 import scipy.version as _v

SimPEG/data/__init__.py

@@ -0,0 +1,20 @@
+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/examples/DC.py

@@ -168,7 +168,6 @@ def genTxRxmat(nelec, spacelec, surfloc, elecini, mesh):
 
 if __name__ == '__main__':
 
-    from SimPEG.regularization import Regularization
     from SimPEG import inverse
     import matplotlib.pyplot as plt
 
@@ -217,7 +216,7 @@ if __name__ == '__main__':
     m0 = mesh.gridCC[:,0]*0+sig2
 
     opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
-    reg = Regularization(mesh)
+    reg = inverse.Regularization(mesh)
     inv = inverse.Inversion(problem, reg, opt, beta0=1e4)
 
     # Check Derivative

SimPEG/examples/Linear.py

@@ -1,14 +1,13 @@
-import numpy as np
-from SimPEG.mesh import TensorMesh
-from SimPEG.forward import Problem
-from SimPEG.regularization import Regularization
-from SimPEG.inverse import *
+from SimPEG import mesh, forward, inverse, np
 import matplotlib.pyplot as plt
 
 
-class LinearProblem(Problem):
+class LinearProblem(forward.Problem):
     """docstring for LinearProblem"""
 
+    def __init__(self, *args, **kwargs):
+        forward.Problem.__init__(self, *args, **kwargs)
+
     def dpred(self, m, u=None):
         return self.G.dot(m)
 
@@ -21,7 +20,7 @@ class LinearProblem(Problem):
 
 def example(N):
     h = np.ones(N)/N
-    M = TensorMesh([h])
+    M = mesh.TensorMesh([h])
 
     nk = 20
     jk = np.linspace(1.,20.,nk)
@@ -35,35 +34,27 @@ def example(N):
     for i in range(nk):
         G[i,:] = g(i)
 
-
-    m_true = np.zeros(M.nC)
-    m_true[M.vectorCCx > 0.3] = 1.
-    m_true[M.vectorCCx > 0.45] = -0.5
-    m_true[M.vectorCCx > 0.6] = 0
-
-
-    d_true = G.dot(m_true)
-    noise = 0.1 * np.random.rand(d_true.size)
-
-    d_obs = d_true + noise
+    mtrue = np.zeros(M.nC)
+    mtrue[M.vectorCCx > 0.3] = 1.
+    mtrue[M.vectorCCx > 0.45] = -0.5
+    mtrue[M.vectorCCx > 0.6] = 0
 
     prob = LinearProblem(M)
     prob.G = G
-    prob.dobs = d_obs
-    prob.std = np.ones_like(d_obs)*0.1
+    data = prob.createSyntheticData(mtrue, std=0.01)
 
-    return prob, m_true
+    return prob, data
 
 
 if __name__ == '__main__':
 
-    prob, m_true = example(100)
+    prob, data = example(100)
     M = prob.mesh
 
-    reg = Regularization(M)
-    opt = InexactGaussNewton(maxIter=20)
-    inv = Inversion(prob,reg,opt,beta0=1e-4)
-    m0 = np.zeros_like(m_true)
+    reg = inverse.Regularization(M)
+    opt = inverse.InexactGaussNewton(maxIter=20)
+    inv = inverse.Inversion(prob,reg,opt,data)
+    m0 = np.zeros_like(data.mtrue)
 
     mrec = inv.run(m0)
 
@@ -73,9 +64,7 @@ if __name__ == '__main__':
 
     plt.figure(2)
 
-    plt.plot(M.vectorCCx, m_true, 'b-')
+    plt.plot(M.vectorCCx, data.mtrue, 'b-')
     plt.plot(M.vectorCCx, mrec, 'r-')
 
-
-
     plt.show()

SimPEG/examples/__init__.py

@@ -0,0 +1,2 @@
+import DC
+import Linear

SimPEG/forward/Problem.py

@@ -1,6 +1,4 @@
-import numpy as np
-from SimPEG.utils import mkvc, sdiag, count, timeIt
-import scipy.sparse as sp
+from SimPEG import utils, data, np, sp
 norm = np.linalg.norm
 
 
@@ -37,9 +35,13 @@ class Problem(object):
         to (locally) find how model parameters change the data, and optimize!
     """
 
-    counter = None
+    __metaclass__ = utils.Save.Savable
 
-    def __init__(self, mesh):
+    counter = None   #: A SimPEG.utils.Counter object
+
+
+    def __init__(self, mesh, *args, **kwargs):
+        utils.setKwargs(self, **kwargs)
         self.mesh = mesh
 
     @property
@@ -65,27 +67,7 @@ class Problem(object):
     def P(self, value):
         self._P = value
 
-    @property
-    def std(self):
-        """
-            Estimated Standard Deviations.
-        """
-        return self._std
-    @std.setter
-    def std(self, value):
-        self._std = value
-
-    @property
-    def dobs(self):
-        """
-            Observed data.
-        """
-        return self._dobs
-    @dobs.setter
-    def dobs(self, value):
-        self._dobs = value
-
-    @count
+    @utils.count
     def dpred(self, m, u=None):
         """
             Predicted data.
@@ -97,8 +79,8 @@ class Problem(object):
             u = self.field(m)
         return self.P*u
 
-    @count
-    def dataResidual(self, m, u=None):
+    @utils.count
+    def dataResidual(self, m, data, u=None):
         """
             :param numpy.array m: geophysical model
             :param numpy.array u: fields
@@ -115,9 +97,9 @@ class Problem(object):
             u is the field of interest; d_obs is the observed data.
         """
 
-        return self.dpred(m, u=u) - self.dobs
+        return self.dpred(m, u=u) - data.dobs
 
-    @timeIt
+    @utils.timeIt
     def J(self, m, v, u=None):
         """
             :param numpy.array m: model
@@ -147,7 +129,7 @@ class Problem(object):
         """
         raise NotImplementedError('J is not yet implemented.')
 
-    @timeIt
+    @utils.timeIt
     def Jt(self, m, v, u=None):
         """
             :param numpy.array m: model
@@ -161,7 +143,7 @@ class Problem(object):
         raise NotImplementedError('Jt is not yet implemented.')
 
 
-    @timeIt
+    @utils.timeIt
     def J_approx(self, m, v, u=None):
         """
 
@@ -176,7 +158,7 @@ class Problem(object):
         """
         return self.J(m, v, u)
 
-    @timeIt
+    @utils.timeIt
     def Jt_approx(self, m, v, u=None):
         """
             :param numpy.array m: model
@@ -238,12 +220,11 @@ 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.
         """
-        dobs = self.dpred(m,u=u)
-        noise = std*abs(dobs)*np.random.randn(*dobs.shape)
-        dobs = dobs+noise
-        eps = np.linalg.norm(mkvc(dobs),2)*1e-5
-        Wd = 1/(abs(dobs)*std+eps)
-        return dobs, Wd
+        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)
 
 
 

SimPEG/forward/__init__.py

@@ -1,4 +1,2 @@
 from Problem import *
-import DCProblem
-from LinearProblem import LinearProblem
 import ModelTransforms

SimPEG/inverse/Inversion.py

@@ -1,13 +1,14 @@
-import numpy as np
-import scipy.sparse as sp
 import SimPEG
-from SimPEG.utils import sdiag, mkvc, setKwargs, checkStoppers, printStoppers, count, timeIt, callHooks
+from SimPEG import utils, sp, np
 from Optimize import Remember
 from BetaSchedule import Cooling
 from SimPEG.inverse import IterationPrinters, StoppingCriteria
 
 class BaseInversion(object):
-    """docstring for BaseInversion"""
+    """BaseInversion(prob, reg, opt, data, **kwargs)
+    """
+
+    __metaclass__ = utils.Save.Savable
 
     maxIter = 1        #: Maximum number of iterations
     name = 'BaseInversion'
@@ -17,14 +18,15 @@ class BaseInversion(object):
     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  = 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, **kwargs):
-        setKwargs(self, **kwargs)
+    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]
@@ -46,8 +48,8 @@ class BaseInversion(object):
             Standard deviation weighting matrix.
         """
         if getattr(self,'_Wd',None) is None:
-            eps = np.linalg.norm(mkvc(self.prob.dobs),2)*1e-5
-            self._Wd = 1/(abs(self.prob.dobs)*self.prob.std+eps)
+            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):
@@ -63,14 +65,14 @@ class BaseInversion(object):
         Note that we do not set the target if it is None, but we return the default value.
         """
         if getattr(self, '_phi_d_target', None) is None:
-            return self.prob.dobs.size #
+            return self.data.dobs.size #
         return self._phi_d_target
 
     @phi_d_target.setter
     def phi_d_target(self, value):
         self._phi_d_target = value
 
-    @timeIt
+    @utils.timeIt
     def run(self, m0):
         """run(m0)
 
@@ -85,25 +87,19 @@ class BaseInversion(object):
             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.
 
-            If you have things that also need to run on startup, you can create a method::
-
-                def _startup*(self, x0):
-                    pass
-
-            Where the * can be any string. If present, _startup* will be called at the start of the default startup call.
-            You may also completely overwrite this function.
-
             :param numpy.ndarray x0: initial x
             :rtype: None
             :return: None
         """
-        callHooks(self,'startup',m0)
 
         if not hasattr(self.reg, '_mref'):
             print 'Regularization has not set mref. SimPEG will set it to m0.'
@@ -115,43 +111,25 @@ class BaseInversion(object):
         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.
 
-            If you have things that also need to run at the end of every iteration, you can create a method::
-
-                def _doStartIteration*(self):
-                    pass
-
-            Where the * can be any string. If present, _doStartIteration* will be called at the start of the default doStartIteration call.
-            You may also completely overwrite this function.
-
             :rtype: None
             :return: None
         """
-        callHooks(self,'doStartIteration')
-
         self._beta = self.getBeta()
 
 
+    @utils.callHooks('doEndIteration')
     def doEndIteration(self):
         """
             **doEndIteration** is called at the end of each run iteration.
 
-            If you have things that also need to run at the end of every iteration, you can create a method::
-
-                def _doEndIteration*(self):
-                    pass
-
-            Where the * can be any string. If present, _doEndIteration* will be called at the start of the default doEndIteration call.
-            You may also completely overwrite this function.
-
             :rtype: None
             :return: None
         """
-        callHooks(self,'doEndIteration')
-
         # store old values
         self.phi_d_last = self.phi_d
         self.phi_m_last = self.phi_m
@@ -203,7 +181,7 @@ class BaseInversion(object):
 
     def stoppingCriteria(self):
         if self.debug: print 'checking stoppingCriteria'
-        return checkStoppers(self, self.stoppers)
+        return utils.checkStoppers(self, self.stoppers)
 
 
     def printDone(self):
@@ -211,9 +189,17 @@ class BaseInversion(object):
             **printDone** is called at the end of the inversion routine.
 
         """
-        printStoppers(self, self.stoppers)
+        utils.printStoppers(self, self.stoppers)
 
-    @timeIt
+    @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)
 
@@ -223,7 +209,7 @@ class BaseInversion(object):
         u = self.prob.field(m)
 
         if self._iter is 0 and self._beta is None:
-            self._beta = self.beta0 = self.estimateBeta0(u=u)
+            self._beta = self.beta0 = self.estimateBeta0(u=u,ratio=self.beta0_ratio)
 
         phi_d = self.dataObj(m, u)
         phi_m = self.reg.modelObj(m)
@@ -253,7 +239,7 @@ class BaseInversion(object):
             out += (operator,)
         return out if len(out) > 1 else out[0]
 
-    @timeIt
+    @utils.timeIt
     def dataObj(self, m, u=None):
         """dataObj(m, u=None)
 
@@ -272,11 +258,11 @@ class BaseInversion(object):
             u is the field of interest; d_obs is the observed data; and W is the weighting matrix.
         """
         # TODO: ensure that this is a data is vector and Wd is a matrix.
-        R = self.Wd*self.prob.dataResidual(m, u=u)
-        R = mkvc(R)
+        R = self.Wd*self.prob.dataResidual(m, self.data, u=u)
+        R = utils.mkvc(R)
         return 0.5*np.vdot(R, R)
 
-    @timeIt
+    @utils.timeIt
     def dataObjDeriv(self, m, u=None):
         """dataObjDeriv(m, u=None)
 
@@ -312,13 +298,13 @@ class BaseInversion(object):
         if u is None:
             u = self.prob.field(m)
 
-        R = self.Wd*self.prob.dataResidual(m, u=u)
+        R = self.Wd*self.prob.dataResidual(m, self.data, u=u)
 
         dmisfit = self.prob.Jt(m, self.Wd * R, u=u)
 
         return dmisfit
 
-    @timeIt
+    @utils.timeIt
     def dataObj2Deriv(self, m, v, u=None):
         """dataObj2Deriv(m, v, u=None)
 
@@ -357,7 +343,7 @@ class BaseInversion(object):
         if u is None:
             u = self.prob.field(m)
 
-        R = self.Wd*self.prob.dataResidual(m, u=u)
+        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.
@@ -365,13 +351,19 @@ class BaseInversion(object):
 
         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, **kwargs):
-        BaseInversion.__init__(self, prob, reg, opt, **kwargs)
+    def __init__(self, prob, reg, opt, data, **kwargs):
+        BaseInversion.__init__(self, prob, reg, opt, data, **kwargs)
 
         self.stoppers.append(StoppingCriteria.phi_d_target_Inversion)
 
@@ -387,8 +379,8 @@ class TimeSteppingInversion(Remember, BaseInversion):
     maxIter = 1
     name = "Time-Stepping SimPEG Inversion"
 
-    def __init__(self, prob, reg, opt, **kwargs):
-        BaseInversion.__init__(self, prob, reg, opt, **kwargs)
+    def __init__(self, prob, reg, opt, data, **kwargs):
+        BaseInversion.__init__(self, prob, reg, opt, data, **kwargs)
 
         self.stoppers.append(StoppingCriteria.phi_d_target_Inversion)
 

SimPEG/inverse/Optimize.py

@@ -1,9 +1,6 @@
-import numpy as np
+from SimPEG import Solver, utils, sp, np
 import matplotlib.pyplot as plt
-from SimPEG.utils import mkvc, sdiag, setKwargs, printTitles, printLine, printStoppers, checkStoppers, count, timeIt, callHooks
 norm = np.linalg.norm
-import scipy.sparse as sp
-from SimPEG import Solver
 
 
 __all__ = ['Minimize', 'Remember', 'SteepestDescent', 'BFGS', 'GaussNewton', 'InexactGaussNewton', 'ProjectedGradient', 'NewtonRoot', 'StoppingCriteria', 'IterationPrinters']
@@ -85,6 +82,8 @@ class Minimize(object):
         Minimize is a general class for derivative based optimization.
     """
 
+    __metaclass__ = utils.Save.Savable
+
     name = "General Optimization Algorithm"  #: The name of the optimization algorithm
 
     maxIter = 20       #: Maximum number of iterations
@@ -110,9 +109,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]
 
-        setKwargs(self, **kwargs)
+        utils.setKwargs(self, **kwargs)
 
-    @timeIt
+    @utils.timeIt
     def minimize(self, evalFunction, x0):
         """minimize(evalFunction, x0)
 
@@ -155,6 +154,7 @@ class Minimize(object):
                 doEndIteration(xt)
 
             printDone()
+            finish()
             return xc
         """
         self.evalFunction = evalFunction
@@ -175,6 +175,7 @@ class Minimize(object):
             self.doEndIteration(xt)
 
         self.printDone()
+        self.finish()
 
         return self.xc
 
@@ -188,6 +189,7 @@ class Minimize(object):
     def parent(self, value):
         self._parent = value
 
+    @utils.callHooks('startup')
     def startup(self, x0):
         """
             **startup** is called at the start of any new minimize call.
@@ -198,19 +200,10 @@ class Minimize(object):
                 xc = x0
                 _iter = _iterLS = 0
 
-            If you have things that also need to run on startup, you can create a method::
-
-                def _startup*(self, x0):
-                    pass
-
-            Where the * can be any string. If present, _startup* will be called at the start of the default startup call.
-            You may also completely overwrite this function.
-
             :param numpy.ndarray x0: initial x
             :rtype: None
             :return: None
         """
-        callHooks(self,'startup',x0)
 
         self._iter = 0
         self._iterLS = 0
@@ -221,7 +214,8 @@ class Minimize(object):
         self.f_last = np.nan
         self.x_last = x0
 
-    @count
+    @utils.count
+    @utils.callHooks('doStartIteration')
     def doStartIteration(self):
         """doStartIteration()
 
@@ -230,7 +224,8 @@ class Minimize(object):
             :rtype: None
             :return: None
         """
-        callHooks(self,'doStartIteration')
+        pass
+
 
     def printInit(self, inLS=False):
         """
@@ -242,8 +237,9 @@ class Minimize(object):
         """
         pad = ' '*10 if inLS else ''
         name = self.name if not inLS else self.nameLS
-        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')
     def printIter(self, inLS=False):
         """
             **printIter** is called directly after function evaluations.
@@ -252,10 +248,8 @@ class Minimize(object):
             parent.printIter function and call that.
 
         """
-        callHooks(self,'printIter',inLS)
-
         pad = ' '*10 if inLS else ''
-        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):
         """
@@ -268,15 +262,29 @@ 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
-        printStoppers(self, stoppers, pad='', stop=stop, done=done)
+        utils.printStoppers(self, stoppers, pad='', stop=stop, done=done)
+
+
+    @utils.callHooks('finish')
+    def finish(self):
+        """finish()
+
+            **finish** is called at the end of the optimization.
+
+            :rtype: None
+            :return: None
+
+        """
+        pass
 
     def stoppingCriteria(self, inLS=False):
         if self._iter == 0:
             self.f0 = self.f
             self.g0 = self.g
-        return checkStoppers(self, self.stoppers if not inLS else self.stoppersLS)
+        return utils.checkStoppers(self, self.stoppers if not inLS else self.stoppersLS)
 
-    @timeIt
+    @utils.timeIt
+    @utils.callHooks('projection')
     def projection(self, p):
         """projection(p)
 
@@ -288,10 +296,9 @@ class Minimize(object):
             :rtype: numpy.ndarray
             :return: p, projected search direction
         """
-        callHooks(self,'projection',p)
         return p
 
-    @timeIt
+    @utils.timeIt
     def findSearchDirection(self):
         """findSearchDirection()
 
@@ -322,7 +329,7 @@ class Minimize(object):
         """
         return -self.g
 
-    @count
+    @utils.count
     def scaleSearchDirection(self, p):
         """scaleSearchDirection(p)
 
@@ -341,7 +348,7 @@ class Minimize(object):
 
     nameLS = "Armijo linesearch" #: The line-search name
 
-    @timeIt
+    @utils.timeIt
     def modifySearchDirection(self, p):
         """modifySearchDirection(p)
 
@@ -379,7 +386,7 @@ class Minimize(object):
 
         return self._LS_xt, self._iterLS < self.maxIterLS
 
-    @count
+    @utils.count
     def modifySearchDirectionBreak(self, p):
         """modifySearchDirectionBreak(p)
 
@@ -401,7 +408,8 @@ class Minimize(object):
         print 'The linesearch got broken. Boo.'
         return p, False
 
-    @count
+    @utils.count
+    @utils.callHooks('doEndIteration')
     def doEndIteration(self, xt):
         """doEndIteration(xt)
 
@@ -411,21 +419,10 @@ class Minimize(object):
 
             self.xc must be updated in this code.
 
-
-            If you have things that also need to run at the end of every iteration, you can create a method::
-
-                def _doEndIteration*(self, xt):
-                    pass
-
-            Where the * can be any string. If present, _doEndIteration* will be called at the start of the default doEndIteration call.
-            You may also completely overwrite this function.
-
             :param numpy.ndarray xt: tested new iterate that ensures a descent direction.
             :rtype: None
             :return: None
         """
-        callHooks(self,'doEndIteration',xt)
-
         # store old values
         self.f_last = self.f
         self.x_last, self.xc = self.xc, xt
@@ -433,6 +430,19 @@ class Minimize(object):
         if self.debug: self.printDone()
 
 
+    def save(self, group):
+        group.setArray('searchDirection', self.searchDirection)
+
+        if getattr(self,'parent',None) is None:
+            group.setArray('x', self.xc)
+        else: # Assume inversion is the parent
+            group.attrs['phi_d'] = self.parent.phi_d
+            group.attrs['phi_m'] = self.parent.phi_m
+            group.setArray('m', self.xc)
+            group.setArray('dpred', self.parent.dpred)
+
+
+
 class Remember(object):
     """
         This mixin remembers all the things you tend to forget.
@@ -517,7 +527,7 @@ class ProjectedGradient(Minimize, Remember):
 
         self.aSet_prev = self.activeSet(x0)
 
-    @count
+    @utils.count
     def projection(self, x):
         """projection(x)
 
@@ -526,7 +536,7 @@ class ProjectedGradient(Minimize, Remember):
         """
         return np.median(np.c_[self.lower,x,self.upper],axis=1)
 
-    @count
+    @utils.count
     def activeSet(self, x):
         """activeSet(x)
 
@@ -535,7 +545,7 @@ class ProjectedGradient(Minimize, Remember):
         """
         return np.logical_or(x == self.lower, x == self.upper)
 
-    @count
+    @utils.count
     def inactiveSet(self, x):
         """inactiveSet(x)
 
@@ -544,7 +554,7 @@ class ProjectedGradient(Minimize, Remember):
         """
         return np.logical_not(self.activeSet(x))
 
-    @count
+    @utils.count
     def bindingSet(self, x):
         """bindingSet(x)
 
@@ -557,7 +567,7 @@ class ProjectedGradient(Minimize, Remember):
         bind_low = np.logical_and(x == self.upper, self.g <= 0)
         return np.logical_or(bind_up, bind_low)
 
-    @timeIt
+    @utils.timeIt
     def findSearchDirection(self):
         """findSearchDirection()
 
@@ -602,7 +612,7 @@ class ProjectedGradient(Minimize, Remember):
             # aSet_after = self.activeSet(self.xc+p)
         return p
 
-    @timeIt
+    @utils.timeIt
     def _doEndIteration_ProjectedGradient(self, xt):
         """_doEndIteration_ProjectedGradient(xt)"""
         aSet = self.activeSet(xt)
@@ -630,13 +640,16 @@ class ProjectedGradient(Minimize, Remember):
 
         if self.debug: print 'doEndIteration.ProjGrad, f_current_decrease: ', f_current_decrease
         if self.debug: print 'doEndIteration.ProjGrad, f_decrease_max: ', self.f_decrease_max
-        if self.debug: print 'doEndIteration.ProjGrad, stopDoingSD: ', self.stopDoingSD
+        if self.debug: print 'doEndIteration.ProjGrad, stopDoingSD: ', self.stopDoingPG
 
 
 class BFGS(Minimize, Remember):
     name  = 'BFGS'
     nbfgs = 10
 
+    def __init__(self, **kwargs):
+        Minimize.__init__(self, **kwargs)
+
     @property
     def bfgsH0(self):
         """
@@ -702,7 +715,10 @@ class BFGS(Minimize, Remember):
 class GaussNewton(Minimize, Remember):
     name = 'Gauss Newton'
 
-    @timeIt
+    def __init__(self, **kwargs):
+        Minimize.__init__(self, **kwargs)
+
+    @utils.timeIt
     def findSearchDirection(self):
         return Solver(self.H).solve(-self.g)
 
@@ -749,7 +765,7 @@ class InexactGaussNewton(BFGS, Minimize, Remember):
     def approxHinv(self, value):
         self._approxHinv = value
 
-    @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)
@@ -759,7 +775,10 @@ class InexactGaussNewton(BFGS, Minimize, Remember):
 class SteepestDescent(Minimize, Remember):
     name = 'Steepest Descent'
 
-    @timeIt
+    def __init__(self, **kwargs):
+        Minimize.__init__(self, **kwargs)
+
+    @utils.timeIt
     def findSearchDirection(self):
         return -self.g
 
@@ -792,15 +811,26 @@ class NewtonRoot(object):
     doLS     = True
 
     def __init__(self, **kwargs):
-        setKwargs(self, **kwargs)
+        utils.setKwargs(self, **kwargs)
 
     def root(self, fun, x):
+        """root(fun, x)
+
+        Function Should have the form::
+
+            def evalFunction(x, return_g=False):
+                    out = (f,)
+                    if return_g:
+                        out += (g,)
+                    return out if len(out) > 1 else out[0]
+
+        """
         if self.comments: print 'Newton Method:\n'
 
         self._iter = 0
         while True:
 
-            [r,J] = fun(x);
+            r, J = fun(x, return_g=True)
             if self.solveTol == 0:
                 Jinv = Solver(J)
                 dh   = - Jinv.solve(r)
@@ -812,13 +842,12 @@ class NewtonRoot(object):
             muLS = 1.
             LScnt  = 1
             xt = x + dh
-            rt, Jt = fun(xt) # TODO: get rid of Jt
+            rt = fun(xt, return_g=False)
 
-            if self.comments: print '\tLinesearch:\n'
+            if self.comments and self.doLS: print '\tLinesearch:\n'
             # Enter Linesearch
             while True and self.doLS:
-                if self.comments:
-                    print '\t\tResid: %e\n'%norm(rt)
+                if self.comments: print '\t\tResid: %e\n'%norm(rt)
                 if norm(rt) <= norm(r) or norm(rt) < self.tol:
                     break
 
@@ -827,14 +856,16 @@ class NewtonRoot(object):
                 print '.'
                 if LScnt > self.maxLS:
                     print 'Newton Method: Line search break.'
-                    root = NaN
-                    return
+                    return None
                 xt = x + muLS*dh
-                rt, Jt = fun(xt) # TODO: get rid of Jt
+                rt = fun(xt, return_g=False)
 
             x = xt
             self._iter += 1
-            if norm(rt) < self.tol or self._iter > self.maxIter:
+            if norm(rt) < self.tol:
+                break
+            if self._iter > self.maxIter:
+                print 'NewtonRoot stopped by maxIters. norm: %4.4e' % norm(rt)
                 break
 
         return x
@@ -854,7 +885,7 @@ if __name__ == '__main__':
 
 
     print 'test the newtonRoot finding.'
-    fun = lambda x: (np.sin(x), 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=False).root(fun,x)
+    pnt = NewtonRoot(comments=True).root(fun,x)
     print pnt

SimPEG/inverse/Regularization.py

@@ -0,0 +1,219 @@
+from SimPEG import utils, np, sp
+
+class Regularization(object):
+    """**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:
+
+        .. math::
+
+            R(m) = \int_\Omega \\frac{\\alpha_x}{2}\left(\\frac{\partial m}{\partial x}\\right)^2 + \\frac{\\alpha_y}{2}\left(\\frac{\partial m}{\partial y}\\right)^2 \partial v
+
+        Our discrete gradient operator works on cell centers and gives the derivative on the cell faces, which is not where we want to be evaluating this integral. We need to average the values back to the cell-centers before we integrate. To avoid null spaces, we square first and then average. In 2D with ij notation it looks like this:
+
+        .. math::
+
+            R(m) \\approx \sum_{ij} \left[\\frac{\\alpha_x}{2}\left[\left(\\frac{m_{i+1,j} - m_{i,j}}{h}\\right)^2 + \left(\\frac{m_{i,j} - m_{i-1,j}}{h}\\right)^2\\right]
+            + \\frac{\\alpha_y}{2}\left[\left(\\frac{m_{i,j+1} - m_{i,j}}{h}\\right)^2 + \left(\\frac{m_{i,j} - m_{i,j-1}}{h}\\right)^2\\right]
+            \\right]h^2
+
+        If we let D_1 be the derivative matrix in the x direction
+
+        .. math::
+
+            \mathbf{D}_1 = \mathbf{I}_2\otimes\mathbf{d}_1
+
+        .. math::
+
+            \mathbf{D}_2 = \mathbf{d}_2\otimes\mathbf{I}_1
+
+        Where d_1 is the one dimensional derivative:
+
+        .. math::
+
+            \mathbf{d}_1 = \\frac{1}{h} \left[ \\begin{array}{cccc}
+            -1 & 1 & & \\\\
+             & \ddots & \ddots&\\\\
+             &  & -1 & 1\end{array} \\right]
+
+        .. math::
+
+            R(m) \\approx \mathbf{v}^\\top \left[\\frac{\\alpha_x}{2}\mathbf{A}_1 (\mathbf{D}_1 m) \odot (\mathbf{D}_1 m) + \\frac{\\alpha_y}{2}\mathbf{A}_2 (\mathbf{D}_2 m) \odot (\mathbf{D}_2 m) \\right]
+
+        Recall that this is really a just point wise multiplication, or a diagonal matrix times a vector. When we multiply by something in a diagonal we can interchange and it gives the same results (i.e. it is point wise)
+
+        .. math::
+
+            \mathbf{a\odot b} = \\text{diag}(\mathbf{a})\mathbf{b} = \\text{diag}(\mathbf{b})\mathbf{a} = \mathbf{b\odot a}
+
+        and the transpose also is true (but the sizes have to make sense...):
+
+        .. math::
+
+            \mathbf{a}^\\top\\text{diag}(\mathbf{b}) = \mathbf{b}^\\top\\text{diag}(\mathbf{a})
+
+        So R(m) can simplify to:
+
+        .. math::
+
+            R(m) \\approx  \mathbf{m}^\\top \left[\\frac{\\alpha_x}{2}\mathbf{D}_1^\\top \\text{diag}(\mathbf{A}_1^\\top\mathbf{v}) \mathbf{D}_1 +  \\frac{\\alpha_y}{2}\mathbf{D}_2^\\top \\text{diag}(\mathbf{A}_2^\\top \mathbf{v}) \mathbf{D}_2 \\right] \mathbf{m}
+
+        We will define W_x as:
+
+        .. math::
+
+            \mathbf{W}_x = \sqrt{\\alpha_x}\\text{diag}\left(\sqrt{\mathbf{A}_1^\\top\mathbf{v}}\\right) \mathbf{D}_1
+
+
+        And then W as a tall matrix of all of the different regularization terms:
+
+        .. math::
+
+            \mathbf{W} = \left[ \\begin{array}{c}
+            \mathbf{W}_s\\\\
+            \mathbf{W}_x\\\\
+            \mathbf{W}_y\end{array} \\right]
+
+        Then we can write
+
+        .. math::
+
+            R(m) \\approx \\frac{1}{2}\mathbf{m^\\top W^\\top W m}
+
+
+    """
+
+    __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")
+
+    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
+
+    @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)
+        return self._Ws
+
+    @property
+    def Wx(self):
+        """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
+        return self._Wx
+
+    @property
+    def Wy(self):
+        """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
+        return self._Wy
+
+    @property
+    def Wz(self):
+        """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
+        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
+        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
+        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
+        return self._Wzz
+
+
+    @property
+    def W(self):
+        """Full regularization matrix W"""
+        if getattr(self, '_W', None) is None:
+            wlist = (self.Ws, self.Wx, self.Wxx)
+            if self.mesh.dim > 1:
+                wlist += (self.Wy, self.Wyy)
+            if self.mesh.dim > 2:
+                wlist += (self.Wz, self.Wzz)
+            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/inverse/__init__.py

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

SimPEG/mesh/BaseMesh.py

@@ -1,6 +1,5 @@
 import numpy as np
-from SimPEG.utils import mkvc
-
+from SimPEG import utils
 
 
 class BaseMesh(object):
@@ -12,6 +11,7 @@ class BaseMesh(object):
     :param numpy.array,list x0: Origin of the mesh (dim, )
 
     """
+
     def __init__(self, n, x0=None):
 
         # Check inputs
@@ -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] = mkvc(xi)
+                x_array[:, i] = utils.mkvc(xi)
             x = x_array
 
         assert type(x) == np.ndarray, "x must be a numpy array"
@@ -91,7 +91,7 @@ class BaseMesh(object):
             if format == 'M':
                 return xx.reshape(nn, order='F')
             elif format == 'V':
-                return mkvc(xx)
+                return utils.mkvc(xx)
 
         def switchKernal(xx):
             """Switches over the different options."""
@@ -101,7 +101,7 @@ class BaseMesh(object):
                 return outKernal(xx, nn)
             elif xType in ['F', 'E']:
                 # This will only deal with components of fields, not full 'F' or 'E'
-                xx = 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]
 
@@ -308,7 +308,7 @@ class BaseMesh(object):
         """
         fget = lambda self: np.array([x for x in [self.nNx, self.nNy, self.nNz] if not x is None])
         return locals()
-    nNv = property(**nNv())    
+    nNv = property(**nNv())
 
     def nEx():
         doc = """

SimPEG/mesh/InnerProducts.py

@@ -81,9 +81,9 @@ class InnerProducts(object):
     def getFaceInnerProduct(self, mu=None, returnP=False):
         """Wrapper function,
 
-        :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getEdgeInnerProduct`
+        :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getFaceInnerProduct`
 
-        :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getEdgeInnerProduct2D`
+        :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getFaceInnerProduct2D`
         """
         if self.dim == 2:
             return getFaceInnerProduct2D(self, mu, returnP)
@@ -93,9 +93,9 @@ class InnerProducts(object):
     def getEdgeInnerProduct(self, sigma=None, returnP=False):
         """Wrapper function,
 
-        :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getFaceInnerProduct`
+        :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getEdgeInnerProduct`
 
-        :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getFaceInnerProduct2D`
+        :py:func:`SimPEG.mesh.InnerProducts.InnerProducts.getEdgeInnerProduct2D`
         """
         if self.dim == 2:
             return getEdgeInnerProduct2D(self, sigma, returnP)

SimPEG/mesh/LogicallyOrthogonalMesh.py

@@ -1,15 +1,14 @@
-import numpy as np
+from SimPEG import utils, np
 from BaseMesh import BaseMesh
 from DiffOperators import DiffOperators
 from InnerProducts import InnerProducts
 from LomView import LomView
-from SimPEG.utils import mkvc, ndgrid, volTetra, indexCube, faceInfo
 
 # 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)), mkvc(length2D(x), 2))
-normalize3D = lambda x: x/np.kron(np.ones((1, 3)), 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):
@@ -21,6 +20,9 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
     .. plot:: examples/mesh/plot_LogicallyOrthogonalMesh.py
 
     """
+
+    __metaclass__ = utils.Save.Savable
+
     _meshType = 'LOM'
 
     def __init__(self, nodes):
@@ -38,7 +40,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
         # Save nodes to private variable _gridN as vectors
         self._gridN = np.ones((nodes[0].size, self.dim))
         for i, node_i in enumerate(nodes):
-            self._gridN[:, i] = mkvc(node_i.astype(float))
+            self._gridN[:, i] = utils.mkvc(node_i.astype(float))
 
     def gridCC():
         doc = "Cell-centered grid."
@@ -69,10 +71,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 = [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 = [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()
@@ -86,10 +88,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 = [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 = [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()
@@ -102,7 +104,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 = [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()
@@ -116,10 +118,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 = [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 = [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()
@@ -133,10 +135,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 = [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 = [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()
@@ -149,7 +151,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 = [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()
@@ -192,25 +194,25 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
         def fget(self):
             if(self._vol is None):
                 if self.dim == 2:
-                    A, B, C, D = indexCube('ABCD', self.n+1)
-                    normal, area = faceInfo(np.c_[self.gridN, np.zeros((self.nN, 1))], A, B, C, D)
+                    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)
                     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 = indexCube('ABCDEFGH', self.n+1)
+                    A, B, C, D, E, F, G, H = utils.indexCube('ABCDEFGH', self.n+1)
 
-                    vol1 = (volTetra(self.gridN, A, B, D, E) +  # cutted edge top
-                            volTetra(self.gridN, B, E, F, G) +  # cutted edge top
-                            volTetra(self.gridN, B, D, E, G) +  # middle
-                            volTetra(self.gridN, B, C, D, G) +  # cutted edge bottom
-                            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 = (volTetra(self.gridN, A, F, B, C) +  # cutted edge top
-                            volTetra(self.gridN, A, E, F, H) +  # cutted edge top
-                            volTetra(self.gridN, A, H, F, C) +  # middle
-                            volTetra(self.gridN, C, H, D, A) +  # cutted edge bottom
-                            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
@@ -226,30 +228,30 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
                 # Compute areas of cell faces
                 if(self.dim == 2):
                     xy = self.gridN
-                    A, B = indexCube('AB', self.n+1, np.array([self.nNx, self.nCy]))
+                    A, B = utils.indexCube('AB', self.n+1, np.array([self.nNx, self.nCy]))
                     edge1 = xy[B, :] - xy[A, :]
                     normal1 = np.c_[edge1[:, 1], -edge1[:, 0]]
                     area1 = length2D(edge1)
-                    A, D = indexCube('AD', self.n+1, np.array([self.nCx, self.nNy]))
+                    A, D = utils.indexCube('AD', self.n+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_[mkvc(area1), 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 = indexCube('AEFB', self.n+1, np.array([self.nNx, self.nCy, self.nCz]))
-                    normal1, area1 = faceInfo(self.gridN, A, E, F, B, average=False, normalizeNormals=False)
+                    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, D, H, E = indexCube('ADHE', self.n+1, np.array([self.nCx, self.nNy, self.nCz]))
-                    normal2, area2 = faceInfo(self.gridN, A, D, H, E, 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, B, C, D = indexCube('ABCD', self.n+1, np.array([self.nCx, self.nCy, self.nNz]))
-                    normal3, area3 = faceInfo(self.gridN, A, B, C, D, 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)
 
-                    self._area = np.r_[mkvc(area1), mkvc(area2), 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()
@@ -289,21 +291,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 = indexCube('AD', self.n+1, np.array([self.nCx, self.nNy]))
+                    A, D = utils.indexCube('AD', self.n+1, np.array([self.nCx, self.nNy]))
                     edge1 = xy[D, :] - xy[A, :]
-                    A, B = indexCube('AB', self.n+1, np.array([self.nNx, self.nCy]))
+                    A, B = utils.indexCube('AB', self.n+1, np.array([self.nNx, self.nCy]))
                     edge2 = xy[B, :] - xy[A, :]
-                    self._edge = np.r_[mkvc(length2D(edge1)), 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 = indexCube('AD', self.n+1, np.array([self.nCx, self.nNy, self.nNz]))
+                    A, D = utils.indexCube('AD', self.n+1, np.array([self.nCx, self.nNy, self.nNz]))
                     edge1 = xyz[D, :] - xyz[A, :]
-                    A, B = indexCube('AB', self.n+1, np.array([self.nNx, self.nCy, self.nNz]))
+                    A, B = utils.indexCube('AB', self.n+1, np.array([self.nNx, self.nCy, self.nNz]))
                     edge2 = xyz[B, :] - xyz[A, :]
-                    A, E = indexCube('AE', self.n+1, np.array([self.nNx, self.nNy, self.nCz]))
+                    A, E = utils.indexCube('AE', self.n+1, np.array([self.nNx, self.nNy, self.nCz]))
                     edge3 = xyz[E, :] - xyz[A, :]
-                    self._edge = np.r_[mkvc(length3D(edge1)), mkvc(length3D(edge2)), 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()
@@ -329,10 +331,10 @@ if __name__ == '__main__':
     h3 = np.cumsum(np.r_[0, np.ones(nc)/(nc)])
     dee3 = True
     if dee3:
-        X, Y, Z = ndgrid(h1, h2, h3, vector=False)
+        X, Y, Z = utils.ndgrid(h1, h2, h3, vector=False)
         M = LogicallyOrthogonalMesh([X, Y, Z])
     else:
-        X, Y = 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/TensorMesh.py

@@ -1,11 +1,8 @@
-import numpy as np
-import scipy.sparse as sp
+from SimPEG import utils, np, sp
 from BaseMesh import BaseMesh
 from TensorView import TensorView
 from DiffOperators import DiffOperators
 from InnerProducts import InnerProducts
-from SimPEG.utils import ndgrid, mkvc, spzeros, interpmat
-
 
 class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
     """
@@ -24,13 +21,20 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
 
     Example of a padded tensor mesh:
 
-    .. plot:: examples/mesh/plot_TensorMesh.py
+    .. plot::
+
+        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])
 
     """
+
+    __metaclass__ = utils.Save.Savable
+
     _meshType = 'TENSOR'
 
     def __init__(self, h_in, x0=None):
@@ -48,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 = [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)
@@ -166,7 +170,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
 
         def fget(self):
             if self._gridCC is None:
-                self._gridCC = ndgrid(self.getTensor('CC'))
+                self._gridCC = utils.ndgrid(self.getTensor('CC'))
             return self._gridCC
         return locals()
     _gridCC = None  # Store grid by default
@@ -177,7 +181,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
 
         def fget(self):
             if self._gridN is None:
-                self._gridN = ndgrid(self.getTensor('N'))
+                self._gridN = utils.ndgrid(self.getTensor('N'))
             return self._gridN
         return locals()
     _gridN = None  # Store grid by default
@@ -188,7 +192,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
 
         def fget(self):
             if self._gridFx is None:
-                self._gridFx = ndgrid(self.getTensor('Fx'))
+                self._gridFx = utils.ndgrid(self.getTensor('Fx'))
             return self._gridFx
         return locals()
     _gridFx = None  # Store grid by default
@@ -199,7 +203,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
 
         def fget(self):
             if self._gridFy is None and self.dim > 1:
-                self._gridFy = ndgrid(self.getTensor('Fy'))
+                self._gridFy = utils.ndgrid(self.getTensor('Fy'))
             return self._gridFy
         return locals()
     _gridFy = None  # Store grid by default
@@ -210,7 +214,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
 
         def fget(self):
             if self._gridFz is None and self.dim > 2:
-                self._gridFz = ndgrid(self.getTensor('Fz'))
+                self._gridFz = utils.ndgrid(self.getTensor('Fz'))
             return self._gridFz
         return locals()
     _gridFz = None  # Store grid by default
@@ -221,7 +225,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
 
         def fget(self):
             if self._gridEx is None:
-                self._gridEx = ndgrid(self.getTensor('Ex'))
+                self._gridEx = utils.ndgrid(self.getTensor('Ex'))
             return self._gridEx
         return locals()
     _gridEx = None  # Store grid by default
@@ -232,7 +236,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
 
         def fget(self):
             if self._gridEy is None and self.dim > 1:
-                self._gridEy = ndgrid(self.getTensor('Ey'))
+                self._gridEy = utils.ndgrid(self.getTensor('Ey'))
             return self._gridEy
         return locals()
     _gridEy = None  # Store grid by default
@@ -243,7 +247,7 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
 
         def fget(self):
             if self._gridEz is None and self.dim > 2:
-                self._gridEz = ndgrid(self.getTensor('Ez'))
+                self._gridEz = utils.ndgrid(self.getTensor('Ez'))
             return self._gridEz
         return locals()
     _gridEz = None  # Store grid by default
@@ -258,13 +262,13 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
                 vh = self.h
                 # Compute cell volumes
                 if(self.dim == 1):
-                    self._vol = mkvc(vh[0])
+                    self._vol = utils.mkvc(vh[0])
                 elif(self.dim == 2):
                     # Cell sizes in each direction
-                    self._vol = 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 = mkvc(np.outer(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
@@ -285,12 +289,12 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
                 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_[mkvc(area1), mkvc(area2)]
+                    self._area = np.r_[utils.mkvc(area1), utils.mkvc(area2)]
                 elif(self.dim == 3):
-                    area1 = np.outer(np.ones(n[0]+1), mkvc(np.outer(vh[1], vh[2])))
-                    area2 = np.outer(vh[0], mkvc(np.outer(np.ones(n[1]+1), vh[2])))
-                    area3 = np.outer(vh[0], mkvc(np.outer(vh[1], np.ones(n[2]+1))))
-                    self._area = np.r_[mkvc(area1), mkvc(area2), 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
@@ -307,16 +311,16 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
                 n = self.n
                 # Compute edge lengths
                 if(self.dim == 1):
-                    self._edge = 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_[mkvc(l1), mkvc(l2)]
+                    self._edge = np.r_[utils.mkvc(l1), utils.mkvc(l2)]
                 elif(self.dim == 3):
-                    l1 = np.outer(vh[0], mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1))))
-                    l2 = np.outer(np.ones(n[0]+1), mkvc(np.outer(vh[1], np.ones(n[2]+1))))
-                    l3 = np.outer(np.ones(n[0]+1), mkvc(np.outer(np.ones(n[1]+1), vh[2])))
-                    self._edge = np.r_[mkvc(l1), mkvc(l2), 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
@@ -406,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 = [spzeros(loc.shape[0], n) for n in nF_nE]
-            components[ind] = 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 = 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

@@ -372,7 +372,7 @@ class TensorView(object):
 
         return animate(fig, animateFrame, frames=mesh.nCv['xyz'.index(normal)])
 
-    def video(mesh,var,function,figsize=(10,8),colorbar=True):
+    def video(mesh, var, function, figsize=(10, 8), colorbar=True, skip=1):
         """
         Call a function for a list of models to create a video.
 
@@ -393,9 +393,11 @@ class TensorView(object):
         if colorbar:
             plt.colorbar(function(var[0],ax,clim,tlt,0))
 
-        def animateFrame(i):
+        frames = np.arange(0,len(var),skip)
+        def animateFrame(j):
+            i = frames[j]
             function(var[i],ax,clim,tlt,i)
 
-        return animate(fig, animateFrame, frames=len(var))
+        return animate(fig, animateFrame, frames=len(frames))
 
 

SimPEG/regularization/Regularization.py

@@ -1,119 +0,0 @@
-from SimPEG.utils import sdiag, count, timeIt, setKwargs
-import numpy as np
-
-class Regularization(object):
-    """docstring for Regularization"""
-
-    @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
-
-    @property
-    def Ws(self):
-        if getattr(self,'_Ws', None) is None:
-            self._Ws = sdiag(self.mesh.vol)
-        return self._Ws
-
-    @property
-    def Wx(self):
-        if getattr(self, '_Wx', None) is None:
-            self._Wx = self.mesh.cellGradx*sdiag(self.mesh.vol)
-        return self._Wx
-
-    @property
-    def Wy(self):
-        if getattr(self, '_Wy', None) is None:
-            self._Wy = self.mesh.cellGrady*sdiag(self.mesh.vol)
-        return self._Wy
-
-    @property
-    def Wz(self):
-        if getattr(self, '_Wz', None) is None:
-            self._Wz = self.mesh.cellGradz*sdiag(self.mesh.vol)
-        return self._Wz
-
-    alpha_s = 1e-6
-    alpha_x = 1.0
-    alpha_y = 1.0
-    alpha_z = 1.0
-
-    counter = None
-
-    def __init__(self, mesh, **kwargs):
-        setKwargs(self, **kwargs)
-        self.mesh = mesh
-
-
-    def pnorm(self, r):
-        return 0.5*r.dot(r)
-
-    @timeIt
-    def modelObj(self, m):
-        mresid = m - self.mref
-
-        mobj = self.alpha_s * self.pnorm( self.Ws * mresid )
-
-        mobj += self.alpha_x * self.pnorm( self.Wx * mresid )
-
-        if self.mesh.dim > 1:
-            mobj += self.alpha_y * self.pnorm( self.Wy * mresid )
-        if self.mesh.dim > 2:
-            mobj += self.alpha_z * self.pnorm( self.Wz * mresid )
-
-        return mobj
-
-    @timeIt
-    def modelObjDeriv(self, m):
-        """
-
-        In 1D:
-
-        .. math::
-
-            m_{\\text{obj}} = {1 \over 2}\\alpha_s  \left\| W_s  (m- m_{\\text{ref}})\\right\|^2_2
-                            + {1 \over 2}\\alpha_x  \left\| W_x  (m- m_{\\text{ref}})\\right\|^2_2
-
-            \\frac{ \partial m_{\\text{obj}} }{\partial m} =
-                            \\alpha_s  W_s^{\\top} W_s  (m - m_{\\text{ref}}) +
-                            \\alpha_x  W_x^{\\top} W_x  (m - m_{\\text{ref}})
-
-
-            \\frac{ \partial^2 m_{\\text{obj}} }{\partial m^2} =
-                            \\alpha_s  W_s^{\\top} W_s +
-                            \\alpha_x  W_x^{\\top} W_x
-
-        """
-
-        mresid = m - self.mref
-
-        mobjDeriv = self.alpha_s * self.Ws.T * ( self.Ws * mresid)
-
-        mobjDeriv = mobjDeriv + self.alpha_x * self.Wx.T * ( self.Wx * mresid)
-
-        if self.mesh.dim > 1:
-            mobjDeriv = mobjDeriv + self.alpha_y * self.Wy.T * ( self.Wy * mresid)
-        if self.mesh.dim > 2:
-            mobjDeriv = mobjDeriv + self.alpha_z * self.Wz.T * ( self.Wz * mresid)
-
-        return mobjDeriv
-
-
-    @timeIt
-    def modelObj2Deriv(self):
-
-        mobj2Deriv = self.alpha_s * self.Ws.T * self.Ws
-
-        mobj2Deriv = mobj2Deriv + self.alpha_x * self.Wx.T * self.Wx
-
-        if self.mesh.dim > 1:
-            mobj2Deriv = mobj2Deriv + self.alpha_y * self.Wy.T * self.Wy
-        if self.mesh.dim > 2:
-            mobj2Deriv = mobj2Deriv + self.alpha_z * self.Wz.T * self.Wz
-
-        return mobj2Deriv
-

SimPEG/regularization/__init__.py

@@ -1 +0,0 @@
-from Regularization import Regularization

SimPEG/tests/test_forward_DCproblem.py

@@ -3,10 +3,9 @@ import unittest
 from SimPEG.mesh import TensorMesh
 from SimPEG.utils import ModelBuilder, sdiag
 from SimPEG.forward import Problem
-from SimPEG.forward.DCProblem import *
+from SimPEG.examples.DC import *
 from TestUtils import checkDerivative
 from scipy.sparse.linalg import dsolve
-from SimPEG.regularization import Regularization
 from SimPEG import inverse
 
 
@@ -44,23 +43,19 @@ class DCProblemTests(unittest.TestCase):
         problem = DCProblem(mesh)
         problem.P = P
         problem.RHS = q
-        dobs, Wd = problem.createSyntheticData(mSynth, std=0.05)
+        data = problem.createSyntheticData(mSynth, std=0.05)
 
         # Now set up the problem to do some minimization
-        problem.W = Wd
-        problem.dobs = dobs
-        problem.std = dobs*0 + 0.05
-
         opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
-        reg = Regularization(mesh)
-        inv = inverse.Inversion(problem, reg, opt, beta0=1e4)
+        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.dobs = dobs
+        self.data = data
 
     def test_misfit(self):
         derChk = lambda m: [self.p.dpred(m), lambda mx: self.p.J(self.m0, mx)]
@@ -71,7 +66,7 @@ class DCProblemTests(unittest.TestCase):
         # 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.dobs.shape[0])
+        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

SimPEG/tests/test_forward_problem.py

@@ -1,8 +1,6 @@
 import numpy as np
 import unittest
-from SimPEG.mesh import TensorMesh
-from SimPEG.forward import Problem
-from SimPEG.regularization import Regularization
+from SimPEG import mesh, forward, inverse
 from TestUtils import checkDerivative
 from scipy.sparse.linalg import dsolve
 
@@ -14,9 +12,9 @@ class ProblemTests(unittest.TestCase):
         a = np.array([1, 1, 1])
         b = np.array([1, 2])
         c = np.array([1, 4])
-        self.mesh2 = TensorMesh([a, b], np.array([3, 5]))
-        self.p2 = Problem(self.mesh2)
-        self.reg = Regularization(self.mesh2)
+        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'

SimPEG/tests/test_optimizers.py

@@ -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()
+        PG = inverse.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.])
@@ -50,5 +50,15 @@ class TestOptimizers(unittest.TestCase):
         print 'x_true: ', x_true
         self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True)
 
+    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)
+        x_true = np.array([np.pi,np.pi,0])
+        print 'Newton Root Finding'
+        print 'xopt: ', xopt
+        print 'x_true: ', x_true
+        self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True)
+
 if __name__ == '__main__':
     unittest.main()

SimPEG/utils/ModelBuilder.py

@@ -172,7 +172,7 @@ def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]):
         if len(shape) is 1:
             smth = np.array([1,10.,1],dtype=float)
         elif len(shape) is 2:
-            smth = np.array([[1,2,1],[7,10,7],[1,2,1]],dtype=float)
+            smth = np.array([[1,7,1],[2,10,2],[1,7,1]],dtype=float)
         elif len(shape) is 3:
             kernal = np.array([1,4,1], dtype=float).reshape((1,3))
             smth = np.array(sp.kron(sp.kron(kernal,kernal.T).todense()[:],kernal).todense()).reshape((3,3,3))

SimPEG/utils/Save.py

@@ -0,0 +1,348 @@
+import numpy as np
+import time
+import re
+
+try:
+    import h5py
+except Exception, e:
+    print 'Warning: SimPEG.utils.Save needs h5py to be installed.'
+
+
+SAVEABLES = {}
+
+def natural_keys(text):
+    '''
+    alist.sort(key=natural_keys) sorts in human order
+    http://nedbatchelder.com/blog/200712/human_sorting.html
+    (See Toothy's implementation in the comments)
+    '''
+    atoi = lambda text: int(text) if text.isdigit() else text
+    return [ atoi(c) for c in re.split('(\d+)', text) ]
+
+
+def preIteration(group):
+    group.attrs['complete'] = False
+    group.attrs['time'] = time.time()
+
+def postIteration(group):
+    group.attrs['time'] = time.time() - group.attrs['time']
+    group.attrs['date'] = time.ctime()
+    group.attrs['complete'] = True
+
+class SimPEGTable:
+    """
+        This is a wrapper class on the HDF5 file.
+    """
+    def __init__(self, name, mode='a'):
+        if '.hdf5' not in name:
+            name += '.hdf5'
+        self.f = h5py.File(name,  mode)
+        self.root = hdf5Group(self,self.f)
+
+        self.inversions = hdf5InversionGroup(self,self.root.addGroup('inversions',soft=True))
+
+    def show(self): self.root.show()
+
+    def saveInversion(self, invObj):
+
+        # Create a new inversion anytime this is run.
+        def _startup_hdf5_inv(invObj, m0):
+            node = self.inversions.addGroup('%d'%self.inversions.numChildren)
+            saveSavable(invObj,node.addGroup('rebuild'))
+            results = node.addGroup('results')
+            preIteration(results)
+            invObj._invNode = results
+            self.f.flush()
+        invObj.hook(_startup_hdf5_inv, overwrite=True)
+
+        # 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))
+            preIteration(invObj._invNodeIt)
+        invObj.hook(_doStartIteration_hdf5_inv, overwrite=True)
+
+        def _doEndIteration_hdf5_inv(invObj):
+            invObj.save(invObj._invNodeIt)
+            postIteration(invObj._invNodeIt)
+            self.f.flush()
+        invObj.hook(_doEndIteration_hdf5_inv, overwrite=True)
+
+        # Delete all iterates that did not finish.
+        def _finish_hdf5_inv(invObj):
+            postIteration(invObj._invNode)
+            for it in invObj._invNode:
+                if not it.attrs['complete']:
+                    del self.f[it.path]
+            del invObj._invNode
+            self.f.flush()
+        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))
+            preIteration(optObj._optNodeIt)
+        invObj.opt.hook(_doStartIteration_hdf5_opt, overwrite=True)
+
+        def _doEndIteration_hdf5_opt(optObj, xt):
+            optObj.save(optObj._optNodeIt)
+            postIteration(optObj._optNodeIt)
+            self.f.flush()
+        invObj.opt.hook(_doEndIteration_hdf5_opt, overwrite=True)
+
+
+
+class hdf5Group(object):
+    """Has some low level support for wrapping the native HDF5-Group class"""
+
+    def __init__(self, T, groupNode):
+        self.T = T
+        # check if you are inputing a hdf5Group rather than a raw node, and act accordingly
+        if issubclass(groupNode.__class__, hdf5Group):
+            self.node = groupNode.node
+        else:
+            self.node = groupNode
+
+        self.childClass = hdf5Group
+        self.parentClass = hdf5Group
+
+    @property
+    def children(self):
+        """Children names in a list
+
+            Use obj[name] to return the actual node.
+        """
+        myChildren = [c for c in self.node]
+        myChildren.sort(key=natural_keys)
+        return myChildren
+
+    @property
+    def numChildren(self):
+        """Returns the len(children)"""
+        return len(self.children)
+
+    @property
+    def parent(self):
+        """Returns parent node"""
+        return self.parentClass(self.T, self.node.parent)
+
+    @property
+    def name(self):
+        return self.path.split('/')[-1]
+
+    @property
+    def path(self):
+        """Returns the root path of the group"""
+        return self.node.name
+
+    @property
+    def attrs(self):
+        """Returns a list of attributes in the group"""
+        return self.node.attrs
+
+    def addGroup(self, name, soft=False):
+        """Adds a child group to the current node."""
+        if name in self.children and soft:
+            return self[name]
+        assert name not in self.children, 'Already a child called: '+self.path+'/'+name
+        return self.childClass(self.T, self.node.create_group(name))
+
+    def setArray(self, name, data):
+        a = self.node.create_dataset(name, data.shape)
+        a[...] = data
+        return a
+
+    def __getitem__(self, val):
+        if type(val) is int:
+            val = self.children[val]
+        child = self.node[val]
+        if type(child) is h5py.Group:
+            child = self.childClass(self.T, child)
+        return child
+
+    def __contains__(self, key):
+        return key in self.children
+
+    def show(self, pad='', maxDepth=1, depth=0):
+        """
+            Recursively show the structure of the database.
+
+            For example::
+
+                <hdf5InversionGroup group "/inversions" (1 member)>
+                    - <hdf5Inversion group "/inversions/0" (4 members)>
+                        - <hdf5InversionIteration group "/inversions/0/0.0" (3 members)>
+                        - <hdf5InversionIteration group "/inversions/0/0.1" (3 members)>
+                        - <hdf5InversionIteration group "/inversions/0/0.2" (3 members)>
+                        - <hdf5InversionIteration group "/inversions/0/0.3" (3 members)>
+        """
+        s = self.__str__()
+        pad += ' '*4
+        if maxDepth <= 0: print s
+        if depth >= maxDepth: return s
+
+        for c in self.children:
+            if issubclass(self[c].__class__, hdf5Group):
+                s += '\n%s- %s' % (pad, self[c].show(pad=pad,depth=depth+1,maxDepth=maxDepth))
+            else:
+                s += '\n%s- %s' % (pad, self[c].__str__())
+        if depth is 0:
+            print s
+        else:
+            return s
+
+    def __str__(self):
+        return '<%s "%s" (%d member%s, attrs=[%s])>' % (self.__class__.__name__, self.path, self.numChildren, '' if self.numChildren == 1 else 's',', '.join([a for a in self.attrs]))
+
+
+
+class hdf5InversionGroup(hdf5Group):
+    def __init__(self, T, groupNode):
+        hdf5Group.__init__(self, T, groupNode)
+        self.childClass = hdf5Inversion
+
+class hdf5Inversion(hdf5Group):
+    def __init__(self, T, groupNode):
+        hdf5Group.__init__(self, T, groupNode)
+        self.parentClass = hdf5InversionGroup
+        self.childClass = hdf5InversionResults
+
+    def rebuild(self):
+        return loadSavable(self['rebuild'])
+
+class hdf5InversionResults(hdf5Group):
+    def __init__(self, T, groupNode):
+        hdf5Group.__init__(self, T, groupNode)
+        self.parentClass = hdf5Inversion
+        self.childClass = hdf5InversionIteration
+
+class hdf5InversionIteration(hdf5Group):
+    def __init__(self, T, groupNode):
+        hdf5Group.__init__(self, T, groupNode)
+        self.parentClass = hdf5InversionResults
+
+
+
+class Savable(type):
+    def __new__(cls, name, bases, attrs):
+        __init__ = attrs['__init__']
+        def init_wrapper(self, *args, **kwargs):
+            self._args_init = args
+            self._kwargs_init = kwargs
+            return __init__(self, *args,**kwargs)
+        attrs['__init__'] = init_wrapper
+
+        newClass = super(Savable, cls).__new__(cls, name, bases, attrs)
+        SAVEABLES[name] = newClass
+        return newClass
+
+
+def saveSavable(obj, group, debug=False):
+    """
+        This creates softlinks if _savable exists in children object.
+
+        The first object is always created.
+    """
+    assert type(obj.__class__) is Savable, 'Can only save objects that are Savable objects.'
+
+    def doSave(grp, name, val):
+        if debug: print name, val
+        if type(val.__class__) is Savable:
+            link = getattr(val,'_savable',None)
+            if link is not None:
+                group.node[name] = h5py.SoftLink(link.path)
+                if debug: 'Created a softlink path to %s' % link.path
+            else:
+                subgrp = grp.addGroup(name)
+                saveSavable(val, subgrp, debug=debug)
+        elif type(val) in [list, tuple]:
+            # Split up, and save each element
+            for i, v in enumerate(val):
+                doSave(grp, name + '[%d]'%i, v)
+        elif type(val) is np.ndarray:
+            grp.setArray(name, val)
+        elif val is None:
+            grp.attrs[name] = 'None'
+        else:
+            # just try saving it as an attr
+            try:
+                grp.attrs[name] = val
+            except Exception, e:
+                print 'Warning: Could not save %s, problems may arise in loading.' % name
+
+    group.attrs['__class__'] = obj.__class__.__name__
+    for arg in obj._kwargs_init:
+        doSave(group, '_kwarg_'+arg, obj._kwargs_init[arg])
+    for i, arg in enumerate(obj._args_init):
+        doSave(group, '_arg%d'%i, arg)
+    obj._savable = group
+
+
+def loadSavable(node, pointers=None):
+    """
+        pointers allow things that point to the same node in the h5py file to
+        be returned as the same object, if they have already been created.
+    """
+
+    if pointers is None: pointers = []
+    for pointer in pointers:
+        if pointer._savable.node == node.node: return pointer
+
+    args = ([a for a in node.attrs if '_arg' in a] + [a for a in node.children if '_arg' in a])
+    kwargs = ([a for a in node.attrs if '_kwarg' in a] + [a for a in node.children if '_kwarg' in a])
+    args.sort(key=natural_keys)
+    kwargs.sort(key=natural_keys)
+
+    def get(node,key):
+        if key in node.children: return node[key]
+        elif key in node.attrs: return node.attrs[key]
+
+    ARGS = []
+    for name in args:
+        val = get(node, name)
+        if val.__class__ is h5py.Dataset: val = val[:]
+        if val is 'None': val = None
+        if '[' in name:  # We are reloading a list
+            ind = int(name[4:name.index('[')])
+            if len(ARGS) is ind: # Create the list
+                ARGS.append([val])
+            else:
+                ARGS[ind].append(val)
+        elif issubclass(val.__class__,hdf5Group):
+            ARGS.append(loadSavable(val,pointers=pointers))
+        else:
+            ind = int(name[4:])
+            ARGS.append(val)
+
+    KWARGS = {}
+    for name in kwargs:
+        val = get(node, name)
+        if val.__class__ is h5py.Dataset: val = val[:]
+        if val is 'None': val = None
+        if '[' in name:  # We are reloading a list
+            key = name[7:name.index('[')]
+            if key not in KWARGS: # Create the list
+                KWARGS[key] = [val]
+            else:
+                KWARGS[key].append(val)
+        elif issubclass(val.__class__,hdf5Group):
+            key = name[7:]
+            KWARGS[key] = loadSavable(val,pointers=pointers)
+        else:
+            key = name[7:]
+            KWARGS[key] = val
+
+    cls = get(node, '__class__')
+    if cls in SAVEABLES:
+        try:
+            out = SAVEABLES[cls](*ARGS, **KWARGS)
+            out._savable = node
+            pointers.append(out)  # Because this is recursive.
+            return out
+        except Exception, e:
+            print 'Warning: %s Class could not be initiated.' % cls
+            print 'ARGS: ', ARGS
+            print 'KWARGS: ', KWARGS
+            return (cls, ARGS, KWARGS, node)
+    else:
+        print 'Warning: %s Class not found in SimPEG.utils.Save.SAVABLES' % cls
+        return (cls, ARGS, KWARGS, node)
+

SimPEG/utils/Solver.py

@@ -73,11 +73,9 @@ class Solver(object):
                 Jacobi = sdiag(1.0/M[1].diagonal())
                 options['M'] = Jacobi
             elif M[0] is 'GS':
-                LL = sp.tril(M[1])
-                UU = sp.triu(M[1])
                 DD = sdiag(M[1].diagonal())
-                Uinv = Solver(UU, flag='U')
-                Linv = Solver(LL, flag='L')
+                Uinv = Solver(M[1], flag='U')
+                Linv = Solver(M[1], flag='L')
                 def GS(f):
                     return Uinv.solve(DD*Linv.solve(f))
                 options['M'] = sp.linalg.LinearOperator( A.shape, GS, dtype=A.dtype )

SimPEG/utils/__init__.py

@@ -3,13 +3,16 @@ 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 lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube, exampleLomGird
+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
+import Save
 import Geophysics
 
 import types
@@ -23,7 +26,7 @@ def hook(obj, method, name=None, overwrite=False, silent=False):
 
         If name is None, the name of the method is used.
     """
-    if name is None: 
+    if name is None:
         name = method.__name__
         if name == '<lambda>':
             raise Exception('Must provide name to hook lambda functions.')
@@ -42,7 +45,7 @@ def setKwargs(obj, **kwargs):
             setattr(obj, attr, kwargs[attr])
         else:
             raise Exception('%s attr is not recognized' % attr)
-    hook(obj,callHooks, silent=True)
+
     hook(obj,hook, silent=True)
     hook(obj,setKwargs, silent=True)
 
@@ -87,11 +90,48 @@ def printStoppers(obj, stoppers, pad='', stop='STOP!', done='DONE!'):
         print pad + stopper['str'] % (l<=r,l,r)
     print pad + "%s%s%s" % ('-'*25,done,'-'*25)
 
-def callHooks(obj, match, *args, **kwargs):
-    for method in [posible for posible in dir(obj) if ('_'+match) in posible]:
-            if getattr(obj,'debug',False): print (match+' is calling self.'+method)
-            getattr(obj,method)(*args, **kwargs)
+def callHooks(match):
+    """
+    Use this to wrap a funciton::
 
+        @callHooks('doEndIteration')
+        def doEndIteration(self):
+            pass
+
+    This will call everything named _doEndIteration* at the beginning of the function call.
+    """
+    def callHooksWrap(f):
+        @wraps(f)
+        def wrapper(self,*args,**kwargs):
+
+            for method in [posible for posible in dir(self) if ('_'+match) in posible]:
+                if getattr(self,'debug',False): print (match+' is calling self.'+method)
+                getattr(self,method)(*args, **kwargs)
+
+            return f(self,*args,**kwargs)
+
+        extra = """
+            If you have things that also need to run in the method %s, you can create a method::
+
+                def _%s*(self, ... ):
+                    pass
+
+            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)
+
+        wrapper.__doc__ += extra
+        return wrapper
+    return callHooksWrap
+
+def dependentProperty(name, value, children, doc):
+    def fget(self): return getattr(self,name,value)
+    def fset(self, val):
+        for child in children:
+            if hasattr(self, child):
+                delattr(self, child)
+        setattr(self, name, val)
+    return property(fget=fget, fset=fset, doc=doc)
 
 
 class Counter(object):

SimPEG/utils/lomutils.py

@@ -4,29 +4,6 @@ from matutils import mkvc, ndgrid, sub2ind
 from sputils import sdiag
 
 
-def exampleLomGird(nC, exType):
-    assert type(nC) == list, "nC must be a list containing the number of nodes"
-    assert len(nC) == 2 or len(nC) == 3, "nC must either two or three dimensions"
-    exType = exType.lower()
-
-    possibleTypes = ['rect', 'rotate']
-    assert exType in possibleTypes, "Not a possible example type."
-
-    if exType == 'rect':
-        return ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)
-    elif exType == 'rotate':
-        if len(nC) == 2:
-            X, Y = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)
-            amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2)
-            amt[amt < 0] = 0
-            return X + (-(Y - 0.5))*amt, Y + (+(X - 0.5))*amt
-        elif len(nC) == 3:
-            X, Y, Z = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)
-            amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2 + (Z - 0.5)**2)
-            amt[amt < 0] = 0
-            return X + (-(Y - 0.5))*amt, Y + (-(Z - 0.5))*amt, Z + (-(X - 0.5))*amt
-
-
 def volTetra(xyz, A, B, C, D):
     """
     Returns the volume for tetrahedras volume specified by the indexes A to D.

SimPEG/utils/meshutils.py

@@ -0,0 +1,58 @@
+import numpy as np
+from scipy import sparse as sp
+from matutils import mkvc, ndgrid, sub2ind
+from sputils import sdiag
+
+def exampleLomGird(nC, exType):
+    assert type(nC) == list, "nC must be a list containing the number of nodes"
+    assert len(nC) == 2 or len(nC) == 3, "nC must either two or three dimensions"
+    exType = exType.lower()
+
+    possibleTypes = ['rect', 'rotate']
+    assert exType in possibleTypes, "Not a possible example type."
+
+    if exType == 'rect':
+        return ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)
+    elif exType == 'rotate':
+        if len(nC) == 2:
+            X, Y = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)
+            amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2)
+            amt[amt < 0] = 0
+            return X + (-(Y - 0.5))*amt, Y + (+(X - 0.5))*amt
+        elif len(nC) == 3:
+            X, Y, Z = ndgrid([np.cumsum(np.r_[0, np.ones(nx)/nx]) for nx in nC], vector=False)
+            amt = 0.5-np.sqrt((X - 0.5)**2 + (Y - 0.5)**2 + (Z - 0.5)**2)
+            amt[amt < 0] = 0
+            return X + (-(Y - 0.5))*amt, Y + (-(Z - 0.5))*amt, Z + (-(X - 0.5))*amt
+
+
+def meshTensors(*args):
+    """
+        **meshTensors** takes any number of tuples that have the form::
+
+            h1 = ( (numPad, sizeStart [, increaseFactor]), (numCore, sizeCode), (numPad, sizeStart [, increaseFactor]) )
+
+        .. plot::
+
+            from SimPEG import mesh, utils
+            M = mesh.TensorMesh(utils.meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0))))
+            M.plotGrid()
+
+    """
+    def padding(num, start, factor=1.3, reverse=False):
+        pad = ((np.ones(num)*factor)**np.arange(num))*start
+        if reverse: pad = pad[::-1]
+        return pad
+    tensors = tuple()
+    for i, arg in enumerate(args):
+        tensors += (np.r_[padding(*arg[0],reverse=True),np.ones(arg[1][0])*arg[1][1],padding(*arg[2])],)
+
+    return list(tensors) if len(tensors) > 1 else tensors[0]
+
+if __name__ == '__main__':
+    from SimPEG import mesh
+    import matplotlib.pyplot as plt
+    M = mesh.TensorMesh(meshTensors(((10,10),(40,10),(10,10)), ((10,10),(20,10),(0,0))))
+    M.plotGrid()
+    plt.gca().axis('tight')
+    plt.show()

SimPEG/visualize/vtk/vtkTools.py

@@ -46,7 +46,7 @@ class vtkTools(object):
 			vZ = mesh.vectorNz
 			zD = mesh.nNz
 		# Use rectilinear VTK grid.
-		# Asaign the spatial information.
+		# Assign the spatial information.
 		vtkObj = vtk.vtkRectilinearGrid()
 		vtkObj.SetDimensions(xD,yD,zD)
 		vtkObj.SetXCoordinates(npsup.numpy_to_vtk(vX,deep=1))
@@ -61,6 +61,7 @@ class vtkTools(object):
 			vtkObj.GetCellData().AddArray(vtkDoubleArr)
 
 		vtkObj.GetCellData().SetActiveScalars(model.keys()[0])
+		vtkObj.Update()
 		return vtkObj
 
 	@staticmethod
@@ -104,7 +105,7 @@ class vtkTools(object):
 		# Cells -cell array
 		FCellArr = vtk.vtkCellArray()
 		FCellArr.SetNumberOfCells(mesh.nF)
-		FCellArr.SetCells(mesh.nF*5,npsup.numpy_to_vtkIdTypeArray(np.vstack([FxCellBlock,FyCellBlock,FzCellBlock]),deep=1))
+		FCellArr.SetCells(mesh.nF,npsup.numpy_to_vtkIdTypeArray(np.vstack([FxCellBlock,FyCellBlock,FzCellBlock]),deep=1))
 		# Cell type
 		FCellType = npsup.numpy_to_vtk(vtk.VTK_QUAD*np.ones(mesh.nF,dtype='uint8'),deep=1)
 		# Cell location
@@ -166,7 +167,7 @@ class vtkTools(object):
 		# Cells -cell array
 		ECellArr = vtk.vtkCellArray()
 		ECellArr.SetNumberOfCells(mesh.nE)
-		ECellArr.SetCells(mesh.nE*3,npsup.numpy_to_vtkIdTypeArray(np.vstack([ExCellBlock,EyCellBlock,EzCellBlock]),deep=1))
+		ECellArr.SetCells(mesh.nE,npsup.numpy_to_vtkIdTypeArray(np.vstack([ExCellBlock,EyCellBlock,EzCellBlock]),deep=1))
 		# Cell type
 		ECellType = npsup.numpy_to_vtk(vtk.VTK_LINE*np.ones(mesh.nE,dtype='uint8'),deep=1)
 		# Cell location
@@ -186,6 +187,7 @@ class vtkTools(object):
 			vtkObj.GetCellData().AddArray(vtkDoubleArr)
 
 		vtkObj.GetCellData().SetActiveScalars(model.keys()[0])
+		vtkObj.Update()
 		return vtkObj
 
 	@staticmethod
@@ -255,8 +257,7 @@ class vtkTools(object):
 		cellThres = vtk.vtkThreshold()
 		cellThres.AllScalarsOn()
 		cellThres.SetInputConnection(cellCore.GetOutputPort())
-		cellThres.ThresholdByUpper(limits[0])
-		cellThres.ThresholdByLower(limits[1])
+		cellThres.ThresholdBetween(limits[0],limits[1])
 		cellThres.Update()
 		return cellThres.GetOutput(), cellCore.GetOutput()
 
@@ -271,8 +272,7 @@ class vtkTools(object):
 		cellThres = vtk.vtkThreshold()
 		cellThres.AllScalarsOn()
 		cellThres.SetInputConnection(cellCore.GetOutputPort())
-		cellThres.ThresholdByUpper(limits[0])
-		cellThres.ThresholdByLower(limits[1])
+		cellThres.ThresholdBetween(limits[0],limits[1])
 		cellThres.Update()
 		return cellThres.GetOutput(), cellCore.GetOutput()
 

SimPEG/visualize/vtk/vtkView.py

@@ -1,6 +1,6 @@
-import numpy as np
+import numpy as np, matplotlib as mpl
 try:
-	import vtk
+	import vtk, vtk.util.numpy_support as npsup
 	#import SimPEG.visualize.vtk.vtkTools as vtkSP # Always get an error for this import
 except Exception, e:
 	print 'VTK import error. Please ensure you have VTK installed to use this visualization package.'
@@ -8,16 +8,16 @@ import SimPEG as simpeg
 
 class vtkView(object):
 	"""
-	Class for storing and view of SimPEG models in VTK (visulization toolkit).
+	Class for storing and view of SimPEG models in VTK (visualization toolkit).
 
 	Inputs:
 	:param mesh, SimPEG mesh.
 	:param propdict, dictionary of property models.
 		Can have these dictionary names:
-		'cell' - cell model; 'face' - face model; 'edge' - edge model
-		The dictionary properties are given as dictionaries with:
+		'C' - cell model; 'F' - face model; 'E' - edge model; ('V' - vector field : NOT SUPPORTED)
+		The dictionary values are given as dictionaries with:
 		{'NameOfThePropertyModel': np.array of the properties}.
-		The property array has to be ordered in compliance with SimPEG standards.
+		The property np.array has to be ordered in compliance with SimPEG standards.
 
 	::
 	Example of usages.
@@ -30,22 +30,7 @@ class vtkView(object):
 		"""
 		"""
 
-		# ToDo: Set the properties up so that there are set/get methods
-		self.name = 'VTK figure of SimPEG model'
-		self.extent = [0,mesh.nCx-1,0,mesh.nCy-1,0,mesh.nCz-1]
-		self.limits = [0, 1e12]
-		self.viewprop = {'cell':0} # Name of the tyep and Int order of the array or name of the vector.
-		self._mesh = mesh
-
-
-		# Set vtk object containers
-		self._cell = None
-		self._faces = None
-		self._edges = None
-
-		self._readPropertyDictionary(propdict)
-
-		# Setup hidden properties
+		# Setup hidden properties, used for the visualization
 		self._ren = None
 		self._iren = None
 		self._renwin = None
@@ -56,6 +41,164 @@ class vtkView(object):
 		self._widget = None
 		self._actor = None
 		self._lut = None
+		# Set vtk object containers
+		self._cells = None
+		self._faces = None
+		self._edges = None
+		self._vectors = None # Not implemented
+		# Set default values
+		self.name = 'VTK figure of SimPEG model'
+		
+
+
+		# Error check the input mesh
+		if type(mesh).__name__ != 'TensorMesh':
+			raise Exception('The input {:s} to vtkView has to be a TensorMesh object'.format(mesh))
+		# Set the mesh
+		self._mesh = mesh
+
+		# Read the property dictionary 
+		self._readPropertyDictionary(propdict)
+		
+
+		
+		
+	# Set/Get properties
+	@property 
+	def cmap(self):
+		''' Colormap to use in vtkView. Colormap is a matplotlib cmap(cm) array, has to be uint8(use flag bytes=True during cmap generation).'''
+		if getattr(self,'_cmap',None) is None:
+			# Set default
+			self._cmap = mpl.cm.hsv(np.arange(0.,1.,0.05),bytes=True)
+		return self._cmap
+	@cmap.setter
+	def cmap(self,value):
+		if value.min() > 0 or value.max() < 255 or value.shape[1] != 4 or value.dtype != np.uint8:
+			raise Exception('Input not an allowed array.\n Use matplotlib.cm to generate an array of size [nrColors,4] and dtype = uint8(flag bytes=True).')
+		self._cmap = value
+	
+	@property 
+	def range(self):
+		''' Range of the colors in vtkView.'''
+		if getattr(self,'_range',None) is None:
+			self._range = np.array(self._getActiveVTKobj().GetArray(self.viewprop.values()[0]).GetRange())
+		return self._range
+	@range.setter
+	def range(self,value):
+		if type(value) not in [tuple, list, np.ndarray] or len(value) != 2 or np.array(value).dtype is not np.dtype('float'):
+			raise Exception('Input not in correct format. \n Has to be a list, tuple or np.arry of 2 floats.')
+		self._range = np.array(value)
+
+	@property 
+	def extent(self):
+		''' Extent of the sub-domain of the model to view'''
+		if getattr(self,'_extent',None) is None:
+			self._extent = [0,self._mesh.nCx-1,0,self._mesh.nCy-1,0,self._mesh.nCz-1]
+		return self._extent
+	@extent.setter
+	def extent(self,value):
+
+		import warnings
+		# Error check
+		valnp = np.array(value,dtype=int)
+		if valnp.dtype != int or len(valnp) != 6:
+			raise Exception('.extent has to be list or nparray of 6 integers.')
+		# Test the range of the values
+		loB = np.zeros(3,dtype=int)
+		upB = np.array(self._mesh.nCv - np.ones(3),dtype=int)
+		# Test the bounds
+		change = 0
+		# Test for lower bounds, can't be smaller the 0
+		tlb = valnp[::2] < loB
+		if tlb.any(): 
+			valnp[::2][tlb] = loB[tlb] 
+			change = 1
+			warnings.warn('Lower bounds smaller then 0')
+		# Test for lower bounds, can't be larger then upB
+		tlub = valnp[::2] > upB
+		if tlub.any(): 
+			valnp[::2][tlub] = upB[tlub] - 1  
+			change = 1
+			warnings.warn('Lower bounds larger then uppermost bounds')
+		# Test for upper bounds, can't be larger the extent of the mesh
+		tub = valnp[1::2] > upB
+		if tub.any(): 
+			valnp[1::2][tub] = upB[tub] 
+			change = 1
+			warnings.warn('Upper bounds greater then number of cells')
+		# Test if lower is smaller the upper
+		tgt = valnp[::2] > valnp[1::2]
+		if tgt.any():
+			valnp[1::2][tgt] = valnp[::2][tgt] + 1
+			change = 1
+			warnings.warn('Lower bounds greater the Upper bounds')
+		# Print a warning
+		if change:
+			warnings.warn('Changed given extent from {:s} to {:s}'.format(value,valnp.tolist()))
+		
+		# Set extent
+		self._extent = valnp
+
+	@property
+	def limits(self):
+		''' Lower and upper limits (cutoffs) of the values to view. '''
+		return getattr(self,'_limits',None)
+	@limits.setter
+	def limits(self,value):
+		if value is None:
+			self._limits = None
+		else:
+			valnp = np.array(value)
+			if valnp.dtype != float or len(valnp) != 2:
+				raise Exception('.limits has to be list or numpy array of 2 floats.')
+			self._limits = valnp
+
+
+	@property 
+	def viewprop(self):
+		''' Controls the property that will be viewed.'''
+
+		if getattr(self,'_viewprop',None) is None:
+			self._viewprop = {'C':0} # Name of the type and Int order of the array or name of the vector.
+		return self._viewprop
+	@viewprop.setter
+	def viewprop(self,value):
+		if type(value) != dict:
+			raise Exception('{:s} has to be a python dictionary containing property type and name index. ')
+		if len(value) > 1:
+			raise Exception('Too many input items in the viewprop dictionary')
+		if value.keys()[0] not in ['C','F','E']:
+			raise Exception('\"{:s}\" is not allowed as a dictionary key. Can be \'C\',\'F\',\'E\'.'.format(propitem[0]))
+		if not(type(self.viewprop.values()[0]) is int or type(self.viewprop.values()[0]) is str):
+			raise Exception('The vtkView.viewprop.values()[0] has the wrong format. Has to be integer or a string with the index.')
+		
+
+		self._viewprop = value
+
+	def _getActiveVTKobj(self):
+		"""
+		Finds the active VTK object.
+		"""
+
+		if self.viewprop.keys()[0] is 'C':
+			vtkCellData = self._cells.GetCellData()
+		elif self.viewprop.keys()[0] is 'F':
+			vtkCellData = self._faces.GetCellData()
+		elif self.viewprop.keys()[0] is 'E':
+			vtkCellData = self._edges.GetCellData()
+
+		return vtkCellData
+
+	def _getActiveArrayName(self):
+		"""
+		Finds the name of the active array.
+		"""
+		actArr = self.viewprop.values()[0]
+		if type(actArr) is str:
+			activeName = actArr
+		elif type(actArr) is int:
+			activeName = self._getActiveVTKobj().GetArrayName(actArr)
+		return activeName
 
 	def _readPropertyDictionary(self,propdict):
 		"""
@@ -64,18 +207,20 @@ class vtkView(object):
 		import SimPEG.visualize.vtk.vtkTools as vtkSP
 
 		# Test the property dictionary
-		if len(propdict) > 3:
-			raise(Exception,'Too many input items in the property dictionary')
+		if type(propdict) != dict:
+			raise Exception('{:s} has to be a python dictionary containing property models. ')
+		if len(propdict) > 4:
+			raise Exception('Too many input items in the property dictionary')
 		for propitem in propdict.iteritems():
-			if propitem[0] in ['cell','face','edge']:
-				if propitem[0] == 'cell':
-					self._cell = vtkSP.makeCellVTKObject(self._mesh,propitem[1])
-				if propitem[0] == 'face':
-					self._face = vtkSP.makeFaceVTKObject(self._mesh,propitem[1])
-				if propitem[0] == 'edge':
-					self._edge = vtkSP.makeEdgeVTKObject(self._mesh,propitem[1])
+			if propitem[0] in ['C','F','E']:
+				if propitem[0] == 'C':
+					self._cells = vtkSP.makeCellVTKObject(self._mesh,propitem[1])
+				if propitem[0] == 'F':
+					self._faces = vtkSP.makeFaceVTKObject(self._mesh,propitem[1])
+				if propitem[0] == 'E':
+					self._edges = vtkSP.makeEdgeVTKObject(self._mesh,propitem[1])
 			else:
-				raise(Exception,'{:s} is not allowed as a dictonary key. Can be \'cell\',\'face\',\'edge\'.'.format(propitem[0]))
+				raise Exception('\"{:s}\" is not allowed as a dictionary key. Can be \'C\',\'F\',\'E\'.'.format(propitem[0]))
 
 	def Show(self):
 		"""
@@ -89,27 +234,34 @@ class vtkView(object):
 		# Make renderwindow. Returns the interactor.
 		self._iren, self._renwin = vtkSP.makeRenderWindow(self._ren)
 
-		imageType = self.viewprop.keys()[0]
-		# Sort out the actor
-		if imageType == 'cell':
-			self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._cell,self.extent,self.limits)
-		elif imageType == 'face':
-			extent = [self._mesh.vectorNx[self.extent[0]], self._mesh.vectorNx[self.extent[1]], self._mesh.vectorNy[self.extent[2]], self._mesh.vectorNy[self.extent[3]], self._mesh.vectorNz[self.extent[4]], self._mesh.vectorNz[self.extent[5]] ]
-			self._vtkobj, self._core = vtkSP.makeUnstructVTKVOIThres(self._face,extent,self.limits)
-		elif imageType == 'edge':
-			extent = [self._mesh.vectorNx[self.extent[0]], self._mesh.vectorNx[self.extent[1]], self._mesh.vectorNy[self.extent[2]], self._mesh.vectorNy[self.extent[3]], self._mesh.vectorNz[self.extent[4]], self._mesh.vectorNz[self.extent[5]] ]
-			self._vtkobj, self._core = vtkSP.makeUnstructVTKVOIThres(self._edge,extent,self.limits)
-		else:
-			raise Exception("{:s} is not a vailid imageType. Has to be 'cell':'face':'edge'".format(imageType))
-
+		
 		# Set the active scalar.
 		if type(self.viewprop.values()[0]) == int:
-			actScalar = self._vtkobj.GetCellData().GetArrayName(self.viewprop.values()[0])
+			actScalar = self._getActiveVTKobj().GetArrayName(self.viewprop.values()[0])
 		elif type(self.viewprop.values()[0]) == str:
 			actScalar = self.viewprop.values()[0]
 		else :
 			raise Exception('The vtkView.viewprop.values()[0] has the wrong format. Has to be interger or a string.')
-		self._vtkobj.GetCellData().SetActiveScalars(actScalar)
+		self._getActiveVTKobj().SetActiveScalars(actScalar)
+		# Sort out the actor
+		imageType = self.viewprop.keys()[0]
+		if imageType == 'C':
+			if self.limits is None:
+				self.limits = self._cells.GetCellData().GetArray(self.viewprop.values()[0]).GetRange()
+			self._vtkobj, self._core = vtkSP.makeRectiVTKVOIThres(self._cells,self.extent,self.limits)
+		elif imageType == 'F':
+			if self.limits is None:
+				self.limits = self._faces.GetCellData().GetArray(self.viewprop.values()[0]).GetRange()
+			extent = [self._mesh.vectorNx[self.extent[0]], self._mesh.vectorNx[self.extent[1]], self._mesh.vectorNy[self.extent[2]], self._mesh.vectorNy[self.extent[3]], self._mesh.vectorNz[self.extent[4]], self._mesh.vectorNz[self.extent[5]] ]
+			self._vtkobj, self._core = vtkSP.makeUnstructVTKVOIThres(self._faces,extent,self.limits)
+		elif imageType == 'E':
+			if self.limits is None:
+				self.limits = self._edges.GetCellData().GetArray(self.viewprop.values()[0]).GetRange()
+			extent = [self._mesh.vectorNx[self.extent[0]], self._mesh.vectorNx[self.extent[1]], self._mesh.vectorNy[self.extent[2]], self._mesh.vectorNy[self.extent[3]], self._mesh.vectorNz[self.extent[4]], self._mesh.vectorNz[self.extent[5]] ]
+			self._vtkobj, self._core = vtkSP.makeUnstructVTKVOIThres(self._edges,extent,self.limits)
+		else:
+			raise Exception("{:s} is not a valid viewprop. Has to be 'C':'F':'E'".format(imageType))
+		#self._vtkobj.GetCellData().SetActiveScalars(actScalar)
 		# Set up the plane, clipper and the user interaction.
 		global intPlane, intActor
 		self._clipper, intPlane = vtkSP.makePlaneClipper(self._vtkobj)
@@ -125,15 +277,27 @@ class vtkView(object):
 
 		self._widget.AddObserver("InteractionEvent",movePlane)
 		lut = vtk.vtkLookupTable()
-		lut.SetNumberOfColors(256)
-		lut.SetHueRange(0,0.66667)
+		lut.SetNumberOfColors(len(self.cmap))
+		lut.SetTable(npsup.numpy_to_vtk(self.cmap))
 		lut.Build()
 		self._lut = lut
+		scalarBar = vtk.vtkScalarBarActor()
+  		scalarBar.SetLookupTable(lut)
+ 		scalarBar.SetTitle(self._getActiveArrayName())
+  		scalarBar.GetPositionCoordinate().SetCoordinateSystemToNormalizedViewport()
+  		scalarBar.GetPositionCoordinate().SetValue(0.1,0.01)
+  		scalarBar.SetOrientationToHorizontal()
+  		scalarBar.SetWidth(0.8)
+  		scalarBar.SetHeight(0.17)
+
+		self._actor.GetMapper().SetScalarRange(self.range)
 		self._actor.GetMapper().SetLookupTable(lut)
 
 		# Set renderer options
 		self._ren.SetBackground(.5,.5,.5)
 		self._ren.AddActor(self._actor)
+		self._ren.AddActor2D(scalarBar)
+  		self._renwin.SetSize(450,450)
 
 		# Start the render Window
 		vtkSP.startRenderWindow(self._iren)
@@ -144,17 +308,43 @@ class vtkView(object):
 
 
 if __name__ == '__main__':
+	
+
 	#Make a mesh and model
 	x0 = np.zeros(3)
-	h1 = np.ones(20)*50
-	h2 = np.ones(10)*100
-	h3 = np.ones(5)*200
+	h1 = np.ones(60)*50
+	h2 = np.ones(60)*100
+	h3 = np.ones(50)*200
 
 	mesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)
 
 	# Make a models that correspond to the cells, faces and edges.
-	models = {'cell':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'face':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'edge':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}}
+	t = np.ones(mesh.nC)
+	t[10000:50000] = 100
+	t[100000:120000] = 100
+	t[100000:120000] = 50
+	models = {'C':{'Test':np.arange(0,mesh.nC),'Model':t, 'AllOnce':np.ones(mesh.nC)},'F':{'Test':np.arange(0,mesh.nF),'AllOnce':np.ones(mesh.nF)},'E':{'Test':np.arange(0,mesh.nE),'AllOnce':np.ones(mesh.nE)}}
 	# Make the vtk viewer object.
 	vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models)
+	# Set the .viewprop for which model to view
+	vtkViewer.viewprop = {'F':'Test'}
 	# Show the image
 	vtkViewer.Show()
+
+	# Set subset of the mesh to view (remove padding)
+	vtkViewer.extent = [4,14,0,7,0,3]
+	vtkViewer.Show()
+
+	# Change viewing property 
+	vtkViewer.viewprop = {'C':'Model'}
+	# Set the color range
+	# Reset extent.
+	vtkViewer.extent = [-1,1000,-1,1000,-1,1000]
+	vtkViewer.range = [0.,100.]
+	vtkViewer.Show()
+	# Change color scale, has to be set to bytes=True.
+	vtkViewer.cmap = mpl.cm.copper(np.arange(0.,1.,0.01),bytes=True)
+	vtkViewer.Show()
+	# Set limits of values to view 
+	vtkViewer.limits = [5.0,100.0]
+	vtkViewer.Show()

docs/api_Optimize.rst

@@ -24,3 +24,12 @@ Beta Schedule
 .. automodule:: SimPEG.inverse.BetaSchedule
     :members:
     :undoc-members:
+
+
+Regularization
+**************
+
+.. automodule:: SimPEG.inverse.Regularization
+    :show-inheritance:
+    :members:
+    :undoc-members:

docs/api_Problem.rst

@@ -14,7 +14,7 @@ Problem
 DCProblem
 *********
 
-.. automodule:: SimPEG.forward.DCProblem
+.. automodule:: SimPEG.examples.DC
     :show-inheritance:
     :members:
     :undoc-members:
@@ -25,7 +25,7 @@ DCProblem
 Linear Problem
 **************
 
-.. automodule:: SimPEG.forward.LinearProblem
+.. automodule:: SimPEG.examples.Linear
     :show-inheritance:
     :members:
     :undoc-members:

docs/api_Utils.rst

@@ -37,6 +37,13 @@ LOM Utilities
     :members:
     :undoc-members:
 
+Mesh Utilities
+**************
+
+.. automodule:: SimPEG.utils.meshutils
+    :members:
+    :undoc-members:
+
 Model Builder Utilities
 ***********************
 

docs/examples/mesh/plot_TensorMesh.py

@@ -1,21 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-from SimPEG.mesh import TensorMesh
-
-pad = 7
-padfactor = 1.4
-xpad = (np.ones(pad)*padfactor)**np.arange(pad)
-ypad = (np.ones(pad)*padfactor)**np.arange(pad)
-
-core = 15
-xcore = np.ones(core)
-ycore = np.ones(core)
-
-h1 = np.r_[xpad[::-1],xcore,xpad]
-h2 = np.r_[ypad[::-1],ycore,ypad]
-
-mesh = TensorMesh([h1, h2])
-mesh.plotGrid()
-plt.axis('tight')
-
-plt.show()

docs/index.rst

@@ -3,14 +3,23 @@
    :alt: SimPEG
    :align: center
 
-SimPEG (Simulation and Parameter Estimation in Geophysics) is a python package for simulation and gradient based parameter estimation in the context of geophysical applications.
+SimPEG (Simulation and Parameter Estimation in Geoscience) is a python
+package for simulation and gradient based parameter estimation in the
+context of geoscience 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:
+The vision is to create a package for finite volume simulation and parameter estimation with
+applications to geophysical imaging and subsurface flow. To enable
+these goals, this package has the following features:
 
-* modular with respect to the spacial discretization
-* is built with the inverse problem in mind
-* supports different hyperbolic solvers (Euler, Semi-Lagrangian, Lagrangian)
-* supports 2D and 3D problems
+* is modular with respect to discretization, physics, optimization, and regularization
+* is built with the (large-scale) inverse problem in mind
+* 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>
 
 
 Meshing & Operators

notebooks/3DRenderingWithvtkTools.ipynb

@@ -1,70 +0,0 @@
-{
- "metadata": {
-  "name": ""
- },
- "nbformat": 3,
- "nbformat_minor": 0,
- "worksheets": [
-  {
-   "cells": [
-    {
-     "cell_type": "code",
-     "collapsed": false,
-     "input": [
-      "import numpy as np, vtk\n",
-      "import SimPEG as simpeg"
-     ],
-     "language": "python",
-     "metadata": {},
-     "outputs": [],
-     "prompt_number": 5
-    },
-    {
-     "cell_type": "code",
-     "collapsed": false,
-     "input": [
-      "#Make a mesh and model\n",
-      "x0 = np.zeros(3)\n",
-      "h1 = np.ones(20)*50\n",
-      "h2 = np.ones(10)*100\n",
-      "h3 = np.ones(5)*200\n",
-      "\n",
-      "mesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n"
-     ],
-     "language": "python",
-     "metadata": {},
-     "outputs": [],
-     "prompt_number": 6
-    },
-    {
-     "cell_type": "code",
-     "collapsed": false,
-     "input": [
-      "# Make a models that correspond to the cells, faces and edges.\n",
-      "models = {'cell':{'Test':np.arange(0,mesh.nC),'AllOnce':np.ones(mesh.nC)},'face':{'Test':np.arange(0,np.sum(mesh.nF)),'AllOnce':np.ones(np.sum(mesh.nF))},'edge':{'Test':np.arange(0,np.sum(mesh.nE)),'AllOnce':np.ones(np.sum(mesh.nE))}}\n",
-      "# Make the vtk viewer object.\n",
-      "vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models)                            \n",
-      "                                                            "
-     ],
-     "language": "python",
-     "metadata": {},
-     "outputs": [],
-     "prompt_number": 7
-    },
-    {
-     "cell_type": "code",
-     "collapsed": false,
-     "input": [
-      "# Show the image \n",
-      "vtkViewer.Show()\n"
-     ],
-     "language": "python",
-     "metadata": {},
-     "outputs": [],
-     "prompt_number": "*"
-    }
-   ],
-   "metadata": {}
-  }
- ]
-}

notebooks/VisualizeWithvtkView-updated.ipynb

@@ -0,0 +1,142 @@
+{
+ "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": {}
+  }
+ ]
+}