Merge branch 'master' of https://github.com/simpeg/simpeg into cylClean

SimPEG/Mesh/View.py

@@ -455,7 +455,6 @@ class TensorView(object):
             ax.set_zlabel('x3')
 
         ax.grid(True)
-        ax.hold(False)
         if showIt: plt.show()
 
     def slicer(mesh, var, imageType='CC', normal='z', index=0, ax=None, clim=None):

SimPEG/Model.py

@@ -58,11 +58,13 @@ class BaseModel(object):
     def example(self):
         return np.random.rand(self.nP)
 
-    def test(self, m=None):
+    def test(self, m=None, **kwargs):
         print 'Testing the %s Class!' % self.__class__.__name__
         if m is None:
             m = self.example()
-        return checkDerivative(lambda m : [self.transform(m), self.transformDeriv(m)], m, plotIt=False)
+            if 'plotIt' not in kwargs:
+                kwargs['plotIt'] = False
+        return checkDerivative(lambda m : [self.transform(m), self.transformDeriv(m)], m, **kwargs)
 
 class BaseNonLinearModel(object):
     """
@@ -308,7 +310,7 @@ class ActiveModel(BaseModel):
             indActive = z
         self.indActive = indActive
         self.indInactive = np.logical_not(indActive)
-        if type(valInactive) in [float, int, long]:
+        if Utils.isScalar(valInactive):
             valInactive = np.ones(self.nC)*float(valInactive)
 
         valInactive[self.indActive] = 0

SimPEG/Solver.py

@@ -1,7 +1,7 @@
 import numpy as np
 import scipy.sparse as sp
 import scipy.sparse.linalg as linalg
-from Utils.matutils import mkvc, sdiag
+from Utils import mkvc, sdiag
 import warnings
 
 DEFAULTS = {'direct':'scipy', 'iter':'scipy', 'triangular':'fortran', 'diagonal':'python'}
@@ -214,12 +214,18 @@ class Solver(object):
     def solveIter(self, b, backend=None, M=None, iterSolver='CG', tol=1e-6, maxIter=50):
         if backend is None: backend = DEFAULTS['iter']
 
-        algorithms = {'CG':sp.linalg.cg}
+        algorithms = {'CG':sp.linalg.cg, 'QMR':sp.linalg.qmr}
         assert iterSolver in algorithms, "iterSolver must be 'CG', or implement it yourself and add it here!"
         alg = algorithms[iterSolver]
 
+        if iterSolver == 'CG':
+            opts = {'M':M}
+        elif iterSolver == 'QMR':
+            #TODO: make preconditioner better.
+            opts = {'M1':np.sqrt(M), 'M2':np.sqrt(M)}
+
         if len(b.shape) == 1 or b.shape[1] == 1:
-            x, self.info = alg(self.A, b, M=M, tol=tol, maxiter=maxIter)
+            x, self.info = alg(self.A, b,  tol=tol, maxiter=maxIter)
         else:
             x = np.empty_like(b)
             for i in range(b.shape[1]):
@@ -381,3 +387,17 @@ if __name__ == '__main__':
     toc = time() - tic
     print x
     print 'Error CG  = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc, 'Info: ', iSolve.info
+
+
+
+    A = -D*D.T
+    A[0,0] *= 10 # remove the constant null space from the matrix
+    e = np.ones(M.nC)
+    b = A.dot(e)
+
+    iSolve = Solver(A, doDirect=False, options={'iterSolver': 'QMR', 'M':'J'})
+    tic = time()
+    x = iSolve.solve(b)
+    toc = time() - tic
+    print x
+    print 'Error QMR  = ', np.linalg.norm(x-e, np.inf), 'Time: ', toc, 'Info: ', iSolve.info

SimPEG/Survey.py

@@ -159,7 +159,7 @@ class BaseSurvey(object):
         return self.mtrue is not None
 
 
-    #TODO: Move this to the data class?
+    #TODO: Move this to the survey class?
     # @property
     # def phi_d_target(self):
     #     """
@@ -178,8 +178,8 @@ class BaseSurvey(object):
     #     self._phi_d_target = value
 
 
-class BaseRxList(object):
-    """SimPEG Receiver List Object"""
+class BaseRx(object):
+    """SimPEG Receiver Object"""
 
     locs = None   #: Locations (nRx x 3)
 
@@ -207,12 +207,15 @@ class BaseTx(object):
     loc    = None #: Location [x,y,z]
 
     rxList = None #: SimPEG Receiver List
-    rxListPair = BaseRxList
+    rxPair = BaseRx
 
     knownTxTypes = None #: Set this to a list of strings to ensure that txType is known
 
     def __init__(self, loc, txType, rxList, **kwargs):
-        assert isinstance(rxList, self.rxListPair), 'rxList must be a %s'%self.rxListPair.__name__
+        assert type(rxList) is list, 'rxList must be a list'
+        for rx in rxList:
+            assert isinstance(rx, self.rxPair), 'rxList must be a %s'%self.rxListPair.__name__
+        assert len(set(rxList)) == len(rxList), 'The rxList must be unique'
 
         self.loc    = loc
         self.txType = txType

SimPEG/Utils/SolverUtils.py

@@ -0,0 +1,77 @@
+import numpy as np
+from matutils import mkvc
+import warnings
+
+def DSolverWrap(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6):
+
+    def __init__(self, A, **kwargs):
+        self.A = A.tocsc()
+        self.kwargs = kwargs
+        if factorize:
+            self.solver = fun(self.A, **kwargs)
+
+    def solve(self, b):
+        if len(b.shape) == 1 or b.shape[1] == 1:
+            b = b.flatten()
+            # Just one RHS
+            if factorize:
+                X = self.solver.solve(b, **self.kwargs)
+            else:
+                X = fun(self.A, b, **self.kwargs)
+        else: # Multiple RHSs
+            X = np.empty_like(b)
+            for i in range(b.shape[1]):
+                if factorize:
+                    X[:,i] = self.solver.solve(b[:,i])
+                else:
+                    X[:,i] = fun(self.A, b[:,i], **self.kwargs)
+
+        if checkAccuracy:
+            nrm = np.linalg.norm(mkvc(self.A*X - b)) / np.linalg.norm(mkvc(b))
+            if nrm > accuracyTol:
+                msg = '### SolverWarning ###: Accuracy on solve is above tolerance: %e > %e' % (nrm, accuracyTol)
+                print msg
+                warnings.warn(msg, RuntimeWarning)
+        return X
+
+    return type(fun.__name__, (object,), {"__init__": __init__, "solve": solve})
+
+
+
+def ISolverWrap(fun, checkAccuracy=True, accuracyTol=1e-5):
+
+    def __init__(self, A, **kwargs):
+        self.A = A.tocsc()
+        self.kwargs = kwargs
+
+    def solve(self, b):
+        if len(b.shape) == 1 or b.shape[1] == 1:
+            b = b.flatten()
+            # Just one RHS
+            out = fun(self.A, b, **self.kwargs)
+            if type(out) is tuple and len(out) == 2:
+                # We are dealing with scipy output with an info!
+                X = out[0]
+                self.info = out[1]
+            else:
+                X = out
+        else: # Multiple RHSs
+            X = np.empty_like(b)
+            for i in range(b.shape[1]):
+                out = fun(self.A, b[:,i], **self.kwargs)
+                if type(out) is tuple and len(out) == 2:
+                    # We are dealing with scipy output with an info!
+                    X[:,i] = out[0]
+                    self.info = out[1]
+                else:
+                    X[:,i] = out
+
+        if checkAccuracy:
+            nrm = np.linalg.norm(mkvc(self.A*X - b)) / np.linalg.norm(mkvc(b))
+            if nrm > accuracyTol:
+                msg = '### SolverWarning ###: Accuracy on solve is above tolerance: %e > %e' % (nrm, accuracyTol)
+                print msg
+                warnings.warn(msg, RuntimeWarning)
+        return X
+
+    return type(fun.__name__, (object,), {"__init__": __init__, "solve": solve})

SimPEG/Utils/__init__.py

@@ -1,9 +1,10 @@
 from matutils import *
-from meshutils import exampleLrmGrid, meshTensors
+from meshutils import exampleLrmGrid, meshTensors, points2nodes
 from lrmutils import volTetra, faceInfo, indexCube
 from interputils import interpmat
 from ipythonutils import easyAnimate as animate
 import ModelBuilder
+import SolverUtils
 
 import types
 import time
@@ -135,7 +136,8 @@ def callHooks(match, mainFirst=False):
 def dependentProperty(name, value, children, doc):
     def fget(self): return getattr(self,name,value)
     def fset(self, val):
-        if getattr(self,name,value) == val: return # it is the same!
+        if isScalar(val) and getattr(self,name,value) == val:
+            return # it is the same!
         for child in children:
             if hasattr(self, child):
                 delattr(self, child)

SimPEG/Utils/meshutils.py

@@ -53,6 +53,27 @@ def meshTensors(*args):
 
     return list(tensors) if len(tensors) > 1 else tensors[0]
 
+def points2nodes(mesh, pts):
+    """
+        Move a list of the nearest nodes to a set of points
+
+        :param simpeg.Mesh.TensorMesh mesh: The mesh
+        :param numpy.ndarray pts: Points to move}
+        :rtype: numpy.ndarray
+        :return: nodeInds
+    """
+
+    pts = np.atleast_2d(pts)
+
+    assert mesh._meshType == 'TENSOR'
+    assert pts.shape[1] == mesh.dim
+
+    nodeInds = np.empty(pts.shape[0], dtype=int)
+
+    for i, pt in enumerate(pts):
+        nodeInds[i] = ((np.tile(pt, (mesh.nN,1)) - mesh.gridN)**2).sum(axis=1).argmin()
+
+    return nodeInds
 if __name__ == '__main__':
     from SimPEG import mesh
     import matplotlib.pyplot as plt

docs/api_Utils.rst

@@ -15,6 +15,14 @@ Matrix Utilities
     :members:
     :undoc-members:
 
+
+Solver Utilities
+================
+
+.. automodule:: SimPEG.Utils.SolverUtils
+    :members:
+    :undoc-members:
+
 LRM Utilities
 =============
 

docs/index.rst

@@ -3,9 +3,11 @@
    :alt: SimPEG
    :align: center
 
-SimPEG (Simulation and Parameter Estimation in Geophysics) is a python
-package for simulation and gradient based parameter estimation in the
-context of geoscience applications.
+Simulation and Parameter Estimation in Geophysics
+*************************************************
+
+SimPEG is a python package for simulation and gradient
+based parameter estimation in the context of geoscience applications.
 
 The vision is to create a package for finite volume simulation and parameter estimation with
 applications to geophysical imaging and subsurface flow. To enable