diff --git a/SimPEG/Model.py b/SimPEG/Maps.py similarity index 87% rename from SimPEG/Model.py rename to SimPEG/Maps.py index 00f8f302..37b210d9 100644 --- a/SimPEG/Model.py +++ b/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)) diff --git a/SimPEG/Mesh/BaseMesh.py b/SimPEG/Mesh/BaseMesh.py index ebd07c2c..9083e59d 100644 --- a/SimPEG/Mesh/BaseMesh.py +++ b/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 diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 4d1b16a2..a8e29768 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/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: diff --git a/SimPEG/Mesh/__init__.py b/SimPEG/Mesh/__init__.py index a79e93dc..1bfa4c6f 100644 --- a/SimPEG/Mesh/__init__.py +++ b/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 diff --git a/SimPEG/Parameters.py b/SimPEG/Parameters.py index fa15ff67..e138ee17 100644 --- a/SimPEG/Parameters.py +++ b/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 diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index 477ba604..0d502e90 100644 --- a/SimPEG/Problem.py +++ b/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 diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 8f7b9ede..5318b7da 100644 --- a/SimPEG/Regularization.py +++ b/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): diff --git a/SimPEG/Tests/test_maps.py b/SimPEG/Tests/test_maps.py new file mode 100644 index 00000000..e689c597 --- /dev/null +++ b/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() diff --git a/SimPEG/Tests/test_model.py b/SimPEG/Tests/test_model.py deleted file mode 100644 index cd4a0e07..00000000 --- a/SimPEG/Tests/test_model.py +++ /dev/null @@ -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() diff --git a/SimPEG/Tests/test_regularization.py b/SimPEG/Tests/test_regularization.py index 3df9180c..71265456 100644 --- a/SimPEG/Tests/test_regularization.py +++ b/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) diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 874f4ffe..49061e09 100644 --- a/SimPEG/__init__.py +++ b/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 diff --git a/Tutorials/Linear.py b/Tutorials/Linear.py index 1bfa7339..1f4bf94e 100644 --- a/Tutorials/Linear.py +++ b/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) diff --git a/docs/api_Model.rst b/docs/api_Maps.rst similarity index 53% rename from docs/api_Model.rst rename to docs/api_Maps.rst index 05963ff0..7974c272 100644 --- a/docs/api_Model.rst +++ b/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: diff --git a/docs/api_Mesh.rst b/docs/api_Mesh.rst index 065dc58e..45172365 100644 --- a/docs/api_Mesh.rst +++ b/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: diff --git a/docs/api_MeshCode.rst b/docs/api_MeshCode.rst index 65820fc8..9bf3ae4f 100644 --- a/docs/api_MeshCode.rst +++ b/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: diff --git a/docs/index.rst b/docs/index.rst index 98f4efb7..fea8b612 100644 --- a/docs/index.rst +++ b/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 ********* diff --git a/docs/simpeg-framework.png b/docs/simpeg-framework.png index f90f8260..0d020d63 100644 Binary files a/docs/simpeg-framework.png and b/docs/simpeg-framework.png differ