Merge pull request #65 from simpeg/Map

Map Updates.

SimPEG/Maps.py

@@ -1,9 +1,9 @@
 import Utils, Parameters, numpy as np, scipy.sparse as sp
 from Tests import checkDerivative
 
-class BaseModel(object):
+class IdentityMap(object):
     """
-    SimPEG Model
+    SimPEG Map
 
     """
 
@@ -66,9 +66,14 @@ class BaseModel(object):
                 kwargs['plotIt'] = False
         return checkDerivative(lambda m : [self.transform(m), self.transformDeriv(m)], m, **kwargs)
 
-class BaseNonLinearModel(object):
+    def _assertMatchesPair(self, pair):
+        assert (isinstance(self, pair) or
+            isinstance(self, ComboMap) and isinstance(self.maps[0], pair)
+            ), "Mapping object must be an instance of a %s class."%(pair.__name__)
+
+class NonLinearMap(object):
     """
-    SimPEG BaseNonLinearModel
+    SimPEG NonLinearMap
 
     """
 
@@ -129,11 +134,11 @@ class BaseNonLinearModel(object):
         raise NotImplementedError('The test is not implemented.')
 
 
-class LogModel(BaseModel):
-    """SimPEG LogModel"""
+class ExpMap(IdentityMap):
+    """SimPEG ExpMap"""
 
     def __init__(self, mesh, **kwargs):
-        BaseModel.__init__(self, mesh, **kwargs)
+        IdentityMap.__init__(self, mesh, **kwargs)
 
     def transform(self, m):
         """
@@ -197,8 +202,8 @@ class LogModel(BaseModel):
         """
         return Utils.sdiag(np.exp(Utils.mkvc(m)))
 
-class Vertical1DModel(BaseModel):
-    """Vertical1DModel
+class Vertical1DMap(IdentityMap):
+    """Vertical1DMap
 
         Given a 1D vector through the last dimension
         of the mesh, this will extend to the full
@@ -206,7 +211,7 @@ class Vertical1DModel(BaseModel):
     """
 
     def __init__(self, mesh, **kwargs):
-        BaseModel.__init__(self, mesh, **kwargs)
+        IdentityMap.__init__(self, mesh, **kwargs)
 
     @property
     def nP(self):
@@ -238,7 +243,7 @@ class Vertical1DModel(BaseModel):
                     ), shape=(repNum, 1))
         return sp.kron(sp.identity(self.nP), repVec)
 
-class Mesh2Mesh(BaseModel):
+class Mesh2Mesh(IdentityMap):
     """
         Takes a model on one mesh are translates it to another mesh.
 
@@ -251,8 +256,8 @@ class Mesh2Mesh(BaseModel):
             M2 = Mesh.TensorMesh([h1,h1])
             V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50)
             v = Utils.mkvc(V)
-            modh = Model.Mesh2Mesh([M,M2])
-            modH = Model.Mesh2Mesh([M2,M])
+            modh = Maps.Mesh2Mesh([M,M2])
+            modH = Maps.Mesh2Mesh([M2,M])
             H = modH.transform(v)
             h = modh.transform(H)
             ax = plt.subplot(131)
@@ -289,7 +294,7 @@ class Mesh2Mesh(BaseModel):
         return self.P
 
 
-class ActiveModel(BaseModel):
+class ActiveCells(IdentityMap):
     """
         Active model parameters.
 
@@ -329,18 +334,18 @@ class ActiveModel(BaseModel):
     def transformDeriv(self, m):
         return self.P
 
-class ComboModel(BaseModel):
-    """Combination of various models."""
+class ComboMap(IdentityMap):
+    """Combination of various maps."""
 
-    def __init__(self, mesh, models, **kwargs):
-        BaseModel.__init__(self, mesh, **kwargs)
+    def __init__(self, mesh, maps, **kwargs):
+        IdentityMap.__init__(self, mesh, **kwargs)
 
-        self.models = []
-        for m in models:
-            if not isinstance(m, BaseModel):
-                self.models += [m(mesh, **kwargs)]
+        self.maps = []
+        for m in maps:
+            if not isinstance(m, IdentityMap):
+                self.maps += [m(mesh, **kwargs)]
             else:
-                self.models += [m]
+                self.maps += [m]
 
     @property
     def nP(self):
@@ -348,25 +353,25 @@ class ComboModel(BaseModel):
 
            The number of cells in the
            last dimension of the mesh."""
-        return self.models[-1].nP
+        return self.maps[-1].nP
 
     def transform(self, m):
-        for model in reversed(self.models):
-            m = model.transform(m)
+        for map_i in reversed(self.maps):
+            m = map_i.transform(m)
         return m
 
     def transformDeriv(self, m):
         deriv = 1
         mi = m
-        for model in reversed(self.models):
-            deriv = model.transformDeriv(mi) * deriv
-            mi = model.transform(mi)
+        for map_i in reversed(self.maps):
+            deriv = map_i.transformDeriv(mi) * deriv
+            mi = map_i.transform(mi)
         return deriv
 
 if __name__ == '__main__':
     from SimPEG import *
     mesh = Mesh.TensorMesh([10,8])
-    combo = ComboModel(mesh, [LogModel, Vertical1DModel])
+    combo = ComboMap(mesh, [ExpMap, Vertical1DMap])
     m = combo.example()
     print m.shape
     print combo.test(np.arange(8))

SimPEG/Mesh/BaseMesh.py

@@ -51,7 +51,7 @@ class BaseMesh(object):
     @property
     def nC(self):
         """
-        Total number of cells in the model.
+        Total number of cells in the mesh.
 
         :rtype: int
         :return: nC

SimPEG/Mesh/TensorMesh.py

@@ -220,8 +220,12 @@ class BaseTensorMesh(BaseRectangularMesh):
         if type(loc) == list:
             loc = np.array(loc)
         assert type(loc) == np.ndarray, "must be a numpy array"
+
+        if loc.size == self.dim:
+            loc = np.atleast_2d(loc)
+
         if self.dim > 1:
-            assert loc.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)"
+            assert loc.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim) not (%d,%d)" % loc.shape
         elif len(loc.shape) == 1:
             loc = loc[:,np.newaxis]
         else:

SimPEG/Mesh/__init__.py

@@ -3,3 +3,4 @@ from CylMesh import CylMesh
 from Cyl1DMesh import Cyl1DMesh
 from LogicallyRectMesh import LogicallyRectMesh
 from TreeMesh import TreeMesh
+from BaseMesh import BaseMesh

SimPEG/Parameters.py

@@ -45,7 +45,7 @@ class Parameter(object):
     @property
     def prob(self): return self.parent.prob
     @property
-    def model(self): return self.parent.model
+    def mapping(self): return self.parent.mapping
     @property
     def mesh(self): return self.parent.mesh
 

SimPEG/Problem.py

@@ -1,5 +1,5 @@
 import Utils, Survey, numpy as np, scipy.sparse as sp
-import Model
+import Maps, Mesh
 
 class BaseProblem(object):
     """
@@ -11,19 +11,17 @@ class BaseProblem(object):
     counter = None   #: A SimPEG.Utils.Counter object
 
     surveyPair = Survey.BaseSurvey   #: A SimPEG.Survey Class
-    modelPair = Model.BaseModel  #: A SimPEG.Model Class
+    mapPair    = Maps.IdentityMap    #: A SimPEG.Map Class
 
-    def __init__(self, model, **kwargs):
+    mapping = None    #: A SimPEG.Map instance.
+    mesh    = None    #: A SimPEG.Mesh instance.
+
+    def __init__(self, mesh, mapping=None, **kwargs):
         Utils.setKwargs(self, **kwargs)
-        assert (isinstance(model, self.modelPair) or
-            isinstance(model, Model.ComboModel) and isinstance(model.models[0], self.modelPair)
-            ), "Model object must be an instance of a %s class."%(self.modelPair.__name__)
-        self.model = model
-
-    @property
-    def mesh(self):
-        """SimPEG mesh that is associated with the model provided."""
-        return self.model.mesh
+        assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
+        self.mesh = mesh
+        self.mapping = mapping or Maps.IdentityMap(mesh)
+        self.mapping._assertMatchesPair(self.mapPair)
 
     @property
     def survey(self):
@@ -53,7 +51,8 @@ class BaseProblem(object):
 
     @Utils.timeIt
     def Jvec(self, m, v, u=None):
-        """
+        """Jvec(m, v, u=None)
+
             Effect of J(m) on a vector v.
 
             :param numpy.array m: model
@@ -66,7 +65,8 @@ class BaseProblem(object):
 
     @Utils.timeIt
     def Jtvec(self, m, v, u=None):
-        """
+        """Jtvec(m, v, u=None)
+
             Effect of transpose of J(m) on a vector v.
 
             :param numpy.array m: model
@@ -80,7 +80,8 @@ class BaseProblem(object):
 
     @Utils.timeIt
     def Jvec_approx(self, m, v, u=None):
-        """
+        """Jvec_approx(m, v, u=None)
+
             Approximate effect of J(m) on a vector v
 
             :param numpy.array m: model
@@ -93,7 +94,8 @@ class BaseProblem(object):
 
     @Utils.timeIt
     def Jtvec_approx(self, m, v, u=None):
-        """
+        """Jtvec_approx(m, v, u=None)
+
             Approximate effect of transpose of J(m) on a vector v.
 
             :param numpy.array m: model

SimPEG/Regularization.py

@@ -1,4 +1,4 @@
-import Utils, Model, Parameters, numpy as np, scipy.sparse as sp
+import Utils, Maps, Mesh, Parameters, numpy as np, scipy.sparse as sp
 
 class BaseRegularization(object):
     """
@@ -6,22 +6,27 @@ class BaseRegularization(object):
 
     This is used to regularize the model space::
 
-        reg = Regularization(mesh, model)
+        reg = Regularization(mesh)
 
     """
 
     __metaclass__ = Utils.SimPEGMetaClass
 
-    modelPair = Model.BaseModel  #: Some regularizations only work on specific models
 
-    model = None    #: A SimPEG.Model instance.
 
     counter = None
 
-    def __init__(self, model, **kwargs):
+    mapPair = Maps.IdentityMap    #: A SimPEG.Map Class
+
+    mapping = None    #: A SimPEG.Map instance.
+    mesh    = None    #: A SimPEG.Mesh instance.
+
+    def __init__(self, mesh, mapping=None, **kwargs):
         Utils.setKwargs(self, **kwargs)
-        assert isinstance(model, self.modelPair), "Incorrect model for this regularization"
-        self.model = model
+        self.mesh = mesh
+        assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
+        self.mapping = mapping or Maps.IdentityMap(mesh)
+        self.mapping._assertMatchesPair(self.mapPair)
 
     mref = Parameters.ParameterProperty('mref', default=None, doc='Reference model.')
 
@@ -47,19 +52,17 @@ class BaseRegularization(object):
     def prob(self): return self.parent.prob
     @property
     def survey(self): return self.parent.survey
-    @property
-    def mesh(self): return self.model.mesh
 
 
     @property
     def W(self):
         """Full regularization weighting matrix W."""
-        return sp.identity(self.model.nP)
+        return sp.identity(self.mapping.nP)
 
 
     @Utils.timeIt
     def modelObj(self, m):
-        r = self.W * self.model.transform(m - self.mref)
+        r = self.W * self.mapping.transform(m - self.mref)
         return 0.5*r.dot(r)
 
     @Utils.timeIt
@@ -79,8 +82,8 @@ class BaseRegularization(object):
             R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
 
         """
-        mTd = self.model.transformDeriv(m - self.mref)
-        return mTd.T * ( self.W.T * ( self.W * self.model.transform(m - self.mref) ) )
+        mTd = self.mapping.transformDeriv(m - self.mref)
+        return mTd.T * ( self.W.T * ( self.W * self.mapping.transform(m - self.mref) ) )
 
     @Utils.timeIt
     def modelObj2Deriv(self, m, v=None):
@@ -104,7 +107,7 @@ class BaseRegularization(object):
             R(m) = \mathbf{W^\\top W}
 
         """
-        mTd = self.model.transformDeriv(m - self.mref)
+        mTd = self.mapping.transformDeriv(m - self.mref)
         if v is None:
             return mTd.T * self.W.T * self.W * mTd
 
@@ -204,8 +207,8 @@ class Tikhonov(BaseRegularization):
     alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction")
     alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction")
 
-    def __init__(self, model, **kwargs):
-        BaseRegularization.__init__(self, model, **kwargs)
+    def __init__(self, mesh, mapping=None, **kwargs):
+        BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs)
 
     @property
     def Ws(self):

SimPEG/Tests/test_maps.py

@@ -0,0 +1,38 @@
+import numpy as np
+import unittest
+from SimPEG import *
+from TestUtils import checkDerivative
+from scipy.sparse.linalg import dsolve
+
+
+class MapTests(unittest.TestCase):
+
+    def setUp(self):
+
+        a = np.array([1, 1, 1])
+        b = np.array([1, 2])
+        self.mesh2 = Mesh.TensorMesh([a, b], x0=np.array([3, 5]))
+        self.mesh22 = Mesh.TensorMesh([b, a], x0=np.array([3, 5]))
+
+    def test_transforms(self):
+        for M in dir(Maps):
+            try:
+                maps = getattr(Maps, M)(self.mesh2)
+                assert isinstance(maps, Maps.IdentityMap)
+            except Exception, e:
+                continue
+            self.assertTrue(maps.test())
+
+    def test_Mesh2MeshMap(self):
+        maps = Maps.Mesh2Mesh([self.mesh22, self.mesh2])
+        self.assertTrue(maps.test())
+
+    def test_comboMaps(self):
+        combos = [(Maps.ExpMap, Maps.Vertical1DMap)]
+        for combo in combos:
+            maps = Maps.ComboMap(self.mesh2, combo)
+            self.assertTrue(maps.test())
+
+
+if __name__ == '__main__':
+    unittest.main()

SimPEG/Tests/test_model.py

@@ -1,38 +0,0 @@
-import numpy as np
-import unittest
-from SimPEG import *
-from TestUtils import checkDerivative
-from scipy.sparse.linalg import dsolve
-
-
-class ModelTests(unittest.TestCase):
-
-    def setUp(self):
-
-        a = np.array([1, 1, 1])
-        b = np.array([1, 2])
-        self.mesh2 = Mesh.TensorMesh([a, b], x0=np.array([3, 5]))
-        self.mesh22 = Mesh.TensorMesh([b, a], x0=np.array([3, 5]))
-
-    def test_modelTransforms(self):
-        for M in dir(Model):
-            try:
-                model = getattr(Model, M)(self.mesh2)
-                assert isinstance(model, Model.BaseModel)
-            except Exception, e:
-                continue
-            self.assertTrue(model.test())
-
-    def test_Mesh2MeshModel(self):
-        model = Model.Mesh2Mesh([self.mesh22, self.mesh2])
-        self.assertTrue(model.test())
-
-    def test_comboModels(self):
-        combos = [(Model.LogModel, Model.Vertical1DModel)]
-        for combo in combos:
-            model = Model.ComboModel(self.mesh2, combo)
-            self.assertTrue(model.test())
-
-
-if __name__ == '__main__':
-    unittest.main()

SimPEG/Tests/test_regularization.py

@@ -19,10 +19,10 @@ class RegularizationTests(unittest.TestCase):
                 continue
             # if 'Regularization' not in R: continue
             print 'Check:', R
-            model = r.modelPair(self.mesh2)
-            reg = r(model)
-            m = model.example()
-            reg.mref = model.example()*0
+            mapping = r.mapPair(self.mesh2)
+            reg = r(self.mesh2, mapping=mapping)
+            m = mapping.example()
+            reg.mref = mapping.example()*0
             passed = checkDerivative(lambda m : [reg.modelObj(m), reg.modelObjDeriv(m)], m, plotIt=False)
             self.assertTrue(passed)
 

SimPEG/__init__.py

@@ -3,7 +3,7 @@ import scipy.sparse as sp
 import Utils
 from Solver import Solver
 import Mesh
-import Model
+import Maps
 import Problem
 import Survey
 import Regularization

Tutorials/Linear.py

@@ -24,42 +24,41 @@ class LinearProblem(Problem.BaseProblem):
 
 
 def example(N):
-    M = Mesh.TensorMesh([N])
+    mesh = Mesh.TensorMesh([N])
 
     nk = 20
     jk = np.linspace(1.,20.,nk)
     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 = lambda k: np.exp(p*jk[k]*mesh.vectorCCx)*np.cos(2*np.pi*q*jk[k]*mesh.vectorCCx)
 
-    G = np.empty((nk, M.nC))
+    G = np.empty((nk, mesh.nC))
 
     for i in range(nk):
         G[i,:] = g(i)
 
-    mtrue = np.zeros(M.nC)
-    mtrue[M.vectorCCx > 0.3] = 1.
-    mtrue[M.vectorCCx > 0.45] = -0.5
-    mtrue[M.vectorCCx > 0.6] = 0
+    mtrue = np.zeros(mesh.nC)
+    mtrue[mesh.vectorCCx > 0.3] = 1.
+    mtrue[mesh.vectorCCx > 0.45] = -0.5
+    mtrue[mesh.vectorCCx > 0.6] = 0
 
 
 
-    model = Model.BaseModel(M)
-    prob = LinearProblem(model, G)
+    prob = LinearProblem(mesh, G)
     survey = prob.createSyntheticSurvey(mtrue, std=0.01)
 
-    return prob, survey, model
+    return prob, survey, mesh
 
 
 if __name__ == '__main__':
 
     import matplotlib.pyplot as plt
 
-    prob, survey, model = example(100)
+    prob, survey, mesh = example(100)
     M = prob.mesh
 
-    reg = Regularization.Tikhonov(model)
+    reg = Regularization.Tikhonov(mesh)
     beta = Parameters.BetaSchedule()
     objFunc = ObjFunction.BaseObjFunction(survey, reg, beta=beta)
     opt = Optimization.InexactGaussNewton(maxIter=20)

docs/api_Maps.rst

@@ -1,45 +1,45 @@
-.. _api_Model:
+.. _api_Maps:
 
 
-Model
-*****
+Maps
+****
 
-A SimPEG model operates on a vector and transforms it to another space.
+A SimPEG Map operates on a vector and transforms it to another space.
 We will use an example commonly applied in electromagnetics (EM) of the
-log-conductivity model (:class:`SimPEG.Model.LogModel`).
+log-conductivity model.
 Electrical conductivity varies over many orders of magnitude, so it is a common
 technique when solving the inverse problem to parameterize and optimize in terms
 of log conductivity. This makes sense not only because it ensures all conductivities
 will be positive, but because this is fundamentally the space where conductivity
-lives (i.e. it varies logarithmically). In SimPEG, we use the term Model to
-describe how to get between these two spaces.
+lives (i.e. it varies logarithmically). In SimPEG, we use a (:class:`SimPEG.Maps.ExpMap`) to
+describe how to map back to conductivity.
 
 The API
 =======
 
-.. autoclass:: SimPEG.Model.BaseModel
+.. autoclass:: SimPEG.Maps.IdentityMap
     :members:
     :undoc-members:
 
-.. autoclass:: SimPEG.Model.BaseNonLinearModel
+.. autoclass:: SimPEG.Maps.NonLinearMap
     :members:
     :undoc-members:
 
-.. autoclass:: SimPEG.Model.ComboModel
+.. autoclass:: SimPEG.Maps.ComboMap
     :members:
     :undoc-members:
 
-Common Models
-=============
+Common Maps
+===========
 
-.. autoclass:: SimPEG.Model.LogModel
+.. autoclass:: SimPEG.Maps.ExpMap
     :members:
     :undoc-members:
 
-.. autoclass:: SimPEG.Model.Vertical1DModel
+.. autoclass:: SimPEG.Maps.Vertical1DMap
     :members:
     :undoc-members:
 
-.. autoclass:: SimPEG.Model.Mesh2Mesh
+.. autoclass:: SimPEG.Maps.Mesh2Mesh
     :members:
     :undoc-members:

docs/api_Mesh.rst

@@ -197,7 +197,6 @@ other types of meshes in this SimPEG framework.
 The API
 =======
 
-.. toctree::
-   :maxdepth: 2
-
-   api_MeshCode
+.. automodule:: SimPEG.Mesh.BaseMesh
+    :members:
+    :undoc-members:

docs/api_MeshCode.rst

@@ -9,6 +9,15 @@ Tensor Mesh
     :undoc-members:
 
 
+Tree Mesh
+=========
+
+.. autoclass:: SimPEG.Mesh.TreeMesh.TreeMesh
+    :show-inheritance:
+    :members:
+    :undoc-members:
+
+
 Logically Rectangular Mesh
 ==========================
 
@@ -25,11 +34,3 @@ Cylindrical 1D Mesh
     :show-inheritance:
     :members:
     :undoc-members:
-
-
-Base Mesh
-=========
-
-.. automodule:: SimPEG.Mesh.BaseMesh
-    :members:
-    :undoc-members:

docs/index.rst

@@ -27,7 +27,7 @@ Meshing & Operators
 *******************
 
 .. toctree::
-   :maxdepth: 2
+   :maxdepth: 3
 
    api_Mesh
    api_DiffOps
@@ -39,9 +39,9 @@ Forward Problems
 .. toctree::
    :maxdepth: 2
 
-   api_Model
-   api_Survey
    api_Problem
+   api_Survey
+   api_Maps
 
 Inversion
 *********