Merge branch 'origin/master'

.gitignore

@@ -1,4 +1,6 @@
 *.pyc
-SimPEG.sublime-project
-SimPEG.sublime-workspace
-docs/_build/
+*.so
+*.sublime-project
+*.sublime-workspace
+docs/_build/
+myNotebooks/*

SimPEG/__init__.py

@@ -5,3 +5,4 @@ import inverse
 import visulize
 import forward
 import regularization
+import examples

SimPEG/forward/DCProblem.py

@@ -5,7 +5,6 @@ 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):
@@ -16,7 +15,7 @@ class DCProblem(ModelTransforms.LogModel, Problem):
 
     """
     def __init__(self, mesh):
-        super(DCProblem, self).__init__(mesh)
+        Problem.__init__(self, mesh)
         self.mesh.setCellGradBC('neumann')
 
     def reshapeFields(self, u):

SimPEG/forward/LinearProblem.py

@@ -13,13 +13,13 @@ class LinearProblem(Problem):
         return self.G.dot(m)
 
     def J(self, m, v, u=None):
-        return G.dot(v)
+        return self.G.dot(v)
 
     def Jt(self, m, v, u=None):
-        return G.T.dot(v)
+        return self.G.T.dot(v)
 
-if __name__ == '__main__':
-    N = 100
+
+def example(N):
     h = np.ones(N)/N
     M = TensorMesh([h])
 
@@ -28,8 +28,6 @@ if __name__ == '__main__':
     p = -0.25
     q = 0.25
 
-
-
     g = lambda k: np.exp(p*jk[k]*M.vectorCCx)*np.cos(2*np.pi*q*jk[k]*M.vectorCCx)
 
     G = np.empty((nk, M.nC))
@@ -38,12 +36,6 @@ if __name__ == '__main__':
         G[i,:] = g(i)
 
 
-
-    plt.figure(1)
-    for i in range(nk):
-        plt.plot(G[i,:])
-
-
     m_true = np.zeros(M.nC)
     m_true[M.vectorCCx > 0.3] = 1.
     m_true[M.vectorCCx > 0.45] = -0.5
@@ -55,29 +47,29 @@ if __name__ == '__main__':
 
     d_obs = d_true + noise
 
-    # plt.figure(3)
-    # plt.plot(d_true,'-o')
-    # plt.plot(d_obs,'r-o')
-
-
-
-
-
     prob = LinearProblem(M)
     prob.G = G
     prob.dobs = d_obs
     prob.std = np.ones_like(d_obs)*0.1
 
+    return prob, m_true
+
+
+if __name__ == '__main__':
+
+    prob, m_true = example(100)
+    M = prob.mesh
+
     reg = Regularization(M)
-
     opt = InexactGaussNewton(maxIter=20)
-
     inv = Inversion(prob,reg,opt,beta0=1e-4)
-
     m0 = np.zeros_like(m_true)
 
     mrec = inv.run(m0)
 
+    plt.figure(1)
+    for i in range(prob.G.shape[0]):
+        plt.plot(prob.G[i,:])
 
     plt.figure(2)
 

SimPEG/forward/Problem.py

@@ -140,7 +140,7 @@ class Problem(object):
             This can often be computed given a vector (i.e. J(v)) rather than stored, as J is a large dense matrix.
 
         """
-        pass
+        raise NotImplementedError('J is not yet implemented.')
 
     def Jt(self, m, v, u=None):
         """
@@ -152,7 +152,7 @@ class Problem(object):
 
             Effect of transpose of J on a vector v.
         """
-        pass
+        raise NotImplementedError('Jt is not yet implemented.')
 
 
     def J_approx(self, m, v, u=None):

SimPEG/inverse/BetaSchedule.py

@@ -8,5 +8,5 @@ class Cooling(object):
 
     def getBeta(self):
         if self._beta is None:
-            return beta0
-        return self._beta / beta_coolingFactor
+            return self.beta0
+        return self._beta / self.beta_coolingFactor

SimPEG/inverse/Inversion.py

@@ -1,35 +1,33 @@
 import numpy as np
 import scipy.sparse as sp
-from SimPEG.utils import sdiag, mkvc
+import SimPEG
+from SimPEG.utils import sdiag, mkvc, setKwargs, checkStoppers, printStoppers
+from Optimize import Remember
+from BetaSchedule import Cooling
 
-class Inversion(object):
-    """docstring for Inversion"""
+class BaseInversion(object):
+    """docstring for BaseInversion"""
 
-    maxIter = 10
-    name = 'SimPEG Inversion'
+    maxIter = 1
+    name = 'BaseInversion'
+    debug = False
+    beta0 = 1e4
 
     def __init__(self, prob, reg, opt, **kwargs):
+        setKwargs(self, **kwargs)
         self.prob = prob
         self.reg = reg
         self.opt = opt
         self.opt.parent = self
-        self.setKwargs(**kwargs)
 
-    def setKwargs(self, **kwargs):
-        """Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist."""
-        for attr in kwargs:
-            if hasattr(self, attr):
-                setattr(self, attr, kwargs[attr])
-            else:
-                raise Exception('%s attr is not recognized' % attr)
+        self.stoppers = [SimPEG.inverse.StoppingCriteria.iteration, SimPEG.inverse.StoppingCriteria.phi_d_target_Inversion]
 
-    def printInit(self):
-        print "%s %s %s" % ('='*22, self.name, '='*22)
-        print "  #    beta     phi_d      phi_m       f     norm(dJ)   #LS"
-        print "%s" % '-'*62
-
-    def printIter(self):
-        print "%3d  %1.2e  %1.2e  %1.2e  %1.2e  %1.2e  %3d" % (self.opt._iter, self._beta, self._phi_d_last, self._phi_m_last, self.opt.f, np.linalg.norm(self.opt.g), self.opt._iterLS)
+        # 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)
 
     @property
     def Wd(self):
@@ -53,34 +51,85 @@ class Inversion(object):
         if getattr(self, '_phi_d_target', None) is None:
             return self.prob.dobs.size #
         return self._phi_d_target
+
     @phi_d_target.setter
     def phi_d_target(self, value):
         self._phi_d_target = value
 
     def run(self, m0):
-        m = m0
-        self._iter = 0
-        self._beta = None
+        self.startup(m0)
         while True:
             self._beta = self.getBeta()
-            m = self.opt.minimize(self.evalFunction,m)
+            self.m = self.opt.minimize(self.evalFunction, self.m)
+            self.doEndIteration()
             if self.stoppingCriteria(): break
-            self._iter += 1
-        return m
 
-    beta0 = 1.e2
-    beta_coolingFactor = 5.
+        self.printDone()
+        return self.m
+
+    def startup(self, m0):
+        """
+            **startup** is called at the start of any new run call.
+
+            If you have things that also need to run on startup, you can create a method::
+
+                def _startup*(self, x0):
+                    pass
+
+            Where the * can be any string. If present, _startup* will be called at the start of the default startup call.
+            You may also completely overwrite this function.
+
+            :param numpy.ndarray x0: initial x
+            :rtype: None
+            :return: None
+        """
+        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)
+
+        self.m = m0
+        self._iter = 0
+        self._beta = None
+
+    def doEndIteration(self):
+        """
+            **doEndIteration** is called at the end of each run iteration.
+
+            If you have things that also need to run at the end of every iteration, you can create a method::
+
+                def _doEndIteration*(self, xt):
+                    pass
+
+            Where the * can be any string. If present, _doEndIteration* will be called at the start of the default doEndIteration call.
+            You may also completely overwrite this function.
+
+            :param numpy.ndarray xt: tested new iterate that ensures a descent direction.
+            :rtype: None
+            :return: None
+        """
+        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)()
+
+        # store old values
+        self.phi_d_last = self.phi_d
+        self.phi_m_last = self.phi_m
+        self._iter += 1
 
     def getBeta(self):
-        if self._beta is None:
-            return self.beta0
-        return self._beta / self.beta_coolingFactor
+        return self.beta0
 
     def stoppingCriteria(self):
-        self._STOP = np.zeros(2,dtype=bool)
-        self._STOP[0] = self._iter >= self.maxIter
-        self._STOP[1] = self._phi_d_last <= self.phi_d_target
-        return np.any(self._STOP)
+        if self.debug: print 'checking stoppingCriteria'
+        return checkStoppers(self, self.stoppers)
+
+
+    def printDone(self):
+        """
+            **printDone** is called at the end of the inversion routine.
+
+        """
+        printStoppers(self, self.stoppers)
 
 
     def evalFunction(self, m, return_g=True, return_H=True):
@@ -89,8 +138,8 @@ class Inversion(object):
         phi_d = self.dataObj(m, u)
         phi_m = self.reg.modelObj(m)
 
-        self._phi_d_last = phi_d
-        self._phi_m_last = phi_m
+        self.phi_d = phi_d
+        self.phi_m = phi_m
 
         f = phi_d + self._beta * phi_m
 
@@ -111,7 +160,7 @@ class Inversion(object):
 
             operator = sp.linalg.LinearOperator( (m.size, m.size), H_fun, dtype=float )
             out += (operator,)
-        return out
+        return out if len(out) > 1 else out[0]
 
 
     def dataObj(self, m, u=None):
@@ -219,3 +268,10 @@ class Inversion(object):
 
         return dmisfit
 
+class Inversion(Cooling, Remember, BaseInversion):
+
+    maxIter = 10
+    name = "SimPEG Inversion"
+
+    def __init__(self, prob, reg, opt, **kwargs):
+        BaseInversion.__init__(self, prob, reg, opt, **kwargs)

SimPEG/inverse/Optimize.py

@@ -1,6 +1,6 @@
 import numpy as np
 import matplotlib.pyplot as plt
-from SimPEG.utils import mkvc, sdiag
+from SimPEG.utils import mkvc, sdiag, setKwargs, printTitles, printLine, printStoppers, checkStoppers
 norm = np.linalg.norm
 import scipy.sparse as sp
 from SimPEG import Solver
@@ -12,6 +12,75 @@ except Exception, e:
     print 'Warning: you may not have the required pubsub installed, use pypubsub. You will not be able to listen to events.'
     doPub = False
 
+class StoppingCriteria(object):
+    """docstring for StoppingCriteria"""
+
+    iteration   = { "str": "%d : maxIter   =     %3d    <= iter          =    %3d",
+                    "left": lambda M: M.maxIter, "right": lambda M: M._iter,
+                    "stopType": "critical"}
+
+    iterationLS = { "str": "%d : maxIterLS =     %3d    <= iterLS          =    %3d",
+                    "left": lambda M: M.maxIterLS, "right": lambda M: M._iterLS,
+                    "stopType": "critical"}
+
+    armijoGoldstein = { "str": "%d :    ft     = %1.4e <= alp*descent     = %1.4e",
+                        "left": lambda M: M._LS_ft, "right": lambda M: M.f + M.LSreduction * M._LS_descent,
+                        "stopType": "optimal"}
+
+    tolerance_f = { "str": "%d : |fc-fOld| = %1.4e <= tolF*(1+|f0|) = %1.4e",
+                    "left":     lambda M: 1 if M._iter==0 else abs(M.f-M.f_last), "right":    lambda M: 0 if M._iter==0 else M.tolF*(1+abs(M.f0)),
+                    "stopType": "optimal"}
+
+    moving_x = { "str": "%d : |xc-x_last| = %1.4e <= tolX*(1+|x0|) = %1.4e",
+                 "left": lambda M: 1 if M._iter==0 else norm(M.xc-M.x_last), "right": lambda M: 0 if M._iter==0 else M.tolX*(1+norm(M.x0)),
+                 "stopType": "optimal"}
+
+    tolerance_g = { "str": "%d : |proj(x-g)-x|    = %1.4e <= tolG          = %1.4e",
+                   "left": lambda M: norm(M.projection(M.xc - M.g) - M.xc), "right": lambda M: M.tolG,
+                   "stopType": "optimal"}
+
+    norm_g = { "str": "%d : |proj(x-g)-x|    = %1.4e <= 1e3*eps       = %1.4e",
+               "left": lambda M: norm(M.projection(M.xc - M.g) - M.xc), "right": lambda M: 1e3*M.eps,
+               "stopType": "critical"}
+
+    bindingSet = { "str": "%d : probSize  =    %3d   <= bindingSet      =    %3d",
+                   "left": lambda M: M.xc.size, "right": lambda M: np.sum(M.bindingSet(M.xc)),
+                   "stopType": "critical"}
+
+    bindingSet_LS = { "str": "%d : probSize  =    %3d   <= bindingSet      =    %3d",
+                      "left": lambda M: M._LS_xt.size, "right": lambda M: np.sum(M.bindingSet(M._LS_xt)),
+                      "stopType": "critical"}
+
+    phi_d_target_Minimize = { "str": "%d : phi_d  = %1.4e <= phi_d_target  = %1.4e ",
+                              "left": lambda M: M.parent.phi_d, "right": lambda M: M.parent.phi_d_target,
+                              "stopType": "critical"}
+
+    phi_d_target_Inversion = { "str": "%d : phi_d  = %1.4e <= phi_d_target  = %1.4e ",
+                               "left": lambda I: I.phi_d, "right": lambda I: I.phi_d_target,
+                               "stopType": "critical"}
+
+
+class IterationPrinters(object):
+    """docstring for IterationPrinters"""
+
+    iteration = {"title": "#", "value": lambda M: M._iter, "width": 5, "format": "%3d"}
+    f = {"title": "f", "value": lambda M: M.f, "width": 10, "format": "%1.2e"}
+    norm_g = {"title": "|proj(x-g)-x|", "value": lambda M: norm(M.projection(M.xc - M.g) - M.xc), "width": 15, "format": "%1.2e"}
+    totalLS = {"title": "LS", "value": lambda M: M._iterLS, "width": 5, "format": "%d"}
+
+    iterationLS = {"title": "#", "value": lambda M: (M._iter, M._iterLS), "width": 5, "format": "%3d.%d"}
+    LS_ft = {"title": "ft", "value": lambda M: M._LS_ft, "width": 10, "format": "%1.2e"}
+    LS_t = {"title": "t", "value": lambda M: M._LS_t, "width": 10, "format": "%0.5f"}
+    LS_armijoGoldstein = {"title": "f + alp*g.T*p", "value": lambda M: M.f + M.LSreduction*M._LS_descent, "width": 16, "format": "%1.2e"}
+
+    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"}
+
+    beta = {"title": "beta", "value": lambda M: M.parent._beta, "width": 10, "format":   "%1.2e"}
+    phi_d = {"title": "phi_d", "value": lambda M: M.parent.phi_d, "width": 10, "format":   "%1.2e"}
+    phi_m = {"title": "phi_m", "value": lambda M: M.parent.phi_m, "width": 10, "format":   "%1.2e"}
 
 
 class Minimize(object):
@@ -22,7 +91,7 @@ class Minimize(object):
 
     """
 
-    name = "GeneralOptimizationAlgorithm"
+    name = "General Optimization Algorithm"
 
     maxIter = 20
     maxIterLS = 10
@@ -34,17 +103,18 @@ class Minimize(object):
     tolG = 1e-1
     eps = 1e-5
 
+    debug   = False
+    debugLS = False
+
     def __init__(self, **kwargs):
         self._id = int(np.random.rand()*1e6)  # create a unique identifier to this program to be used in pubsub
-        self.setKwargs(**kwargs)
+        self.stoppers = [StoppingCriteria.tolerance_f, StoppingCriteria.moving_x, StoppingCriteria.tolerance_g, StoppingCriteria.norm_g, StoppingCriteria.iteration]
+        self.stoppersLS = [StoppingCriteria.armijoGoldstein, StoppingCriteria.iterationLS]
 
-    def setKwargs(self, **kwargs):
-        """Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist."""
-        for attr in kwargs:
-            if hasattr(self, attr):
-                setattr(self, attr, kwargs[attr])
-            else:
-                raise Exception('%s attr is not recognized' % attr)
+        self.printers = [IterationPrinters.iteration, IterationPrinters.f, IterationPrinters.norm_g, IterationPrinters.totalLS]
+        self.printersLS = [IterationPrinters.iterationLS, IterationPrinters.LS_ft, IterationPrinters.LS_t, IterationPrinters.LS_armijoGoldstein]
+
+        setKwargs(self, **kwargs)
 
     def minimize(self, evalFunction, x0):
         """
@@ -59,6 +129,14 @@ class Minimize(object):
 
             (f[, g][, H]) = evalFunction(x, return_g=False, return_H=False )
 
+            def evalFunction(x, return_g=False, return_H=False):
+                out = (f,)
+                if return_g:
+                    out += (g,)
+                if return_H:
+                    out += (H,)
+                return out if len(out) > 1 else out[0]
+
 
         Events are fired with the following inputs via pypubsub::
 
@@ -146,19 +224,33 @@ class Minimize(object):
                 xc = x0
                 _iter = _iterLS = 0
 
+            If you have things that also need to run on startup, you can create a method::
+
+                def _startup*(self, x0):
+                    pass
+
+            Where the * can be any string. If present, _startup* will be called at the start of the default startup call.
+            You may also completely overwrite this function.
+
             :param numpy.ndarray x0: initial x
             :rtype: None
             :return: None
         """
+        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)
+
         self._iter = 0
         self._iterLS = 0
-        self._STOP = np.zeros((5,1),dtype=bool)
 
+        x0 = self.projection(x0)  # ensure that we start of feasible.
         self.x0 = x0
         self.xc = x0
-        self.xOld = x0
+        self.f_last = np.nan
+        self.x_last = x0
 
-    def printInit(self):
+
+    def printInit(self, inLS=False):
         """
             **printInit** is called at the beginning of the optimization routine.
 
@@ -166,15 +258,12 @@ class Minimize(object):
             parent.printInit function and call that.
 
         """
-        if doPub: pub.sendMessage('Minimize.printInit', minimize=self)
-        if self.parent is not None and hasattr(self.parent, 'printInit'):
-            self.parent.printInit()
-            return
-        print "%s %s %s" % ('='*22, self.name, '='*22)
-        print "iter\tJc\t\tnorm(dJ)\tLS"
-        print "%s" % '-'*57
+        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)
 
-    def printIter(self):
+    def printIter(self, inLS=False):
         """
             **printIter** is called directly after function evaluations.
 
@@ -182,13 +271,11 @@ class Minimize(object):
             parent.printIter function and call that.
 
         """
-        if doPub: pub.sendMessage('Minimize.printIter', minimize=self)
-        if self.parent is not None and hasattr(self.parent, 'printIter'):
-            self.parent.printIter()
-            return
-        print "%3d\t%1.2e\t%1.2e\t%d" % (self._iter, self.f, norm(self.g), self._iterLS)
+        if doPub and not inLS: pub.sendMessage('Minimize.printIter', minimize=self)
+        pad = ' '*10 if inLS else ''
+        printLine(self, self.printers if not inLS else self.printersLS, pad=pad)
 
-    def printDone(self):
+    def printDone(self, inLS=False):
         """
             **printDone** is called at the end of the optimization routine.
 
@@ -196,31 +283,19 @@ class Minimize(object):
             parent.printDone function and call that.
 
         """
-        if doPub: pub.sendMessage('Minimize.printDone', minimize=self)
-        if self.parent is not None and hasattr(self.parent, 'printDone'):
-            self.parent.printDone()
-            return
-        print "%s STOP! %s" % ('-'*25,'-'*25)
-        # TODO: put controls on gradient value, min model update, and function value
-        if self._iter > 0:
-            print "%d : |fc-fOld| = %1.4e <= tolF*(1+|fStop|) = %1.4e"  % (self._STOP[0], abs(self.f-self.fOld), self.tolF*(1+abs(self.fStop)))
-            print "%d : |xc-xOld| = %1.4e <= tolX*(1+|x0|)    = %1.4e"  % (self._STOP[1], norm(self.xc-self.xOld), self.tolX*(1+norm(self.x0)))
-        print "%d : |g|       = %1.4e <= tolG*(1+|fStop|) = %1.4e"  % (self._STOP[2], norm(self.g), self.tolG*(1+abs(self.fStop)))
-        print "%d : |g|       = %1.4e <= 1e3*eps          = %1.4e"  % (self._STOP[3], norm(self.g), 1e3*self.eps)
-        print "%d : iter      = %3d\t <= maxIter\t       = %3d"     % (self._STOP[4], self._iter, self.maxIter)
-        print "%s DONE! %s\n" % ('='*25,'='*25)
+        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):
+
+    def stoppingCriteria(self, inLS=False):
         if self._iter == 0:
-            self.fStop = self.f  # Save this for stopping criteria
+            self.f0 = self.f
+            self.g0 = self.g
+        return checkStoppers(self, self.stoppers if not inLS else self.stoppersLS)
 
-        # check stopping rules
-        self._STOP[0] = self._iter > 0 and (abs(self.f-self.fOld) <= self.tolF*(1+abs(self.fStop)))
-        self._STOP[1] = self._iter > 0 and (norm(self.xc-self.xOld) <= self.tolX*(1+norm(self.x0)))
-        self._STOP[2] = norm(self.g) <= self.tolG*(1+abs(self.fStop))
-        self._STOP[3] = norm(self.g) <= 1e3*self.eps
-        self._STOP[4] = self._iter >= self.maxIter
-        return all(self._STOP[0:3]) | any(self._STOP[3:])
 
     def projection(self, p):
         """
@@ -278,6 +353,8 @@ class Minimize(object):
             p = self.maxStep*p/np.abs(p.max())
         return p
 
+    nameLS = "Armijo linesearch"
+
     def modifySearchDirection(self, p):
         """
             **modifySearchDirection** changes the search direction based on some sort of linesearch or trust-region criteria.
@@ -296,20 +373,23 @@ class Minimize(object):
             :rtype: numpy.ndarray,bool
             :return: (xt, passLS)
         """
-        # Armijo linesearch
-        descent = np.inner(self.g, p)
-        t = 1
-        iterLS = 0
-        while iterLS < self.maxIterLS:
-            xt = self.projection(self.xc + t*p)
-            ft = self.evalFunction(xt, return_g=False, return_H=False)
-            if ft < self.f + t*self.LSreduction*descent:
-                break
-            iterLS += 1
-            t = self.LSshorten*t
+        # Projected Armijo linesearch
+        self._LS_t = 1
+        self._iterLS = 0
+        while self._iterLS < self.maxIterLS:
+            self._LS_xt      = self.projection(self.xc + self._LS_t*p)
+            self._LS_ft      = self.evalFunction(self._LS_xt, return_g=False, return_H=False)
+            self._LS_descent = np.inner(self.g, self._LS_xt - self.xc)  # this takes into account multiplying by t, but is important for projection.
+            if self.stoppingCriteria(inLS=True): break
+            self._iterLS += 1
+            self._LS_t = self.LSshorten*self._LS_t
+            if self.debugLS:
+                if self._iterLS == 1: self.printInit(inLS=True)
+                self.printIter(inLS=True)
 
-        self._iterLS = iterLS
-        return xt, iterLS < self.maxIterLS
+        if self.debugLS and self._iterLS > 0: self.printDone(inLS=True)
+
+        return self._LS_xt, self._iterLS < self.maxIterLS
 
     def modifySearchDirectionBreak(self, p):
         """
@@ -327,35 +407,227 @@ class Minimize(object):
             :rtype: numpy.ndarray,bool
             :return: (xt, breakCaught)
         """
+        self.printDone(inLS=True)
         print 'The linesearch got broken. Boo.'
         return p, False
 
     def doEndIteration(self, xt):
         """
-        **doEndIteration** is called at the end of each minimize iteration.
+            **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.
+            By default, function values and x locations are shuffled to store 1 past iteration in memory.
 
-        self.xc must be updated in this code.
+            self.xc must be updated in this code.
 
-        :param numpy.ndarray xt: tested new iterate that ensures a descent direction.
-        :rtype: None
-        :return: None
+
+            If you have things that also need to run at the end of every iteration, you can create a method::
+
+                def _doEndIteration*(self, xt):
+                    pass
+
+            Where the * can be any string. If present, _doEndIteration* will be called at the start of the default doEndIteration call.
+            You may also completely overwrite this function.
+
+            :param numpy.ndarray xt: tested new iterate that ensures a descent direction.
+            :rtype: None
+            :return: None
         """
+        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)
+
         # store old values
-        self.fOld = self.f
-        self.xOld, self.xc = self.xc, xt
+        self.f_last = self.f
+        self.x_last, self.xc = self.xc, xt
         self._iter += 1
+        if self.debug: self.printDone()
 
 
-class GaussNewton(Minimize):
-    name = 'GaussNewton'
+
+class Remember(object):
+    """
+        This mixin remembers all the things you tend to forget.
+
+        You can remember parameters directly, naming the str in Minimize,
+        or pass a tuple with the name and the function that takes Minimize.
+
+        For Example::
+
+            opt.remember('f',('norm_g', lambda M: np.linalg.norm(M.g)))
+
+            opt.minimize(evalFunction, x0)
+
+            opt.recall('f')
+
+        The param name (str) can also be located in the parent (if no conflicts),
+        and it will be looked up by default.
+    """
+
+    _rememberThese = []
+
+    def remember(self, *args):
+        self._rememberThese = args
+
+    def recall(self, param):
+        assert param in self._rememberList, "You didn't tell me to remember "+param+", you gotta tell me what to remember!"
+        return self._rememberList[param]
+
+    def _startupRemember(self, x0):
+        self._rememberList = {}
+        for param in self._rememberThese:
+            if type(param) is str:
+                self._rememberList[param] = []
+            elif type(param) is tuple:
+                self._rememberList[param[0]] = []
+
+    def _doEndIterationRemember(self, *args):
+        for param in self._rememberThese:
+            if type(param) is str:
+                if self.debug: print 'Remember is remembering: ' + param
+                val = getattr(self, param, None)
+                if val is None and getattr(self, 'parent', None) is not None:
+                    # Look to the parent for the param if not found here.
+                    val = getattr(self.parent, param, None)
+                self._rememberList[param].append( val )
+            elif type(param) is tuple:
+                if self.debug: print 'Remember is remembering: ' + param[0]
+                self._rememberList[param[0]].append( param[1](self) )
+
+
+
+class ProjectedGradient(Minimize, Remember):
+    name = 'Projected Gradient'
+
+    maxIterCG = 10
+    tolCG = 1e-3
+
+    lower = -np.inf
+    upper = np.inf
+
+    def __init__(self,**kwargs):
+        super(ProjectedGradient, self).__init__(**kwargs)
+
+        self.stoppers.append(StoppingCriteria.bindingSet)
+        self.stoppersLS.append(StoppingCriteria.bindingSet_LS)
+
+        self.printers.extend([ IterationPrinters.itType, IterationPrinters.aSet, IterationPrinters.bSet, IterationPrinters.comment ])
+
+
+    def _startup(self, x0):
+        # ensure bound vectors are the same size as the model
+        if type(self.lower) is not np.ndarray:
+            self.lower = np.ones_like(x0)*self.lower
+        if type(self.upper) is not np.ndarray:
+            self.upper = np.ones_like(x0)*self.upper
+
+        self.explorePG = True
+        self.exploreCG = False
+        self.stopDoingPG = False
+
+        self._itType = 'SD'
+        self.projComment = ''
+
+        self.aSet_prev = self.activeSet(x0)
+
+    def projection(self, x):
+        """Make sure we are feasible."""
+        return np.median(np.c_[self.lower,x,self.upper],axis=1)
+
+    def activeSet(self, x):
+        """If we are on a bound"""
+        return np.logical_or(x == self.lower, x == self.upper)
+
+    def inactiveSet(self, x):
+        """The free variables."""
+        return np.logical_not(self.activeSet(x))
+
+    def bindingSet(self, x):
+        """
+            If we are on a bound and the negative gradient points away from the feasible set.
+
+            Optimality condition. (Satisfies Kuhn-Tucker) MoreToraldo91
+
+        """
+        bind_up  = np.logical_and(x == self.lower, self.g >= 0)
+        bind_low = np.logical_and(x == self.upper, self.g <= 0)
+        return np.logical_or(bind_up, bind_low)
+
+    def findSearchDirection(self):
+        self.aSet_prev = self.activeSet(self.xc)
+        allBoundsAreActive = sum(self.aSet_prev) == self.xc.size
+
+        if self.debug: print 'findSearchDirection: stopDoingPG: ', self.stopDoingPG
+        if self.debug: print 'findSearchDirection: explorePG: ', self.explorePG
+        if self.debug: print 'findSearchDirection: exploreCG: ', self.exploreCG
+        if self.debug: print 'findSearchDirection: aSet', np.sum(self.activeSet(self.xc))
+        if self.debug: print 'findSearchDirection: bSet', np.sum(self.bindingSet(self.xc))
+        if self.debug: print 'findSearchDirection: allBoundsAreActive: ', allBoundsAreActive
+
+        if self.explorePG or not self.exploreCG or allBoundsAreActive:
+            if self.debug: print 'findSearchDirection.PG: doingPG'
+            self._itType = 'SD'
+            p = -self.g
+        else:
+            if self.debug: print 'findSearchDirection.CG: doingCG'
+            # Reset the max decrease each time you do a CG iteration
+            self.f_decrease_max = -np.inf
+
+            self._itType = '.CG.'
+
+            iSet  = self.inactiveSet(self.xc)  # The inactive set (free variables)
+            bSet = self.bindingSet(self.xc)
+            shape = (self.xc.size, np.sum(iSet))
+            v = np.ones(shape[1])
+            i = np.where(iSet)[0]
+            j = np.arange(shape[1])
+            if self.debug: print 'findSearchDirection.CG: Z.shape', shape
+            Z = sp.csr_matrix((v, (i, j)), shape=shape)
+
+            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 )
+            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
+
+    def _doEndIteration_ProjectedGradient(self, xt):
+        aSet = self.activeSet(xt)
+        bSet = self.bindingSet(xt)
+
+        self.explorePG = not np.all(aSet == self.aSet_prev) # explore proximal gradient
+        self.exploreCG = np.all(aSet == bSet) # explore conjugate gradient
+
+        f_current_decrease = self.f_last - self.f
+        self.projComment = ''
+        if self._iter < 1:
+            # Note that this is reset on every CG iteration.
+            self.f_decrease_max = -np.inf
+        else:
+            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.explorePG = False
+                self.exploreCG = True
+        # implement 3.8, MoreToraldo91
+        #self.eta_2 * max_decrease where max decrease
+        # if true go to CG
+        # don't do too many steps of PG in a row.
+
+        if self.debug: print 'doEndIteration.ProjGrad, f_current_decrease: ', f_current_decrease
+        if self.debug: print 'doEndIteration.ProjGrad, f_decrease_max: ', self.f_decrease_max
+        if self.debug: print 'doEndIteration.ProjGrad, stopDoingSD: ', self.stopDoingSD
+
+class GaussNewton(Minimize, Remember):
+    name = 'Gauss Newton'
     def findSearchDirection(self):
         return Solver(self.H).solve(-self.g)
 
 
-class InexactGaussNewton(Minimize):
-    name = 'InexactGaussNewton'
+class InexactGaussNewton(Minimize, Remember):
+    name = 'Inexact Gauss Newton'
 
     maxIterCG = 10
     tolCG = 1e-5
@@ -366,8 +638,8 @@ class InexactGaussNewton(Minimize):
         return p
 
 
-class SteepestDescent(Minimize):
-    name = 'SteepestDescent'
+class SteepestDescent(Minimize, Remember):
+    name = 'Steepest Descent'
     def findSearchDirection(self):
         return -self.g
 
@@ -377,9 +649,9 @@ if __name__ == '__main__':
     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')
+    # 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])

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]
@@ -331,7 +331,7 @@ class BaseMesh(object):
         return locals()
     nEz = property(**nEz())
 
-    def nE():
+    def nEv():
         doc = """
         Total number of edges in each direction
 
@@ -346,6 +346,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 +403,7 @@ class BaseMesh(object):
         return locals()
     nFz = property(**nFz())
 
-    def nF():
+    def nFv():
         doc = """
         Total number of faces in each direction
 
@@ -406,6 +418,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 +443,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 +464,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

@@ -110,6 +110,12 @@ class Cyl1DMesh(object):
         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

SimPEG/mesh/DiffOperators.py

@@ -1,11 +1,6 @@
 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
 
 
 def checkBC(bc):
@@ -39,11 +34,6 @@ def ddxCellGrad(n, bc):
     return D
 
 
-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")
-
-
 class DiffOperators(object):
     """
         Class creates the differential operators that you need!
@@ -107,6 +97,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.
@@ -182,7 +204,6 @@ class DiffOperators(object):
         return locals()
     cellGradx = property(**cellGradx())
 
-
     def cellGrady():
         doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
         def fget(self):
@@ -203,8 +224,6 @@ class DiffOperators(object):
         return locals()
     cellGrady = property(**cellGrady())
 
-
-
     def cellGradz():
         doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
         def fget(self):
@@ -222,7 +241,6 @@ class DiffOperators(object):
         return locals()
     cellGradz = property(**cellGradz())
 
-
     def edgeCurl():
         doc = "Construct the 3D curl operator."
 
@@ -265,6 +283,8 @@ class DiffOperators(object):
     _edgeCurl = None
     edgeCurl = property(**edgeCurl())
 
+    # --------------- Averaging ---------------------
+
     def aveF2CC():
         doc = "Construct the averaging operator on cell faces to cell centers."
 
@@ -274,10 +294,10 @@ 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])),
+                    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])),
+                    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
@@ -295,10 +315,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 +336,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

@@ -396,7 +396,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/setup.py

@@ -0,0 +1,8 @@
+import os
+print 'Compiling TriSolve.'
+os.system('f2py -c utils/TriSolve.f -m TriSolve')
+print 'TriSolve Compiled! yay.'
+print 'Moving TriSolve into Utils.'
+os.system('mv TriSolve.so utils/TriSolve.so')
+print 'Thats it. Well Done Computer.'
+

SimPEG/tests/TestUtils.py

@@ -5,11 +5,17 @@ from SimPEG.utils import mkvc, sdiag
 from SimPEG import utils
 from SimPEG.mesh import TensorMesh, LogicallyOrthogonalMesh
 import numpy as np
+import scipy.sparse as sp
 import unittest
 import inspect
 
-happiness = ['The test be workin!', 'You get a gold star!', 'Yay passed!', 'Happy little convergence test!', 'That was easy!', 'Testing is important.', 'You are awesome.', 'Go Test Go!', 'Once upon a time, a happy little test passed.', 'And then everyone was happy.']
-sadness = ['No gold star for you.','Try again soon.','Thankfully, persistence is a great substitute for talent.','It might be easier to call this a feature...','Coffee break?', 'Boooooooo  :(', 'Testing is important. Do it again.']
+try:
+    import getpass
+    name = getpass.getuser()[0].upper() + getpass.getuser()[1:]
+except Exception, e:
+    name = 'You'
+happiness = ['The test be workin!', 'You get a gold star!', 'Yay passed!', 'Happy little convergence test!', 'That was easy!', 'Testing is important.', 'You are awesome.', 'Go Test Go!', 'Once upon a time, a happy little test passed.', 'And then everyone was happy.','Not just a pretty face '+name,'You deserve a pat on the back!','Well done '+name+'!', 'Awesome, '+name+', just awesome.']
+sadness = ['No gold star for you.','Try again soon.','Thankfully, persistence is a great substitute for talent.','It might be easier to call this a feature...','Coffee break?', 'Boooooooo  :(', 'Testing is important. Do it again.',"Did you put your clever trousers on today?",'Just think about a dancing dinosaur and life will get better!','You had so much promise '+name+', oh well...', name.upper()+' ERROR!','Get on it '+name+'!', 'You break it, you fix it.']
 
 class OrderTest(unittest.TestCase):
     """
@@ -174,14 +180,14 @@ def Rosenbrock(x, return_g=True, return_H=True):
 
     f = 100*(x[1]-x[0]**2)**2+(1-x[0])**2
     g = np.array([2*(200*x[0]**3-200*x[0]*x[1]+x[0]-1), 200*(x[1]-x[0]**2)])
-    H = np.array([[-400*x[1]+1200*x[0]**2+2, -400*x[0]], [-400*x[0], 200]])
+    H = sp.csr_matrix(np.array([[-400*x[1]+1200*x[0]**2+2, -400*x[0]], [-400*x[0], 200]]))
 
     out = (f,)
     if return_g:
         out += (g,)
     if return_H:
         out += (H,)
-    return out
+    return out if len(out) > 1 else out[0]
 
 def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None):
     """
@@ -269,6 +275,28 @@ def checkDerivative(fctn, x0, num=7, plotIt=True, dx=None):
     return passTest
 
 
+
+def getQuadratic(A, b):
+    """
+        Given A and b, this returns a quadratic, Q
+
+        .. math::
+
+            \mathbf{Q( x ) = 0.5 x A x + b x}
+    """
+    def Quadratic(x, return_g=True, return_H=True):
+        f = 0.5 * x.dot( A.dot(x)) + b.dot( x )
+        out = (f,)
+        if return_g:
+            g = A.dot(x) + b
+            out += (g,)
+        if return_H:
+            H = A
+            out += (H,)
+        return out if len(out) > 1 else out[0]
+    return Quadratic
+
+
 if __name__ == '__main__':
 
     def simplePass(x):

SimPEG/tests/__init__.py

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

SimPEG/tests/api_TestResults.rst

@@ -1,355 +0,0 @@
-.. _api_TestResults:
-
-.. raw:: html
-    <style type="text/css" media="screen">
-    body        { font-family: verdana, arial, helvetica, sans-serif; font-size: 80%; }
-    table       { font-size: 100%; }
-    pre         { }
-    
-    /* -- heading ---------------------------------------------------------------------- */
-    h1 {
-        font-size: 16pt;
-        color: gray;
-    }
-    .heading {
-        margin-top: 0ex;
-        margin-bottom: 1ex;
-    }
-    
-    .heading .attribute {
-        margin-top: 1ex;
-        margin-bottom: 0;
-    }
-    
-    .heading .description {
-        margin-top: 4ex;
-        margin-bottom: 6ex;
-    }
-    
-    /* -- css div popup ------------------------------------------------------------------------ */
-    a.popup_link {
-    }
-    
-    a.popup_link:hover {
-        color: red;
-    }
-    
-    .popup_window {
-        display: none;
-        position: relative;
-        left: 0px;
-        top: 0px;
-        /*border: solid #627173 1px; */
-        padding: 10px;
-        background-color: #E6E6D6;
-        font-family: "Lucida Console", "Courier New", Courier, monospace;
-        text-align: left;
-        font-size: 8pt;
-        width: 500px;
-    }
-    
-    }
-    /* -- report ------------------------------------------------------------------------ */
-    #show_detail_line {
-        margin-top: 3ex;
-        margin-bottom: 1ex;
-    }
-    #result_table {
-        width: 80%;
-        border-collapse: collapse;
-        border: 1px solid #777;
-    }
-    #header_row {
-        font-weight: bold;
-        color: white;
-        background-color: #777;
-    }
-    #result_table td {
-        border: 1px solid #777;
-        padding: 2px;
-    }
-    #total_row  { font-weight: bold; }
-    .passClass  { background-color: #6c6; }
-    .failClass  { background-color: #c60; }
-    .errorClass { background-color: #c00; }
-    .passCase   { color: #6c6; }
-    .failCase   { color: #c60; font-weight: bold; }
-    .errorCase  { color: #c00; font-weight: bold; }
-    .hiddenRow  { display: none; }
-    .testcase   { margin-left: 2em; }
-    
-    
-    /* -- ending ---------------------------------------------------------------------- */
-    #ending {
-    }
-    
-    </style>
-    
-    <body>
-    <script language="javascript" type="text/javascript"><!--
-    output_list = Array();
-    
-    /* level - 0:Summary; 1:Failed; 2:All */
-    function showCase(level) {
-        trs = document.getElementsByTagName("tr");
-        for (var i = 0; i < trs.length; i++) {
-            tr = trs[i];
-            id = tr.id;
-            if (id.substr(0,2) == 'ft') {
-                if (level < 1) {
-                    tr.className = 'hiddenRow';
-                }
-                else {
-                    tr.className = '';
-                }
-            }
-            if (id.substr(0,2) == 'pt') {
-                if (level > 1) {
-                    tr.className = '';
-                }
-                else {
-                    tr.className = 'hiddenRow';
-                }
-            }
-        }
-    }
-    
-    
-    function showClassDetail(cid, count) {
-        var id_list = Array(count);
-        var toHide = 1;
-        for (var i = 0; i < count; i++) {
-            tid0 = 't' + cid.substr(1) + '.' + (i+1);
-            tid = 'f' + tid0;
-            tr = document.getElementById(tid);
-            if (!tr) {
-                tid = 'p' + tid0;
-                tr = document.getElementById(tid);
-            }
-            id_list[i] = tid;
-            if (tr.className) {
-                toHide = 0;
-            }
-        }
-        for (var i = 0; i < count; i++) {
-            tid = id_list[i];
-            if (toHide) {
-                document.getElementById('div_'+tid).style.display = 'none'
-                document.getElementById(tid).className = 'hiddenRow';
-            }
-            else {
-                document.getElementById(tid).className = '';
-            }
-        }
-    }
-    
-    
-    function showTestDetail(div_id){
-        var details_div = document.getElementById(div_id)
-        var displayState = details_div.style.display
-        // alert(displayState)
-        if (displayState != 'block' ) {
-            displayState = 'block'
-            details_div.style.display = 'block'
-        }
-        else {
-            details_div.style.display = 'none'
-        }
-    }
-    
-    
-    function html_escape(s) {
-        s = s.replace(/&/g,'&amp;');
-        s = s.replace(/</g,'&lt;');
-        s = s.replace(/>/g,'&gt;');
-        return s;
-    }
-    
-    /* obsoleted by detail in <div>
-    function showOutput(id, name) {
-        var w = window.open("", //url
-                        name,
-                        "resizable,scrollbars,status,width=800,height=450");
-        d = w.document;
-        d.write("<pre>");
-        d.write(html_escape(output_list[id]));
-        d.write("\n");
-        d.write("<a href='javascript:window.close()'>close</a>\n");
-        d.write("</pre>\n");
-        d.close();
-    }
-    */
-    --></script>
-    
-    <div class='heading'>
-    <h1>Test Report</h1>
-    <p class='attribute'><strong>Start Time:</strong> 2013-11-05 15:24:44</p>
-    <p class='attribute'><strong>Duration:</strong> 0:00:00.007500</p>
-    <p class='attribute'><strong>Status:</strong> Pass 22</p>
-    
-    <p class='description'>This demonstrates the report output by Prasanna.Yelsangikar.</p>
-    </div>
-    
-    
-    
-    <p id='show_detail_line'>Show
-    <a href='javascript:showCase(0)'>Summary</a>
-    <a href='javascript:showCase(1)'>Failed</a>
-    <a href='javascript:showCase(2)'>All</a>
-    </p>
-    <table id='result_table'>
-    <colgroup>
-    <col align='left' />
-    <col align='right' />
-    <col align='right' />
-    <col align='right' />
-    <col align='right' />
-    <col align='right' />
-    </colgroup>
-    <tr id='header_row'>
-        <td>Test Group/Test case</td>
-        <td>Count</td>
-        <td>Pass</td>
-        <td>Fail</td>
-        <td>Error</td>
-        <td>View</td>
-    </tr>
-    
-    <tr class='passClass'>
-        <td>test_basemesh.TestBaseMesh</td>
-        <td>11</td>
-        <td>11</td>
-        <td>0</td>
-        <td>0</td>
-        <td><a href="javascript:showClassDetail('c1',11)">Detail</a></td>
-    </tr>
-    
-    <tr id='pt1.1' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_meshDimensions</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt1.2' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_nc</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt1.3' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_nc_xyz</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt1.4' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_ne</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt1.5' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_nf</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt1.6' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_numbers</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt1.7' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_r_CC_M</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt1.8' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_r_E_M</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt1.9' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_r_E_V</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt1.10' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_r_F_M</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt1.11' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_r_F_V</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr class='passClass'>
-        <td>test_basemesh.TestMeshNumbers2D</td>
-        <td>11</td>
-        <td>11</td>
-        <td>0</td>
-        <td>0</td>
-        <td><a href="javascript:showClassDetail('c2',11)">Detail</a></td>
-    </tr>
-    
-    <tr id='pt2.1' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_meshDimensions</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt2.2' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_nc</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt2.3' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_nc_xyz</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt2.4' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_ne</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt2.5' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_nf</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt2.6' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_numbers</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt2.7' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_r_CC_M</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt2.8' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_r_E_M</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt2.9' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_r_E_V</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt2.10' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_r_F_M</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='pt2.11' class='hiddenRow'>
-        <td class='none'><div class='testcase'>test_mesh_r_F_V</div></td>
-        <td colspan='5' align='center'>pass</td>
-    </tr>
-    
-    <tr id='total_row'>
-        <td>Total</td>
-        <td>22</td>
-        <td>22</td>
-        <td>0</td>
-        <td>0</td>
-        <td>&nbsp;</td>
-    </tr>
-    </table>
-    

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_Solver.py

@@ -58,23 +58,39 @@ class TestSolver(unittest.TestCase):
         x = solve.solve(rhs)
         self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
 
-    def test_directLower_1(self):
+    def test_directLower_1_python(self):
         AL = sparse.tril(self.A)
-        solve = Solver(AL, doDirect=True, flag='L', options={})
+        solve = Solver(AL, doDirect=True, flag='L', options={'backend':'python'})
         e = np.ones(self.M.nC)
         rhs = AL.dot(e)
         x = solve.solve(rhs)
         self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
 
-    def test_directLower_M(self):
+    def test_directLower_M_python(self):
         AL = sparse.tril(self.A)
-        solve = Solver(AL, doDirect=True, flag='L', options={})
+        solve = Solver(AL, doDirect=True, flag='L', options={'backend':'python'})
         e = np.ones((self.M.nC,numRHS))
         rhs = AL.dot(e)
         x = solve.solve(rhs)
+
+    def test_directLower_1_fortran(self):
+        AL = sparse.tril(self.A)
+        solve = Solver(AL, doDirect=True, flag='L', options={'backend':'fortran'})
+        e = np.ones(self.M.nC)
+        rhs = AL.dot(e)
+        x = solve.solve(rhs)
         self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
 
-    def test_directUpper_1(self):
+    def test_directLower_M_fortran(self):
+        AL = sparse.tril(self.A)
+        solve = Solver(AL, doDirect=True, flag='L', options={'backend':'fortran'})
+        e = np.ones((self.M.nC,numRHS))
+        rhs = AL.dot(e)
+        x = solve.solve(rhs)
+        self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
+        self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
+
+    def test_directUpper_1_python(self):
         AU = sparse.triu(self.A)
         solve = Solver(AU, doDirect=True, flag='U', options={})
         e = np.ones(self.M.nC)
@@ -82,7 +98,7 @@ class TestSolver(unittest.TestCase):
         x = solve.solve(rhs)
         self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
 
-    def test_directUpper_M(self):
+    def test_directUpper_M_python(self):
         AU = sparse.triu(self.A)
         solve = Solver(AU, doDirect=True, flag='U', options={})
         e = np.ones((self.M.nC,numRHS))
@@ -90,6 +106,23 @@ class TestSolver(unittest.TestCase):
         x = solve.solve(rhs)
         self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
 
+
+    def test_directUpper_1_fortran(self):
+        AU = sparse.triu(self.A)
+        solve = Solver(AU, doDirect=True, flag='U', options={'backend':'fortran'})
+        e = np.ones(self.M.nC)
+        rhs = AU.dot(e)
+        x = solve.solve(rhs)
+        self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
+
+    def test_directUpper_M_fortran(self):
+        AU = sparse.triu(self.A)
+        solve = Solver(AU, doDirect=True, flag='U', options={'backend':'fortran'})
+        e = np.ones((self.M.nC,numRHS))
+        rhs = AU.dot(e)
+        x = solve.solve(rhs)
+        self.assertTrue(np.linalg.norm(e-x,np.inf) < TOL, True)
+
     def test_directDiagonal_1(self):
         AD = sdiag(self.A.diagonal())
         solve = Solver(AD, doDirect=True, flag='D', options={})

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_operators.py

@@ -155,6 +155,109 @@ 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_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()
+
 
 if __name__ == '__main__':
     unittest.main()

SimPEG/tests/test_optimizers.py

@@ -0,0 +1,54 @@
+import unittest
+from SimPEG import Solver
+from SimPEG.mesh import TensorMesh
+from SimPEG.utils import sdiag
+import numpy as np
+import scipy.sparse as sp
+from SimPEG import inverse
+from SimPEG.tests import getQuadratic, Rosenbrock
+
+TOL = 1e-2
+
+class TestOptimizers(unittest.TestCase):
+
+    def setUp(self):
+        self.A = sp.identity(2).tocsr()
+        self.b = np.array([-5,-5])
+
+    def test_GN_Rosenbrock(self):
+        GN = inverse.GaussNewton()
+        xopt = GN.minimize(Rosenbrock,np.array([0,0]))
+        x_true = np.array([1.,1.])
+        print 'xopt: ', xopt
+        print 'x_true: ', x_true
+        self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True)
+
+    def test_GN_quadratic(self):
+        GN = inverse.GaussNewton()
+        xopt = GN.minimize(getQuadratic(self.A,self.b),np.array([0,0]))
+        x_true = np.array([5.,5.])
+        print 'xopt: ', xopt
+        print 'x_true: ', x_true
+        self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True)
+
+    def test_ProjGradient_quadraticBounded(self):
+        PG = inverse.ProjectedGradient()
+        PG.lower, PG.upper = -2, 2
+        xopt = PG.minimize(getQuadratic(self.A,self.b),np.array([0,0]))
+        x_true = np.array([2.,2.])
+        print 'xopt: ', xopt
+        print 'x_true: ', x_true
+        self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True)
+
+    def test_ProjGradient_quadratic1Bound(self):
+        myB = np.array([-5,1])
+        PG = inverse.ProjectedGradient()
+        PG.lower, PG.upper = -2, 2
+        xopt = PG.minimize(getQuadratic(self.A,myB),np.array([0,0]))
+        x_true = np.array([2.,-1.])
+        print 'xopt: ', xopt
+        print 'x_true: ', x_true
+        self.assertTrue(np.linalg.norm(xopt-x_true,2) < TOL, True)
+
+if __name__ == '__main__':
+    unittest.main()

SimPEG/utils/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,39 @@
+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)
+
+    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/Solver.py

@@ -1,7 +1,22 @@
 import numpy as np
 import scipy.sparse as sparse
 import scipy.sparse.linalg as linalg
+from SimPEG.utils import mkvc
 
+DEFAULTS = {'direct':'scipy', 'forward':'fortran', 'backward':'fortran', 'diagonal':'python'}
+
+try:
+    import TriSolve
+except Exception, e:
+    try:
+        import os
+        # Note: this may not work from SublimeText, if that is the case, just run the command in your shell.
+        os.system('f2py -c TriSolve.f -m TriSolve')
+        import TriSolve
+    except Exception, e:
+        print 'Warning: Python backend is being used for solver. Run setup.py from the command line.'
+        DEFAULTS['forward'] = 'python'
+        DEFAULTS['backward'] = 'python'
 
 class Solver(object):
     """
@@ -55,11 +70,11 @@ class Solver(object):
         elif self.flag is None and not self.doDirect:
             return self.solveIter(b, **self.options)
         elif self.flag == 'U':
-            return self.solveBackward(b)
+            return self.solveBackward(b, **self.options)
         elif self.flag == 'L':
-            return self.solveForward(b)
+            return self.solveForward(b, **self.options)
         elif self.flag == 'D':
-            return self.solveDiagonal(b)
+            return self.solveDiagonal(b, **self.options)
         else:
             raise Exception('Unknown flag.')
         pass
@@ -69,7 +84,7 @@ class Solver(object):
         del self.dsolve
         self.dsolve = None
 
-    def solveDirect(self, b, factorize=False, backend='scipy'):
+    def solveDirect(self, b, factorize=False, backend=None):
         """
             Use solve instead of this interface.
 
@@ -78,6 +93,8 @@ class Solver(object):
             :rtype: numpy.ndarray
             :return: x
         """
+        if backend is None: backend = DEFAULTS['direct']
+
         assert np.shape(self.A)[1] == np.shape(b)[0], 'Dimension mismatch'
 
         if factorize and self.dsolve is None:
@@ -104,7 +121,7 @@ class Solver(object):
     def solveIter(self, b, M=None, iterSolver='CG'):
         pass
 
-    def solveBackward(self, b, backend='python'):
+    def solveBackward(self, b, backend=None):
         """
             Use solve instead of this interface.
 
@@ -114,21 +131,29 @@ class Solver(object):
             :rtype: numpy.ndarray
             :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)
         vals = self.A.data
         rowptr = self.A.indptr
         colind = self.A.indices
-        x = np.empty_like(b)   # empty() is faster than zeros().
-        for i in reversed(xrange(self.A.shape[0])):
-            ith_row = vals[rowptr[i] : rowptr[i+1]]
-            cols = colind[rowptr[i] : rowptr[i+1]]
-            x_vals = x[cols]
-            x[i] = (b[i] - np.dot(ith_row[1:], x_vals[1:])) / ith_row[0]
+        if backend == 'fortran':
+            if len(b.shape) == 1 or b.shape[1] == 1:
+                x = TriSolve.backward(vals, rowptr, colind, b, self.A.data.size, b.size, 1)
+                x = mkvc(x)
+            else:
+                x = TriSolve.backward(vals, rowptr, colind, b, self.A.data.size, b.shape[0], b.shape[1])
+        elif backend == 'python':
+            x = np.empty_like(b)   # empty() is faster than zeros().
+            for i in reversed(xrange(self.A.shape[0])):
+                ith_row = vals[rowptr[i] : rowptr[i+1]]
+                cols = colind[rowptr[i] : rowptr[i+1]]
+                x_vals = x[cols]
+                x[i] = (b[i] - np.dot(ith_row[1:], x_vals[1:])) / ith_row[0]
         return x
 
-    def solveForward(self, b, backend='python'):
+    def solveForward(self, b, backend=None):
         """
             Use solve instead of this interface.
 
@@ -138,21 +163,29 @@ class Solver(object):
             :rtype: numpy.ndarray
             :return: x
         """
+        if backend is None: backend = DEFAULTS['forward']
         if type(self.A) is not sparse.csr.csr_matrix:
             from scipy.sparse import csr_matrix
             self.A = csr_matrix(self.A)
         vals = self.A.data
         rowptr = self.A.indptr
         colind = self.A.indices
-        x = np.empty_like(b)   # empty() is faster than zeros().
-        for i in xrange(self.A.shape[0]):
-            ith_row = vals[rowptr[i] : rowptr[i+1]]
-            cols = colind[rowptr[i] : rowptr[i+1]]
-            x_vals = x[cols]
-            x[i] = (b[i] - np.dot(ith_row[:-1], x_vals[:-1])) / ith_row[-1]
+        if backend == 'fortran':
+            if len(b.shape) == 1 or b.shape[1] == 1:
+                x = TriSolve.forward(vals, rowptr, colind, b, self.A.data.size, b.size, 1)
+                x = mkvc(x)
+            else:
+                x = TriSolve.forward(vals, rowptr, colind, b, self.A.data.size, b.shape[0], b.shape[1])
+        elif backend == 'python':
+            x = np.empty_like(b)   # empty() is faster than zeros().
+            for i in xrange(self.A.shape[0]):
+                ith_row = vals[rowptr[i] : rowptr[i+1]]
+                cols = colind[rowptr[i] : rowptr[i+1]]
+                x_vals = x[cols]
+                x[i] = (b[i] - np.dot(ith_row[:-1], x_vals[:-1])) / ith_row[-1]
         return x
 
-    def solveDiagonal(self, b, backend='python'):
+    def solveDiagonal(self, b, backend=None):
         """
             Use solve instead of this interface.
 
@@ -162,6 +195,8 @@ class Solver(object):
             :rtype: numpy.ndarray
             :return: x
         """
+        if backend is None: backend = DEFAULTS['diagonal']
+
         diagA = self.A.diagonal()
         if len(b.shape) == 1 or b.shape[1] == 1:
             # Just one RHS
@@ -205,3 +240,20 @@ if __name__ == '__main__':
     print np.linalg.norm(e-x,np.inf)
 
 
+    n = 6000
+    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)
+    pSolve = Solver(A,flag='L',options={'backend':'python'});
+    fSolve = Solver(A,flag='L',options={'backend':'fortran'})
+    tic = time()
+    x = pSolve.solve(b)
+    toc = time() - tic
+    print 'Error Forward Python  = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc
+    tic = time()
+    x = fSolve.solve(b)
+    toc = time() - tic
+    print 'Error Forward Fortran = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc

SimPEG/utils/TriSolve.f

@@ -0,0 +1,64 @@
+c File TriSolve.f
+       subroutine forward(al, ial, jal, b, nv, n, nRHS, x)
+       double precision al(nv)
+       integer ial(n+1)
+       integer jal(nv)
+       double precision b(n,nRHS)
+       double precision x(n,nRHS)
+       integer nv
+       integer n
+       integer nRHS
+       integer rhs
+cf2py  intent(in) :: al
+cf2py  intent(in) :: ial
+cf2py  intent(in) :: jal
+cf2py  intent(in) :: b
+cf2py  intent(in) :: nv
+cf2py  intent(in) :: n
+cf2py  intent(in) :: nRHS
+cf2py  intent(out) :: x
+       real ( kind = 8 ) t
+
+       do rhs = 1, nRHS
+          do k = 1, n
+             t = b(k,rhs)
+             do j = ial(k)+1, ial(k+1)
+                t = t - al(j) * x(jal(j)+1,rhs)
+             end do
+             x(k,rhs) = t/al(ial(k+1))
+          end do
+       end do
+       end subroutine forward
+
+
+       subroutine backward(au,iau, jau, b, nv, n, nRHS, x)
+       double precision au(nv)
+       integer iau(n+1)
+       integer jau(nv)
+       double precision b(n,nRHS)
+       double precision x(n,nRHS)
+       integer nv
+       integer n
+       integer nRHS
+       integer rhs
+cf2py  intent(in) :: au
+cf2py  intent(in) :: iau
+cf2py  intent(in) :: jau
+cf2py  intent(in) :: b
+cf2py  intent(in) :: nv
+cf2py  intent(in) :: n
+cf2py  intent(in) :: nRHS
+cf2py  intent(out) :: x
+       real ( kind = 8 ) t
+
+       do rhs = 1, nRHS
+          do k = n, 1, -1
+             t = b(k,rhs)
+             do j = iau(k)+1, iau(k+1)
+                t = t - au(j) * x(jau(j)+1,rhs)
+             end do
+             x(k,rhs) = t/au(iau(k)+1)
+          end do
+       end do
+
+       end subroutine backward

SimPEG/utils/__init__.py

@@ -3,9 +3,60 @@ import sputils
 import lomutils
 import interputils
 import ModelBuilder
-import Solver
-from Solver import Solver
 from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind
-from sputils import spzeros, kron3, speye, sdiag
+from sputils import spzeros, kron3, speye, sdiag, ddx, av
 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
+
+def setKwargs(obj, **kwargs):
+    """Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist."""
+    for attr in kwargs:
+        if hasattr(obj, attr):
+            setattr(obj, attr, kwargs[attr])
+        else:
+            raise Exception('%s attr is not recognized' % attr)
+
+def printTitles(obj, printers, name='Print Titles', pad=''):
+    titles = ''
+    widths = 0
+    for printer in printers:
+        titles += ('{:^%i}'%printer['width']).format(printer['title']) + ''
+        widths += printer['width']
+    print pad + "{0} {1} {0}".format('='*((widths-1-len(name))/2), name)
+    print pad + titles
+    print pad + "%s" % '-'*widths
+
+def printLine(obj, printers, pad=''):
+    values = ''
+    for printer in printers:
+        values += ('{:^%i}'%printer['width']).format(printer['format'] % printer['value'](obj))
+    print pad + values
+
+def checkStoppers(obj, stoppers):
+    # check stopping rules
+    optimal = []
+    critical = []
+    for stopper in stoppers:
+        l = stopper['left'](obj)
+        r = stopper['right'](obj)
+        if stopper['stopType'] == 'optimal':
+            optimal.append(l <= r)
+        if stopper['stopType'] == 'critical':
+            critical.append(l <= r)
+
+    if obj.debug: print 'checkStoppers.optimal: ', optimal
+    if obj.debug: print 'checkStoppers.critical: ', critical
+
+    return (len(optimal)>0 and all(optimal)) | (len(critical)>0 and any(critical))
+
+def printStoppers(obj, stoppers, pad='', stop='STOP!', done='DONE!'):
+    print pad + "%s%s%s" % ('-'*25,stop,'-'*25)
+    for stopper in stoppers:
+        l = stopper['left'](obj)
+        r = stopper['right'](obj)
+        print pad + stopper['str'] % (l<=r,l,r)
+    print pad + "%s%s%s" % ('-'*25,done,'-'*25)

SimPEG/utils/ipythonUtils.py

@@ -0,0 +1,28 @@
+from tempfile import NamedTemporaryFile
+import matplotlib.pyplot as plt
+from matplotlib import animation
+
+# http://jakevdp.github.io/blog/2013/05/12/embedding-matplotlib-animations/
+# http://www.renevolution.com/how-to-install-ffmpeg-on-mac-os-x/
+
+VIDEO_TAG = """<video controls loop>
+ <source src="data:video/x-m4v;base64,{0}" type="video/mp4">
+ Your browser does not support the video tag.
+</video>"""
+
+def anim_to_html(anim):
+    if not hasattr(anim, '_encoded_video'):
+        with NamedTemporaryFile(suffix='.mp4') as f:
+            anim.save(f.name, fps=20, extra_args=['-vcodec', 'libx264', '-pix_fmt', 'yuv420p'])
+            video = open(f.name, "rb").read()
+        anim._encoded_video = video.encode("base64")
+
+    return VIDEO_TAG.format(anim._encoded_video)
+
+def display_animation(anim):
+    plt.close(anim._fig)
+    return anim_to_html(anim)
+
+animation.Animation._repr_html_ = display_animation
+
+easyAnimate = animation.FuncAnimation

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,13 @@ 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 cell-centres to nodes."""
+    return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr")

docs/requirements.txt

@@ -1,2 +1,3 @@
 numpy
 pypubsub
+ipython