diff --git a/SimPEG/InvProblem.py b/SimPEG/InvProblem.py index 5f3d1abb..0b68a121 100644 --- a/SimPEG/InvProblem.py +++ b/SimPEG/InvProblem.py @@ -14,8 +14,8 @@ class BaseInvProblem(object): debug = False #: Print debugging information counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things - reg = None #: Regularization dmisfit = None #: DataMisfit + reg = None #: Regularization opt = None #: Optimization program u_current = None #: The most current evaluated field @@ -42,7 +42,7 @@ class BaseInvProblem(object): if self.debug: print 'Calling InvProblem.startup' if self.reg.mref is None: - print 'Regularization has not set mref. SimPEG.InvProblem will set it to m0.' + print 'SimPEG.InvProblem will set Regularization.mref to m0.' self.reg.mref = m0 self.phi_d = np.nan @@ -50,7 +50,7 @@ class BaseInvProblem(object): self.m_current = m0 - print 'Setting bfgsH0 to the inverse of the modelObj2Deriv. Done using direct methods.' + print 'SimPEG.InvProblem is setting bfgsH0 to the inverse of the modelObj2Deriv. \n ***Done using direct methods***' self.opt.bfgsH0 = Solver(self.reg.modelObj2Deriv(self.m_current)) @Utils.timeIt diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 59572bbd..4333a617 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -1,39 +1,6 @@ import Utils, numpy as np, scipy.sparse as sp from Tests import checkDerivative -class Model(np.ndarray): - - def __new__(cls, input_array, mapping=None): - assert isinstance(mapping, IdentityMap), 'mapping must be a SimPEG.Mapping' - obj = np.asarray(input_array).view(cls) - obj._mapping = mapping - if not obj.size == mapping.nP: - raise Exception('Incorrect size for array.') - return obj - - def __array_finalize__(self, obj): - if obj is None: return - self._mapping = getattr(obj, '_mapping', None) - - @property - def mapping(self): - return self._mapping - - @property - def transform(self): - if getattr(self, '_transform', None) is None: - self._transform = self.mapping.transform(self.view(np.ndarray)) - return self._transform - - @property - def transformDeriv(self): - if getattr(self, '_transformDeriv', None) is None: - self._transformDeriv = self.mapping.transformDeriv(self.view(np.ndarray)) - return self._transformDeriv - - def test(self, **kwargs): - return self.mapping.test(self.view(np.ndarray),**kwargs) - class IdentityMap(object): """ @@ -66,7 +33,7 @@ class IdentityMap(object): """ return (self.mesh.nC, self.nP) - def transform(self, m): + def _transform(self, m): """ Changes the model into the physical property. @@ -81,7 +48,7 @@ class IdentityMap(object): """ return m - def transformInverse(self, D): + def inverse(self, D): """ Changes the physical property into the model. @@ -96,7 +63,7 @@ class IdentityMap(object): """ raise NotImplementedError('The transformInverse is not implemented.') - def transformDeriv(self, m): + def deriv(self, m): """ The derivative of the transformation. @@ -105,7 +72,7 @@ class IdentityMap(object): :return: derivative of transformed model """ - return sp.identity(m.size) + return sp.identity(self.nP) def test(self, m=None, **kwargs): """Test the derivative of the mapping. @@ -116,12 +83,12 @@ class IdentityMap(object): :return: passed the test? """ - print 'Testing the %s Class!' % self.__class__.__name__ + print 'Testing %s' % str(self) if m is None: m = np.random.rand(self.nP) if 'plotIt' not in kwargs: kwargs['plotIt'] = False - return checkDerivative(lambda m : [self.transform(m), self.transformDeriv(m)], m, **kwargs) + return checkDerivative(lambda m : [self * m, self.deriv(m)], m, **kwargs) def _assertMatchesPair(self, pair): assert (isinstance(self, pair) or @@ -129,104 +96,90 @@ class IdentityMap(object): ), "Mapping object must be an instance of a %s class."%(pair.__name__) def __mul__(self, val): - if isinstance(val, ComboMap): - return ComboMap(self.mesh, [self] + val.maps) - elif isinstance(val, IdentityMap): + if isinstance(val, IdentityMap): + if not self.shape[1] == val.shape[0]: + raise ValueError('Dimension mismatch in %s and %s.' % (str(self), str(val))) return ComboMap(self.mesh, [self, val]) elif isinstance(val, np.ndarray): - return self.transform(val) + if not self.shape[1] == val.shape[0]: + raise ValueError('Dimension mismatch in %s and np.ndarray%s.' % (str(self), str(val.shape))) + return self._transform(val) + raise Exception('Unrecognized data type to multiply. Try a map or a numpy.ndarray!') -class NonLinearMap(object): - """ - SimPEG NonLinearMap + def __str__(self): + return "%s(%d,%d)" % (self.__class__.__name__, self.shape[0], self.shape[1]) - """ +class ComboMap(IdentityMap): + """Combination of various maps.""" - __metaclass__ = Utils.SimPEGMetaClass + def __init__(self, mesh, maps, **kwargs): + IdentityMap.__init__(self, mesh, **kwargs) - counter = None #: A SimPEG.Utils.Counter object - mesh = None #: A SimPEG Mesh + self.maps = [] + for ii, m in enumerate(maps): + assert isinstance(m, IdentityMap), 'Unrecognized data type, inherit from an IdentityMap or ComboMap!' + if ii > 0 and not self.shape[1] == m.shape[0]: + prev = self.maps[-1] + errArgs = (prev.__name__, prev.shape[0], prev.shape[1], m.__name__, m.shape[0], m.shape[1]) + raise ValueError('Dimension mismatch in map[%s] (%i, %i) and map[%s] (%i, %i).' % errArgs) - def __init__(self, mesh): - self.mesh = mesh + if isinstance(m, ComboMap): + self.maps += m.maps + elif isinstance(m, IdentityMap): + self.maps += [m] - def transform(self, u, m): - """ - :param numpy.array u: fields - :param numpy.array m: model - :rtype: numpy.array - :return: transformed model - - The *transform* changes the model into the physical property. - - """ - return m - - def transformDerivU(self, u, m): - """ - :param numpy.array u: fields - :param numpy.array m: model - :rtype: scipy.csr_matrix - :return: derivative of transformed model - - The *transform* changes the model into the physical property. - The *transformDerivU* provides the derivative of the *transform* with respect to the fields. - """ - raise NotImplementedError('The transformDerivU is not implemented.') - - - def transformDerivM(self, u, m): - """ - :param numpy.array u: fields - :param numpy.array m: model - :rtype: scipy.csr_matrix - :return: derivative of transformed model - - The *transform* changes the model into the physical property. - The *transformDerivU* provides the derivative of the *transform* with respect to the model. - """ - raise NotImplementedError('The transformDerivM is not implemented.') + @property + def shape(self): + return (self.maps[0].shape[0], self.maps[-1].shape[1]) @property def nP(self): - """Number of parameters in the model.""" - return self.mesh.nC + """Number of model properties. - def example(self): - raise NotImplementedError('The example is not implemented.') + The number of cells in the + last dimension of the mesh.""" + return self.maps[-1].nP - def test(self, m=None): - raise NotImplementedError('The test is not implemented.') + def _transform(self, m): + for map_i in reversed(self.maps): + m = map_i * m + return m + + def deriv(self, m): + deriv = 1 + mi = m + for map_i in reversed(self.maps): + deriv = map_i.deriv(mi) * deriv + mi = map_i * mi + return deriv + + def __str__(self): + return 'ComboMap[%s]%s' % (' * '.join([m.__str__() for m in self.maps]), str(self.shape)) class ExpMap(IdentityMap): - """SimPEG ExpMap""" + """ + + Changes the model into the physical property. + + A common example of this is to invert for electrical conductivity + in log space. In this case, your model will be log(sigma) and to + get back to sigma, you can take the exponential: + + .. math:: + + m = \log{\sigma} + + \exp{m} = \exp{\log{\sigma}} = \sigma + """ def __init__(self, mesh, **kwargs): IdentityMap.__init__(self, mesh, **kwargs) - def transform(self, m): - """ - :param numpy.array m: model - :rtype: numpy.array - :return: transformed model - - The *transform* changes the model into the physical property. - - A common example of this is to invert for electrical conductivity - in log space. In this case, your model will be log(sigma) and to - get back to sigma, you can take the exponential: - - .. math:: - - m = \log{\sigma} - - \exp{m} = \exp{\log{\sigma}} = \sigma - """ + def _transform(self, m): return np.exp(Utils.mkvc(m)) - - def transformInverse(self, D): + def inverse(self, D): """ :param numpy.array D: physical property :rtype: numpy.array @@ -242,7 +195,7 @@ class ExpMap(IdentityMap): return np.log(Utils.mkvc(D)) - def transformDeriv(self, m): + def deriv(self, m): """ :param numpy.array m: model :rtype: scipy.csr_matrix @@ -286,7 +239,7 @@ class Vertical1DMap(IdentityMap): last dimension of the mesh.""" return self.mesh.vnC[self.mesh.dim-1] - def transform(self, m): + def _transform(self, m): """ :param numpy.array m: model :rtype: numpy.array @@ -295,7 +248,7 @@ class Vertical1DMap(IdentityMap): repNum = self.mesh.vnC[:self.mesh.dim-1].prod() return Utils.mkvc(m).repeat(repNum) - def transformDeriv(self, m): + def deriv(self, m): """ :param numpy.array m: model :rtype: scipy.csr_matrix @@ -326,13 +279,18 @@ class Mesh2Mesh(IdentityMap): self.P = self.mesh2.getInterpolationMat(self.mesh.gridCC,'CC',zerosOutside=True) + @property + def shape(self): + """Number of parameters in the model.""" + return (self.mesh.nC, self.mesh2.nC) + @property def nP(self): """Number of parameters in the model.""" return self.mesh2.nC - def transform(self, m): + def _transform(self, m): return self.P*m - def transformDeriv(self, m): + def deriv(self, m): return self.P @@ -366,57 +324,20 @@ class ActiveCells(IdentityMap): inds = np.nonzero(self.indActive)[0] self.P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(self.nC, self.nP)) + @property + def shape(self): + return (self.nC, self.nP) + @property def nP(self): """Number of parameters in the model.""" return self.indActive.sum() - def transform(self, m): + def _transform(self, m): return self.P*m + self.valInactive - def transformDeriv(self, m): + def deriv(self, m): return self.P -class ComboMap(IdentityMap): - """Combination of various maps.""" - - def __init__(self, mesh, maps, **kwargs): - IdentityMap.__init__(self, mesh, **kwargs) - - self.maps = [] - for m in maps: - if not isinstance(m, IdentityMap): - self.maps += [m(mesh, **kwargs)] - else: - self.maps += [m] - - @property - def nP(self): - """Number of model properties. - - The number of cells in the - last dimension of the mesh.""" - return self.maps[-1].nP - - def transform(self, m): - for map_i in reversed(self.maps): - m = map_i.transform(m) - return m - - def transformDeriv(self, m): - deriv = 1 - mi = m - for map_i in reversed(self.maps): - deriv = map_i.transformDeriv(mi) * deriv - mi = map_i.transform(mi) - return deriv - - def __mul__(self, val): - if isinstance(val, ComboMap): - return ComboMap(self.mesh, self.maps + val.maps) - elif isinstance(val, IdentityMap): - return ComboMap(self.mesh, self.maps + [val]) - elif isinstance(val, np.ndarray): - return self.transform(val) class ComplexMap(IdentityMap): """ComplexMap @@ -438,11 +359,11 @@ class ComplexMap(IdentityMap): def shape(self): return (self.nP/2,self.nP) - def transform(self, m): + def _transform(self, m): nC = self.mesh.nC return m[:nC] + m[nC:]*1j - def transformDeriv(self, m): + def deriv(self, m): nC = self.nP/2 shp = (nC, nC*2) def fwd(v): @@ -451,5 +372,5 @@ class ComplexMap(IdentityMap): return np.r_[v.real,v.imag] return Utils.SimPEGLinearOperator(shp,fwd,adj) - transformInverse = transformDeriv + inverse = deriv diff --git a/SimPEG/Models.py b/SimPEG/Models.py new file mode 100644 index 00000000..fdefc31d --- /dev/null +++ b/SimPEG/Models.py @@ -0,0 +1,36 @@ +import numpy as np +from Maps import * + + +class Model(np.ndarray): + + def __new__(cls, input_array, mapping=None): + assert isinstance(mapping, IdentityMap), 'mapping must be a SimPEG.Mapping' + obj = np.asarray(input_array).view(cls) + obj._mapping = mapping + if not obj.size == mapping.nP: + raise Exception('Incorrect size for array.') + return obj + + def __array_finalize__(self, obj): + if obj is None: return + self._mapping = getattr(obj, '_mapping', None) + + @property + def mapping(self): + return self._mapping + + @property + def transform(self): + if getattr(self, '_transform', None) is None: + self._transform = self.mapping * self.view(np.ndarray) + return self._transform + + @property + def transformDeriv(self): + if getattr(self, '_transformDeriv', None) is None: + self.deriv = self.mapping.deriv(self.view(np.ndarray)) + return self.deriv + + def test(self, **kwargs): + return self.mapping.test(self.view(np.ndarray),**kwargs) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 90dfe4d5..cd103357 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -61,7 +61,7 @@ class BaseRegularization(object): @Utils.timeIt def modelObj(self, m): - r = self.W * self.mapping.transform(m - self.mref) + r = self.W * ( self.mapping * (m - self.mref) ) return 0.5*r.dot(r) @Utils.timeIt @@ -81,8 +81,9 @@ class BaseRegularization(object): R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} """ - mTd = self.mapping.transformDeriv(m - self.mref) - return mTd.T * ( self.W.T * ( self.W * self.mapping.transform(m - self.mref) ) ) + mD = self.mapping.deriv(m - self.mref) + r = self.W * ( self.mapping * (m - self.mref) ) + return mD.T * ( self.W.T * r ) @Utils.timeIt def modelObj2Deriv(self, m, v=None): @@ -106,11 +107,11 @@ class BaseRegularization(object): R(m) = \mathbf{W^\\top W} """ - mTd = self.mapping.transformDeriv(m - self.mref) + mD = self.mapping.deriv(m - self.mref) if v is None: - return mTd.T * self.W.T * self.W * mTd + return mD.T * self.W.T * self.W * mD - return mTd.T * ( self.W.T * ( self.W * ( mTd * v) ) ) + return mD.T * ( self.W.T * ( self.W * ( mD * v) ) ) diff --git a/SimPEG/Tests/test_maps.py b/SimPEG/Tests/test_maps.py index e518f182..04623fd6 100644 --- a/SimPEG/Tests/test_maps.py +++ b/SimPEG/Tests/test_maps.py @@ -28,23 +28,19 @@ class MapTests(unittest.TestCase): 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()) - def test_mapMultiplication(self): M = Mesh.TensorMesh([2,3]) expMap = Maps.ExpMap(M) vertMap = Maps.Vertical1DMap(M) combo = expMap*vertMap m = np.arange(3.0) - t = combo * m t_true = np.exp(np.r_[0,0,1,1,2,2.]) - self.assertLess(np.linalg.norm(t-t_true,np.inf),TOL) + self.assertLess(np.linalg.norm((combo * m)-t_true,np.inf),TOL) + self.assertLess(np.linalg.norm((expMap * vertMap * m)-t_true,np.inf),TOL) + self.assertLess(np.linalg.norm(expMap * (vertMap * m)-t_true,np.inf),TOL) + self.assertLess(np.linalg.norm((expMap * vertMap) * m-t_true,np.inf),TOL) #Try making a model - mod = Maps.Model(m,mapping=combo) + mod = Models.Model(m, mapping=combo) # print mod.transform # import matplotlib.pyplot as plt # plt.colorbar(M.plotImage(mod.transform)[0]) @@ -53,13 +49,24 @@ class MapTests(unittest.TestCase): self.assertTrue(mod.test(plotIt=False)) - self.assertRaises(Exception,Maps.Model,np.r_[1.0],mapping=combo) + self.assertRaises(Exception,Models.Model,np.r_[1.0],mapping=combo) + + self.assertRaises(ValueError, lambda: combo * (vertMap * expMap)) + self.assertRaises(ValueError, lambda: (combo * vertMap) * expMap) + self.assertRaises(ValueError, lambda: vertMap * expMap) + self.assertRaises(ValueError, lambda: expMap * np.ones(100)) + self.assertRaises(ValueError, lambda: expMap * np.ones((100.0,1))) + self.assertRaises(ValueError, lambda: expMap * np.ones((100.0,5))) + self.assertRaises(ValueError, lambda: combo * np.ones(100)) + self.assertRaises(ValueError, lambda: combo * np.ones((100.0,1))) + self.assertRaises(ValueError, lambda: combo * np.ones((100.0,5))) def test_activeCells(self): M = Mesh.TensorMesh([2,4],'0C') + expMap = Maps.ExpMap(M) actMap = Maps.ActiveCells(M, M.vectorCCy <=0, 10, nC=M.nCy) vertMap = Maps.Vertical1DMap(M) - mod = Maps.Model(np.r_[1,2.],vertMap * actMap) + mod = Models.Model(np.r_[1,2.],vertMap * actMap) # import matplotlib.pyplot as plt # plt.colorbar(M.plotImage(mod.transform)[0]) # plt.show() @@ -67,5 +74,22 @@ class MapTests(unittest.TestCase): self.assertTrue(mod.test()) + def test_tripleMultiply(self): + M = Mesh.TensorMesh([2,4],'0C') + expMap = Maps.ExpMap(M) + vertMap = Maps.Vertical1DMap(M) + actMap = Maps.ActiveCells(M, M.vectorCCy <=0, 10, nC=M.nCy) + m = np.r_[1,2.] + t_true = np.exp(np.r_[1,1,2,2,10,10,10,10.]) + self.assertLess(np.linalg.norm((expMap * vertMap * actMap * m)-t_true,np.inf),TOL) + self.assertLess(np.linalg.norm(((expMap * vertMap * actMap) * m)-t_true,np.inf),TOL) + self.assertLess(np.linalg.norm((expMap * vertMap * (actMap * m))-t_true,np.inf),TOL) + self.assertLess(np.linalg.norm((expMap * (vertMap * actMap) * m)-t_true,np.inf),TOL) + self.assertLess(np.linalg.norm(((expMap * vertMap) * actMap * m)-t_true,np.inf),TOL) + + self.assertRaises(ValueError, lambda: expMap * actMap * vertMap ) + self.assertRaises(ValueError, lambda: actMap * vertMap * expMap ) + + if __name__ == '__main__': unittest.main() diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index 5ba6ea55..10251c76 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -4,6 +4,7 @@ import Utils from Utils.SolverUtils import * import Mesh import Maps +import Models import Problem import Survey import Regularization