Merge branch 'develop' into visulization

Conflicts:
	notebooks/3DRenderingWithvtkTools.ipynb

SimPEG/GCEtools/gceStartup.txt

@@ -0,0 +1,11 @@
+# Check project status
+gcutil getproject --project=<ProjectName> --cache_flag_values
+
+# Start an instance
+gcutil addinstance <instanceName>
+
+# Log in 
+gcutil ssh <instanceName>
+
+# Shut down
+gcutil deleteinstance <instanceName>

SimPEG/GCEtools/startup.sh

@@ -0,0 +1,22 @@
+#! /bin/bash
+sudo aptitude -y update
+sudo aptitude -y upgrade
+sudo aptitude -y install gcc gfortran git libopenmpi-dev python-pip python-dev
+sudo aptitude -y install ipython python-scipy python-numpy python-nose python-pip python-matplotlib
+sudo aptitude -y install libmumps-ptscotch-4.10.0 libmumps-ptscotch-dev
+sudo aptitude -y install libblas-dev liblapack-dev
+
+sudo pip install mpi4py
+sudo pip install pymumps
+
+sudo pip install scipy --upgrade
+sudo pip install numpy --upgrade
+sudo pip install ipython --upgrade
+
+git clone https://bitbucket.org/rcockett/simpeg.git
+cd simpeg/SimPEG/
+python setup.py
+cd ~
+
+echo export PYTHONPATH=/home/$USER/simpeg/ >> .bashrc
+source .bashrc

SimPEG/__init__.py

@@ -2,5 +2,11 @@ import utils
 from utils import Solver
 import mesh
 import inverse
+import visualize
 import forward
 import regularization
+import examples
+
+import scipy.version as _v
+if _v.version < '0.13.0':
+    print 'Warning: upgrade your scipy to 0.13.0'

SimPEG/forward/DCProblem.py

@@ -1,11 +1,10 @@
 from SimPEG.mesh import TensorMesh
-from SimPEG.forward import Problem, SyntheticProblem, ModelTransforms
+from SimPEG.forward import Problem, ModelTransforms
 from SimPEG.tests import checkDerivative
 from SimPEG.utils import ModelBuilder, sdiag, mkvc
 from SimPEG import Solver
 import numpy as np
 import scipy.sparse as sp
-import scipy.sparse.linalg as linalg
 
 
 class DCProblem(ModelTransforms.LogModel, Problem):
@@ -202,23 +201,17 @@ if __name__ == '__main__':
     P = Q.T
 
     # Create some data
-    class syntheticDCProblem(DCProblem, SyntheticProblem):
-        pass
+    problem = DCProblem(mesh)
+    problem.P = P
+    problem.RHS = q
+    dobs, Wd = problem.createSyntheticData(mSynth, std=0.05)
 
-    synthetic = syntheticDCProblem(mesh);
-    synthetic.P = P
-    synthetic.RHS = q
-    dobs, Wd = synthetic.createData(mSynth, std=0.05)
-
-    u = synthetic.field(mSynth)
-    u = synthetic.reshapeFields(u)
+    u = problem.field(mSynth)
+    u = problem.reshapeFields(u)
     mesh.plotImage(u[:,10])
     # plt.show()
 
     # Now set up the problem to do some minimization
-    problem = DCProblem(mesh)
-    problem.P = P
-    problem.RHS = q
     problem.dobs = dobs
     problem.std = dobs*0 + 0.05
     m0 = mesh.gridCC[:,0]*0+sig2

SimPEG/forward/Problem.py

@@ -1,5 +1,5 @@
 import numpy as np
-from SimPEG.utils import mkvc, sdiag
+from SimPEG.utils import mkvc, sdiag, count, timeIt
 import scipy.sparse as sp
 norm = np.linalg.norm
 
@@ -37,6 +37,8 @@ class Problem(object):
         to (locally) find how model parameters change the data, and optimize!
     """
 
+    counter = None
+
     def __init__(self, mesh):
         self.mesh = mesh
 
@@ -83,6 +85,7 @@ class Problem(object):
     def dobs(self, value):
         self._dobs = value
 
+    @count
     def dpred(self, m, u=None):
         """
             Predicted data.
@@ -94,6 +97,7 @@ class Problem(object):
             u = self.field(m)
         return self.P*u
 
+    @count
     def dataResidual(self, m, u=None):
         """
             :param numpy.array m: geophysical model
@@ -113,6 +117,7 @@ class Problem(object):
 
         return self.dpred(m, u=u) - self.dobs
 
+    @timeIt
     def J(self, m, v, u=None):
         """
             :param numpy.array m: model
@@ -142,6 +147,7 @@ class Problem(object):
         """
         raise NotImplementedError('J is not yet implemented.')
 
+    @timeIt
     def Jt(self, m, v, u=None):
         """
             :param numpy.array m: model
@@ -155,6 +161,7 @@ class Problem(object):
         raise NotImplementedError('Jt is not yet implemented.')
 
 
+    @timeIt
     def J_approx(self, m, v, u=None):
         """
 
@@ -169,6 +176,7 @@ class Problem(object):
         """
         return self.J(m, v, u)
 
+    @timeIt
     def Jt_approx(self, m, v, u=None):
         """
             :param numpy.array m: model
@@ -218,32 +226,19 @@ class Problem(object):
         """
         return sp.eye(m.size)
 
-
-
-
-class SyntheticProblem(object):
-    """
-        Has helpful functions when dealing with synthetic problems
-
-        To use this class, inherit to your problem::
-
-            class mySyntheticExample(Problem, SyntheticProblem):
-                pass
-    """
-    def createData(self, m, std=0.05):
+    def createSyntheticData(self, m, std=0.05, u=None):
         """
+            Create synthetic data given a model, and a standard deviation.
+
             :param numpy.array m: geophysical model
             :param numpy.array std: standard deviation
             :rtype: numpy.array, numpy.array
             :return: dobs, Wd
 
-            Create synthetic data given a model, and a standard deviation.
-
             Returns the observed data with random Gaussian noise
             and Wd which is the same size as data, and can be used to weight the inversion.
         """
-        dobs = self.dpred(m)
-        dobs = dobs
+        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

SimPEG/inverse/Inversion.py

@@ -1,17 +1,22 @@
 import numpy as np
 import scipy.sparse as sp
 import SimPEG
-from SimPEG.utils import sdiag, mkvc, setKwargs, checkStoppers, printStoppers
+from SimPEG.utils import sdiag, mkvc, setKwargs, checkStoppers, printStoppers, count, timeIt, callHooks
 from Optimize import Remember
 from BetaSchedule import Cooling
+from SimPEG.inverse import IterationPrinters, StoppingCriteria
 
 class BaseInversion(object):
     """docstring for BaseInversion"""
 
-    maxIter = 1
+    maxIter = 1        #: Maximum number of iterations
     name = 'BaseInversion'
-    debug = False
-    beta0 = 1e4
+
+    debug   = False    #: Print debugging information
+
+    comment = ''       #: Used by some functions to indicate what is going on in the algorithm
+    counter = None     #: Set this to a SimPEG.utils.Counter() if you want to count things
+
 
     def __init__(self, prob, reg, opt, **kwargs):
         setKwargs(self, **kwargs)
@@ -20,14 +25,18 @@ class BaseInversion(object):
         self.opt = opt
         self.opt.parent = self
 
-        self.stoppers = [SimPEG.inverse.StoppingCriteria.iteration, SimPEG.inverse.StoppingCriteria.phi_d_target_Inversion]
+        self.stoppers = [StoppingCriteria.iteration]
 
         # Check if we have inserted printers into the optimization
-        if not np.any([p is SimPEG.inverse.IterationPrinters.phi_d for p in self.opt.printers]):
-            self.opt.printers.insert(1,SimPEG.inverse.IterationPrinters.beta)
-            self.opt.printers.insert(2,SimPEG.inverse.IterationPrinters.phi_d)
-            self.opt.printers.insert(3,SimPEG.inverse.IterationPrinters.phi_m)
-            self.opt.stoppers.append(SimPEG.inverse.StoppingCriteria.phi_d_target_Minimize)
+        if IterationPrinters.phi_d not in self.opt.printers:
+            self.opt.printers.insert(1,IterationPrinters.beta)
+            self.opt.printers.insert(2,IterationPrinters.phi_d)
+            self.opt.printers.insert(3,IterationPrinters.phi_m)
+
+        if not hasattr(opt, '_bfgsH0') and hasattr(opt, 'bfgsH0'): # Check if it has been set by the user and the default is not being used.
+            print 'Setting bfgsH0 to the inverse of the modelObj2Deriv. Done using direct methods.'
+            opt.bfgsH0 = SimPEG.Solver(reg.modelObj2Deriv())
+
 
     @property
     def Wd(self):
@@ -38,6 +47,9 @@ class BaseInversion(object):
             eps = np.linalg.norm(mkvc(self.prob.dobs),2)*1e-5
             self._Wd = 1/(abs(self.prob.dobs)*self.prob.std+eps)
         return self._Wd
+    @Wd.setter
+    def Wd(self, value):
+        self._Wd = value
 
     @property
     def phi_d_target(self):
@@ -56,7 +68,13 @@ class BaseInversion(object):
     def phi_d_target(self, value):
         self._phi_d_target = value
 
+    @timeIt
     def run(self, m0):
+        """run(m0)
+
+            Runs the inversion!
+
+        """
         self.startup(m0)
         while True:
             self._beta = self.getBeta()
@@ -83,13 +101,17 @@ class BaseInversion(object):
             :rtype: None
             :return: None
         """
-        for method in [posible for posible in dir(self) if '_startup' in posible]:
-            if self.debug: print 'startup is calling self.'+method
-            getattr(self,method)(m0)
+        callHooks(self,'startup',m0)
+
+        if not hasattr(self.reg, '_mref'):
+            print 'Regularization has not set mref. SimPEG will set it to m0.'
+            self.reg.mref = m0
 
         self.m = m0
         self._iter = 0
         self._beta = None
+        self.phi_d_last = np.nan
+        self.phi_m_last = np.nan
 
     def doEndIteration(self):
         """
@@ -97,28 +119,75 @@ class BaseInversion(object):
 
             If you have things that also need to run at the end of every iteration, you can create a method::
 
-                def _doEndIteration*(self, xt):
+                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.
 
-            :param numpy.ndarray xt: tested new iterate that ensures a descent direction.
             :rtype: None
             :return: None
         """
-        for method in [posible for posible in dir(self) if '_doEndIteration' in posible]:
-            if self.debug: print 'doEndIteration is calling self.'+method
-            getattr(self,method)()
+        callHooks(self,'doEndIteration')
 
         # store old values
         self.phi_d_last = self.phi_d
         self.phi_m_last = self.phi_m
         self._iter += 1
 
+    @property
+    def beta0(self):
+        if getattr(self,'_beta0',None) is None:
+            self._beta0 = self.estimateBeta0()
+        return self._beta0
+    @beta0.setter
+    def beta0(self, value):
+        self._beta0 = value
+
     def getBeta(self):
         return self.beta0
 
+    def estimateBeta0(self, u=None, ratio=0.1):
+        """estimateBeta0(u=None, ratio=0.1)
+
+            The initial beta is calculated by comparing the estimated
+            eigenvalues of JtJ and WtW.
+
+            To estimate the eigenvector of **A**, we will use one iteration
+            of the *Power Method*:
+
+            .. math::
+
+                \mathbf{x_1 = A x_0}
+
+            Given this (very course) approximation of the eigenvector,
+            we can use the *Rayleigh quotient* to approximate the largest eigenvalue.
+
+            .. math::
+
+                \lambda_0 = \\frac{\mathbf{x^\\top A x}}{\mathbf{x^\\top x}}
+
+            We will approximate the largest eigenvalue for both JtJ and WtW, and
+            use some ratio of the quotient to estimate beta0.
+
+            .. math::
+
+                \\beta_0 = \gamma \\frac{\mathbf{x^\\top J^\\top J x}}{\mathbf{x^\\top W^\\top W x}}
+
+
+            :param numpy.array u: fields
+            :param float ratio: desired ratio of the eigenvalues, default is 0.1
+            :rtype: float
+            :return: beta0
+        """
+        if u is None:
+            u = self.prob.field(self.m)
+
+        x0 = np.random.rand(*self.m.shape)
+        t = x0.dot(self.dataObj2Deriv(self.m,x0,u=u))
+        b = x0.dot(self.reg.modelObj2Deriv()*x0)
+        return ratio*(t/b)
+
     def stoppingCriteria(self):
         if self.debug: print 'checking stoppingCriteria'
         return checkStoppers(self, self.stoppers)
@@ -131,8 +200,12 @@ class BaseInversion(object):
         """
         printStoppers(self, self.stoppers)
 
-
+    @timeIt
     def evalFunction(self, m, return_g=True, return_H=True):
+        """evalFunction(m, return_g=True, return_H=True)
+
+
+        """
 
         u = self.prob.field(m)
         phi_d = self.dataObj(m, u)
@@ -154,17 +227,18 @@ class BaseInversion(object):
         if return_H:
             def H_fun(v):
                 phi_d2Deriv = self.dataObj2Deriv(m, v, u=u)
-                phi_m2Deriv = self.reg.modelObj2Deriv(m)*v
+                phi_m2Deriv = self.reg.modelObj2Deriv()*v
 
                 return phi_d2Deriv + self._beta * phi_m2Deriv
 
-            operator = sp.linalg.LinearOperator( (m.size, m.size), H_fun, dtype=float )
+            operator = sp.linalg.LinearOperator( (m.size, m.size), H_fun, dtype=m.dtype )
             out += (operator,)
         return out if len(out) > 1 else out[0]
 
-
+    @timeIt
     def dataObj(self, m, u=None):
-        """
+        """dataObj(m, u=None)
+
             :param numpy.array m: geophysical model
             :param numpy.array u: fields
             :rtype: float
@@ -184,8 +258,10 @@ class BaseInversion(object):
         R = mkvc(R)
         return 0.5*np.vdot(R, R)
 
+    @timeIt
     def dataObjDeriv(self, m, u=None):
-        """
+        """dataObjDeriv(m, u=None)
+
             :param numpy.array m: geophysical model
             :param numpy.array u: fields
             :rtype: numpy.array
@@ -224,9 +300,12 @@ class BaseInversion(object):
 
         return dmisfit
 
+    @timeIt
     def dataObj2Deriv(self, m, v, u=None):
-        """
+        """dataObj2Deriv(m, v, u=None)
+
             :param numpy.array m: geophysical model
+            :param numpy.array v: vector to multiply
             :param numpy.array u: fields
             :rtype: numpy.array
             :return: data misfit derivative
@@ -263,7 +342,7 @@ class BaseInversion(object):
         R = self.Wd*self.prob.dataResidual(m, u=u)
 
         # TODO: abstract to different norms a little cleaner.
-        #                                 \/ it goes here. in l2 it is the identity.
+        #                                        \/ it goes here. in l2 it is the identity.
         dmisfit = self.prob.Jt_approx(m, self.Wd * self.Wd * self.prob.J_approx(m, v, u=u), u=u)
 
         return dmisfit
@@ -275,3 +354,33 @@ class Inversion(Cooling, Remember, BaseInversion):
 
     def __init__(self, prob, reg, opt, **kwargs):
         BaseInversion.__init__(self, prob, reg, opt, **kwargs)
+
+        self.stoppers.append(StoppingCriteria.phi_d_target_Inversion)
+
+        if StoppingCriteria.phi_d_target_Minimize not in self.opt.stoppers:
+            self.opt.stoppers.append(StoppingCriteria.phi_d_target_Minimize)
+
+class TimeSteppingInversion(Remember, BaseInversion):
+    """
+        A slightly different view on regularization parameters,
+        let Beta be viewed as 1/dt, and timestep by updating the
+        reference model every optimization iteration.
+    """
+    maxIter = 1
+    name = "Time-Stepping SimPEG Inversion"
+
+    def __init__(self, prob, reg, opt, **kwargs):
+        BaseInversion.__init__(self, prob, reg, opt, **kwargs)
+
+        self.stoppers.append(StoppingCriteria.phi_d_target_Inversion)
+
+        if StoppingCriteria.phi_d_target_Minimize not in self.opt.stoppers:
+            self.opt.stoppers.append(StoppingCriteria.phi_d_target_Minimize)
+
+    def _startup_TimeSteppingInversion(self, m0):
+
+        def _doEndIteration_updateMref(self, xt):
+            if self.debug: 'Updating the reference model.'
+            self.parent.reg.mref = self.xc
+
+        self.opt.hook(_doEndIteration_updateMref, overwrite=True)

SimPEG/inverse/Optimize.py

@@ -1,16 +1,13 @@
 import numpy as np
 import matplotlib.pyplot as plt
-from SimPEG.utils import mkvc, sdiag, setKwargs, printTitles, printLine, printStoppers, checkStoppers
+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
 
-try:
-    from pubsub import pub
-    doPub = True
-except Exception, e:
-    print 'Warning: you may not have the required pubsub installed, use pypubsub. You will not be able to listen to events.'
-    doPub = False
+
+__all__ = ['Minimize', 'Remember', 'SteepestDescent', 'BFGS', 'GaussNewton', 'InexactGaussNewton', 'ProjectedGradient', 'NewtonRoot', 'StoppingCriteria', 'IterationPrinters']
+
 
 class StoppingCriteria(object):
     """docstring for StoppingCriteria"""
@@ -76,7 +73,7 @@ class IterationPrinters(object):
     itType = {"title": "itType", "value": lambda M: M._itType, "width": 8, "format": "%s"}
     aSet = {"title": "aSet", "value": lambda M: np.sum(M.activeSet(M.xc)), "width": 8, "format": "%d"}
     bSet = {"title": "bSet", "value": lambda M: np.sum(M.bindingSet(M.xc)), "width": 8, "format": "%d"}
-    comment = {"title": "Comment", "value": lambda M: M.projComment, "width": 7, "format": "%s"}
+    comment = {"title": "Comment", "value": lambda M: M.comment, "width": 12, "format": "%s"}
 
     beta = {"title": "beta", "value": lambda M: M.parent._beta, "width": 10, "format":   "%1.2e"}
     phi_d = {"title": "phi_d", "value": lambda M: M.parent.phi_d, "width": 10, "format":   "%1.2e"}
@@ -85,29 +82,28 @@ class IterationPrinters(object):
 
 class Minimize(object):
     """
-
         Minimize is a general class for derivative based optimization.
-
-
     """
 
-    name = "General Optimization Algorithm"
+    name = "General Optimization Algorithm"  #: The name of the optimization algorithm
 
-    maxIter = 20
-    maxIterLS = 10
-    maxStep = np.inf
-    LSreduction = 1e-4
-    LSshorten = 0.5
-    tolF = 1e-1
-    tolX = 1e-1
-    tolG = 1e-1
-    eps = 1e-5
+    maxIter = 20       #: Maximum number of iterations
+    maxIterLS = 10     #: Maximum number of iterations for the line-search
+    maxStep = np.inf   #: Maximum step possible, used in scaling before the line-search.
+    LSreduction = 1e-4 #: Expected decrease in the line-search
+    LSshorten = 0.5    #: Line-search step is shortened by this amount each time.
+    tolF = 1e-1        #: Tolerance on function value decrease
+    tolX = 1e-1        #: Tolerance on norm(x) movement
+    tolG = 1e-1        #: Tolerance on gradient norm
+    eps = 1e-5         #: Small value
 
-    debug   = False
-    debugLS = False
+    debug   = False    #: Print debugging information
+    debugLS = False    #: Print debugging information for the line-search
+
+    comment = ''       #: Used by some functions to indicate what is going on in the algorithm
+    counter = None     #: Set this to a SimPEG.utils.Counter() if you want to count things
 
     def __init__(self, **kwargs):
-        self._id = int(np.random.rand()*1e6)  # create a unique identifier to this program to be used in pubsub
         self.stoppers = [StoppingCriteria.tolerance_f, StoppingCriteria.moving_x, StoppingCriteria.tolerance_g, StoppingCriteria.norm_g, StoppingCriteria.iteration]
         self.stoppersLS = [StoppingCriteria.armijoGoldstein, StoppingCriteria.iterationLS]
 
@@ -116,8 +112,10 @@ class Minimize(object):
 
         setKwargs(self, **kwargs)
 
+    @timeIt
     def minimize(self, evalFunction, x0):
-        """
+        """minimize(evalFunction, x0)
+
         Minimizes the function (evalFunction) starting at the location x0.
 
         :param def evalFunction: function handle that evaluates: f, g, H = F(x)
@@ -138,27 +136,6 @@ class Minimize(object):
                 return out if len(out) > 1 else out[0]
 
 
-        Events are fired with the following inputs via pypubsub::
-
-            Minimize.printInit              (minimize)
-            Minimize.evalFunction           (minimize, f, g, H)
-            Minimize.printIter              (minimize)
-            Minimize.searchDirection        (minimize, p)
-            Minimize.scaleSearchDirection   (minimize, p)
-            Minimize.modifySearchDirection  (minimize, xt, passLS)
-            Minimize.endIteration           (minimize, xt)
-            Minimize.printDone              (minimize)
-
-        To hook into one of these events (must have pypubsub installed)::
-
-            from pubsub import pub
-            def listener(minimize,p):
-                print 'The search direction is:  ', p
-            pub.subscribe(listener, 'Minimize.searchDirection')
-
-        You can use pubsub communication to debug your code, it is not used internally.
-
-
         The algorithm for general minimization is as follows::
 
             startup(x0)
@@ -185,20 +162,15 @@ class Minimize(object):
 
         while True:
             self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True)
-            if doPub: pub.sendMessage('Minimize.evalFunction', minimize=self, f=self.f, g=self.g, H=self.H)
             self.printIter()
             if self.stoppingCriteria(): break
             p = self.findSearchDirection()
-            if doPub: pub.sendMessage('Minimize.searchDirection', minimize=self, p=p)
             p = self.scaleSearchDirection(p)
-            if doPub: pub.sendMessage('Minimize.scaleSearchDirection', minimize=self, p=p)
             xt, passLS = self.modifySearchDirection(p)
-            if doPub: pub.sendMessage('Minimize.modifySearchDirection', minimize=self, xt=xt, passLS=passLS)
             if not passLS:
                 xt, caught = self.modifySearchDirectionBreak(p)
                 if not caught: return self.xc
             self.doEndIteration(xt)
-            if doPub: pub.sendMessage('Minimize.endIteration', minimize=self, xt=xt)
 
         self.printDone()
 
@@ -236,9 +208,7 @@ class Minimize(object):
             :rtype: None
             :return: None
         """
-        for method in [posible for posible in dir(self) if '_startup' in posible]:
-            if self.debug: print 'startup is calling self.'+method
-            getattr(self,method)(x0)
+        callHooks(self,'startup',x0)
 
         self._iter = 0
         self._iterLS = 0
@@ -258,7 +228,6 @@ class Minimize(object):
             parent.printInit function and call that.
 
         """
-        if doPub and not inLS: pub.sendMessage('Minimize.printInit', minimize=self)
         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)
@@ -271,7 +240,8 @@ class Minimize(object):
             parent.printIter function and call that.
 
         """
-        if doPub and not inLS: pub.sendMessage('Minimize.printIter', minimize=self)
+        callHooks(self,'printIter',inLS)
+
         pad = ' '*10 if inLS else ''
         printLine(self, self.printers if not inLS else self.printersLS, pad=pad)
 
@@ -283,22 +253,21 @@ class Minimize(object):
             parent.printDone function and call that.
 
         """
-        if doPub and not inLS: pub.sendMessage('Minimize.printDone', minimize=self)
         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)
 
-
     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)
 
-
+    @timeIt
     def projection(self, p):
-        """
+        """projection(p)
+
             projects the search direction.
 
             by default, no projection is applied.
@@ -307,10 +276,13 @@ class Minimize(object):
             :rtype: numpy.ndarray
             :return: p, projected search direction
         """
+        callHooks(self,'projection',p)
         return p
 
+    @timeIt
     def findSearchDirection(self):
-        """
+        """findSearchDirection()
+
             **findSearchDirection** should return an approximation of:
 
             .. math::
@@ -338,8 +310,10 @@ class Minimize(object):
         """
         return -self.g
 
+    @count
     def scaleSearchDirection(self, p):
-        """
+        """scaleSearchDirection(p)
+
             **scaleSearchDirection** should scale the search direction if appropriate.
 
             Set the parameter **maxStep** in the minimize object, to scale back the gradient to a maximum size.
@@ -353,10 +327,12 @@ class Minimize(object):
             p = self.maxStep*p/np.abs(p.max())
         return p
 
-    nameLS = "Armijo linesearch"
+    nameLS = "Armijo linesearch" #: The line-search name
 
+    @timeIt
     def modifySearchDirection(self, p):
-        """
+        """modifySearchDirection(p)
+
             **modifySearchDirection** changes the search direction based on some sort of linesearch or trust-region criteria.
 
             By default, an Armijo backtracking linesearch is preformed with the following parameters:
@@ -391,8 +367,10 @@ class Minimize(object):
 
         return self._LS_xt, self._iterLS < self.maxIterLS
 
+    @count
     def modifySearchDirectionBreak(self, p):
-        """
+        """modifySearchDirectionBreak(p)
+
             Code is called if modifySearchDirection fails
             to find a descent direction.
 
@@ -411,8 +389,10 @@ class Minimize(object):
         print 'The linesearch got broken. Boo.'
         return p, False
 
+    @count
     def doEndIteration(self, xt):
-        """
+        """doEndIteration(xt)
+
             **doEndIteration** is called at the end of each minimize iteration.
 
             By default, function values and x locations are shuffled to store 1 past iteration in memory.
@@ -432,9 +412,7 @@ class Minimize(object):
             :rtype: None
             :return: None
         """
-        for method in [posible for posible in dir(self) if '_doEndIteration' in posible]:
-            if self.debug: print 'doEndIteration is calling self.'+method
-            getattr(self,method)(xt)
+        callHooks(self,'doEndIteration',xt)
 
         # store old values
         self.f_last = self.f
@@ -443,7 +421,6 @@ class Minimize(object):
         if self.debug: self.printDone()
 
 
-
 class Remember(object):
     """
         This mixin remembers all the things you tend to forget.
@@ -494,12 +471,11 @@ class Remember(object):
                 self._rememberList[param[0]].append( param[1](self) )
 
 
-
 class ProjectedGradient(Minimize, Remember):
     name = 'Projected Gradient'
 
-    maxIterCG = 10
-    tolCG = 1e-3
+    maxIterCG = 5
+    tolCG = 1e-1
 
     lower = -np.inf
     upper = np.inf
@@ -525,24 +501,41 @@ class ProjectedGradient(Minimize, Remember):
         self.stopDoingPG = False
 
         self._itType = 'SD'
-        self.projComment = ''
+        self.comment = ''
 
         self.aSet_prev = self.activeSet(x0)
 
+    @count
     def projection(self, x):
-        """Make sure we are feasible."""
+        """projection(x)
+
+            Make sure we are feasible.
+
+        """
         return np.median(np.c_[self.lower,x,self.upper],axis=1)
 
+    @count
     def activeSet(self, x):
-        """If we are on a bound"""
+        """activeSet(x)
+
+            If we are on a bound
+
+        """
         return np.logical_or(x == self.lower, x == self.upper)
 
+    @count
     def inactiveSet(self, x):
-        """The free variables."""
+        """inactiveSet(x)
+
+            The free variables.
+
+        """
         return np.logical_not(self.activeSet(x))
 
+    @count
     def bindingSet(self, x):
-        """
+        """bindingSet(x)
+
             If we are on a bound and the negative gradient points away from the feasible set.
 
             Optimality condition. (Satisfies Kuhn-Tucker) MoreToraldo91
@@ -552,7 +545,12 @@ class ProjectedGradient(Minimize, Remember):
         bind_low = np.logical_and(x == self.upper, self.g <= 0)
         return np.logical_or(bind_up, bind_low)
 
+    @timeIt
     def findSearchDirection(self):
+        """findSearchDirection()
+
+            Finds the search direction based on either CG or steepest descent.
+        """
         self.aSet_prev = self.activeSet(self.xc)
         allBoundsAreActive = sum(self.aSet_prev) == self.xc.size
 
@@ -586,13 +584,15 @@ class ProjectedGradient(Minimize, Remember):
             def reduceHess(v):
                 # Z is tall and skinny
                 return Z.T*(self.H*(Z*v))
-            operator = sp.linalg.LinearOperator( (shape[1], shape[1]), reduceHess, dtype=float )
+            operator = sp.linalg.LinearOperator( (shape[1], shape[1]), reduceHess, dtype=self.xc.dtype )
             p, info  = sp.linalg.cg(operator, -Z.T*self.g, tol=self.tolCG, maxiter=self.maxIterCG)
             p = Z*p  #  bring up to full size
             # aSet_after = self.activeSet(self.xc+p)
         return p
 
+    @timeIt
     def _doEndIteration_ProjectedGradient(self, xt):
+        """_doEndIteration_ProjectedGradient(xt)"""
         aSet = self.activeSet(xt)
         bSet = self.bindingSet(xt)
 
@@ -600,7 +600,7 @@ class ProjectedGradient(Minimize, Remember):
         self.exploreCG = np.all(aSet == bSet) # explore conjugate gradient
 
         f_current_decrease = self.f_last - self.f
-        self.projComment = ''
+        self.comment = ''
         if self._iter < 1:
             # Note that this is reset on every CG iteration.
             self.f_decrease_max = -np.inf
@@ -608,7 +608,7 @@ class ProjectedGradient(Minimize, Remember):
             self.f_decrease_max = max(self.f_decrease_max, f_current_decrease)
             self.stopDoingPG = f_current_decrease < 0.25 * self.f_decrease_max
             if self.stopDoingPG:
-                self.projComment = 'Stop SD'
+                self.comment = 'Stop SD'
                 self.explorePG = False
                 self.exploreCG = True
         # implement 3.8, MoreToraldo91
@@ -620,40 +620,229 @@ class ProjectedGradient(Minimize, Remember):
         if self.debug: print 'doEndIteration.ProjGrad, f_decrease_max: ', self.f_decrease_max
         if self.debug: print 'doEndIteration.ProjGrad, stopDoingSD: ', self.stopDoingSD
 
+
+class BFGS(Minimize, Remember):
+    name  = 'BFGS'
+    nbfgs = 10
+
+    @property
+    def bfgsH0(self):
+        """
+            Approximate Hessian used in preconditioning the problem.
+
+            Must be a SimPEG.Solver
+        """
+        _bfgsH0 = getattr(self,'_bfgsH0',None)
+        if _bfgsH0 is None:
+            return Solver(sp.identity(self.xc.size).tocsc(), flag='D')
+        return _bfgsH0
+    @bfgsH0.setter
+    def bfgsH0(self, value):
+        assert type(value) is Solver, 'bfgsH0 must be a SimPEG.Solver'
+        self._bfgsH0 = value
+
+    def _startup_BFGS(self,x0):
+        self._bfgscnt = -1
+        self._bfgsY   = np.zeros((x0.size, self.nbfgs))
+        self._bfgsS   = np.zeros((x0.size, self.nbfgs))
+        if not np.any([p is IterationPrinters.comment for p in self.printers]):
+            self.printers.append(IterationPrinters.comment)
+
+    def bfgs(self, d):
+        n  = self._bfgscnt
+        nn = ktop = min(self._bfgsS.shape[1],n)
+        return self.bfgsrec(ktop,n,nn,self._bfgsS,self._bfgsY,d)
+
+    def bfgsrec(self,k,n,nn,S,Y,d):
+        """BFGS recursion"""
+        if k < 0:
+            d = self.bfgsH0.solve(d)
+        else:
+            khat    = 0 if nn is 0 else np.mod(n-nn+k,nn)
+            gamma   = np.vdot(S[:,khat],d)/np.vdot(Y[:,khat],S[:,khat])
+            d       = d - gamma*Y[:,khat]
+            d       = self.bfgsrec(k-1,n,nn,S,Y,d)
+            d       = d + (gamma - np.vdot(Y[:,khat],d)/np.vdot(Y[:,khat],S[:,khat]))*S[:,khat]
+        return d
+
+    def findSearchDirection(self):
+        return self.bfgs(-self.g)
+
+    def _doEndIteration_BFGS(self, xt):
+        if self._iter is 0:
+            self.g_last = self.g
+            return
+
+        yy = self.g - self.g_last;
+        ss = self.xc - xt;
+        self.g_last = self.g
+
+        if yy.dot(ss) > 0:
+            self._bfgscnt += 1
+            ktop = np.mod(self._bfgscnt,self.nbfgs)
+            self._bfgsY[:,ktop] = yy
+            self._bfgsS[:,ktop] = ss
+            self.comment = ''
+        else:
+            self.comment = 'Skip BFGS'
+
+
 class GaussNewton(Minimize, Remember):
     name = 'Gauss Newton'
+
+    @timeIt
     def findSearchDirection(self):
         return Solver(self.H).solve(-self.g)
 
 
-class InexactGaussNewton(Minimize, Remember):
+class InexactGaussNewton(BFGS, Minimize, Remember):
+    """
+        Minimizes using CG as the inexact solver of
+
+        .. math::
+
+            \mathbf{H p = -g}
+
+        By default BFGS is used as the preconditioner.
+
+        Use *nbfgs* to set the memory limitation of BFGS.
+
+        To set the initial H0 to be used in BFGS, set *bfgsH0* to be a SimPEG.Solver
+
+    """
+
+    def __init__(self, **kwargs):
+        Minimize.__init__(self, **kwargs)
+
     name = 'Inexact Gauss Newton'
 
-    maxIterCG = 10
-    tolCG = 1e-5
+    maxIterCG = 5
+    tolCG = 1e-1
 
+    @property
+    def approxHinv(self):
+        """
+            The approximate Hessian inverse is used to precondition CG.
+
+            Default uses BFGS, with an initial H0 of *bfgsH0*.
+
+            Must be a scipy.sparse.linalg.LinearOperator
+        """
+        _approxHinv = getattr(self,'_approxHinv',None)
+        if _approxHinv is None:
+            M = sp.linalg.LinearOperator( (self.xc.size, self.xc.size), self.bfgs, dtype=self.xc.dtype )
+            return M
+        return _approxHinv
+    @approxHinv.setter
+    def approxHinv(self, value):
+        self._approxHinv = value
+
+    @timeIt
     def findSearchDirection(self):
-        # TODO: use BFGS as a preconditioner or gauss sidel of the WtW or solve WtW directly
-        p, info = sp.linalg.cg(self.H, -self.g, tol=self.tolCG, maxiter=self.maxIterCG)
+        Hinv = Solver(self.H, doDirect=False, options={'iterSolver': 'CG', 'M': self.approxHinv, 'tol': self.tolCG, 'maxIter': self.maxIterCG})
+        p = Hinv.solve(-self.g)
         return p
 
 
 class SteepestDescent(Minimize, Remember):
     name = 'Steepest Descent'
+
+    @timeIt
     def findSearchDirection(self):
         return -self.g
 
+
+class NewtonRoot(object):
+    """
+        Newton Method - Root Finding
+
+        root = newtonRoot(fun,x);
+
+        Where fun is the function that returns the function value as well as the
+        gradient.
+
+        For iterative solving of dh = -J\\r, use O.solveTol = TOL. For direct
+        solves, use SOLVETOL = 0 (default)
+
+        Rowan Cockett
+        16-May-2013 16:29:51
+        University of British Columbia
+        rcockett@eos.ubc.ca
+
+    """
+
+    tol      = 1.000e-06
+    solveTol = 0 # Default direct solve.
+    maxIter  = 20
+    stepDcr  = 0.5
+    maxLS    = 30
+    comments = False
+    doLS     = True
+
+    def __init__(self, **kwargs):
+        setKwargs(self, **kwargs)
+
+    def root(self, fun, x):
+        if self.comments: print 'Newton Method:\n'
+
+        self._iter = 0
+        while True:
+
+            [r,J] = fun(x);
+            if self.solveTol == 0:
+                Jinv = Solver(J)
+                dh   = - Jinv.solve(r)
+            else:
+                raise NotImplementedError('Iterative solve on NewtonRoot is not yet implemented.')
+                # M = @(x) tril(J)\(diag(J).*(triu(J)\x));
+                # [dh, ~] = bicgstab(J,-r,O.solveTol,500,M);
+
+            muLS = 1.
+            LScnt  = 1
+            xt = x + dh
+            rt, Jt = fun(xt) # TODO: get rid of Jt
+
+            if self.comments: print '\tLinesearch:\n'
+            # Enter Linesearch
+            while True and self.doLS:
+                if self.comments:
+                    print '\t\tResid: %e\n'%norm(rt)
+                if norm(rt) <= norm(r) or norm(rt) < self.tol:
+                    break
+
+                muLS = muLS*self.stepDcr
+                LScnt = LScnt + 1
+                print '.'
+                if LScnt > self.maxLS:
+                    print 'Newton Method: Line search break.'
+                    root = NaN
+                    return
+                xt = x + muLS*dh
+                rt, Jt = fun(xt) # TODO: get rid of Jt
+
+            x = xt
+            self._iter += 1
+            if norm(rt) < self.tol or self._iter > self.maxIter:
+                break
+
+        return x
+
+
+
 if __name__ == '__main__':
     from SimPEG.tests import Rosenbrock, checkDerivative
     import matplotlib.pyplot as plt
     x0 = np.array([2.6, 3.7])
     checkDerivative(Rosenbrock, x0, plotIt=False)
 
-    # def listener1(minimize,p):
-    #     print 'hi:  ', p
-    # if doPub: pub.subscribe(listener1, 'Minimize.searchDirection')
-
     xOpt = GaussNewton(maxIter=20,tolF=1e-10,tolX=1e-10,tolG=1e-10).minimize(Rosenbrock,x0)
     print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1])
     xOpt = SteepestDescent(maxIter=30, maxIterLS=15,tolF=1e-10,tolX=1e-10,tolG=1e-10).minimize(Rosenbrock, x0)
     print "xOpt=[%f, %f]" % (xOpt[0], xOpt[1])
+
+
+    print 'test the newtonRoot finding.'
+    fun = lambda x: (np.sin(x), sdiag(np.cos(x)))
+    x = np.array([np.pi-0.3, np.pi+0.1, 0])
+    pnt = NewtonRoot(comments=False).root(fun,x)
+    print pnt

SimPEG/mesh/BaseMesh.py

@@ -102,7 +102,7 @@ class BaseMesh(object):
             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
-                nn = self.nF if xType == 'F' else self.nE
+                nn = self.nFv if xType == 'F' else self.nEv
                 nn = np.r_[0, nn]
 
                 nx = [0, 0, 0]
@@ -228,6 +228,17 @@ class BaseMesh(object):
         return locals()
     nC = property(**nC())
 
+    def nCv():
+        doc = """
+        Total number of cells in each direction
+
+        :rtype: numpy.array (dim, )
+        :return: [nCx, nCy, nCz]
+        """
+        fget = lambda self: np.array([x for x in [self.nCx, self.nCy, self.nCz] if not x is None])
+        return locals()
+    nCv = property(**nCv())
+
     def nNx():
         doc = """
         Number of nodes in the x-direction
@@ -288,6 +299,17 @@ class BaseMesh(object):
         return locals()
     nN = property(**nN())
 
+    def nNv():
+        doc = """
+        Total number of nodes in each direction
+
+        :rtype: numpy.array (dim, )
+        :return: [nNx, nNy, nNz]
+        """
+        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())    
+
     def nEx():
         doc = """
         Number of x-edges in each direction
@@ -331,7 +353,7 @@ class BaseMesh(object):
         return locals()
     nEz = property(**nEz())
 
-    def nE():
+    def nEv():
         doc = """
         Total number of edges in each direction
 
@@ -346,6 +368,18 @@ class BaseMesh(object):
         """
         fget = lambda self: np.array([np.prod(x) for x in [self.nEx, self.nEy, self.nEz] if not x is None])
         return locals()
+    nEv = property(**nEv())
+
+    def nE():
+        doc = """
+        Total number of edges.
+
+        :rtype: int
+        :return: sum([prod(nEx), prod(nEy), prod(nEz)])
+
+        """
+        fget = lambda self: np.sum(self.nEv)
+        return locals()
     nE = property(**nE())
 
     def nFx():
@@ -391,7 +425,7 @@ class BaseMesh(object):
         return locals()
     nFz = property(**nFz())
 
-    def nF():
+    def nFv():
         doc = """
         Total number of faces in each direction
 
@@ -406,6 +440,19 @@ class BaseMesh(object):
         """
         fget = lambda self: np.array([np.prod(x) for x in [self.nFx, self.nFy, self.nFz] if not x is None])
         return locals()
+    nFv = property(**nFv())
+
+
+    def nF():
+        doc = """
+        Total number of faces.
+
+        :rtype: int
+        :return: sum([prod(nFx), prod(nFy), prod(nFz)])
+
+        """
+        fget = lambda self: np.sum(self.nFv)
+        return locals()
     nF = property(**nF())
 
     def normals():
@@ -418,13 +465,13 @@ class BaseMesh(object):
 
         def fget(self):
             if self.dim == 2:
-                nX = np.c_[np.ones(self.nF[0]), np.zeros(self.nF[0])]
-                nY = np.c_[np.zeros(self.nF[1]), np.ones(self.nF[1])]
+                nX = np.c_[np.ones(self.nFv[0]), np.zeros(self.nFv[0])]
+                nY = np.c_[np.zeros(self.nFv[1]), np.ones(self.nFv[1])]
                 return np.r_[nX, nY]
             elif self.dim == 3:
-                nX = np.c_[np.ones(self.nF[0]), np.zeros(self.nF[0]), np.zeros(self.nF[0])]
-                nY = np.c_[np.zeros(self.nF[1]), np.ones(self.nF[1]), np.zeros(self.nF[1])]
-                nZ = np.c_[np.zeros(self.nF[2]), np.zeros(self.nF[2]), np.ones(self.nF[2])]
+                nX = np.c_[np.ones(self.nFv[0]), np.zeros(self.nFv[0]), np.zeros(self.nFv[0])]
+                nY = np.c_[np.zeros(self.nFv[1]), np.ones(self.nFv[1]), np.zeros(self.nFv[1])]
+                nZ = np.c_[np.zeros(self.nFv[2]), np.zeros(self.nFv[2]), np.ones(self.nFv[2])]
                 return np.r_[nX, nY, nZ]
         return locals()
     normals = property(**normals())
@@ -439,13 +486,13 @@ class BaseMesh(object):
 
         def fget(self):
             if self.dim == 2:
-                tX = np.c_[np.ones(self.nE[0]), np.zeros(self.nE[0])]
-                tY = np.c_[np.zeros(self.nE[1]), np.ones(self.nE[1])]
+                tX = np.c_[np.ones(self.nEv[0]), np.zeros(self.nEv[0])]
+                tY = np.c_[np.zeros(self.nEv[1]), np.ones(self.nEv[1])]
                 return np.r_[tX, tY]
             elif self.dim == 3:
-                tX = np.c_[np.ones(self.nE[0]), np.zeros(self.nE[0]), np.zeros(self.nE[0])]
-                tY = np.c_[np.zeros(self.nE[1]), np.ones(self.nE[1]), np.zeros(self.nE[1])]
-                tZ = np.c_[np.zeros(self.nE[2]), np.zeros(self.nE[2]), np.ones(self.nE[2])]
+                tX = np.c_[np.ones(self.nEv[0]), np.zeros(self.nEv[0]), np.zeros(self.nEv[0])]
+                tY = np.c_[np.zeros(self.nEv[1]), np.ones(self.nEv[1]), np.zeros(self.nEv[1])]
+                tZ = np.c_[np.zeros(self.nEv[2]), np.zeros(self.nEv[2]), np.ones(self.nEv[2])]
                 return np.r_[tX, tY, tZ]
         return locals()
     tangents = property(**tangents())

SimPEG/mesh/Cyl1DMesh.py

@@ -62,11 +62,11 @@ class Cyl1DMesh(object):
     # Counting
     ####################################################
 
-    def nCr():
+    def nCx():
         doc = "Number of cells in the radial direction"
         fget = lambda self: self.hr.size
         return locals()
-    nCr = property(**nCr())
+    nCx = property(**nCx())
 
     def nCz():
         doc = "Number of cells in the z direction"
@@ -76,10 +76,16 @@ class Cyl1DMesh(object):
 
     def nC():
         doc = "Total number of cells"
-        fget = lambda self: self.nCr * self.nCz
+        fget = lambda self: self.nCx * self.nCz
         return locals()
     nC = property(**nC())
 
+    def nCv():
+        doc = "Total number of cells in each direction"
+        fget = lambda self: np.array([self.nCx, self.nCz])
+        return locals()
+    nCv = property(**nCv())    
+
     def nNr():
         doc = "Number of nodes in the radial direction"
         fget = lambda self: self.hr.size
@@ -106,10 +112,16 @@ class Cyl1DMesh(object):
 
     def nFz():
         doc = "Number of z faces"
-        fget = lambda self: self.nNz * self.nCr
+        fget = lambda self: self.nNz * self.nCx
         return locals()
     nFz = property(**nFz())
 
+    def nFv():
+        doc = "Total number of faces in each direction"
+        fget = lambda self: np.array([self.nFr, self.nFz])
+        return locals()
+    nFv = property(**nFv())
+
     def nF():
         doc = "Total number of faces"
         fget = lambda self: self.nFr + self.nFz
@@ -236,12 +248,12 @@ class Cyl1DMesh(object):
         def fget(self):
             if self._edgeCurl is None:
                 #1D Difference matricies
-                dr = sp.spdiags((np.ones((self.nCr+1, 1))*[-1, 1]).T, [-1,0], self.nCr, self.nCr, format="csr")
+                dr = sp.spdiags((np.ones((self.nCx+1, 1))*[-1, 1]).T, [-1,0], self.nCx, self.nCx, format="csr")
                 dz = sp.spdiags((np.ones((self.nCz+1, 1))*[-1, 1]).T, [0,1], self.nCz, self.nCz+1, format="csr")
 
                 #2D Difference matricies
                 Dr = sp.kron(sp.eye(self.nNz), dr)
-                Dz = -sp.kron(dz, sp.eye(self.nCr))  #Not sure about this negative
+                Dz = -sp.kron(dz, sp.eye(self.nCx))  #Not sure about this negative
 
                 #Edge curl operator
                 self._edgeCurl = sp.diags(1/self.area,0)*sp.vstack((Dz, Dr))*sp.diags(self.edge,0)
@@ -255,7 +267,7 @@ class Cyl1DMesh(object):
         def fget(self):
             if self._aveE2CC is None:
                 az = sp.spdiags(0.5*np.ones((2, self.nNz)), [-1,0], self.nNz, self.nCz, format='csr')
-                ar = sp.spdiags(0.5*np.ones((2, self.nCr)), [0, 1], self.nCr, self.nCr, format='csr')
+                ar = sp.spdiags(0.5*np.ones((2, self.nCx)), [0, 1], self.nCx, self.nCx, format='csr')
                 ar[0,0] = 1
                 self._aveE2CC = sp.kron(az, ar).T
             return self._aveE2CC
@@ -268,10 +280,10 @@ class Cyl1DMesh(object):
         def fget(self):
             if self._aveF2CC is None:
                 az = sp.spdiags(0.5*np.ones((2, self.nNz)), [-1,0], self.nNz, self.nCz, format='csr')
-                ar = sp.spdiags(0.5*np.ones((2, self.nCr)), [0, 1], self.nCr, self.nCr, format='csr')
+                ar = sp.spdiags(0.5*np.ones((2, self.nCx)), [0, 1], self.nCx, self.nCx, format='csr')
                 ar[0,0] = 1
                 Afr = sp.kron(sp.eye(self.nCz),ar)
-                Afz = sp.kron(az,sp.eye(self.nCr))
+                Afz = sp.kron(az,sp.eye(self.nCx))
                 self._aveF2CC = sp.vstack((Afr,Afz)).T
             return self._aveF2CC
         return locals()
@@ -305,7 +317,7 @@ class Cyl1DMesh(object):
         elif type(materialProp) is float:
             materialProp = np.ones(self.nC)*materialProp
         elif materialProp.shape == (self.nCz,):
-            materialProp = materialProp.repeat(self.nCr)
+            materialProp = materialProp.repeat(self.nCx)
         materialProp = mkvc(materialProp)
         assert materialProp.shape == (self.nC,), "materialProp incorrect shape"
 
@@ -377,7 +389,7 @@ class Cyl1DMesh(object):
                     dFz = np.sum(dFz**2, axis=1)
 
                     indBL = np.argmin(dFz)         # Face below and to the left
-                    indAL = indBL + self.nCr        # Face above and to the left
+                    indAL = indBL + self.nCx        # Face above and to the left
 
                     zF_BL = self.gridFz[indBL,:]
                     zF_AL = self.gridFz[indAL,:]

SimPEG/mesh/DiffOperators.py

@@ -1,15 +1,16 @@
 import numpy as np
 from scipy import sparse as sp
-from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros
-
-
-def ddx(n):
-    """Define 1D derivatives, inner, this means we go from n+1 to n+1"""
-    return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr")
+from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros, ddx, av, avExtrap
 
 
 def checkBC(bc):
-    """ Checks if boundary condition 'bc' is valid. """
+    """
+
+        Checks if boundary condition 'bc' is valid.
+
+        Each bc must be either 'dirichlet' or 'neumann'
+
+    """
     if(type(bc) is str):
         bc = [bc, bc]
     assert type(bc) is list, 'bc must be a list'
@@ -22,7 +23,33 @@ def checkBC(bc):
 
 
 def ddxCellGrad(n, bc):
-    """Create 1D derivative operator from cell-centres to nodes this means we go from n to n+1"""
+    """
+        Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1
+
+        For Cell-Centered **Dirichlet**, use a ghost point::
+
+            (u_1 - u_g)/hf = grad
+
+                u_g       u_1      u_2
+                 *    |    *   |    *     ...
+                      ^
+                      0
+
+            u_g = - u_1
+            grad = 2*u1/dx
+            negitive on the other side.
+
+        For Cell-Centered **Neumann**, use a ghost point::
+
+            (u_1 - u_g)/hf = 0
+
+                u_g       u_1      u_2
+                 *    |    *   |    *     ...
+
+            u_g = u_1
+            grad = 0;  put a zero in.
+
+    """
     bc = checkBC(bc)
 
     D = sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [-1, 0], n+1, n, format="csr")
@@ -38,10 +65,55 @@ def ddxCellGrad(n, bc):
         D[-1, -1] = 0
     return D
 
+def ddxCellGradBC(n, bc):
+    """
 
-def av(n):
-    """Define 1D averaging operator from cell-centres to nodes."""
-    return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr")
+        Create 1D derivative operator from cell-centers to nodes this means we go from n to n+1
+
+        For Cell-Centered **Dirichlet**, use a ghost point::
+
+            (u_1 - u_g)/hf = grad
+
+             u_g       u_1      u_2
+              *    |    *   |    *     ...
+                   ^
+                  u_b
+
+        We know the value at the boundary (u_b)::
+
+            (u_g+u_1)/2 = u_b               (the average)
+            u_g = 2*u_b - u_1
+
+            So plug in to gradient:
+
+            (u_1 - (2*u_b - u_1))/hf = grad
+            2*(u_1-u_b)/hf = grad
+
+        Separate, because BC are known (and can move to RHS later)::
+
+            ( 2/hf )*u_1 + ( -2/hf )*u_b = grad
+
+                           (   ^   ) JUST RETURN THIS
+
+
+    """
+    bc = checkBC(bc)
+
+    ij   = (np.array([0, n]),np.array([0, 1]))
+    vals = np.zeros(2)
+
+    # Set the first side
+    if(bc[0] == 'dirichlet'):
+        vals[0] = -2
+    elif(bc[0] == 'neumann'):
+        vals[0] = 0
+    # Set the second side
+    if(bc[1] == 'dirichlet'):
+        vals[1] = 2
+    elif(bc[1] == 'neumann'):
+        vals[1] = 0
+    D = sp.csr_matrix((vals, ij), shape=(n+1,2))
+    return D
 
 
 class DiffOperators(object):
@@ -80,6 +152,73 @@ class DiffOperators(object):
     _faceDiv = None
     faceDiv = property(**faceDiv())
 
+    def faceDivx():
+        doc = "Construct divergence operator in the x component (face-stg to cell-centres)."
+
+        def fget(self):
+            if(self._faceDivx is None):
+                # The number of cell centers in each direction
+                n = self.n
+                # Compute faceDivergence operator on faces
+                if(self.dim == 1):
+                    D1 = ddx(n[0])
+                elif(self.dim == 2):
+                    D1 = sp.kron(speye(n[1]), ddx(n[0]))
+                elif(self.dim == 3):
+                    D1 = kron3(speye(n[2]), speye(n[1]), ddx(n[0]))
+                # Compute areas of cell faces & volumes
+                S = self.r(self.area, 'F','Fx', 'V')
+                V = self.vol
+                self._faceDivx = sdiag(1/V)*D1*sdiag(S)
+
+            return self._faceDivx
+        return locals()
+    _faceDivx = None
+    faceDivx = property(**faceDivx())
+
+    def faceDivy():
+        doc = "Construct divergence operator in the y component (face-stg to cell-centres)."
+
+        def fget(self):
+            if(self.dim < 2): return None
+            if(self._faceDivy is None):
+                # The number of cell centers in each direction
+                n = self.n
+                # Compute faceDivergence operator on faces
+                if(self.dim == 2):
+                    D2 = sp.kron(ddx(n[1]), speye(n[0]))
+                elif(self.dim == 3):
+                    D2 = kron3(speye(n[2]), ddx(n[1]), speye(n[0]))
+                # Compute areas of cell faces & volumes
+                S = self.r(self.area, 'F','Fy', 'V')
+                V = self.vol
+                self._faceDivy = sdiag(1/V)*D2*sdiag(S)
+
+            return self._faceDivy
+        return locals()
+    _faceDivy = None
+    faceDivy = property(**faceDivy())
+
+    def faceDivz():
+        doc = "Construct divergence operator in the z component (face-stg to cell-centres)."
+
+        def fget(self):
+            if(self.dim < 3): return None
+            if(self._faceDivz is None):
+                # The number of cell centers in each direction
+                n = self.n
+                # Compute faceDivergence operator on faces
+                D3 = kron3(ddx(n[2]), speye(n[1]), speye(n[0]))
+                # Compute areas of cell faces & volumes
+                S = self.r(self.area, 'F','Fz', 'V')
+                V = self.vol
+                self._faceDivz = sdiag(1/V)*D3*sdiag(S)
+
+            return self._faceDivz
+        return locals()
+    _faceDivz = None
+    faceDivz = property(**faceDivz())
+
     def nodalGrad():
         doc = "Construct gradient operator (nodes to edges)."
 
@@ -107,6 +246,38 @@ class DiffOperators(object):
     _nodalGrad = None
     nodalGrad = property(**nodalGrad())
 
+    def nodalLaplacian():
+        doc = "Construct laplacian operator (nodes to edges)."
+
+        def fget(self):
+            if(self._nodalLaplacian is None):
+                print 'Warning: Laplacian has not been tested rigorously.'
+                # The number of cell centers in each direction
+                n = self.n
+                # Compute divergence operator on faces
+                if(self.dim == 1):
+                    D1 = sdiag(1./self.hx) * ddx(mesh.nCx)
+                    L  = - D1.T*D1
+                elif(self.dim == 2):
+                    D1 = sdiag(1./self.hx) * ddx(n[0])
+                    D2 = sdiag(1./self.hy) * ddx(n[1])
+                    L1 = sp.kron(speye(n[1]+1), - D1.T * D1)
+                    L2 = sp.kron(- D2.T * D2, speye(n[0]+1))
+                    L  = L1 + L2
+                elif(self.dim == 3):
+                    D1 = sdiag(1./self.hx) * ddx(n[0])
+                    D2 = sdiag(1./self.hy) * ddx(n[1])
+                    D3 = sdiag(1./self.hz) * ddx(n[2])
+                    L1 = kron3(speye(n[2]+1), speye(n[1]+1), - D1.T * D1)
+                    L2 = kron3(speye(n[2]+1), - D2.T * D2, speye(n[0]+1))
+                    L3 = kron3(- D3.T * D3, speye(n[1]+1), speye(n[0]+1))
+                    L  = L1 + L2 + L3
+                self._nodalLaplacian = L
+            return self._nodalLaplacian
+        return locals()
+    _nodalLaplacian = None
+    nodalLaplacian = property(**nodalLaplacian())
+
     def setCellGradBC(self, BC):
         """
         Function that sets the boundary conditions for cell-centred derivative operators.
@@ -129,17 +300,19 @@ class DiffOperators(object):
         for i, bc_i in enumerate(BC):
             BC[i] = checkBC(bc_i)
 
-        self._cellGrad = None  # ensure we create a new gradient next time we call it
-        self._cellGradBC = BC
+        # ensure we create a new gradient next time we call it
+        self._cellGrad        = None
+        self._cellGradBC      = None
+        self._cellGradBC_list = BC
         return BC
-    _cellGradBC = 'neumann'
+    _cellGradBC_list = 'neumann'
 
     def cellGrad():
         doc = "The cell centered Gradient, takes you to cell faces."
 
         def fget(self):
             if(self._cellGrad is None):
-                BC = self.setCellGradBC(self._cellGradBC)
+                BC = self.setCellGradBC(self._cellGradBC_list)
                 n = self.n
                 if(self.dim == 1):
                     G = ddxCellGrad(n[0], BC[0])
@@ -154,13 +327,40 @@ class DiffOperators(object):
                     G = sp.vstack((G1, G2, G3), format="csr")
                 # Compute areas of cell faces & volumes
                 S = self.area
-                V = self.vol
-                self._cellGrad = sdiag(S)*G*sdiag(1/V)
+                V = self.aveCC2F*self.vol  # Average volume between adjacent cells
+                self._cellGrad = sdiag(S/V)*G
             return self._cellGrad
         return locals()
     _cellGrad = None
     cellGrad = property(**cellGrad())
 
+    def cellGradBC():
+        doc = "The cell centered Gradient boundary condition matrix"
+
+        def fget(self):
+            if(self._cellGradBC is None):
+                BC = self.setCellGradBC(self._cellGradBC_list)
+                n = self.n
+                if(self.dim == 1):
+                    G = ddxCellGradBC(n[0], BC[0])
+                elif(self.dim == 2):
+                    G1 = sp.kron(speye(n[1]), ddxCellGradBC(n[0], BC[0]))
+                    G2 = sp.kron(ddxCellGradBC(n[1], BC[1]), speye(n[0]))
+                    G = sp.block_diag((G1, G2), format="csr")
+                elif(self.dim == 3):
+                    G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGradBC(n[0], BC[0]))
+                    G2 = kron3(speye(n[2]), ddxCellGradBC(n[1], BC[1]), speye(n[0]))
+                    G3 = kron3(ddxCellGradBC(n[2], BC[2]), speye(n[1]), speye(n[0]))
+                    G = sp.block_diag((G1, G2, G3), format="csr")
+                # Compute areas of cell faces & volumes
+                S = self.area
+                V = self.aveCC2F*self.vol  # Average volume between adjacent cells
+                self._cellGradBC = sdiag(S/V)*G
+            return self._cellGradBC
+        return locals()
+    _cellGradBC = None
+    cellGradBC = property(**cellGradBC())
+
     def cellGradx():
         doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
 
@@ -175,19 +375,17 @@ class DiffOperators(object):
                 elif(self.dim == 3):
                     G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC))
                 # Compute areas of cell faces & volumes
-                S = self.r(self.area, 'F','Fx', 'V')
-                V = self.vol
-                self._cellGradx = sdiag(S)*G1*sdiag(1/V)
+                V = self.aveCC2F*self.vol
+                L = self.r(self.area/V, 'F','Fx', 'V')
+                self._cellGradx = sdiag(L)*G1
             return self._cellGradx
         return locals()
     cellGradx = property(**cellGradx())
 
-
     def cellGrady():
         doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
         def fget(self):
-            if self.dim < 2:
-                return None
+            if self.dim < 2: return None
             if getattr(self, '_cellGrady', None) is None:
                 BC = ['neumann', 'neumann']
                 n = self.n
@@ -196,33 +394,29 @@ class DiffOperators(object):
                 elif(self.dim == 3):
                     G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0]))
                 # Compute areas of cell faces & volumes
-                S = self.r(self.area, 'F','Fy', 'V')
-                V = self.vol
-                self._cellGrady = sdiag(S)*G2*sdiag(1/V)
+                V = self.aveCC2F*self.vol
+                L = self.r(self.area/V, 'F','Fy', 'V')
+                self._cellGrady = sdiag(L)*G2
             return self._cellGrady
         return locals()
     cellGrady = property(**cellGrady())
 
-
-
     def cellGradz():
         doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
         def fget(self):
-            if self.dim < 3:
-                return None
+            if self.dim < 3: return None
             if getattr(self, '_cellGradz', None) is None:
                 BC = ['neumann', 'neumann']
                 n = self.n
                 G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0]))
                 # Compute areas of cell faces & volumes
-                S = self.r(self.area, 'F','Fz', 'V')
-                V = self.vol
-                self._cellGradz = sdiag(S)*G3*sdiag(1/V)
+                V = self.aveCC2F*self.vol
+                L = self.r(self.area/V, 'F','Fz', 'V')
+                self._cellGradz = sdiag(L)*G3
             return self._cellGradz
         return locals()
     cellGradz = property(**cellGradz())
 
-
     def edgeCurl():
         doc = "Construct the 3D curl operator."
 
@@ -265,6 +459,8 @@ class DiffOperators(object):
     _edgeCurl = None
     edgeCurl = property(**edgeCurl())
 
+    # --------------- Averaging ---------------------
+
     def aveF2CC():
         doc = "Construct the averaging operator on cell faces to cell centers."
 
@@ -274,17 +470,37 @@ class DiffOperators(object):
                 if(self.dim == 1):
                     self._aveF2CC = av(n[0])
                 elif(self.dim == 2):
-                    self._aveF2CC = sp.hstack((sp.kron(speye(n[1]), av(n[0])),
-                                               sp.kron(av(n[1]), speye(n[0]))), format="csr")
+                    self._aveF2CC = (0.5)*sp.hstack((sp.kron(speye(n[1]), av(n[0])),
+                                                     sp.kron(av(n[1]), speye(n[0]))), format="csr")
                 elif(self.dim == 3):
-                    self._aveF2CC = sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])),
-                                               kron3(speye(n[2]), av(n[1]), speye(n[0])),
-                                               kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr")
+                    self._aveF2CC = (1./3.)*sp.hstack((kron3(speye(n[2]), speye(n[1]), av(n[0])),
+                                                       kron3(speye(n[2]), av(n[1]), speye(n[0])),
+                                                       kron3(av(n[2]), speye(n[1]), speye(n[0]))), format="csr")
             return self._aveF2CC
         return locals()
     _aveF2CC = None
     aveF2CC = property(**aveF2CC())
 
+    def aveCC2F():
+        doc = "Construct the averaging operator on cell cell centers to faces."
+
+        def fget(self):
+            if(self._aveCC2F is None):
+                n = self.n
+                if(self.dim == 1):
+                    self._aveCC2F = avExtrap(n[0])
+                elif(self.dim == 2):
+                    self._aveCC2F = sp.vstack((sp.kron(speye(n[1]), avExtrap(n[0])),
+                                               sp.kron(avExtrap(n[1]), speye(n[0]))), format="csr")
+                elif(self.dim == 3):
+                    self._aveCC2F = sp.vstack((kron3(speye(n[2]), speye(n[1]), avExtrap(n[0])),
+                                               kron3(speye(n[2]), avExtrap(n[1]), speye(n[0])),
+                                               kron3(avExtrap(n[2]), speye(n[1]), speye(n[0]))), format="csr")
+            return self._aveCC2F
+        return locals()
+    _aveCC2F = None
+    aveCC2F = property(**aveCC2F())
+
     def aveE2CC():
         doc = "Construct the averaging operator on cell edges to cell centers."
 
@@ -295,10 +511,10 @@ class DiffOperators(object):
                 if(self.dim == 1):
                     raise Exception('Edge Averaging does not make sense in 1D: Use Identity?')
                 elif(self.dim == 2):
-                    self._aveE2CC = sp.hstack((sp.kron(av(n[1]), speye(n[0])),
+                    self._aveE2CC = 0.5*sp.hstack((sp.kron(av(n[1]), speye(n[0])),
                                                sp.kron(speye(n[1]), av(n[0]))), format="csr")
                 elif(self.dim == 3):
-                    self._aveE2CC = sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])),
+                    self._aveE2CC = (1./3)*sp.hstack((kron3(av(n[2]), av(n[1]), speye(n[0])),
                                                kron3(av(n[2]), speye(n[1]), av(n[0])),
                                                kron3(speye(n[2]), av(n[1]), av(n[0]))), format="csr")
             return self._aveE2CC
@@ -316,37 +532,57 @@ class DiffOperators(object):
                 if(self.dim == 1):
                     self._aveN2CC = av(n[0])
                 elif(self.dim == 2):
-                    self._aveN2CC = sp.hstack((sp.kron(av(n[1]), av(n[0])),
-                                               sp.kron(av(n[1]), av(n[0]))), format="csr")
+                    self._aveN2CC = sp.kron(av(n[1]), av(n[0])).tocsr()
                 elif(self.dim == 3):
-                    self._aveN2CC = sp.hstack((kron3(av(n[2]), av(n[1]), av(n[0])),
-                                               kron3(av(n[2]), av(n[1]), av(n[0])),
-                                               kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr")
+                    self._aveN2CC = kron3(av(n[2]), av(n[1]), av(n[0])).tocsr()
             return self._aveN2CC
         return locals()
     _aveN2CC = None
     aveN2CC = property(**aveN2CC())
 
-    def aveN2CCv():
-        doc = "Construct the averaging operator on cell nodes to cell centers, keeping each dimension separate."
+    def aveN2E():
+        doc = "Construct the averaging operator on cell nodes to cell edges, keeping each dimension separate."
 
         def fget(self):
-            if(self._aveN2CCv is None):
+            if(self._aveN2E is None):
                 # The number of cell centers in each direction
                 n = self.n
                 if(self.dim == 1):
-                    self._aveN2CCv = av(n[0])
+                    self._aveN2E = av(n[0])
                 elif(self.dim == 2):
-                    self._aveN2CCv = sp.block_diag((sp.kron(av(n[1]), av(n[0])),
-                                                    sp.kron(av(n[1]), av(n[0]))), format="csr")
+                    self._aveN2E = sp.vstack((sp.kron(speye(n[1]+1), av(n[0])),
+                                              sp.kron(av(n[1]), speye(n[0]+1))), format="csr")
                 elif(self.dim == 3):
-                    self._aveN2CCv = sp.block_diag((kron3(av(n[2]), av(n[1]), av(n[0])),
-                                                    kron3(av(n[2]), av(n[1]), av(n[0])),
-                                                    kron3(av(n[2]), av(n[1]), av(n[0]))), format="csr")
-            return self._aveN2CCv
+                    self._aveN2E = sp.vstack((kron3(speye(n[2]+1), speye(n[1]+1), av(n[0])),
+                                              kron3(speye(n[2]+1), av(n[1]), speye(n[0]+1)),
+                                              kron3(av(n[2]), speye(n[1]+1), speye(n[0]+1))), format="csr")
+            return self._aveN2E
         return locals()
-    _aveN2CCv = None
-    aveN2CCv = property(**aveN2CCv())
+    _aveN2E = None
+    aveN2E = property(**aveN2E())
+
+    def aveN2F():
+        doc = "Construct the averaging operator on cell nodes to cell faces, keeping each dimension separate."
+
+        def fget(self):
+            if(self._aveN2F is None):
+                # The number of cell centers in each direction
+                n = self.n
+                if(self.dim == 1):
+                    self._aveN2F = av(n[0])
+                elif(self.dim == 2):
+                    self._aveN2F = sp.vstack((sp.kron(av(n[1]), speye(n[0]+1)),
+                                              sp.kron(speye(n[1]+1), av(n[0]))), format="csr")
+                elif(self.dim == 3):
+                    self._aveN2F = sp.vstack((kron3(av(n[2]), av(n[1]), speye(n[0]+1)),
+                                              kron3(av(n[2]), speye(n[1]+1), av(n[0])),
+                                              kron3(speye(n[2]+1), av(n[1]), av(n[0]))), format="csr")
+            return self._aveN2F
+        return locals()
+    _aveN2F = None
+    aveN2F = property(**aveN2F())
+
+    # --------------- Methods ---------------------
 
     def getMass(self, materialProp=None, loc='e'):
         """ Produces mass matricies.

SimPEG/mesh/InnerProducts.py

@@ -174,8 +174,8 @@ def getFaceInnerProduct(mesh, mu=None, returnP=False):
 
     def Pxxx(pos):
         ind1 = sub2ind(mesh.nFx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]])
-        ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nF[0]
-        ind3 = sub2ind(mesh.nFz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nF[0] + mesh.nF[1]
+        ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nFv[0]
+        ind3 = sub2ind(mesh.nFz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nFv[0] + mesh.nFv[1]
 
         IND = np.r_[ind1, ind2, ind3].flatten()
 
@@ -286,7 +286,7 @@ def getFaceInnerProduct2D(mesh, mu=None, returnP=False):
 
     def Pxx(pos):
         ind1 = sub2ind(mesh.nFx, np.c_[ii + pos[0][0], jj + pos[0][1]])
-        ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nF[0]
+        ind2 = sub2ind(mesh.nFy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nFv[0]
 
         IND = np.r_[ind1, ind2].flatten()
 
@@ -387,8 +387,8 @@ def getEdgeInnerProduct(mesh, sigma=None, returnP=False):
 
     def Pxxx(pos):
         ind1 = sub2ind(mesh.nEx, np.c_[ii + pos[0][0], jj + pos[0][1], kk + pos[0][2]])
-        ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nE[0]
-        ind3 = sub2ind(mesh.nEz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nE[0] + mesh.nE[1]
+        ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1], kk + pos[1][2]]) + mesh.nEv[0]
+        ind3 = sub2ind(mesh.nEz, np.c_[ii + pos[2][0], jj + pos[2][1], kk + pos[2][2]]) + mesh.nEv[0] + mesh.nEv[1]
 
         IND = np.r_[ind1, ind2, ind3].flatten()
 
@@ -499,7 +499,7 @@ def getEdgeInnerProduct2D(mesh, sigma=None, returnP=False):
 
     def Pxx(pos):
         ind1 = sub2ind(mesh.nEx, np.c_[ii + pos[0][0], jj + pos[0][1]])
-        ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nE[0]
+        ind2 = sub2ind(mesh.nEy, np.c_[ii + pos[1][0], jj + pos[1][1]]) + mesh.nEv[0]
 
         IND = np.r_[ind1, ind2].flatten()
 

SimPEG/mesh/LogicallyOrthogonalMesh.py

@@ -45,8 +45,7 @@ class LogicallyOrthogonalMesh(BaseMesh, DiffOperators, InnerProducts, LomView):
 
         def fget(self):
             if self._gridCC is None:
-                ccV = (self.aveN2CCv*mkvc(self.gridN))
-                self._gridCC = ccV.reshape((-1, self.dim), order='F')
+                self._gridCC = np.concatenate([self.aveN2CC*self.gridN[:,i] for i in range(self.dim)]).reshape((-1,self.dim), order='F')
             return self._gridCC
         return locals()
     _gridCC = None  # Store grid by default

SimPEG/mesh/TensorMesh.py

@@ -26,17 +26,26 @@ class TensorMesh(BaseMesh, TensorView, DiffOperators, InnerProducts):
 
     .. plot:: examples/mesh/plot_TensorMesh.py
 
+    For a quick tensor mesh on a (10x12x15) unit cube::
+
+        mesh = TensorMesh([10, 12, 15])
+
     """
     _meshType = 'TENSOR'
 
-    def __init__(self, h, x0=None):
-        super(TensorMesh, self).__init__(np.array([x.size for x in h]), x0)
-
-        assert len(h) == len(self.x0), "Dimension mismatch. x0 != len(h)"
-
-        for i, h_i in enumerate(h):
+    def __init__(self, h_in, x0=None):
+        assert type(h_in) is list, 'h_in must be a list'
+        h = range(len(h_in))
+        for i, h_i in enumerate(h_in):
+            if type(h_i) in [int, long, float]:
+                # This gives you something over the unit cube.
+                h_i = np.ones(int(h_i))/int(h_i)
             assert type(h_i) == np.ndarray, ("h[%i] is not a numpy array." % i)
             assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i)
+            h[i] = h_i[:] # make a copy.
+
+        BaseMesh.__init__(self, np.array([x.size for x in h]), x0)
+        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]
@@ -396,7 +405,7 @@ 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.nF if 'F' in locType else self.nE
+            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))
             Q = sp.hstack(components)

SimPEG/mesh/TensorView.py

@@ -2,7 +2,7 @@ import numpy as np
 import matplotlib.pyplot as plt
 import matplotlib
 from mpl_toolkits.mplot3d import Axes3D
-from SimPEG.utils import mkvc
+from SimPEG.utils import mkvc, animate
 
 
 class TensorView(object):
@@ -14,7 +14,7 @@ class TensorView(object):
     def __init__(self):
         pass
 
-    def plotImage(self, I, imageType='CC', figNum=1,ax=None,direction='z',numbering=True,annotationColor='w',showIt=False):
+    def plotImage(self, I, imageType='CC', figNum=1,ax=None,direction='z',numbering=True,annotationColor='w',showIt=False,clim=None):
         """
         Mesh.plotImage(I)
 
@@ -89,18 +89,18 @@ class TensorView(object):
             # Determine the subplot number: 131, 121
             numPlots = 130 if plotAll else len(imageType)/2*10+100
             pltNum = 1
-            fx, fy, fz = self.r(I,'F','F','M')
+            fxyz = self.r(I,'F','F','M')
             if plotAll or 'Fx' in imageType:
                 ax_x = plt.subplot(numPlots+pltNum)
-                self.plotImage(fx, imageType='Fx', ax=ax_x, **options)
+                self.plotImage(fxyz[0], imageType='Fx', ax=ax_x, **options)
                 pltNum +=1
             if plotAll or 'Fy' in imageType:
                 ax_y = plt.subplot(numPlots+pltNum)
-                self.plotImage(fy, imageType='Fy', ax=ax_y, **options)
+                self.plotImage(fxyz[1], imageType='Fy', ax=ax_y, **options)
                 pltNum +=1
             if plotAll or 'Fz' in imageType:
                 ax_z = plt.subplot(numPlots+pltNum)
-                self.plotImage(fz, imageType='Fz', ax=ax_z, **options)
+                self.plotImage(fxyz[2], imageType='Fz', ax=ax_z, **options)
                 pltNum +=1
             return
         else:
@@ -141,7 +141,9 @@ class TensorView(object):
                 C = I[:].reshape(self.nEy, order='F')
                 C = 0.5*(C[:-1,:] + C[1:,:] )
 
-            ph = ax.pcolormesh(self.vectorNx, self.vectorNy, C.T)
+            if clim is None:
+                clim = [C.min(),C.max()]
+            ph = ax.pcolormesh(self.vectorNx, self.vectorNy, C.T, vmin=clim[0], vmax=clim[1])
             ax.axis('tight')
             ax.set_xlabel("x")
             ax.set_ylabel("y")
@@ -196,7 +198,10 @@ class TensorView(object):
                 xx = np.r_[0, np.cumsum(np.kron(np.ones((nX, 1)), self.hx).ravel())]
                 yy = np.r_[0, np.cumsum(np.kron(np.ones((nY, 1)), self.hy).ravel())]
                 # Plot the mesh
-                ph = ax.pcolormesh(xx, yy, C.T)
+
+                if clim is None:
+                    clim = [C.min(),C.max()]
+                ph = ax.pcolormesh(xx, yy, C.T, vmin=clim[0], vmax=clim[1])
                 # Plot the lines
                 gx =  np.arange(nX+1)*(self.vectorNx[-1]-self.x0[0])
                 gy =  np.arange(nY+1)*(self.vectorNy[-1]-self.x0[1])
@@ -336,3 +341,60 @@ class TensorView(object):
             ax.set_ylabel('x2')
             ax.set_zlabel('x3')
             if showIt: plt.show()
+
+    def slicer(mesh, var, imageType='CC', normal='z', index=0, ax=None, clim=None):
+        assert normal in 'xyz', 'normal must be x, y, or z'
+        if ax is None: ax = plt.subplot(111)
+        I = mesh.r(var,'CC','CC','M')
+        axes = [p for p in 'xyz' if p not in normal.lower()]
+        if normal is 'x': I = I[index,:,:]
+        if normal is 'y': I = I[:,index,:]
+        if normal is 'z': I = I[:,:,index]
+        if clim is None: clim = [I.min(),I.max()]
+        p = ax.pcolormesh(getattr(mesh,'vectorN'+axes[0]),getattr(mesh,'vectorN'+axes[1]),I.T,vmin=clim[0],vmax=clim[1])
+        ax.axis('tight')
+        ax.set_xlabel(axes[0])
+        ax.set_ylabel(axes[1])
+        return p
+
+    def videoSlicer(mesh,var,imageType='CC',normal='z',figsize=(10,8)):
+        assert mesh.dim > 2, 'This is for 3D meshes only.'
+        # First set up the figure, the axis, and the plot element we want to animate
+        fig = plt.figure(figsize=figsize)
+        ax = plt.axes()
+        clim = [var.min(),var.max()]
+        plt.colorbar(mesh.slicer(var, imageType=imageType, normal=normal, index=0, ax=ax, clim=clim))
+        tlt = plt.title(normal)
+
+        def animateFrame(i):
+            mesh.slicer(var, imageType=imageType, normal=normal, index=i, ax=ax, clim=clim)
+            tlt.set_text(normal.upper()+('-Slice: %d, %4.4f' % (i,getattr(mesh,'vectorCC'+normal)[i])))
+
+        return animate(fig, animateFrame, frames=mesh.nCv['xyz'.index(normal)])
+
+    def video(mesh,var,function,figsize=(10,8)):
+        """
+        Call a function for a list of models to create a video.
+
+        ::
+
+            def function(var, ax, clim, tlt, i):
+                tlt.set_text('%%d'%%i)
+                return mesh.plotImage(var, imageType='CC', ax=ax, clim=clim)
+
+            mesh.video([model1, model2, ..., modeln],function)
+        """
+        # First set up the figure, the axis, and the plot element we want to animate
+        fig = plt.figure(figsize=figsize)
+        ax = plt.axes()
+        VAR = np.concatenate(var)
+        clim = [VAR.min(),VAR.max()]
+        tlt = plt.title('')
+        plt.colorbar(function(var[0],ax,clim,tlt,0))
+
+        def animateFrame(i):
+            function(var[i],ax,clim,tlt,i)
+
+        return animate(fig, animateFrame, frames=len(var))
+
+

SimPEG/regularization/Regularization.py

@@ -1,4 +1,4 @@
-from SimPEG.utils import sdiag
+from SimPEG.utils import sdiag, count, timeIt, setKwargs
 import numpy as np
 
 class Regularization(object):
@@ -7,7 +7,7 @@ class Regularization(object):
     @property
     def mref(self):
         if getattr(self, '_mref', None) is None:
-            self._mref = np.zeros(self.mesh.nC);
+            return np.zeros(self.mesh.nC);
         return self._mref
     @mref.setter
     def mref(self, value):
@@ -40,21 +40,22 @@ class Regularization(object):
             self._Wz = sdiag(a)*self.mesh.cellGradz
         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):
+    def __init__(self, mesh, **kwargs):
+        setKwargs(self, **kwargs)
         self.mesh = mesh
-        self._Wx = None
-        self._Wy = None
-        self._Wz = None
-        self.alpha_s = 1e-6
-        self.alpha_x = 1
-        self.alpha_y = 1
-        self.alpha_z = 1
+
 
     def pnorm(self, r):
         return 0.5*r.dot(r)
 
+    @timeIt
     def modelObj(self, m):
         mresid = m - self.mref
 
@@ -69,6 +70,7 @@ class Regularization(object):
 
         return mobj
 
+    @timeIt
     def modelObjDeriv(self, m):
         """
 
@@ -104,8 +106,8 @@ class Regularization(object):
         return mobjDeriv
 
 
-    def modelObj2Deriv(self, m):
-        mresid = m - self.mref
+    @timeIt
+    def modelObj2Deriv(self):
 
         mobj2Deriv = self.alpha_s * self.Ws.T * self.Ws
 

SimPEG/tests/HTMLTestRunner.py

@@ -241,7 +241,8 @@ function showClassDetail(cid, count) {
     for (var i = 0; i < count; i++) {
         tid = id_list[i];
         if (toHide) {
-            document.getElementById('div_'+tid).style.display = 'none'
+            var divTid = document.getElementById('div_'+tid);
+            if(divTid !== null){divTid.style.display = 'none';}
             document.getElementById(tid).className = 'hiddenRow';
         }
         else {

SimPEG/tests/TestUtils.py

@@ -189,7 +189,7 @@ def Rosenbrock(x, return_g=True, return_H=True):
         out += (H,)
     return out if len(out) > 1 else out[0]
 
-def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None):
+def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None, expectedOrder=2, tolerance=0.85, eps=1e-10):
     """
         Basic derivative check
 
@@ -201,6 +201,9 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None):
         :param int num: number of times to reduce step length, h
         :param bool plotIt: if you would like to plot
         :param numpy.array dx: step direction
+        :param int expectedOrder: The order that you expect the derivative to yield.
+        :param float tolerance: The tolerance on the expected order.
+        :param float eps: What is zero?
         :rtype: bool
         :return: did you pass the test?!
 
@@ -243,9 +246,6 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None):
         order1 = np.log10(E1[:-1]/E1[1:])
         print "%d\t%1.2e\t%1.3e\t\t%1.3e\t\t%1.3f" % (i, t[i], E0[i], E1[i], np.nan if i == 0 else order1[i-1])
 
-    tolerance = 0.9
-    expectedOrder = 2
-    eps = 1e-10
     order0 = order0[E0[1:] > eps]
     order1 = order1[E1[1:] > eps]
     belowTol = order1.size == 0 and order0.size > 0
@@ -276,16 +276,16 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None):
 
 
 
-def getQuadratic(A, b):
+def getQuadratic(A, b, c=0):
     """
-        Given A and b, this returns a quadratic, Q
+        Given A, b and c, this returns a quadratic, Q
 
         .. math::
 
-            \mathbf{Q( x ) = 0.5 x A x + b x}
+            \mathbf{Q( x ) = 0.5 x A x + b x} + c
     """
     def Quadratic(x, return_g=True, return_H=True):
-        f = 0.5 * x.dot( A.dot(x)) + b.dot( x )
+        f = 0.5 * x.dot( A.dot(x)) + b.dot( x ) + c
         out = (f,)
         if return_g:
             g = A.dot(x) + b

SimPEG/tests/runTests.py

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

SimPEG/tests/test_LogicallyOrthogonalMesh.py

@@ -44,50 +44,54 @@ class BasicLOMTests(unittest.TestCase):
 
     def test_tangents(self):
         T = self.LOM2.tangents
-        self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM2.nE[0])))
-        self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM2.nE[0])))
-        self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM2.nE[1])))
-        self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM2.nE[1])))
+        self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM2.nEv[0])))
+        self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM2.nEv[0])))
+        self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM2.nEv[1])))
+        self.assertTrue(np.all(self.LOM2.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM2.nEv[1])))
 
         T = self.LOM3.tangents
-        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM3.nE[0])))
-        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM3.nE[0])))
-        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.LOM3.nE[0])))
+        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[0] == np.ones(self.LOM3.nEv[0])))
+        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[1] == np.zeros(self.LOM3.nEv[0])))
+        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ex', 'V')[2] == np.zeros(self.LOM3.nEv[0])))
 
-        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM3.nE[1])))
-        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM3.nE[1])))
-        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.LOM3.nE[1])))
+        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[0] == np.zeros(self.LOM3.nEv[1])))
+        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[1] == np.ones(self.LOM3.nEv[1])))
+        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ey', 'V')[2] == np.zeros(self.LOM3.nEv[1])))
 
-        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.LOM3.nE[2])))
-        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.LOM3.nE[2])))
-        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.LOM3.nE[2])))
+        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[0] == np.zeros(self.LOM3.nEv[2])))
+        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[1] == np.zeros(self.LOM3.nEv[2])))
+        self.assertTrue(np.all(self.LOM3.r(T, 'E', 'Ez', 'V')[2] == np.ones(self.LOM3.nEv[2])))
 
     def test_normals(self):
         N = self.LOM2.normals
-        self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM2.nF[0])))
-        self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM2.nF[0])))
-        self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM2.nF[1])))
-        self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM2.nF[1])))
+        self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM2.nFv[0])))
+        self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM2.nFv[0])))
+        self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM2.nFv[1])))
+        self.assertTrue(np.all(self.LOM2.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM2.nFv[1])))
 
         N = self.LOM3.normals
-        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM3.nF[0])))
-        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM3.nF[0])))
-        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.LOM3.nF[0])))
+        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[0] == np.ones(self.LOM3.nFv[0])))
+        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[1] == np.zeros(self.LOM3.nFv[0])))
+        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fx', 'V')[2] == np.zeros(self.LOM3.nFv[0])))
 
-        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM3.nF[1])))
-        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM3.nF[1])))
-        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.LOM3.nF[1])))
+        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[0] == np.zeros(self.LOM3.nFv[1])))
+        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[1] == np.ones(self.LOM3.nFv[1])))
+        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fy', 'V')[2] == np.zeros(self.LOM3.nFv[1])))
 
-        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.LOM3.nF[2])))
-        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.LOM3.nF[2])))
-        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.LOM3.nF[2])))
+        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[0] == np.zeros(self.LOM3.nFv[2])))
+        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[1] == np.zeros(self.LOM3.nFv[2])))
+        self.assertTrue(np.all(self.LOM3.r(N, 'F', 'Fz', 'V')[2] == np.ones(self.LOM3.nFv[2])))
 
     def test_grid(self):
+        self.assertTrue(np.all(self.LOM2.gridCC == self.TM2.gridCC))
+        self.assertTrue(np.all(self.LOM2.gridN == self.TM2.gridN))
         self.assertTrue(np.all(self.LOM2.gridFx == self.TM2.gridFx))
         self.assertTrue(np.all(self.LOM2.gridFy == self.TM2.gridFy))
         self.assertTrue(np.all(self.LOM2.gridEx == self.TM2.gridEx))
         self.assertTrue(np.all(self.LOM2.gridEy == self.TM2.gridEy))
 
+        self.assertTrue(np.all(self.LOM3.gridCC == self.TM3.gridCC))
+        self.assertTrue(np.all(self.LOM3.gridN == self.TM3.gridN))
         self.assertTrue(np.all(self.LOM3.gridFx == self.TM3.gridFx))
         self.assertTrue(np.all(self.LOM3.gridFy == self.TM3.gridFy))
         self.assertTrue(np.all(self.LOM3.gridFz == self.TM3.gridFz))

SimPEG/tests/test_basemesh.py

@@ -38,15 +38,17 @@ class TestBaseMesh(unittest.TestCase):
 
     def test_mesh_numbers(self):
         c = self.mesh.nC == 36
-        f = np.all(self.mesh.nF == [42, 54, 48])
-        e = np.all(self.mesh.nE == [72, 56, 63])
+        fv = np.all(self.mesh.nFv == [42, 54, 48])
+        ev = np.all(self.mesh.nEv == [72, 56, 63])
+        f = np.all(self.mesh.nF == np.sum([42, 54, 48]))
+        e = np.all(self.mesh.nE == np.sum([72, 56, 63]))
 
-        self.assertTrue(np.all([c, f, e]))
+        self.assertTrue(np.all([c, fv, ev, f, e]))
 
     def test_mesh_r_E_V(self):
-        ex = np.ones(self.mesh.nE[0])
-        ey = np.ones(self.mesh.nE[1])*2
-        ez = np.ones(self.mesh.nE[2])*3
+        ex = np.ones(self.mesh.nEv[0])
+        ey = np.ones(self.mesh.nEv[1])*2
+        ez = np.ones(self.mesh.nEv[2])*3
         e = np.r_[ex, ey, ez]
         tex = self.mesh.r(e, 'E', 'Ex', 'V')
         tey = self.mesh.r(e, 'E', 'Ey', 'V')
@@ -60,9 +62,9 @@ class TestBaseMesh(unittest.TestCase):
         self.assertTrue(np.all(tez == ez))
 
     def test_mesh_r_F_V(self):
-        fx = np.ones(self.mesh.nF[0])
-        fy = np.ones(self.mesh.nF[1])*2
-        fz = np.ones(self.mesh.nF[2])*3
+        fx = np.ones(self.mesh.nFv[0])
+        fy = np.ones(self.mesh.nFv[1])*2
+        fz = np.ones(self.mesh.nFv[2])*3
         f = np.r_[fx, fy, fz]
         tfx = self.mesh.r(f, 'F', 'Fx', 'V')
         tfy = self.mesh.r(f, 'F', 'Fy', 'V')
@@ -146,14 +148,16 @@ class TestMeshNumbers2D(unittest.TestCase):
 
     def test_mesh_numbers(self):
         c = self.mesh.nC == 12
-        f = np.all(self.mesh.nF == [14, 18])
-        e = np.all(self.mesh.nE == [18, 14])
+        fv = np.all(self.mesh.nFv == [14, 18])
+        ev = np.all(self.mesh.nEv == [18, 14])
+        f = np.all(self.mesh.nF == np.sum([14, 18]))
+        e = np.all(self.mesh.nE == np.sum([18, 14]))
 
-        self.assertTrue(np.all([c, f, e]))
+        self.assertTrue(np.all([c, fv, ev, f, e]))
 
     def test_mesh_r_E_V(self):
-        ex = np.ones(self.mesh.nE[0])
-        ey = np.ones(self.mesh.nE[1])*2
+        ex = np.ones(self.mesh.nEv[0])
+        ey = np.ones(self.mesh.nEv[1])*2
         e = np.r_[ex, ey]
         tex = self.mesh.r(e, 'E', 'Ex', 'V')
         tey = self.mesh.r(e, 'E', 'Ey', 'V')
@@ -165,8 +169,8 @@ class TestMeshNumbers2D(unittest.TestCase):
         self.assertRaises(AssertionError,   self.mesh.r, e, 'E', 'Ez', 'V')
 
     def test_mesh_r_F_V(self):
-        fx = np.ones(self.mesh.nF[0])
-        fy = np.ones(self.mesh.nF[1])*2
+        fx = np.ones(self.mesh.nFv[0])
+        fy = np.ones(self.mesh.nFv[1])*2
         f = np.r_[fx, fy]
         tfx = self.mesh.r(f, 'F', 'Fx', 'V')
         tfy = self.mesh.r(f, 'F', 'Fy', 'V')

SimPEG/tests/test_forward_DCproblem.py

@@ -2,7 +2,7 @@ import numpy as np
 import unittest
 from SimPEG.mesh import TensorMesh
 from SimPEG.utils import ModelBuilder, sdiag
-from SimPEG.forward import Problem, SyntheticProblem
+from SimPEG.forward import Problem
 from SimPEG.forward.DCProblem import *
 from TestUtils import checkDerivative
 from scipy.sparse.linalg import dsolve
@@ -40,18 +40,13 @@ class DCProblemTests(unittest.TestCase):
         P = Q.T
 
         # Create some data
-        class syntheticDCProblem(DCProblem, SyntheticProblem):
-            pass
 
-        synthetic = syntheticDCProblem(mesh);
-        synthetic.P = P
-        synthetic.RHS = q
-        dobs, Wd = synthetic.createData(mSynth, std=0.05)
-
-        # Now set up the problem to do some minimization
         problem = DCProblem(mesh)
         problem.P = P
         problem.RHS = q
+        dobs, Wd = 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

SimPEG/tests/test_operators.py

@@ -1,6 +1,7 @@
 import numpy as np
 import unittest
 from TestUtils import OrderTest
+import matplotlib.pyplot as plt
 
 MESHTYPES = ['uniformTensorMesh', 'uniformLOM', 'rotateLOM']
 call2 = lambda fun, xyz: fun(xyz[:, 0], xyz[:, 1])
@@ -48,8 +49,148 @@ class TestCurl(OrderTest):
         self.orderTest()
 
 
-class TestFaceDiv(OrderTest):
-    name = "Face Divergence"
+class TestCellGrad1D_InhomogeneousDirichlet(OrderTest):
+    name = "Cell Grad 1D - Dirichlet"
+    meshTypes = ['uniformTensorMesh']
+    meshDimension = 1
+    expectedOrders = 1 # because of the averaging involved in the ghost point. u_b = (u_n + u_g)/2
+    meshSizes = [8, 16, 32, 64]
+
+    def getError(self):
+        #Test function
+        fx = lambda x: -2*np.pi*np.sin(2*np.pi*x)
+        sol = lambda x: np.cos(2*np.pi*x)
+
+
+        xc = sol(self.M.gridCC)
+
+        gradX_anal = fx(self.M.gridFx)
+
+        bc = np.array([1,1])
+        self.M.setCellGradBC('dirichlet')
+        gradX = self.M.cellGrad.dot(xc) + self.M.cellGradBC*bc
+
+        err = np.linalg.norm((gradX-gradX_anal), np.inf)
+
+        return err
+
+    def test_order(self):
+        self.orderTest()
+
+class TestCellGrad2D_Dirichlet(OrderTest):
+    name = "Cell Grad 2D - Dirichlet"
+    meshTypes = ['uniformTensorMesh']
+    meshDimension = 2
+    meshSizes = [8, 16, 32, 64]
+
+    def getError(self):
+        #Test function
+        fx = lambda x, y: 2*np.pi*np.cos(2*np.pi*x)*np.sin(2*np.pi*y)
+        fy = lambda x, y: 2*np.pi*np.cos(2*np.pi*y)*np.sin(2*np.pi*x)
+        sol = lambda x, y: np.sin(2*np.pi*x)*np.sin(2*np.pi*y)
+
+        xc = call2(sol, self.M.gridCC)
+
+        Fc = cartF2(self.M, fx, fy)
+        gradX_anal = self.M.projectFaceVector(Fc)
+
+        self.M.setCellGradBC('dirichlet')
+        gradX = self.M.cellGrad.dot(xc)
+
+        err = np.linalg.norm((gradX-gradX_anal), np.inf)
+
+        return err
+
+    def test_order(self):
+        self.orderTest()
+
+
+class TestCellGrad3D_Dirichlet(OrderTest):
+    name = "Cell Grad 3D - Dirichlet"
+    meshTypes = ['uniformTensorMesh']
+    meshDimension = 3
+    meshSizes = [8, 16, 32]
+
+    def getError(self):
+        #Test function
+        fx = lambda x, y, z: 2*np.pi*np.cos(2*np.pi*x)*np.sin(2*np.pi*y)*np.sin(2*np.pi*z)
+        fy = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*x)*np.cos(2*np.pi*y)*np.sin(2*np.pi*z)
+        fz = lambda x, y, z: 2*np.pi*np.sin(2*np.pi*x)*np.sin(2*np.pi*y)*np.cos(2*np.pi*z)
+        sol = lambda x, y, z: np.sin(2*np.pi*x)*np.sin(2*np.pi*y)*np.sin(2*np.pi*z)
+
+        xc = call3(sol, self.M.gridCC)
+
+        Fc = cartF3(self.M, fx, fy, fz)
+        gradX_anal = self.M.projectFaceVector(Fc)
+
+        self.M.setCellGradBC('dirichlet')
+        gradX = self.M.cellGrad.dot(xc)
+
+        err = np.linalg.norm((gradX-gradX_anal), np.inf)
+
+        return err
+
+    def test_order(self):
+        self.orderTest()
+
+class TestCellGrad2D_Neumann(OrderTest):
+    name = "Cell Grad 2D - Neumann"
+    meshTypes = ['uniformTensorMesh']
+    meshDimension = 2
+    meshSizes = [8, 16, 32, 64]
+
+    def getError(self):
+        #Test function
+        fx = lambda x, y: -2*np.pi*np.sin(2*np.pi*x)*np.cos(2*np.pi*y)
+        fy = lambda x, y: -2*np.pi*np.sin(2*np.pi*y)*np.cos(2*np.pi*x)
+        sol = lambda x, y: np.cos(2*np.pi*x)*np.cos(2*np.pi*y)
+
+        xc = call2(sol, self.M.gridCC)
+
+        Fc = cartF2(self.M, fx, fy)
+        gradX_anal = self.M.projectFaceVector(Fc)
+
+        self.M.setCellGradBC('neumann')
+        gradX = self.M.cellGrad.dot(xc)
+
+        err = np.linalg.norm((gradX-gradX_anal), np.inf)
+
+        return err
+
+    def test_order(self):
+        self.orderTest()
+
+
+class TestCellGrad3D_Neumann(OrderTest):
+    name = "Cell Grad 3D - Neumann"
+    meshTypes = ['uniformTensorMesh']
+    meshDimension = 3
+    meshSizes = [8, 16, 32]
+
+    def getError(self):
+        #Test function
+        fx = lambda x, y, z: -2*np.pi*np.sin(2*np.pi*x)*np.cos(2*np.pi*y)*np.cos(2*np.pi*z)
+        fy = lambda x, y, z: -2*np.pi*np.cos(2*np.pi*x)*np.sin(2*np.pi*y)*np.cos(2*np.pi*z)
+        fz = lambda x, y, z: -2*np.pi*np.cos(2*np.pi*x)*np.cos(2*np.pi*y)*np.sin(2*np.pi*z)
+        sol = lambda x, y, z: np.cos(2*np.pi*x)*np.cos(2*np.pi*y)*np.cos(2*np.pi*z)
+
+        xc = call3(sol, self.M.gridCC)
+
+        Fc = cartF3(self.M, fx, fy, fz)
+        gradX_anal = self.M.projectFaceVector(Fc)
+
+        self.M.setCellGradBC('neumann')
+        gradX = self.M.cellGrad.dot(xc)
+
+        err = np.linalg.norm((gradX-gradX_anal), np.inf)
+
+        return err
+
+    def test_order(self):
+        self.orderTest()
+
+class TestFaceDiv3D(OrderTest):
+    name = "Face Divergence 3D"
     meshTypes = MESHTYPES
     meshSizes = [8, 16, 32]
 
@@ -155,6 +296,130 @@ class TestNodalGrad2D(OrderTest):
     def test_order(self):
         self.orderTest()
 
+class TestAveraging2D(OrderTest):
+    name = "Averaging 2D"
+    meshTypes = MESHTYPES
+    meshDimension = 2
+
+    def getError(self):
+        num = self.getAve(self.M) * self.getHere(self.M)
+        err = np.linalg.norm((self.getThere(self.M)-num), np.inf)
+        return err
+
+    def test_orderN2CC(self):
+        self.name = "Averaging 2D: N2CC"
+        fun = lambda x, y: (np.cos(x)+np.sin(y))
+        self.getHere = lambda M: call2(fun, M.gridN)
+        self.getThere = lambda M: call2(fun, M.gridCC)
+        self.getAve = lambda M: M.aveN2CC
+        self.orderTest()
+
+    def test_orderN2F(self):
+        self.name = "Averaging 2D: N2F"
+        fun = lambda x, y: (np.cos(x)+np.sin(y))
+        self.getHere = lambda M: call2(fun, M.gridN)
+        self.getThere = lambda M: np.r_[call2(fun, M.gridFx), call2(fun, M.gridFy)]
+        self.getAve = lambda M: M.aveN2F
+        self.orderTest()
+
+    def test_orderN2E(self):
+        self.name = "Averaging 2D: N2E"
+        fun = lambda x, y: (np.cos(x)+np.sin(y))
+        self.getHere = lambda M: call2(fun, M.gridN)
+        self.getThere = lambda M: np.r_[call2(fun, M.gridEx), call2(fun, M.gridEy)]
+        self.getAve = lambda M: M.aveN2E
+        self.orderTest()
+
+    def test_orderF2CC(self):
+        self.name = "Averaging 2D: F2CC"
+        fun = lambda x, y: (np.cos(x)+np.sin(y))
+        self.getHere = lambda M: np.r_[call2(fun, M.gridFx), call2(fun, M.gridFy)]
+        self.getThere = lambda M: call2(fun, M.gridCC)
+        self.getAve = lambda M: M.aveF2CC
+        self.orderTest()
+
+    def test_orderCC2F(self):
+        self.name = "Averaging 2D: CC2F"
+        fun = lambda x, y: (np.cos(x)+np.sin(y))
+        self.getHere = lambda M: call2(fun, M.gridCC)
+        self.getThere = lambda M: np.r_[call2(fun, M.gridFx), call2(fun, M.gridFy)]
+        self.getAve = lambda M: M.aveCC2F
+        self.expectedOrders = 1
+        self.orderTest()
+        self.expectedOrders = 2
+
+
+    def test_orderE2CC(self):
+        self.name = "Averaging 2D: E2CC"
+        fun = lambda x, y: (np.cos(x)+np.sin(y))
+        self.getHere = lambda M: np.r_[call2(fun, M.gridEx), call2(fun, M.gridEy)]
+        self.getThere = lambda M: call2(fun, M.gridCC)
+        self.getAve = lambda M: M.aveE2CC
+        self.orderTest()
+
+
+class TestAveraging3D(OrderTest):
+    name = "Averaging 3D"
+    meshTypes = MESHTYPES
+    meshDimension = 3
+
+    def getError(self):
+        num = self.getAve(self.M) * self.getHere(self.M)
+        err = np.linalg.norm((self.getThere(self.M)-num), np.inf)
+        return err
+
+    def test_orderN2CC(self):
+        self.name = "Averaging 3D: N2CC"
+        fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
+        self.getHere = lambda M: call3(fun, M.gridN)
+        self.getThere = lambda M: call3(fun, M.gridCC)
+        self.getAve = lambda M: M.aveN2CC
+        self.orderTest()
+
+    def test_orderN2F(self):
+        self.name = "Averaging 3D: N2F"
+        fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
+        self.getHere = lambda M: call3(fun, M.gridN)
+        self.getThere = lambda M: np.r_[call3(fun, M.gridFx), call3(fun, M.gridFy), call3(fun, M.gridFz)]
+        self.getAve = lambda M: M.aveN2F
+        self.orderTest()
+
+    def test_orderN2E(self):
+        self.name = "Averaging 3D: N2E"
+        fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
+        self.getHere = lambda M: call3(fun, M.gridN)
+        self.getThere = lambda M: np.r_[call3(fun, M.gridEx), call3(fun, M.gridEy), call3(fun, M.gridEz)]
+        self.getAve = lambda M: M.aveN2E
+        self.orderTest()
+
+    def test_orderF2CC(self):
+        self.name = "Averaging 3D: F2CC"
+        fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
+        self.getHere = lambda M: np.r_[call3(fun, M.gridFx), call3(fun, M.gridFy), call3(fun, M.gridFz)]
+        self.getThere = lambda M: call3(fun, M.gridCC)
+        self.getAve = lambda M: M.aveF2CC
+        self.orderTest()
+
+
+    def test_orderE2CC(self):
+        self.name = "Averaging 3D: E2CC"
+        fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
+        self.getHere = lambda M: np.r_[call3(fun, M.gridEx), call3(fun, M.gridEy), call3(fun, M.gridEz)]
+        self.getThere = lambda M: call3(fun, M.gridCC)
+        self.getAve = lambda M: M.aveE2CC
+        self.orderTest()
+
+    def test_orderCC2F(self):
+        self.name = "Averaging 3D: CC2F"
+        fun = lambda x, y, z: (np.cos(x)+np.sin(y)+np.exp(z))
+        self.getHere = lambda M: call3(fun, M.gridCC)
+        self.getThere = lambda M: np.r_[call3(fun, M.gridFx), call3(fun, M.gridFy), call3(fun, M.gridFz)]
+        self.getAve = lambda M: M.aveCC2F
+        self.expectedOrders = 1
+        self.orderTest()
+        self.expectedOrders = 2
+
+
 
 if __name__ == '__main__':
     unittest.main()

SimPEG/utils/Geophysics/__init__.py

@@ -0,0 +1 @@
+import emSources

SimPEG/utils/Geophysics/emSources/__init__.py

@@ -0,0 +1 @@
+from emSources import MagneticDipoleVectorPotential

SimPEG/utils/Geophysics/emSources/emSources.py

@@ -0,0 +1,40 @@
+import numpy as np
+from scipy.constants import mu_0, pi
+
+def MagneticDipoleVectorPotential(txLoc, obsLoc, component, dipoleMoment=(0., 0., 1.)):
+    """
+        Calculate the vector potential of a set of magnetic dipoles
+        at given locations 'ref. <http://en.wikipedia.org/wiki/Dipole#Magnetic_vector_potential>'
+
+        :param numpy.ndarray txLoc: Location of the transmitter(s) (x, y, z)
+        :param numpy.ndarray obsLoc: Where the potentials will be calculated (x, y, z)
+        :param str component: The component to calculate - 'x', 'y', or 'z'
+        :param numpy.ndarray dipoleMoment: The vector dipole moment
+        :rtype: numpy.ndarray
+        :return: The vector potential each dipole at each observation location
+    """
+
+    if component=='x':
+        dimInd = 0
+    elif component=='y':
+        dimInd = 1
+    elif component=='z':
+        dimInd = 2
+    else:
+        raise ValueError('Invalid component')
+
+    txLoc = np.atleast_2d(txLoc)
+    obsLoc = np.atleast_2d(obsLoc)
+    dipoleMoment = np.atleast_2d(dipoleMoment)
+
+    nEdges = obsLoc.shape[0]
+    nTx = txLoc.shape[0]
+
+    m = np.array(dipoleMoment).repeat(nEdges, axis=0)
+    A = np.empty((nEdges, nTx))
+    for i in range(nTx):
+        dR = obsLoc - txLoc[i, np.newaxis].repeat(nEdges, axis=0)
+        mCr = np.cross(m, dR)
+        r = np.sqrt((dR**2).sum(axis=1))
+        A[:, i] = -(mu_0/(4*pi)) * mCr[:,dimInd]/(r**3)
+    return A

SimPEG/utils/ModelBuilder.py

@@ -1,4 +1,6 @@
 import numpy as np
+import scipy.ndimage as ndi
+import scipy.sparse as sp
 
 
 def getIndecesBlock(p0,p1,ccMesh):
@@ -130,6 +132,68 @@ def scalarConductivity(ccMesh,pFunction):
 
     return sigma
 
+
+
+def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]):
+    """
+        Create a random model by convolving a kernal with a
+        uniformly distributed model.
+
+        :param int,tuple shape: shape of the model.
+        :param int seed: pick which model to produce, prints the seed if you don't choose.
+        :param numpy.ndarray,list anisotropy: this is the (3 x n) blurring kernal that is used.
+        :param int its: number of smoothing iterations
+        :param list bounds: bounds on the model, len(list) == 2
+        :rtype: numpy.ndarray
+        :return: M, the model
+
+
+        .. plot::
+
+            import matplotlib.pyplot as plt
+            import SimPEG.utils.ModelBuilder as MB
+            plt.colorbar(plt.imshow(MB.randomModel((50,50),bounds=[-4,0])))
+            plt.title('A very cool, yet completely random model.')
+            plt.show()
+
+
+    """
+
+    if seed is None:
+        seed = np.random.randint(1e3)
+        print 'Using a seed of: ', seed
+
+    if type(shape) in [int, long, float]:
+        shape = (shape,) # make it a tuple for consistency
+
+    np.random.seed(seed)
+    mr = np.random.rand(*shape)
+    if anisotropy is None:
+        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)
+        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))
+    else:
+        assert len(anisotropy.shape) is len(shape), 'Anisotropy must be the same shape.'
+        smth = np.array(anisotropy,dtype=float)
+
+    smth = smth/smth.sum() # normalize
+    mi = mr
+    for i in range(its):
+        mi = ndi.convolve(mi, smth)
+
+    # scale the model to live between the bounds.
+    mi = (mi - mi.min())/(mi.max()-mi.min()) # scaled between 0 and 1
+    mi = mi*(bounds[1]-bounds[0])+bounds[0]
+
+
+    return mi
+
+
+
 if __name__ == '__main__':
 
     from SimPEG.mesh import TensorMesh

SimPEG/utils/Solver.py

@@ -1,9 +1,10 @@
 import numpy as np
-import scipy.sparse as sparse
+import scipy.sparse as sp
 import scipy.sparse.linalg as linalg
-from SimPEG.utils import mkvc
+from SimPEG.utils import mkvc, sdiag
+import warnings
 
-DEFAULTS = {'direct':'scipy', 'forward':'fortran', 'backward':'fortran', 'diagonal':'python'}
+DEFAULTS = {'direct':'scipy', 'iter':'scipy', 'forward':'fortran', 'backward':'fortran', 'diagonal':'python'}
 
 try:
     import TriSolve
@@ -18,6 +19,11 @@ except Exception, e:
         DEFAULTS['forward'] = 'python'
         DEFAULTS['backward'] = 'python'
 
+try:
+    import mumps
+except Exception, e:
+    print 'Warning: mumps solver not available.'
+
 class Solver(object):
     """
         Solver is a light wrapper on the various types of
@@ -45,13 +51,44 @@ class Solver(object):
     def __init__(self, A, doDirect=True, flag=None, options={}):
         assert type(doDirect) is bool, 'doDirect must be a boolean'
         assert flag in [None, 'L', 'U', 'D'], "flag must be set to None, 'L', 'U', or 'D'"
-
+        assert type(options) is dict, 'options must be a dictionary object'
         self.A = A
 
         self.dsolve = None
         self.doDirect = doDirect
         self.flag = flag
         self.options = options
+        if doDirect: return
+
+        # Now deal with iterative stuff only
+        if 'M' not in options:
+            warnings.warn("You should provide a preconditioner, M.", UserWarning)
+            return
+        M = options['M']
+        if type(M) is sp.linalg.LinearOperator:
+            return
+        PreconditionerList = ['J','GS']
+        if type(M) is str:
+            assert M in PreconditionerList, "M must be in the known preconditioner list. ['J','GS']"
+            M = (M,A) # use A as the base for the preconditioner.
+        if type(M) is tuple:
+            assert type(M[0]) is str and M[0] in PreconditionerList, "M as a tuple must be (str, Matrix) where str is in ['J','GS']: e.g. ('J', WtW) where J stands for Jacobi, and WtW is a sparse matrix."
+            if M[0] is 'J':
+                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')
+                def GS(f):
+                    return Uinv.solve(DD*Linv.solve(f))
+                options['M'] = sp.linalg.LinearOperator( A.shape, GS, dtype=A.dtype )
+
+        else:
+            raise Exception('M must be a LinearOperator or a tuple')
+
 
     def solve(self, b):
         """
@@ -81,6 +118,9 @@ class Solver(object):
 
     def clean(self):
         """Cleans up the memory"""
+        if self.options.has_key('backend'):
+            if self.options['backend'] == 'mumps':
+                self.mctx.destroy()
         del self.dsolve
         self.dsolve = None
 
@@ -88,6 +128,7 @@ class Solver(object):
         """
             Use solve instead of this interface.
 
+            :param numpy.ndarray b: the right hand side
             :param bool factorize: if you want to factorize and store factors
             :param str backend: which backend to use. Default is scipy
             :rtype: numpy.ndarray
@@ -97,6 +138,22 @@ class Solver(object):
 
         assert np.shape(self.A)[1] == np.shape(b)[0], 'Dimension mismatch'
 
+        if backend == 'scipy':
+            X = self.solveDirect_scipy(b, factorize)
+        elif backend == 'mumps':
+            X = self.solveDirect_mumps(b, factorize)
+
+        return X
+
+    def solveDirect_scipy(self, b, factorize):
+        """
+            Use solve instead of this interface.
+
+            :param numpy.ndarray b: the right hand side
+            :param bool factorize: if you want to factorize and store factors
+            :rtype: numpy.ndarray
+            :return: x
+        """
         if factorize and self.dsolve is None:
             self.A = self.A.tocsc()  # for efficiency
             self.dsolve = linalg.factorized(self.A)
@@ -118,8 +175,62 @@ class Solver(object):
 
         return X
 
-    def solveIter(self, b, M=None, iterSolver='CG'):
-        pass
+    def solveDirect_mumps(self, b, factorize):
+        """
+            Use solve instead of this interface.
+
+            :param numpy.ndarray b: the right hand side
+            :param bool factorize: if you want to factorize and store factors
+            :rtype: numpy.ndarray
+            :return: x
+        """        
+        if factorize and self.dsolve is None:
+            self.mctx = mumps.DMumpsContext()
+            self.mctx.set_icntl(14, 60)
+            # self.mctx.set_silent()
+            self.mctx.set_centralized_sparse(self.A)
+            self.mctx.run(job=4)
+
+            def mdsolve(rhs):
+                x = rhs.copy()
+                self.mctx.set_rhs(x)
+                self.mctx.run(job=3)
+                return x
+
+            self.dsolve = mdsolve
+
+        if len(b.shape) == 1 or b.shape[1] == 1:
+            # Just one RHS
+            if factorize:
+                X = self.dsolve(b)
+            else:
+                X = mumps.spsolve(self.A, b)
+
+        else:
+            # Multiple RHSs
+            X = np.empty_like(b)
+            for i in range(b.shape[1]):
+                if factorize:
+                    X[:,i] = self.dsolve(b[:,i])
+                else:
+                    X[:,i] = mumps.spsolve(self.A,b[:,i])
+    
+        return X
+
+    def solveIter(self, b, backend=None, M=None, iterSolver='CG', tol=1e-6, maxIter=50):
+        if backend is None: backend = DEFAULTS['iter']
+
+        algorithms = {'CG':sp.linalg.cg}
+        assert iterSolver in algorithms, "iterSolver must be 'CG', or implement it yourself and add it here!"
+        alg = algorithms[iterSolver]
+
+        if len(b.shape) == 1 or b.shape[1] == 1:
+            x, self.info = alg(self.A, b, M=M, tol=tol, maxiter=maxIter)
+        else:
+            x = np.empty_like(b)
+            for i in range(b.shape[1]):
+                x[:,i], self.info = alg(self.A, b[:,i], M=M, tol=tol, maxiter=maxIter)
+        return x
 
     def solveBackward(self, b, backend=None):
         """
@@ -132,9 +243,8 @@ class Solver(object):
             :return: x
         """
         if backend is None: backend = DEFAULTS['backward']
-        if type(self.A) is not sparse.csr.csr_matrix:
-            from scipy.sparse import csr_matrix
-            self.A = csr_matrix(self.A)
+        if type(self.A) is not sp.csr.csr_matrix:
+            self.A = sp.csr_matrix(self.A)
         vals = self.A.data
         rowptr = self.A.indptr
         colind = self.A.indices
@@ -164,7 +274,7 @@ class Solver(object):
             :return: x
         """
         if backend is None: backend = DEFAULTS['forward']
-        if type(self.A) is not sparse.csr.csr_matrix:
+        if type(self.A) is not sp.csr.csr_matrix:
             from scipy.sparse import csr_matrix
             self.A = csr_matrix(self.A)
         vals = self.A.data
@@ -240,13 +350,13 @@ if __name__ == '__main__':
     print np.linalg.norm(e-x,np.inf)
 
 
-    n = 6000
+    n = 600
     A_dense = np.random.random((n,n))
     L = np.tril(np.dot(A_dense, A_dense))  # Positive definite is better conditioned.
     e = np.ones(n)
     b = np.dot(L, e)
 
-    A = sparse.csr_matrix(L)
+    A = sp.csr_matrix(L)
     pSolve = Solver(A,flag='L',options={'backend':'python'});
     fSolve = Solver(A,flag='L',options={'backend':'fortran'})
     tic = time()
@@ -257,3 +367,17 @@ if __name__ == '__main__':
     x = fSolve.solve(b)
     toc = time() - tic
     print 'Error Forward Fortran = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc
+
+
+
+    A = -D*D.T
+    A[0,0] *= 10 # remove the constant null space from the matrix
+    e = np.ones(M.nC)
+    b = A.dot(e)
+
+    iSolve = Solver(A, doDirect=False,options={'M':('GS',A)})
+    tic = time()
+    x = iSolve.solve(b)
+    toc = time() - tic
+    print x
+    print 'Error CG  = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc, 'Info: ', iSolve.info

SimPEG/utils/__init__.py

@@ -4,12 +4,36 @@ import lomutils
 import interputils
 import ModelBuilder
 from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind
-from sputils import spzeros, kron3, speye, sdiag
+from sputils import spzeros, kron3, speye, sdiag, ddx, av, avExtrap
 from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube, exampleLomGird
 from interputils import interpmat
 from ipythonUtils import easyAnimate as animate
 import Solver
 from Solver import Solver
+import Geophysics
+
+import types
+import time
+import numpy as np
+from functools import wraps
+
+def hook(obj, method, name=None, overwrite=False, silent=False):
+    """
+        This dynamically binds a method to the instance of the class.
+
+        If name is None, the name of the method is used.
+    """
+    if name is None: 
+        name = method.__name__
+        if name == '<lambda>':
+            raise Exception('Must provide name to hook lambda functions.')
+    if not hasattr(obj,name) or overwrite:
+        setattr(obj, name, types.MethodType( method, obj ))
+        if getattr(obj,'debug',False):
+            print 'Method '+name+' was added to class.'
+    elif not silent or getattr(obj,'debug',False):
+        print 'Method '+name+' was not overwritten.'
+
 
 def setKwargs(obj, **kwargs):
     """Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist."""
@@ -18,6 +42,9 @@ 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)
 
 def printTitles(obj, printers, name='Print Titles', pad=''):
     titles = ''
@@ -59,3 +86,120 @@ def printStoppers(obj, stoppers, pad='', stop='STOP!', done='DONE!'):
         r = stopper['right'](obj)
         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)
+
+
+
+class Counter(object):
+    """
+        Counter allows anything that calls it to record iterations and
+        timings in a simple way.
+
+        Also has plotting functions that allow quick recalls of data.
+
+        If you want to use this, import *count* or *timeIt* and use them as decorators on class methods.
+
+        ::
+
+            class MyClass(object):
+                def __init__(self, url):
+                    self.counter = Counter()
+
+                @count
+                def MyMethod(self):
+                    pass
+
+                @timeIt
+                def MySecondMethod(self):
+                    pass
+
+            c = MyClass('blah')
+            for i in range(100): c.MyMethod()
+            for i in range(300): c.MySecondMethod()
+            c.counter.summary()
+
+    """
+    def __init__(self):
+        self._countList = {}
+        self._timeList = {}
+
+    def count(self, prop):
+        """
+            Increases the count of the property.
+        """
+        assert type(prop) is str, 'The property must be a string.'
+        if prop not in self._countList:
+            self._countList[prop] = 0
+        self._countList[prop] += 1
+
+    def countTic(self, prop):
+        """
+            Times a property call, this is the init call.
+        """
+        assert type(prop) is str, 'The property must be a string.'
+        if prop not in self._timeList:
+            self._timeList[prop] = []
+        self._timeList[prop].append(-time.time())
+
+    def countToc(self, prop):
+        """
+            Times a property call, this is the end call.
+        """
+        assert type(prop) is str, 'The property must be a string.'
+        assert prop in self._timeList, 'The property must already be in the dictionary.'
+        self._timeList[prop][-1] += time.time()
+
+    def summary(self):
+        """
+            Provides a text summary of the current counters and timers.
+        """
+        print 'Counters:'
+        for prop in sorted(self._countList):
+            print "  {0:<40}: {1:8d}".format(prop,self._countList[prop])
+        print '\nTimes:'+' '*40+'mean      sum'
+        for prop in sorted(self._timeList):
+            l = len(self._timeList[prop])
+            a = np.array(self._timeList[prop])
+            print "  {0:<40}: {1:4.2e}, {2:4.2e}, {3:4d}x".format(prop,a.mean(),a.sum(),l)
+
+def count(f):
+    @wraps(f)
+    def wrapper(self,*args,**kwargs):
+        counter = getattr(self,'counter',None)
+        if type(counter) is Counter: counter.count(self.__class__.__name__+'.'+f.__name__)
+        out = f(self,*args,**kwargs)
+        return out
+    return wrapper
+
+def timeIt(f):
+    @wraps(f)
+    def wrapper(self,*args,**kwargs):
+        counter = getattr(self,'counter',None)
+        if type(counter) is Counter: counter.countTic(self.__class__.__name__+'.'+f.__name__)
+        out = f(self,*args,**kwargs)
+        if type(counter) is Counter: counter.countToc(self.__class__.__name__+'.'+f.__name__)
+        return out
+    return wrapper
+
+
+if __name__ == '__main__':
+    class MyClass(object):
+        def __init__(self, url):
+            self.counter = Counter()
+
+        @count
+        def MyMethod(self):
+            pass
+
+        @timeIt
+        def MySecondMethod(self):
+            pass
+
+    c = MyClass('blah')
+    for i in range(100): c.MyMethod()
+    for i in range(300): c.MySecondMethod()
+    c.counter.summary()

SimPEG/utils/sputils.py

@@ -1,5 +1,6 @@
 from scipy import sparse as sp
 from matutils import mkvc
+import numpy as np
 
 
 def sdiag(h):
@@ -20,3 +21,18 @@ def kron3(A, B, C):
 def spzeros(n1, n2):
     """spzeros"""
     return sp.coo_matrix((n1, n2)).tocsr()
+
+
+def ddx(n):
+    """Define 1D derivatives, inner, this means we go from n+1 to n"""
+    return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr")
+
+
+def av(n):
+    """Define 1D averaging operator from nodes to cell-centers."""
+    return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr")
+
+def avExtrap(n):
+    """Define 1D averaging operator from cell-centers to nodes."""
+    Av = sp.spdiags((0.5*np.ones((n, 1))*[1, 1]).T, [-1, 0], n+1, n, format="csr") + sp.csr_matrix(([0.5,0.5],([0,n],[0,n-1])),shape=(n+1,n))
+    return Av

SimPEG/visualize/__init__.py

@@ -0,0 +1,2 @@
+import vtk
+#import mpl

SimPEG/visualize/vtk/__init__.py

@@ -0,0 +1,2 @@
+from vtkTools import vtkTools
+from vtkView import vtkView

SimPEG/visualize/vtk/vtkTools.py

@@ -0,0 +1,385 @@
+import numpy as np
+try:
+	import vtk, vtk.util.numpy_support as npsup, pdb
+except Exception, e:
+	print 'VTK import error. Please ensure you have VTK installed to use this visualization package.'
+from SimPEG.utils import mkvc
+
+
+class vtkTools(object):
+	"""
+	Class that interacts with VTK visulization toolkit.
+
+	"""
+
+	def __init__(self):
+		""" Initializes the VTK vtkTools.
+
+		"""
+
+		pass
+
+	@staticmethod
+	def makeCellVTKObject(mesh,model):
+		"""
+		Make and return a cell based VTK object for a simpeg mesh and model.
+
+		Input:
+		:param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
+		:param model, dictionary of numpy.array - Name('s) and array('s). Match number of cells
+
+		Output:
+        :rtype: vtkRecilinearGrid object
+        :return: vtkObj
+		"""
+
+		# Deal with dimensionalities
+		if mesh.dim >= 1:
+			vX = mesh.vectorNx
+			xD = mesh.nNx
+			yD,zD = 1,1
+			vY, vZ = np.array([0,0])
+		if mesh.dim >= 2:
+			vY = mesh.vectorNy
+			yD = mesh.nNy
+		if mesh.dim == 3:
+			vZ = mesh.vectorNz
+			zD = mesh.nNz
+		# Use rectilinear VTK grid.
+		# Asaign the spatial information.
+		vtkObj = vtk.vtkRectilinearGrid()
+		vtkObj.SetDimensions(xD,yD,zD)
+		vtkObj.SetXCoordinates(npsup.numpy_to_vtk(vX,deep=1))
+		vtkObj.SetYCoordinates(npsup.numpy_to_vtk(vY,deep=1))
+		vtkObj.SetZCoordinates(npsup.numpy_to_vtk(vZ,deep=1))
+
+		# Assign the model('s) to the object
+		for item in model.iteritems():
+			# Convert numpy array
+			vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1)
+			vtkDoubleArr.SetName(item[0])
+			vtkObj.GetCellData().AddArray(vtkDoubleArr)
+
+		vtkObj.GetCellData().SetActiveScalars(model.keys()[0])
+		return vtkObj
+
+	@staticmethod
+	def makeFaceVTKObject(mesh,model):
+		"""
+		Make and return a face based VTK object for a simpeg mesh and model.
+
+		Input:
+		:param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
+		:param model, dictionary of numpy.array - Name('s) and array('s).
+			Property array must be order hstack(Fx,Fy,Fz)
+
+		Output:
+        :rtype: vtkUnstructuredGrid object
+        :return: vtkObj
+		"""
+
+		## Convert simpeg mesh to VTK properties
+		# Convert mesh nodes to vtkPoints
+		vtkPts = vtk.vtkPoints()
+		vtkPts.SetData(npsup.numpy_to_vtk(mesh.gridN,deep=1))
+
+		# Define the face "cells"
+		# Using VTK_QUAD cell for faces (see VTK file format)
+		nodeMat = mesh.r(np.arange(mesh.nN,dtype='int64'),'N','N','M')
+		def faceR(mat,length):
+			return mat.T.reshape((length,1))
+		# First direction
+		nTFx = np.prod(mesh.nFx)
+		FxCellBlock = np.hstack([ 4*np.ones((nTFx,1),dtype='int64'),faceR(nodeMat[:,:-1,:-1],nTFx),faceR(nodeMat[:,1: ,:-1],nTFx),faceR(nodeMat[:,1: ,1: ],nTFx),faceR(nodeMat[:,:-1,1: ],nTFx)] )
+		FyCellBlock = np.array([],dtype='int64')
+		FzCellBlock = np.array([],dtype='int64')
+		# Second direction
+		if mesh.dim >= 2:
+			nTFy = np.prod(mesh.nFy)
+			FyCellBlock = np.hstack([ 4*np.ones((nTFy,1),dtype='int64'),faceR(nodeMat[:-1,:,:-1],nTFy),faceR(nodeMat[1: ,:,:-1],nTFy),faceR(nodeMat[1: ,:,1: ],nTFy),faceR(nodeMat[:-1,:,1: ],nTFy)] )
+		# Third direction
+		if mesh.dim == 3:
+			nTFz = np.prod(mesh.nFz)
+			FzCellBlock = np.hstack([ 4*np.ones((nTFz,1),dtype='int64'),faceR(nodeMat[:-1,:-1,:],nTFz),faceR(nodeMat[1: ,:-1,:],nTFz),faceR(nodeMat[1: ,1: ,:],nTFz),faceR(nodeMat[:-1,1: ,:],nTFz)] )
+		# 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))
+		# Cell type
+		FCellType = npsup.numpy_to_vtk(vtk.VTK_QUAD*np.ones(mesh.nF,dtype='uint8'),deep=1)
+		# Cell location
+		FCellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,mesh.nF*5,5,dtype='int64'),deep=1)
+
+		## Make the object
+		vtkObj = vtk.vtkUnstructuredGrid()
+		# Set the objects properties
+		vtkObj.SetPoints(vtkPts)
+		vtkObj.SetCells(FCellType,FCellLoc,FCellArr)
+
+		# Assign the model('s) to the object
+		for item in model.iteritems():
+			# Convert numpy array
+			vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1)
+			vtkDoubleArr.SetName(item[0])
+			vtkObj.GetCellData().AddArray(vtkDoubleArr)
+
+		vtkObj.GetCellData().SetActiveScalars(model.keys()[0])
+		vtkObj.Update()
+		return vtkObj
+
+	@staticmethod
+	def makeEdgeVTKObject(mesh,model):
+		"""
+		Make and return a edge based VTK object for a simpeg mesh and model.
+
+		Input:
+		:param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
+		:param model, dictionary of numpy.array - Name('s) and array('s).
+			Property array must be order hstack(Ex,Ey,Ez)
+
+		Output:
+        :rtype: vtkUnstructuredGrid object
+        :return: vtkObj
+		"""
+
+		## Convert simpeg mesh to VTK properties
+		# Convert mesh nodes to vtkPoints
+		vtkPts = vtk.vtkPoints()
+		vtkPts.SetData(npsup.numpy_to_vtk(mesh.gridN,deep=1))
+
+		# Define the face "cells"
+		# Using VTK_QUAD cell for faces (see VTK file format)
+		nodeMat = mesh.r(np.arange(mesh.nN,dtype='int64'),'N','N','M')
+		def edgeR(mat,length):
+			return mat.T.reshape((length,1))
+		# First direction
+		nTEx = np.prod(mesh.nEx)
+		ExCellBlock = np.hstack([ 2*np.ones((nTEx,1),dtype='int64'),edgeR(nodeMat[:-1,:,:],nTEx),edgeR(nodeMat[1:,:,:],nTEx)])
+		# Second direction
+		if mesh.dim >= 2:
+			nTEy = np.prod(mesh.nEy)
+			EyCellBlock = np.hstack([ 2*np.ones((nTEy,1),dtype='int64'),edgeR(nodeMat[:,:-1,:],nTEy),edgeR(nodeMat[:,1:,:],nTEy)])
+		# Third direction
+		if mesh.dim == 3:
+			nTEz = np.prod(mesh.nEz)
+			EzCellBlock = np.hstack([ 2*np.ones((nTEz,1),dtype='int64'),edgeR(nodeMat[:,:,:-1],nTEz),edgeR(nodeMat[:,:,1:],nTEz)])
+		# 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))
+		# Cell type
+		ECellType = npsup.numpy_to_vtk(vtk.VTK_LINE*np.ones(mesh.nE,dtype='uint8'),deep=1)
+		# Cell location
+		ECellLoc = npsup.numpy_to_vtkIdTypeArray(np.arange(0,mesh.nE*3,3,dtype='int64'),deep=1)
+
+		## Make the object
+		vtkObj = vtk.vtkUnstructuredGrid()
+		# Set the objects properties
+		vtkObj.SetPoints(vtkPts)
+		vtkObj.SetCells(ECellType,ECellLoc,ECellArr)
+
+		# Assign the model('s) to the object
+		for item in model.iteritems():
+			# Convert numpy array
+			vtkDoubleArr = npsup.numpy_to_vtk(item[1],deep=1)
+			vtkDoubleArr.SetName(item[0])
+			vtkObj.GetCellData().AddArray(vtkDoubleArr)
+
+		vtkObj.GetCellData().SetActiveScalars(model.keys()[0])
+		return vtkObj
+
+	@staticmethod
+	def makeRenderWindow(ren):
+		renwin = vtk.vtkRenderWindow()
+		renwin.AddRenderer(ren)
+		iren = vtk.vtkRenderWindowInteractor()
+		iren.GetInteractorStyle().SetCurrentStyleToTrackballCamera()
+		iren.SetRenderWindow(renwin)
+
+		return iren, renwin
+
+
+	@staticmethod
+	def closeRenderWindow(iren):
+		renwin = iren.GetRenderWindow()
+		renwin.Finalize()
+		iren.TerminateApp()
+
+		del iren, renwin
+
+	@staticmethod
+	def makeVTKActor(vtkObj):
+		""" Makes a vtk mapper and Actor"""
+		mapper = vtk.vtkDataSetMapper()
+		mapper.SetInput(vtkObj)
+		actor = vtk.vtkActor()
+		actor.SetMapper(mapper)
+		actor.GetProperty().SetColor(0,0,0)
+		actor.GetProperty().SetRepresentationToWireframe()
+		return actor
+
+	@staticmethod
+	def makeVTKLODActor(vtkObj,clipper):
+		"""Make LOD vtk Actor"""
+		selectMapper = vtk.vtkDataSetMapper()
+		selectMapper.SetInputConnection(clipper.GetOutputPort())
+		selectMapper.SetScalarVisibility(1)
+		selectMapper.SetColorModeToMapScalars()
+		selectMapper.SetScalarModeToUseCellData()
+		selectMapper.SetScalarRange(clipper.GetInputDataObject(0,0).GetCellData().GetArray(0).GetRange())
+
+		selectActor = vtk.vtkLODActor()
+		selectActor.SetMapper(selectMapper)
+		selectActor.GetProperty().SetEdgeColor(1,0.5,0)
+		selectActor.GetProperty().SetEdgeVisibility(0)
+		selectActor.VisibilityOn()
+		selectActor.SetScale(1.01, 1.01, 1.01)
+		return selectActor
+
+	@staticmethod
+	def setScalar2View(vtkObj,scalarName):
+		""" Sets the sclar to view """
+		useArr = vtkObj.GetCellData().GetArray(scalarName)
+		if useArr == None:
+			raise IOError('Nerty array {:s} in the vtkObject'.format(scalarName))
+		vtkObj.GetCellData().SetActiveScalars(scalarName)
+
+	@staticmethod
+	def makeRectiVTKVOIThres(vtkObj,VOI,limits):
+		"""Make volume of interest and threshold for rectilinear grid."""
+		# Check for the input
+		cellCore = vtk.vtkExtractRectilinearGrid()
+		cellCore.SetVOI(VOI)
+		cellCore.SetInput(vtkObj)
+
+		cellThres = vtk.vtkThreshold()
+		cellThres.AllScalarsOn()
+		cellThres.SetInputConnection(cellCore.GetOutputPort())
+		cellThres.ThresholdByUpper(limits[0])
+		cellThres.ThresholdByLower(limits[1])
+		cellThres.Update()
+		return cellThres.GetOutput(), cellCore.GetOutput()
+
+	@staticmethod
+	def makeUnstructVTKVOIThres(vtkObj,extent,limits):
+		"""Make volume of interest and threshold for rectilinear grid."""
+		# Check for the input
+		cellCore = vtk.vtkExtractUnstructuredGrid()
+		cellCore.SetExtent(extent)
+		cellCore.SetInput(vtkObj)
+
+		cellThres = vtk.vtkThreshold()
+		cellThres.AllScalarsOn()
+		cellThres.SetInputConnection(cellCore.GetOutputPort())
+		cellThres.ThresholdByUpper(limits[0])
+		cellThres.ThresholdByLower(limits[1])
+		cellThres.Update()
+		return cellThres.GetOutput(), cellCore.GetOutput()
+
+	@staticmethod
+	def makePlaneClipper(vtkObj):
+		"""Makes a plane and clipper """
+		plane = vtk.vtkPlane()
+		clipper = vtk.vtkClipDataSet()
+		clipper.SetInputConnection(vtkObj.GetProducerPort())
+		clipper.SetClipFunction(plane)
+		clipper.InsideOutOff()
+		return clipper, plane
+
+	@staticmethod
+	def makePlaneWidget(vtkObj,iren,plane,actor):
+		"""Make an interactive planeWidget"""
+
+		# Callback function
+		def movePlane(obj, events):
+		    obj.GetPlane(intPlane)
+		    intActor.VisibilityOn()
+
+		# Associate the line widget with the interactor
+		planeWidget = vtk.vtkImplicitPlaneWidget()
+		planeWidget.SetInteractor(iren)
+		planeWidget.SetPlaceFactor(1.25)
+		planeWidget.SetInput(vtkObj)
+		planeWidget.PlaceWidget()
+		#planeWidget.AddObserver("InteractionEvent", movePlane)
+		planeWidget.SetScaleEnabled(0)
+		planeWidget.SetEnabled(1)
+		planeWidget.SetOutlineTranslation(0)
+		planeWidget.GetPlaneProperty().SetOpacity(0.1)
+		return planeWidget
+
+
+	@staticmethod
+	def startRenderWindow(iren):
+		""" Start a vtk rendering window"""
+		iren.Initialize()
+		renwin = iren.GetRenderWindow()
+		renwin.Render()
+		iren.Start()
+
+
+	# Simple write/read VTK xml model functions.
+	@staticmethod
+	def writeVTPFile(fileName,vtkPolyObject):
+	    '''Function to write vtk polydata file (vtp).'''
+	    polyWriter = vtk.vtkXMLPolyDataWriter()
+	    polyWriter.SetInput(vtkPolyObject)
+	    polyWriter.SetFileName(fileName)
+	    polyWriter.Update()
+
+	@staticmethod
+	def writeVTUFile(fileName,vtkUnstructuredGrid):
+	    '''Function to write vtk unstructured grid (vtu).'''
+	    Writer = vtk.vtkXMLUnstructuredGridWriter()
+	    Writer.SetInput(vtkUnstructuredGrid)
+	    Writer.SetFileName(fileName)
+	    Writer.Update()
+
+	@staticmethod
+	def writeVTRFile(fileName,vtkRectilinearGrid):
+	    '''Function to write vtk rectilinear grid (vtr).'''
+	    Writer = vtk.vtkXMLRectilinearGridWriter()
+	    Writer.SetInput(vtkRectilinearGrid)
+	    Writer.SetFileName(fileName)
+	    Writer.Update()
+
+	@staticmethod
+	def writeVTSFile(fileName,vtkStructuredGrid):
+	    '''Function to write vtk structured grid (vts).'''
+	    Writer = vtk.vtkXMLStructuredGridWriter()
+	    Writer.SetInput(vtkStructuredGrid)
+	    Writer.SetFileName(fileName)
+	    Writer.Update()
+
+	@staticmethod
+	def readVTSFile(fileName):
+	    '''Function to read vtk structured grid (vts) and return a grid object.'''
+	    Reader = vtk.vtkXMLStructuredGridReader()
+	    Reader.SetFileName(fileName)
+	    Reader.Update()
+	    return Reader.GetOutput()
+
+	@staticmethod
+	def readVTUFile(fileName):
+	    '''Function to read vtk structured grid (vtu) and return a grid object.'''
+	    Reader = vtk.vtkXMLUnstructuredGridReader()
+	    Reader.SetFileName(fileName)
+	    Reader.Update()
+	    return Reader.GetOutput()
+
+	@staticmethod
+	def readVTRFile(fileName):
+	    '''Function to read vtk structured grid (vtr) and return a grid object.'''
+	    Reader = vtk.vtkXMLRectilinearGridReader()
+	    Reader.SetFileName(fileName)
+	    Reader.Update()
+	    return Reader.GetOutput()
+
+	@staticmethod
+	def readVTPFile(fileName):
+	    '''Function to read vtk structured grid (vtp) and return a grid object.'''
+	    Reader = vtk.vtkXMLPolyDataReader()
+	    Reader.SetFileName(fileName)
+	    Reader.Update()
+	    return Reader.GetOutput()
+

SimPEG/visualize/vtk/vtkView.py

@@ -0,0 +1,160 @@
+import numpy as np
+try:
+	import vtk
+	#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.'
+import SimPEG as simpeg
+
+class vtkView(object):
+	"""
+	Class for storing and view of SimPEG models in VTK (visulization 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:
+		{'NameOfThePropertyModel': np.array of the properties}.
+		The property array has to be ordered in compliance with SimPEG standards.
+
+	::
+	Example of usages.
+
+		ToDo
+
+	"""
+
+	def __init__(self,mesh,propdict):
+		"""
+		"""
+
+		# 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
+		self._ren = None
+		self._iren = None
+		self._renwin = None
+		self._core = None
+		self._viewobj = None
+		self._plane = None
+		self._clipper = None
+		self._widget = None
+		self._actor = None
+		self._lut = None
+
+	def _readPropertyDictionary(self,propdict):
+		"""
+		Reads the property and assigns to the 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')
+		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])
+			else:
+				raise(Exception,'{:s} is not allowed as a dictonary key. Can be \'cell\',\'face\',\'edge\'.'.format(propitem[0]))
+
+	def Show(self):
+		"""
+		Open the VTK figure window and show the mesh.
+		"""
+		#vtkSP = simpeg.visualize.vtk.vtkTools
+		import SimPEG.visualize.vtk.vtkTools as vtkSP
+
+		# Make a renderer
+		self._ren = vtk.vtkRenderer()
+		# 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])
+		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)
+		# Set up the plane, clipper and the user interaction.
+		global intPlane, intActor
+		self._clipper, intPlane = vtkSP.makePlaneClipper(self._vtkobj)
+		intActor = vtkSP.makeVTKLODActor(self._vtkobj,self._clipper)
+		self._widget = vtkSP.makePlaneWidget(self._vtkobj,self._iren,self._clipper.GetClipFunction(),self._actor)
+		# Callback function
+		self._plane = intPlane
+		self._actor = intActor
+		def movePlane(obj, events):
+			global intPlane, intActor
+			obj.GetPlane(intPlane)
+			intActor.VisibilityOn()
+
+		self._widget.AddObserver("InteractionEvent",movePlane)
+		lut = vtk.vtkLookupTable()
+		lut.SetNumberOfColors(256)
+		lut.SetHueRange(0,0.66667)
+		lut.Build()
+		self._lut = lut
+		self._actor.GetMapper().SetLookupTable(lut)
+
+		# Set renderer options
+		self._ren.SetBackground(.5,.5,.5)
+		self._ren.AddActor(self._actor)
+
+		# Start the render Window
+		vtkSP.startRenderWindow(self._iren)
+		# Close the window when exited
+		vtkSP.closeRenderWindow(self._iren)
+		del self._iren, self._renwin
+
+
+
+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
+
+	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))}}
+	# Make the vtk viewer object.
+	vtkViewer = simpeg.visualize.vtk.vtkView(mesh,models)
+	# Show the image
+	vtkViewer.Show()

docs/api_LogicallyOrthogonalMesh.rst

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

docs/api_Optimize.rst

@@ -4,6 +4,8 @@ Optimize
 ********
 
 .. automodule:: SimPEG.inverse.Optimize
+    :show-inheritance:
+    :private-members:
     :members:
     :undoc-members:
 
@@ -12,6 +14,7 @@ Inversion
 *********
 
 .. automodule:: SimPEG.inverse.Inversion
+    :show-inheritance:
     :members:
     :undoc-members:
 

docs/api_Problem.rst

@@ -1,21 +1,24 @@
 .. _api_Problem:
 
 
-
 Problem
 *******
 
 .. automodule:: SimPEG.forward.Problem
+    :show-inheritance:
     :members:
     :undoc-members:
+    :inherited-members:
 
 
 DCProblem
 *********
 
 .. automodule:: SimPEG.forward.DCProblem
+    :show-inheritance:
     :members:
     :undoc-members:
+    :inherited-members:
 
 
 
@@ -23,6 +26,7 @@ Linear Problem
 **************
 
 .. automodule:: SimPEG.forward.LinearProblem
+    :show-inheritance:
     :members:
     :undoc-members:
-
+    :inherited-members:

docs/api_Solver.rst

@@ -1,9 +0,0 @@
-.. _api_Solver:
-
-Solver
-******
-
-.. automodule:: SimPEG.utils.Solver
-    :members:
-    :undoc-members:
-

docs/api_TensorMesh.rst

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

docs/api_Utils.rst

@@ -1,24 +1,52 @@
 .. _api_Utils:
 
+
+Solver
+******
+
+.. automodule:: SimPEG.utils.Solver
+    :members:
+    :undoc-members:
+
+
 Utilities
 *********
 
+.. automodule:: SimPEG.utils
+    :members:
+    :undoc-members:
+
+Matrix Utilities
+****************
+
 .. automodule:: SimPEG.utils.matutils
     :members:
     :undoc-members:
 
+Sparse Utilities
+****************
+
 .. automodule:: SimPEG.utils.sputils
     :members:
     :undoc-members:
 
+LOM Utilities
+*************
+
 .. automodule:: SimPEG.utils.lomutils
     :members:
     :undoc-members:
 
+Model Builder Utilities
+***********************
+
 .. automodule:: SimPEG.utils.ModelBuilder
     :members:
     :undoc-members:
 
+Interpolation Utilities
+***********************
+
 .. automodule:: SimPEG.utils.interputils
     :members:
     :undoc-members:

docs/index.rst

@@ -1,5 +1,7 @@
-SimPEG
-======
+.. image:: simpeg-logo.png
+   :width: 300 px
+   :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.
 
@@ -56,7 +58,6 @@ Utility Codes
 .. toctree::
    :maxdepth: 2
 
-   api_Solver
    api_Utils
 
 

notebooks/3DRenderingWithvtkTools.ipynb

@@ -1,6 +1,10 @@
 {
  "metadata": {
+<<<<<<< HEAD
   "name": "3D rendering with vtkTools"
+=======
+  "name": ""
+>>>>>>> develop
  },
  "nbformat": 3,
  "nbformat_minor": 0,
@@ -10,33 +14,75 @@
     {
      "cell_type": "code",
      "collapsed": false,
+<<<<<<< HEAD
      "input": "import numpy as np, vtk\nimport SimPEG as simpeg",
      "language": "python",
      "metadata": {},
      "outputs": [],
      "prompt_number": 1
+=======
+     "input": [
+      "import numpy as np, vtk\n",
+      "import SimPEG as simpeg"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [],
+     "prompt_number": 5
+>>>>>>> develop
     },
     {
      "cell_type": "code",
      "collapsed": false,
+<<<<<<< HEAD
      "input": "#Make a mesh and model\nx0 = np.zeros(3)\nh1 = np.ones(20)*5\nh2 = np.ones(10)*10\nh3 = np.ones(5)*20\n\nmesh = simpeg.mesh.TensorMesh([h1,h2,h3],x0)\n",
      "language": "python",
      "metadata": {},
      "outputs": [],
      "prompt_number": 2
+=======
+     "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
+>>>>>>> develop
     },
     {
      "cell_type": "code",
      "collapsed": false,
+<<<<<<< HEAD
      "input": "# Make a models that correspond to the cells, faces and edges.\nmodels = {'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.\nvtkViewer = simpeg.visulize.vtk.vtkView(mesh,models)                            \n                                                            ",
      "language": "python",
      "metadata": {},
      "outputs": [],
      "prompt_number": 3
+=======
+     "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
+>>>>>>> develop
     },
     {
      "cell_type": "code",
      "collapsed": false,
+<<<<<<< HEAD
      "input": "# Show the image \nvtkViewer.Show(imageType='cell')\n",
      "language": "python",
      "metadata": {},
@@ -51,6 +97,16 @@
      "metadata": {},
      "outputs": [],
      "prompt_number": 4
+=======
+     "input": [
+      "# Show the image \n",
+      "vtkViewer.Show()\n"
+     ],
+     "language": "python",
+     "metadata": {},
+     "outputs": [],
+     "prompt_number": "*"
+>>>>>>> develop
     }
    ],
    "metadata": {}

notebooks/BFGS.ipynb

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