From 83cb5ce46a9dd3e908e95935856f7910dc8bd65a Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sat, 28 Nov 2015 12:55:24 -0800 Subject: [PATCH] - Meshlesses Identity Map (takes nP instead of a mesh) - Tikhonov regularization if active cells are used (don't take derivs across interfaces between active cells and not) - testing improvements: test 1D, 2D, 3D on a random tensor mesh , also test that for a constant mref, phi_m(ref) = 0 --- SimPEG/Maps.py | 72 +++++++++++++++++++++++++++++ SimPEG/Regularization.py | 49 +++++++++++++++++--- tests/base/test_regularization.py | 75 ++++++++++++++++++++++++++----- 3 files changed, 179 insertions(+), 17 deletions(-) diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 5b4782ac..7c46fa1d 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -118,6 +118,78 @@ class IdentityMap(object): def __str__(self): return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1]) + +class IdentityMap_Meshless(IdentityMap): + + def __init__(self, nP=None, **kwargs): + IdentityMap.__init__(self, None, **kwargs) + self._nP = nP + + @property + def nP(self): + """ + :rtype: int + :return: number of parameters in the model + """ + if self._nP is None: + return '*' + return self._nP + + @property + def shape(self): + """ + The default shape is (mesh.nC, nP). + + :rtype: (int,int) + :return: shape of the operator as a tuple + """ + if self._nP is None: + return ('*', '*') + return (self.nP, self.nP) + + + def _transform(self, m): + """ + Changes the model into the physical property. + + .. note:: + + This can be called by the __mul__ property against a numpy.ndarray. + + :param numpy.array m: model + :rtype: numpy.array + :return: transformed model + + """ + return m + + def inverse(self, D): + """ + Changes the physical property into the model. + + .. note:: + + The *transformInverse* may not be easy to create in general. + + :param numpy.array D: physical property + :rtype: numpy.array + :return: model + + """ + raise NotImplementedError('The transformInverse is not implemented.') + + def deriv(self, m): + """ + The derivative of the transformation. + + :param numpy.array m: model + :rtype: scipy.csr_matrix + :return: derivative of transformed model + + """ + return sp.identity(self.nP) + + class ComboMap(IdentityMap): """Combination of various maps.""" diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index e3506290..7e98f073 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -20,12 +20,13 @@ class BaseRegularization(object): mesh = None #: A SimPEG.Mesh instance. mref = None #: Reference model. - def __init__(self, mesh, mapping=None, **kwargs): + def __init__(self, mesh, mapping=None, indActive=None, **kwargs): Utils.setKwargs(self, **kwargs) 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) + self.indActive = indActive @property def parent(self): @@ -112,8 +113,6 @@ class BaseRegularization(object): return mD.T * ( self.W.T * ( self.W * ( mD * v) ) ) - - class Tikhonov(BaseRegularization): """**Tikhonov Regularization** @@ -205,14 +204,18 @@ 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, mesh, mapping=None, **kwargs): + def __init__(self, mesh, mapping=None, indActive = None, **kwargs): BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs) + self.indActive = indActive @property def Ws(self): """Regularization matrix Ws""" if getattr(self,'_Ws', None) is None: - self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5) + self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5) + if self.indActive is not None: + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + self._Ws = Pac.T * self._Ws * Pac return self._Ws @property @@ -221,6 +224,13 @@ class Tikhonov(BaseRegularization): if getattr(self, '_Wx', None) is None: Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFx].T*self.mesh.vol self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx + + if self.indActive is not None: + indActive_Fx = (self.mesh.aveFx2CC.T * self.indActive) == 1 + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + Pafx = Utils.speye(self.mesh.nFx)[:,indActive_Fx] + self._Wx = Pafx.T*self._Wx*Pac + return self._Wx @property @@ -229,6 +239,13 @@ class Tikhonov(BaseRegularization): if getattr(self, '_Wy', None) is None: Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFx:np.sum(self.mesh.vnF[:2])].T*self.mesh.vol self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady + + if self.indActive is not None: + indActive_Fy = (self.mesh.aveFy2CC.T * self.indActive) == 1 + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + Pafy = Utils.speye(self.mesh.nFy)[:,indActive_Fy] + self._Wy = Pafy.T*self._Wy*Pac + return self._Wy @property @@ -237,6 +254,13 @@ class Tikhonov(BaseRegularization): if getattr(self, '_Wz', None) is None: Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.vnF[:2]):].T*self.mesh.vol self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz + + if self.indActive is not None: + indActive_Fz = (self.mesh.aveFz2CC.T * self.indActive) == 1 + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + Pafz = Utils.speye(self.mesh.nFz)[:,indActive_Fz] + self._Wz = Pafz.T*self._Wz*Pac + return self._Wz @property @@ -244,6 +268,11 @@ class Tikhonov(BaseRegularization): """Regularization matrix Wxx""" if getattr(self, '_Wxx', None) is None: self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx + + if self.indActive is not None: + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + self._Wxx = Pac.T*self._Wxx*Pac + return self._Wxx @property @@ -251,6 +280,11 @@ class Tikhonov(BaseRegularization): """Regularization matrix Wyy""" if getattr(self, '_Wyy', None) is None: self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady + + if self.indActive is not None: + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + self._Wyy = Pac.T*self._Wyy*Pac + return self._Wyy @property @@ -258,6 +292,11 @@ class Tikhonov(BaseRegularization): """Regularization matrix Wzz""" if getattr(self, '_Wzz', None) is None: self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz + + if self.indActive is not None: + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + self._Wzz = Pac.T*self._Wzz*Pac + return self._Wzz @property diff --git a/tests/base/test_regularization.py b/tests/base/test_regularization.py index af7da692..3243e51d 100644 --- a/tests/base/test_regularization.py +++ b/tests/base/test_regularization.py @@ -4,11 +4,17 @@ from SimPEG import * from scipy.sparse.linalg import dsolve import inspect +TOL = 1e-20 class RegularizationTests(unittest.TestCase): def setUp(self): - self.mesh2 = Mesh.TensorMesh([3, 2]) + hx, hy, hz = np.random.rand(10), np.random.rand(9), np.random.rand(8) + hx, hy, hz = hx/hx.sum(), hy/hy.sum(), hz/hz.sum() + mesh1 = Mesh.TensorMesh([hx]) + mesh2 = Mesh.TensorMesh([hx, hy]) + mesh3 = Mesh.TensorMesh([hx, hy, hz]) + self.meshlist = [mesh1,mesh2, mesh3] def test_regularization(self): for R in dir(Regularization): @@ -16,18 +22,63 @@ class RegularizationTests(unittest.TestCase): if not inspect.isclass(r): continue if not issubclass(r, Regularization.BaseRegularization): continue - # if 'Regularization' not in R: continue - mapping = r.mapPair(self.mesh2) - reg = r(self.mesh2, mapping=mapping) - m = np.random.rand(mapping.nP) - reg.mref = m[:]*np.mean(m) + + for i, mesh in enumerate(self.meshlist): - print 'Check:', R - passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False) - self.assertTrue(passed) - print 'Check 2 Deriv:', R - passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False) - self.assertTrue(passed) + print 'Testing %iD'%mesh.dim + + mapping = r.mapPair(mesh) + reg = r(mesh, mapping=mapping) + m = np.random.rand(mapping.nP) + reg.mref = np.ones_like(m)*np.mean(m) + + print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref) + passed = reg.eval(reg.mref) < TOL + self.assertTrue(passed) + + print 'Check:', R + passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False) + self.assertTrue(passed) + + print 'Check 2 Deriv:', R + passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False) + self.assertTrue(passed) + + def test_regularization_ActiveCells(self): + for R in dir(Regularization): + r = getattr(Regularization, R) + if not inspect.isclass(r): continue + if not issubclass(r, Regularization.BaseRegularization): + continue + + for i, mesh in enumerate(self.meshlist): + + print 'Testing Active Cells %iD'%(mesh.dim) + + if mesh.dim == 1: + indAct = Utils.mkvc(mesh.gridCC <= 0.8) + elif mesh.dim == 2: + indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5) + elif mesh.dim == 3: + indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5) + + mapping = Maps.IdentityMap_Meshless(nP=indAct.nonzero()[0].size) + + reg = r(mesh, mapping=mapping, indActive=indAct) + m = np.random.rand(mesh.nC)[indAct] + reg.mref = np.ones_like(m)*np.mean(m) + + print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref) + passed = reg.eval(reg.mref) < TOL + self.assertTrue(passed) + + print 'Check:', R + passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False) + self.assertTrue(passed) + + print 'Check 2 Deriv:', R + passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False) + self.assertTrue(passed) if __name__ == '__main__':