From 7d34f596ff03078356570fe100dd5ef264d45da5 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 17 Jun 2014 09:47:03 -0600 Subject: [PATCH] More tests on more meshes, IPDeriv returns a function. --- SimPEG/Mesh/InnerProducts.py | 12 +- SimPEG/Mesh/TensorMesh.py | 6 +- SimPEG/Tests/test_massMatricesDerivs.py | 222 +++++++++++++++++++----- SimPEG/Utils/matutils.py | 7 - 4 files changed, 192 insertions(+), 55 deletions(-) diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 4ea9a3d2..f133410b 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -114,8 +114,14 @@ class InnerProducts(object): :param TensorType tensorType: type of the tensor: TensorType(mesh, sigma) :param list P: list of projection matrices :param bool doFast: do a faster implementation if available. + :rtype: function + :return: dMdmu(u), the derivative of the inner product matrix (u) + + Given u, dMdmu returns (nF, nC*nA) + + :param np.ndarray u: vector that multiplies dMdmu :rtype: scipy.csr_matrix - :return: dMdm, the derivative of the inner product matrix (nF, nC*nA) + :return: dMdmu, the derivative of the inner product matrix for a certain u """ assert isinstance(tensorType, TensorType), 'tensorType must be an instance of TensorType.' fast = None @@ -131,7 +137,7 @@ class InnerProducts(object): def innerProductDeriv(v): return self._getInnerProductDeriv(tensorType, P, self.nF, v) - return DerivOperator(innerProductDeriv) + return innerProductDeriv def getEdgeInnerProductDeriv(self, tensorType, P=None, doFast=True): @@ -155,7 +161,7 @@ class InnerProducts(object): P = self._getInnerProductProjectionMatrices('E', tensorType=tensorType) def innerProductDeriv(v): return self._getInnerProductDeriv(tensorType, P, self.nE, v) - return DerivOperator(innerProductDeriv) + return innerProductDeriv def _getInnerProductDeriv(self, tensorType, P, n, v): """ diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index faf45cf2..96de4be5 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -319,21 +319,21 @@ class BaseTensorMesh(BaseRectangularMesh): # return self.dim * Av.T * V * ones def scalarInnerProductDeriv(v): return Utils.sdiag(v) * self.dim * Av.T * V * ones - return Utils.DerivOperator(scalarInnerProductDeriv) + return scalarInnerProductDeriv if tensorType == 1: Av = getattr(self, 'ave'+projType+'2CC') V = Utils.sdiag(self.vol) def isotropicInnerProductDeriv(v): return Utils.sdiag(v) * self.dim * Av.T * V - return Utils.DerivOperator(isotropicInnerProductDeriv) + return isotropicInnerProductDeriv if tensorType == 2: # anisotropic Av = getattr(self, 'ave'+projType+'2CCV') V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) def anisotropicInnerProductDeriv(v): return Utils.sdiag(v) * Av.T * V - return Utils.DerivOperator(anisotropicInnerProductDeriv) + return anisotropicInnerProductDeriv diff --git a/SimPEG/Tests/test_massMatricesDerivs.py b/SimPEG/Tests/test_massMatricesDerivs.py index 72e36f3e..dd208c68 100644 --- a/SimPEG/Tests/test_massMatricesDerivs.py +++ b/SimPEG/Tests/test_massMatricesDerivs.py @@ -6,102 +6,240 @@ from TestUtils import checkDerivative class TestInnerProductsDerivs(unittest.TestCase): - def doTestFace(self, h, rep, fast): - mesh = Mesh.TensorMesh(h) + def doTestFace(self, h, rep, fast, meshType): + if meshType == 'LRM': + hRect = Utils.exampleLrmGrid(h,'rotate') + mesh = Mesh.LogicallyRectMesh(hRect) + elif meshType == 'Tree': + mesh = Mesh.TreeMesh(h) + elif meshType == 'Tensor': + mesh = Mesh.TensorMesh(h) v = np.random.rand(mesh.nF) sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep) Md = mesh.getFaceInnerProductDeriv(Utils.TensorType(mesh, sig), doFast=fast) def fun(sig): M = mesh.getFaceInnerProduct(sig) - return M*v, Md*v + return M*v, Md(v) + print meshType, 'Face', h, rep, fast return checkDerivative(fun, sig, num=5, plotIt=False) - def doTestEdge(self, h, rep, fast): - mesh = Mesh.TensorMesh(h) + def doTestEdge(self, h, rep, fast, meshType): + if meshType == 'LRM': + hRect = Utils.exampleLrmGrid(h,'rotate') + mesh = Mesh.LogicallyRectMesh(hRect) + elif meshType == 'Tree': + mesh = Mesh.TreeMesh(h) + elif meshType == 'Tensor': + mesh = Mesh.TensorMesh(h) v = np.random.rand(mesh.nE) sig = np.random.rand(1) if rep is 0 else np.random.rand(mesh.nC*rep) Md = mesh.getEdgeInnerProductDeriv(Utils.TensorType(mesh, sig), doFast=fast) def fun(sig): M = mesh.getEdgeInnerProduct(sig) - return M*v, Md*v + return M*v, Md(v) + print meshType, 'Edge', h, rep, fast return checkDerivative(fun, sig, num=5, plotIt=False) def test_FaceIP_1D_float(self): - self.assertTrue(self.doTestFace([10],0, False)) + self.assertTrue(self.doTestFace([10],0, False, 'Tensor')) def test_FaceIP_2D_float(self): - self.assertTrue(self.doTestFace([10, 4],0, False)) + self.assertTrue(self.doTestFace([10, 4],0, False, 'Tensor')) def test_FaceIP_3D_float(self): - self.assertTrue(self.doTestFace([10, 4, 5],0, False)) + self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'Tensor')) def test_FaceIP_1D_isotropic(self): - self.assertTrue(self.doTestFace([10],1, False)) + self.assertTrue(self.doTestFace([10],1, False, 'Tensor')) def test_FaceIP_2D_isotropic(self): - self.assertTrue(self.doTestFace([10, 4],1, False)) + self.assertTrue(self.doTestFace([10, 4],1, False, 'Tensor')) def test_FaceIP_3D_isotropic(self): - self.assertTrue(self.doTestFace([10, 4, 5],1, False)) + self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'Tensor')) def test_FaceIP_2D_anisotropic(self): - self.assertTrue(self.doTestFace([10, 4],2, False)) + self.assertTrue(self.doTestFace([10, 4],2, False, 'Tensor')) def test_FaceIP_3D_anisotropic(self): - self.assertTrue(self.doTestFace([10, 4, 5],3, False)) + self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'Tensor')) def test_FaceIP_2D_tensor(self): - self.assertTrue(self.doTestFace([10, 4],3, False)) + self.assertTrue(self.doTestFace([10, 4],3, False, 'Tensor')) def test_FaceIP_3D_tensor(self): - self.assertTrue(self.doTestFace([10, 4, 5],6, False)) + self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'Tensor')) def test_FaceIP_1D_float_fast(self): - self.assertTrue(self.doTestFace([10],0, True)) + self.assertTrue(self.doTestFace([10],0, True, 'Tensor')) def test_FaceIP_2D_float_fast(self): - self.assertTrue(self.doTestFace([10, 4],0, True)) + self.assertTrue(self.doTestFace([10, 4],0, True, 'Tensor')) def test_FaceIP_3D_float_fast(self): - self.assertTrue(self.doTestFace([10, 4, 5],0, True)) + self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Tensor')) def test_FaceIP_1D_isotropic_fast(self): - self.assertTrue(self.doTestFace([10],1, True)) + self.assertTrue(self.doTestFace([10],1, True, 'Tensor')) def test_FaceIP_2D_isotropic_fast(self): - self.assertTrue(self.doTestFace([10, 4],1, True)) + self.assertTrue(self.doTestFace([10, 4],1, True, 'Tensor')) def test_FaceIP_3D_isotropic_fast(self): - self.assertTrue(self.doTestFace([10, 4, 5],1, True)) + self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'Tensor')) def test_FaceIP_2D_anisotropic_fast(self): - self.assertTrue(self.doTestFace([10, 4],2, True)) + self.assertTrue(self.doTestFace([10, 4],2, True, 'Tensor')) def test_FaceIP_3D_anisotropic_fast(self): - self.assertTrue(self.doTestFace([10, 4, 5],3, True)) + self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Tensor')) def test_EdgeIP_1D_float(self): - self.assertTrue(self.doTestEdge([10],0, False)) + self.assertTrue(self.doTestEdge([10],0, False, 'Tensor')) def test_EdgeIP_2D_float(self): - self.assertTrue(self.doTestEdge([10, 4],0, False)) + self.assertTrue(self.doTestEdge([10, 4],0, False, 'Tensor')) def test_EdgeIP_3D_float(self): - self.assertTrue(self.doTestEdge([10, 4, 5],0, False)) + self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'Tensor')) def test_EdgeIP_1D_isotropic(self): - self.assertTrue(self.doTestEdge([10],1, False)) + self.assertTrue(self.doTestEdge([10],1, False, 'Tensor')) def test_EdgeIP_2D_isotropic(self): - self.assertTrue(self.doTestEdge([10, 4],1, False)) + self.assertTrue(self.doTestEdge([10, 4],1, False, 'Tensor')) def test_EdgeIP_3D_isotropic(self): - self.assertTrue(self.doTestEdge([10, 4, 5],1, False)) + self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'Tensor')) def test_EdgeIP_2D_anisotropic(self): - self.assertTrue(self.doTestEdge([10, 4],2, False)) + self.assertTrue(self.doTestEdge([10, 4],2, False, 'Tensor')) def test_EdgeIP_3D_anisotropic(self): - self.assertTrue(self.doTestEdge([10, 4, 5],3, False)) + self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'Tensor')) def test_EdgeIP_2D_tensor(self): - self.assertTrue(self.doTestEdge([10, 4],3, False)) + self.assertTrue(self.doTestEdge([10, 4],3, False, 'Tensor')) def test_EdgeIP_3D_tensor(self): - self.assertTrue(self.doTestEdge([10, 4, 5],6, False)) + self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'Tensor')) def test_EdgeIP_1D_float_fast(self): - self.assertTrue(self.doTestEdge([10],0, True)) + self.assertTrue(self.doTestEdge([10],0, True, 'Tensor')) def test_EdgeIP_2D_float_fast(self): - self.assertTrue(self.doTestEdge([10, 4],0, True)) + self.assertTrue(self.doTestEdge([10, 4],0, True, 'Tensor')) def test_EdgeIP_3D_float_fast(self): - self.assertTrue(self.doTestEdge([10, 4, 5],0, True)) + self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'Tensor')) def test_EdgeIP_1D_isotropic_fast(self): - self.assertTrue(self.doTestEdge([10],1, True)) + self.assertTrue(self.doTestEdge([10],1, True, 'Tensor')) def test_EdgeIP_2D_isotropic_fast(self): - self.assertTrue(self.doTestEdge([10, 4],1, True)) + self.assertTrue(self.doTestEdge([10, 4],1, True, 'Tensor')) def test_EdgeIP_3D_isotropic_fast(self): - self.assertTrue(self.doTestEdge([10, 4, 5],1, True)) + self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'Tensor')) def test_EdgeIP_2D_anisotropic_fast(self): - self.assertTrue(self.doTestEdge([10, 4],2, True)) + self.assertTrue(self.doTestEdge([10, 4],2, True, 'Tensor')) def test_EdgeIP_3D_anisotropic_fast(self): - self.assertTrue(self.doTestEdge([10, 4, 5],3, True)) + self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'Tensor')) + + def test_FaceIP_2D_float_LRM(self): + self.assertTrue(self.doTestFace([10, 4],0, False, 'LRM')) + def test_FaceIP_3D_float_LRM(self): + self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'LRM')) + def test_FaceIP_2D_isotropic_LRM(self): + self.assertTrue(self.doTestFace([10, 4],1, False, 'LRM')) + def test_FaceIP_3D_isotropic_LRM(self): + self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'LRM')) + def test_FaceIP_2D_anisotropic_LRM(self): + self.assertTrue(self.doTestFace([10, 4],2, False, 'LRM')) + def test_FaceIP_3D_anisotropic_LRM(self): + self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'LRM')) + def test_FaceIP_2D_tensor_LRM(self): + self.assertTrue(self.doTestFace([10, 4],3, False, 'LRM')) + def test_FaceIP_3D_tensor_LRM(self): + self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'LRM')) + + def test_FaceIP_2D_float_fast_LRM(self): + self.assertTrue(self.doTestFace([10, 4],0, True, 'LRM')) + def test_FaceIP_3D_float_fast_LRM(self): + self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'LRM')) + def test_FaceIP_2D_isotropic_fast_LRM(self): + self.assertTrue(self.doTestFace([10, 4],1, True, 'LRM')) + def test_FaceIP_3D_isotropic_fast_LRM(self): + self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'LRM')) + def test_FaceIP_2D_anisotropic_fast_LRM(self): + self.assertTrue(self.doTestFace([10, 4],2, True, 'LRM')) + def test_FaceIP_3D_anisotropic_fast_LRM(self): + self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'LRM')) + + def test_EdgeIP_2D_float_LRM(self): + self.assertTrue(self.doTestEdge([10, 4],0, False, 'LRM')) + def test_EdgeIP_3D_float_LRM(self): + self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'LRM')) + def test_EdgeIP_2D_isotropic_LRM(self): + self.assertTrue(self.doTestEdge([10, 4],1, False, 'LRM')) + def test_EdgeIP_3D_isotropic_LRM(self): + self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'LRM')) + def test_EdgeIP_2D_anisotropic_LRM(self): + self.assertTrue(self.doTestEdge([10, 4],2, False, 'LRM')) + def test_EdgeIP_3D_anisotropic_LRM(self): + self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'LRM')) + def test_EdgeIP_2D_tensor_LRM(self): + self.assertTrue(self.doTestEdge([10, 4],3, False, 'LRM')) + def test_EdgeIP_3D_tensor_LRM(self): + self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'LRM')) + + def test_EdgeIP_2D_float_fast_LRM(self): + self.assertTrue(self.doTestEdge([10, 4],0, True, 'LRM')) + def test_EdgeIP_3D_float_fast_LRM(self): + self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'LRM')) + def test_EdgeIP_2D_isotropic_fast_LRM(self): + self.assertTrue(self.doTestEdge([10, 4],1, True, 'LRM')) + def test_EdgeIP_3D_isotropic_fast_LRM(self): + self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'LRM')) + def test_EdgeIP_2D_anisotropic_fast_LRM(self): + self.assertTrue(self.doTestEdge([10, 4],2, True, 'LRM')) + def test_EdgeIP_3D_anisotropic_fast_LRM(self): + self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'LRM')) + + + + + def test_FaceIP_2D_float_Tree(self): + self.assertTrue(self.doTestFace([10, 4],0, False, 'Tree')) + def test_FaceIP_3D_float_Tree(self): + self.assertTrue(self.doTestFace([10, 4, 5],0, False, 'Tree')) + def test_FaceIP_2D_isotropic_Tree(self): + self.assertTrue(self.doTestFace([10, 4],1, False, 'Tree')) + def test_FaceIP_3D_isotropic_Tree(self): + self.assertTrue(self.doTestFace([10, 4, 5],1, False, 'Tree')) + def test_FaceIP_2D_anisotropic_Tree(self): + self.assertTrue(self.doTestFace([10, 4],2, False, 'Tree')) + def test_FaceIP_3D_anisotropic_Tree(self): + self.assertTrue(self.doTestFace([10, 4, 5],3, False, 'Tree')) + def test_FaceIP_2D_tensor_Tree(self): + self.assertTrue(self.doTestFace([10, 4],3, False, 'Tree')) + def test_FaceIP_3D_tensor_Tree(self): + self.assertTrue(self.doTestFace([10, 4, 5],6, False, 'Tree')) + + def test_FaceIP_2D_float_fast_Tree(self): + self.assertTrue(self.doTestFace([10, 4],0, True, 'Tree')) + def test_FaceIP_3D_float_fast_Tree(self): + self.assertTrue(self.doTestFace([10, 4, 5],0, True, 'Tree')) + def test_FaceIP_2D_isotropic_fast_Tree(self): + self.assertTrue(self.doTestFace([10, 4],1, True, 'Tree')) + def test_FaceIP_3D_isotropic_fast_Tree(self): + self.assertTrue(self.doTestFace([10, 4, 5],1, True, 'Tree')) + def test_FaceIP_2D_anisotropic_fast_Tree(self): + self.assertTrue(self.doTestFace([10, 4],2, True, 'Tree')) + def test_FaceIP_3D_anisotropic_fast_Tree(self): + self.assertTrue(self.doTestFace([10, 4, 5],3, True, 'Tree')) + + def test_EdgeIP_2D_float_Tree(self): + self.assertTrue(self.doTestEdge([10, 4],0, False, 'Tree')) + def test_EdgeIP_3D_float_Tree(self): + self.assertTrue(self.doTestEdge([10, 4, 5],0, False, 'Tree')) + def test_EdgeIP_2D_isotropic_Tree(self): + self.assertTrue(self.doTestEdge([10, 4],1, False, 'Tree')) + def test_EdgeIP_3D_isotropic_Tree(self): + self.assertTrue(self.doTestEdge([10, 4, 5],1, False, 'Tree')) + def test_EdgeIP_2D_anisotropic_Tree(self): + self.assertTrue(self.doTestEdge([10, 4],2, False, 'Tree')) + def test_EdgeIP_3D_anisotropic_Tree(self): + self.assertTrue(self.doTestEdge([10, 4, 5],3, False, 'Tree')) + def test_EdgeIP_2D_tensor_Tree(self): + self.assertTrue(self.doTestEdge([10, 4],3, False, 'Tree')) + def test_EdgeIP_3D_tensor_Tree(self): + self.assertTrue(self.doTestEdge([10, 4, 5],6, False, 'Tree')) + + def test_EdgeIP_2D_float_fast_Tree(self): + self.assertTrue(self.doTestEdge([10, 4],0, True, 'Tree')) + def test_EdgeIP_3D_float_fast_Tree(self): + self.assertTrue(self.doTestEdge([10, 4, 5],0, True, 'Tree')) + def test_EdgeIP_2D_isotropic_fast_Tree(self): + self.assertTrue(self.doTestEdge([10, 4],1, True, 'Tree')) + def test_EdgeIP_3D_isotropic_fast_Tree(self): + self.assertTrue(self.doTestEdge([10, 4, 5],1, True, 'Tree')) + def test_EdgeIP_2D_anisotropic_fast_Tree(self): + self.assertTrue(self.doTestEdge([10, 4],2, True, 'Tree')) + def test_EdgeIP_3D_anisotropic_fast_Tree(self): + self.assertTrue(self.doTestEdge([10, 4, 5],3, True, 'Tree')) + if __name__ == '__main__': unittest.main() diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py index 1a4c4a0e..691aaabb 100644 --- a/SimPEG/Utils/matutils.py +++ b/SimPEG/Utils/matutils.py @@ -349,10 +349,3 @@ class SimPEGLinearOperator(LinearOperator): @property def T(self): return self.__class__((self.shape[1],self.shape[0]),self.rmatvec,rmatvec=self.matvec,matmat=self.matmat) - - -class DerivOperator(object): - def __init__(self, f): - self.f = f - def __mul__(self, v): - return self.f(v)