From fdfc1c600e73cddf754f8af8e675f802a219405e Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 19 Feb 2014 18:49:30 -0800 Subject: [PATCH 01/18] fixed sub2ind and ind2sub (thin wrappers on numpy) as per Issue #58 --- SimPEG/Tests/test_utils.py | 15 ++++++++++++++- SimPEG/Utils/matutils.py | 39 +++++++++++++------------------------- 2 files changed, 27 insertions(+), 27 deletions(-) diff --git a/SimPEG/Tests/test_utils.py b/SimPEG/Tests/test_utils.py index fea231f2..12d79fb4 100644 --- a/SimPEG/Tests/test_utils.py +++ b/SimPEG/Tests/test_utils.py @@ -1,6 +1,6 @@ import numpy as np import unittest -from SimPEG.Utils import mkvc, ndgrid, indexCube, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal +from SimPEG.Utils import mkvc, ndgrid, indexCube, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal,sub2ind,ind2sub from SimPEG.Tests import checkDerivative @@ -64,6 +64,19 @@ class TestSequenceFunctions(unittest.TestCase): self.assertTrue(np.all(XYZ[:, 1] == X2_test)) self.assertTrue(np.all(XYZ[:, 2] == X3_test)) + def test_sub2ind(self): + x = np.ones((5,2)) + self.assertTrue(np.all(sub2ind(x.shape, [0,0]) == [0])) + self.assertTrue(np.all(sub2ind(x.shape, [4,0]) == [4])) + self.assertTrue(np.all(sub2ind(x.shape, [0,1]) == [5])) + self.assertTrue(np.all(sub2ind(x.shape, [4,1]) == [9])) + self.assertTrue(np.all(sub2ind(x.shape, [[0,0],[4,0],[0,1],[4,1]]) == [0,4,5,9])) + + def test_ind2sub(self): + x = np.ones((5,2)) + self.assertTrue(np.all(ind2sub(x.shape, [0,4,5,9])[0] == [0,4,0,4])) + self.assertTrue(np.all(ind2sub(x.shape, [0,4,5,9])[1] == [0,0,1,1])) + def test_indexCube_2D(self): nN = np.array([3, 3]) self.assertTrue(np.all(indexCube('A', nN) == np.array([0, 1, 3, 4]))) diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py index 24286cdd..b2e1daf2 100644 --- a/SimPEG/Utils/matutils.py +++ b/SimPEG/Utils/matutils.py @@ -97,35 +97,22 @@ def ndgrid(*args, **kwargs): else: return XYZ[2], XYZ[1], XYZ[0] - -def ind2sub(shape, ind): - """From the given shape, returns the subscrips of the given index""" - revshp = [] - revshp.extend(shape) - mult = [1] - for i in range(0, len(revshp)-1): - mult.extend([mult[i]*revshp[i]]) - mult = np.array(mult).reshape(len(mult)) - - sub = [] - - for i in range(0, len(shape)): - sub.extend([np.math.floor(ind / mult[i])]) - ind = ind - (np.math.floor(ind/mult[i]) * mult[i]) - return sub - +def ind2sub(shape, inds): + """From the given shape, returns the subscripts of the given index""" + if type(inds) is not np.ndarray: + inds = np.array(inds) + assert len(inds.shape) == 1, 'Indexing must be done as a 1D row vector, e.g. [3,6,6,...]' + return np.unravel_index(inds, shape, order='F') def sub2ind(shape, subs): """From the given shape, returns the index of the given subscript""" - revshp = list(shape) - mult = [1] - for i in range(0, len(revshp)-1): - mult.extend([mult[i]*revshp[i]]) - mult = np.array(mult).reshape(len(mult), 1) - - idx = np.dot((subs), (mult)) - return idx - + if type(subs) is not np.ndarray: + subs = np.array(subs) + if subs.size == len(shape): + subs = subs[np.newaxis,:] + assert subs.shape[1] == len(shape), 'Indexing must be done as a column vectors. e.g. [[3,6],[6,2],...]' + inds = np.ravel_multi_index(subs.T, shape, order='F') + return mkvc(inds) def getSubArray(A, ind): """subArray""" From 8b295ed7eb3bc3519584aa74d913dc282f69375d Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Thu, 20 Feb 2014 00:09:46 -0800 Subject: [PATCH 02/18] interp updates --- SimPEG/Mesh/TensorMesh.py | 43 ++++++++++++++++++++++++++++--------- SimPEG/Utils/interputils.py | 8 ++++--- 2 files changed, 38 insertions(+), 13 deletions(-) diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 12253bb5..21662d8f 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -367,7 +367,7 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): return [t for t in ten if t is not None] - def isInside(self, pts): + def isInside(self, pts, locType='N'): """ Determines if a set of points are inside a mesh. @@ -376,15 +376,23 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): :return inside, numpy array of booleans """ - pts = np.atleast_2d(pts) - inside = (pts[:,0] >= self.vectorNx.min()) & (pts[:,0] <= self.vectorNx.max()) + tensors = self.getTensor(locType) + if type(pts) == list: + pts = np.array(pts) + assert type(pts) == np.ndarray, "must be a numpy array" if self.dim > 1: - inside = inside & ((pts[:,1] >= self.vectorNy.min()) & (pts[:,1] <= self.vectorNy.max())) - if self.dim > 2: - inside = inside & ((pts[:,2] >= self.vectorNz.min()) & (pts[:,2] <= self.vectorNz.max())) + assert pts.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" + elif len(pts.shape) == 1: + pts = pts[:,np.newaxis] + else: + assert pts.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" + + inside = np.ones(pts.shape[0],dtype=bool) + for i, tensor in enumerate(tensors): + inside = inside & (pts[:,i] >= tensor.min()) & (pts[:,i] <= tensor.max()) return inside - def getInterpolationMat(self, loc, locType): + def getInterpolationMat(self, loc, locType, zerosOutside=False): """ Produces interpolation matrix :param numpy.ndarray loc: Location of points to interpolate to @@ -404,8 +412,21 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): 'CC' -> scalar field defined on cell centers """ - loc = np.atleast_2d(loc) - assert np.all(self.isInside(loc)), "Points outside of mesh" + if type(loc) == list: + loc = np.array(loc) + assert type(loc) == np.ndarray, "must be a numpy array" + if self.dim > 1: + assert loc.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" + elif len(loc.shape) == 1: + loc = loc[:,np.newaxis] + else: + assert loc.shape[1] == self.dim, "must be a column vector of shape (nPts, mesh.dim)" + + if zerosOutside is False: + assert np.all(self.isInside(loc)), "Points outside of mesh" + else: + indZeros = np.logical_not(self.isInside(loc)) + loc[indZeros, :] = np.array([v.mean() for v in self.getTensor('CC')]) ind = 0 if 'x' in locType else 1 if 'y' in locType else 2 if 'z' in locType else -1 if locType in ['Fx','Fy','Fz','Ex','Ey','Ez'] and self.dim >= ind: @@ -417,7 +438,9 @@ class TensorMesh(BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): Q = Utils.interpmat(loc, *self.getTensor(locType)) else: raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim)) - return Q + if zerosOutside: + Q[indZeros, :] = 0 + return Q.tocsr() if __name__ == '__main__': print('Welcome to tensor mesh!') diff --git a/SimPEG/Utils/interputils.py b/SimPEG/Utils/interputils.py index 38308492..907b9456 100644 --- a/SimPEG/Utils/interputils.py +++ b/SimPEG/Utils/interputils.py @@ -20,6 +20,8 @@ def _interp_point_1D(x, xr_i): elif xr_i - x[im] < 0: # Point on the right ind_x1 = im-1 ind_x2 = im + ind_x1 = max(min(ind_x1, x.size-1), 0) + ind_x2 = max(min(ind_x2, x.size-1), 0) dx1 = xr_i - x[ind_x1] dx2 = x[ind_x2] - xr_i return ind_x1, ind_x2, dx1, dx2 @@ -77,7 +79,7 @@ def _interpmat1D(locs, x): inds = [ind_x1, ind_x2] vals = [(1-dx1/Dx),(1-dx2/Dx)] Q[i, inds] = vals - return Q.tocsr() + return Q @@ -114,7 +116,7 @@ def _interpmat2D(locs, x, y): Q[i, mkvc(inds)] = vals - return Q.tocsr() + return Q @@ -162,7 +164,7 @@ def _interpmat3D(locs, x, y, z): Q[i, mkvc(inds)] = vals - return Q.tocsr() + return Q if __name__ == '__main__': From 9f2c976be1af8b51de9deb36cd7b1cebdd57a086 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Thu, 20 Feb 2014 00:38:13 -0800 Subject: [PATCH 03/18] interp bug fixes --- SimPEG/Utils/interputils.py | 79 +++++++++++++++++++++---------------- 1 file changed, 46 insertions(+), 33 deletions(-) diff --git a/SimPEG/Utils/interputils.py b/SimPEG/Utils/interputils.py index 907b9456..04aed35d 100644 --- a/SimPEG/Utils/interputils.py +++ b/SimPEG/Utils/interputils.py @@ -24,7 +24,13 @@ def _interp_point_1D(x, xr_i): ind_x2 = max(min(ind_x2, x.size-1), 0) dx1 = xr_i - x[ind_x1] dx2 = x[ind_x2] - xr_i - return ind_x1, ind_x2, dx1, dx2 + + Dx = x[ind_x2] - x[ind_x1] + if ind_x1 == ind_x2: + dx1 = 0.5 + dx2 = 0.5 + Dx = 1 + return ind_x1, ind_x2, dx1, dx2, Dx def interpmat(locs, x, y=None, z=None): @@ -72,13 +78,16 @@ def _interpmat1D(locs, x): Q = sp.lil_matrix((npts, nx)) for i in range(npts): - ind_x1, ind_x2, dx1, dx2 = _interp_point_1D(x, locs[i]) - dv = (x[ind_x2] - x[ind_x1]) - Dx = x[ind_x2] - x[ind_x1] + ind_x1, ind_x2, dx1, dx2, Dx = _interp_point_1D(x, locs[i]) + # dv = (x[ind_x2] - x[ind_x1]) # Get the row in the matrix inds = [ind_x1, ind_x2] vals = [(1-dx1/Dx),(1-dx2/Dx)] - Q[i, inds] = vals + + for I, V in zip(inds, vals): + Q[i, I] += V + # Q[i, mkvc(inds)] = vals + return Q @@ -93,13 +102,10 @@ def _interpmat2D(locs, x, y): for i in range(npts): - ind_x1, ind_x2, dx1, dx2 = _interp_point_1D(x, locs[i, 0]) - ind_y1, ind_y2, dy1, dy2 = _interp_point_1D(y, locs[i, 1]) + ind_x1, ind_x2, dx1, dx2, Dx = _interp_point_1D(x, locs[i, 0]) + ind_y1, ind_y2, dy1, dy2, Dy = _interp_point_1D(y, locs[i, 1]) - dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) - - Dx = x[ind_x2] - x[ind_x1] - Dy = y[ind_y2] - y[ind_y1] + # dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) # Get the row in the matrix @@ -114,7 +120,11 @@ def _interpmat2D(locs, x, y): (1-dx2/Dx)*(1-dy1/Dy), (1-dx2/Dx)*(1-dy2/Dy)] - Q[i, mkvc(inds)] = vals + + for I, V in zip(mkvc(inds), vals): + Q[i, I] += V + # Q[i, mkvc(inds)] = vals + return Q @@ -131,15 +141,11 @@ def _interpmat3D(locs, x, y, z): for i in range(npts): - ind_x1, ind_x2, dx1, dx2 = _interp_point_1D(x, locs[i, 0]) - ind_y1, ind_y2, dy1, dy2 = _interp_point_1D(y, locs[i, 1]) - ind_z1, ind_z2, dz1, dz2 = _interp_point_1D(z, locs[i, 2]) + ind_x1, ind_x2, dx1, dx2, Dx = _interp_point_1D(x, locs[i, 0]) + ind_y1, ind_y2, dy1, dy2, Dy = _interp_point_1D(y, locs[i, 1]) + ind_z1, ind_z2, dz1, dz2, Dz = _interp_point_1D(z, locs[i, 2]) - dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) *(z[ind_z2] - z[ind_z1]) - - Dx = x[ind_x2] - x[ind_x1] - Dy = y[ind_y2] - y[ind_y1] - Dz = z[ind_z2] - z[ind_z1] + # dv = (x[ind_x2] - x[ind_x1]) * (y[ind_y2] - y[ind_y1]) *(z[ind_z2] - z[ind_z1]) # Get the row in the matrix @@ -162,22 +168,29 @@ def _interpmat3D(locs, x, y, z): (1-dx2/Dx)*(1-dy1/Dy)*(1-dz2/Dz), (1-dx2/Dx)*(1-dy2/Dy)*(1-dz2/Dz)] - Q[i, mkvc(inds)] = vals + for I, V in zip(mkvc(inds), vals): + Q[i, I] += V + # Q[i, mkvc(inds)] = vals return Q if __name__ == '__main__': - import SimPEG - import numpy as np + from SimPEG import * import matplotlib.pyplot as plt - locs = np.random.rand(50)*0.8+0.1 - x = np.linspace(0,1,7) - dense = np.linspace(0,1,200) - fun = lambda x: np.cos(2*np.pi*x) - Q = SimPEG.Utils.interpmat(locs, x) - plt.plot(x, fun(x), 'bs-') - plt.plot(dense, fun(dense), 'y:') - plt.plot(locs, Q*fun(x), 'mo') - plt.plot(locs, fun(locs), 'rx') - plt.show() + + M = Mesh.TensorMesh([4]) + M.vectorCCx + Q = M.getInterpolationMat(np.array([[0],[0.126],[0.127]]),'CC',zerosOutside=True) + print Q.todense() + + # locs = np.random.rand(50)*0.8+0.1 + # x = np.linspace(0,1,7) + # dense = np.linspace(0,1,200) + # fun = lambda x: np.cos(2*np.pi*x) + # Q = Utils.interpmat(locs, x) + # plt.plot(x, fun(x), 'bs-') + # plt.plot(dense, fun(dense), 'y:') + # plt.plot(locs, Q*fun(x), 'mo') + # plt.plot(locs, fun(locs), 'rx') + # plt.show() From 361fb719f0ce5415b0cedaee6f799e3a6ed95a52 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Thu, 20 Feb 2014 10:06:10 -0800 Subject: [PATCH 04/18] Fix and test interpolation. Issue #43 --- SimPEG/Tests/test_interpolation.py | 16 ++++++++++++++++ SimPEG/Utils/interputils.py | 26 ++++++++++---------------- 2 files changed, 26 insertions(+), 16 deletions(-) diff --git a/SimPEG/Tests/test_interpolation.py b/SimPEG/Tests/test_interpolation.py index 018f9091..e2634cb7 100644 --- a/SimPEG/Tests/test_interpolation.py +++ b/SimPEG/Tests/test_interpolation.py @@ -2,6 +2,9 @@ import numpy as np import unittest from TestUtils import OrderTest from SimPEG.Utils import mkvc +from SimPEG import Mesh +import unittest + MESHTYPES = ['uniformTensorMesh', 'randomTensorMesh'] TOLERANCES = [0.9, 0.5] @@ -50,6 +53,19 @@ class TestInterpolation1D(OrderTest): self.name = 'Interpolation 1D: N' self.orderTest() +class TestOutliersInterp1D(unittest.TestCase): + + def setUp(self): + pass + + def test_outliers(self): + M = Mesh.TensorMesh([4]) + Q = M.getInterpolationMat(np.array([[0],[0.126],[0.127]]),'CC',zerosOutside=True) + x = np.arange(4)+1 + self.assertTrue(np.all(Q*x == [1,1.004,1.008])) + Q = M.getInterpolationMat(np.array([[-1],[0.126],[0.127]]),'CC',zerosOutside=True) + self.assertTrue(np.all(Q*x == [0,1.004,1.008])) + class TestInterpolation2d(OrderTest): name = "Interpolation 2D" LOCS = np.random.rand(50,2)*0.6+0.2 diff --git a/SimPEG/Utils/interputils.py b/SimPEG/Utils/interputils.py index 04aed35d..3f9b7e71 100644 --- a/SimPEG/Utils/interputils.py +++ b/SimPEG/Utils/interputils.py @@ -178,19 +178,13 @@ def _interpmat3D(locs, x, y, z): if __name__ == '__main__': from SimPEG import * import matplotlib.pyplot as plt - - M = Mesh.TensorMesh([4]) - M.vectorCCx - Q = M.getInterpolationMat(np.array([[0],[0.126],[0.127]]),'CC',zerosOutside=True) - print Q.todense() - - # locs = np.random.rand(50)*0.8+0.1 - # x = np.linspace(0,1,7) - # dense = np.linspace(0,1,200) - # fun = lambda x: np.cos(2*np.pi*x) - # Q = Utils.interpmat(locs, x) - # plt.plot(x, fun(x), 'bs-') - # plt.plot(dense, fun(dense), 'y:') - # plt.plot(locs, Q*fun(x), 'mo') - # plt.plot(locs, fun(locs), 'rx') - # plt.show() + locs = np.random.rand(50)*0.8+0.1 + x = np.linspace(0,1,7) + dense = np.linspace(0,1,200) + fun = lambda x: np.cos(2*np.pi*x) + Q = Utils.interpmat(locs, x) + plt.plot(x, fun(x), 'bs-') + plt.plot(dense, fun(dense), 'y:') + plt.plot(locs, Q*fun(x), 'mo') + plt.plot(locs, fun(locs), 'rx') + plt.show() From 089ac427c4b9a43b14f629022cb8ee5c413c2c5e Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Thu, 20 Feb 2014 10:22:12 -0800 Subject: [PATCH 05/18] Model that allows interpolation from one mesh to another. --- SimPEG/Model.py | 50 ++++++++++++++++++++++++++++++++++++++ SimPEG/Tests/test_model.py | 7 +++++- 2 files changed, 56 insertions(+), 1 deletion(-) diff --git a/SimPEG/Model.py b/SimPEG/Model.py index 5ea87720..1fcd2f0b 100644 --- a/SimPEG/Model.py +++ b/SimPEG/Model.py @@ -174,6 +174,56 @@ class Vertical1DModel(BaseModel): ), shape=(repNum, 1)) return sp.kron(sp.identity(self.nP), repVec) +class Mesh2Mesh(BaseModel): + """ + Takes a model on one mesh are translates it to another mesh. + + .. plot:: + + from SimPEG import * + M = Mesh.TensorMesh([100,100]) + h1 = Utils.meshTensors(((7,6,1.5),(10,6),(7,6,1.5))) + h1 = h1/h1.sum() + 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]) + H = modH.transform(v) + h = modh.transform(H) + ax = plt.subplot(131) + M.plotImage(v, ax=ax) + ax.set_title('Fine Mesh (Original)') + ax = plt.subplot(132) + M2.plotImage(H,clim=[0,1],ax=ax) + ax.set_title('Course Mesh') + ax = plt.subplot(133) + M.plotImage(h,clim=[0,1],ax=ax) + ax.set_title('Fine Mesh (Interpolated)') + + """ + + def __init__(self, meshes, **kwargs): + Utils.setKwargs(self, **kwargs) + + assert type(meshes) is list, "meshes must be a list of two meshes" + assert len(meshes) == 2, "meshes must be a list of two meshes" + assert meshes[0].dim == meshes[1].dim, """The two meshes must be the same dimension""" + + self.mesh = meshes[0] + self.mesh2 = meshes[1] + + self.P = self.mesh2.getInterpolationMat(self.mesh.gridCC,'CC',zerosOutside=True) + + @property + def nP(self): + """Number of parameters in the model.""" + return self.mesh2.nC + def transform(self, m): + return self.P*m + def transformDeriv(self, m): + return self.P + class ComboModel(BaseModel): """Combination of various models.""" diff --git a/SimPEG/Tests/test_model.py b/SimPEG/Tests/test_model.py index e6a8844b..cd4a0e07 100644 --- a/SimPEG/Tests/test_model.py +++ b/SimPEG/Tests/test_model.py @@ -11,7 +11,8 @@ class ModelTests(unittest.TestCase): a = np.array([1, 1, 1]) b = np.array([1, 2]) - self.mesh2 = Mesh.TensorMesh([a, b], np.array([3, 5])) + 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): @@ -22,6 +23,10 @@ class ModelTests(unittest.TestCase): 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: From a649129acf85bb1a614ab8d06d4bf25956962dbb Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Thu, 20 Feb 2014 16:04:11 -0800 Subject: [PATCH 06/18] first attempt at activeCell model --- SimPEG/Model.py | 41 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 41 insertions(+) diff --git a/SimPEG/Model.py b/SimPEG/Model.py index 1fcd2f0b..5cc26dd1 100644 --- a/SimPEG/Model.py +++ b/SimPEG/Model.py @@ -224,6 +224,47 @@ class Mesh2Mesh(BaseModel): def transformDeriv(self, m): return self.P + +class ActiveModel(BaseModel): + """ + Active model parameters. + + """ + + indActive = None #: Active Cells + valInactive = None #: Values of inactive Cells + nC = None #: Number of cells in the full model + + def __init__(self, mesh, indActive, valInactive, nC=None): + self.mesh = mesh + + self.nC = nC or mesh.nC + + if indActive.dtype is not bool: + z = np.zeros(mesh.nC,dtype=bool) + z[indActive] = True + indActive = z + self.indActive = indActive + self.indInactive = np.logical_not(indActive) + if type(valInactive) in [float, int, long]: + valInactive = np.ones(mesh.nC)*float(valInactive) + + valInactive[self.indActive] = 0 + self.valInactive = valInactive + + 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 nP(self): + """Number of parameters in the model.""" + return self.indActive.sum() + + def transform(self, m): + return self.P*m + self.valInactive + def transformDeriv(self, m): + return self.P + class ComboModel(BaseModel): """Combination of various models.""" From 6490922dc99a473b33d6465a779323060d9c984b Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Thu, 20 Feb 2014 16:13:16 -0800 Subject: [PATCH 07/18] updates to combo model --- SimPEG/Model.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/SimPEG/Model.py b/SimPEG/Model.py index 5cc26dd1..f6d4d6d4 100644 --- a/SimPEG/Model.py +++ b/SimPEG/Model.py @@ -241,13 +241,13 @@ class ActiveModel(BaseModel): self.nC = nC or mesh.nC if indActive.dtype is not bool: - z = np.zeros(mesh.nC,dtype=bool) + z = np.zeros(self.nC,dtype=bool) z[indActive] = True indActive = z self.indActive = indActive self.indInactive = np.logical_not(indActive) if type(valInactive) in [float, int, long]: - valInactive = np.ones(mesh.nC)*float(valInactive) + valInactive = np.ones(self.nC)*float(valInactive) valInactive[self.indActive] = 0 self.valInactive = valInactive @@ -270,7 +270,13 @@ class ComboModel(BaseModel): def __init__(self, mesh, models, **kwargs): BaseModel.__init__(self, mesh, **kwargs) - self.models = [m(mesh, **kwargs) for m in models] + + self.models = [] + for m in models: + if not isinstance(m, BaseModel): + self.models += [m(mesh, **kwargs)] + else: + self.models += [m] @property def nP(self): From 544d65bea10158bded18db706180e386179df3bc Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Thu, 20 Feb 2014 16:43:12 -0800 Subject: [PATCH 08/18] regularization updates. --- SimPEG/Regularization.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index ed3d03d1..8a2e9831 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -14,14 +14,16 @@ class BaseRegularization(object): modelPair = Model.BaseModel #: Some regularizations only work on specific models + mesh = None #: A SimPEG.Mesh instance. model = None #: A SimPEG.Model instance. counter = None - def __init__(self, model, **kwargs): + def __init__(self, mesh, model, **kwargs): Utils.setKwargs(self, **kwargs) assert isinstance(model, self.modelPair), "Incorrect model for this regularization" self.model = model + self.mesh = mesh mref = Parameters.ParameterProperty('mref', default=None, doc='Reference model.') @@ -59,7 +61,7 @@ class BaseRegularization(object): @Utils.timeIt def modelObj(self, m): - r = self.W * (m - self.mref) + r = self.W * self.model.transform(m - self.mref) return 0.5*r.dot(r) @Utils.timeIt @@ -79,7 +81,7 @@ class BaseRegularization(object): R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} """ - return self.W.T * ( self.W * (m - self.mref) ) + return self.W.T * ( self.W * self.model.transform(m - self.mref) ) @Utils.timeIt def modelObj2Deriv(self): @@ -193,8 +195,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, model, **kwargs): + BaseRegularization.__init__(self, mesh, model, **kwargs) @property def Ws(self): From 2c233c55249b5ad2534c34ad1fe4855ef98fabff Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Thu, 20 Feb 2014 16:53:05 -0800 Subject: [PATCH 09/18] reg updates --- SimPEG/Regularization.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 8a2e9831..ce2440f0 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -19,11 +19,11 @@ class BaseRegularization(object): counter = None - def __init__(self, mesh, model, **kwargs): + def __init__(self, model, **kwargs): Utils.setKwargs(self, **kwargs) assert isinstance(model, self.modelPair), "Incorrect model for this regularization" self.model = model - self.mesh = mesh + self.mesh = model.mesh mref = Parameters.ParameterProperty('mref', default=None, doc='Reference model.') @@ -195,8 +195,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, mesh, model, **kwargs): - BaseRegularization.__init__(self, mesh, model, **kwargs) + def __init__(self, model, **kwargs): + BaseRegularization.__init__(self, model, **kwargs) @property def Ws(self): From 7e465ef3968259e7aea59255027114e97813fc3d Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Thu, 20 Feb 2014 19:57:44 -0800 Subject: [PATCH 10/18] fix regularization --- SimPEG/Regularization.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index ce2440f0..4fbf5d4b 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -14,7 +14,6 @@ class BaseRegularization(object): modelPair = Model.BaseModel #: Some regularizations only work on specific models - mesh = None #: A SimPEG.Mesh instance. model = None #: A SimPEG.Model instance. counter = None @@ -23,7 +22,6 @@ class BaseRegularization(object): Utils.setKwargs(self, **kwargs) assert isinstance(model, self.modelPair), "Incorrect model for this regularization" self.model = model - self.mesh = model.mesh mref = Parameters.ParameterProperty('mref', default=None, doc='Reference model.') From a0ddf267442f6ea865131c18b2139f2f84594377 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 21 Feb 2014 10:24:29 -0800 Subject: [PATCH 11/18] updates to Travis. --- .travis.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.travis.yml b/.travis.yml index 1094912b..de86afc0 100644 --- a/.travis.yml +++ b/.travis.yml @@ -5,8 +5,8 @@ virtualenv: system_site_packages: true before_install: - sudo apt-get install -qq gcc gfortran libblas-dev liblapack-dev python-numpy python-scipy python-matplotlib python-pip - - sudo pip install scipy --upgrade - - sudo pip install numpy --upgrade + - sudo pip install scipy --upgrade > tempSP.txt + - sudo pip install numpy --upgrade > tempNP.txt - cd SimPEG - python setup.py - cd ../ From 12ee598a09209e39f798857638c7dec7eedc4c4b Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 21 Feb 2014 10:48:26 -0800 Subject: [PATCH 12/18] update travis...annoying output will have to stay :( --- .travis.yml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/.travis.yml b/.travis.yml index de86afc0..1094912b 100644 --- a/.travis.yml +++ b/.travis.yml @@ -5,8 +5,8 @@ virtualenv: system_site_packages: true before_install: - sudo apt-get install -qq gcc gfortran libblas-dev liblapack-dev python-numpy python-scipy python-matplotlib python-pip - - sudo pip install scipy --upgrade > tempSP.txt - - sudo pip install numpy --upgrade > tempNP.txt + - sudo pip install scipy --upgrade + - sudo pip install numpy --upgrade - cd SimPEG - python setup.py - cd ../ From 872a8a178b4719939be7eb0c0f5fb10c071d7581 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 21 Feb 2014 17:06:43 -0800 Subject: [PATCH 13/18] fixed float issue in innerproducts --- SimPEG/Mesh/InnerProducts.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 7a1a7869..7c00967b 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -317,6 +317,9 @@ def _makeTensor(M, sigma): if sigma is None: # default is ones sigma = np.ones((M.nC, 1)) + if type(sigma) is float: + sigma = np.ones((M.nC, 1))*sigma + if M.dim == 1: if sigma.size == M.nC: # Isotropic! sigma = mkvc(sigma) # ensure it is a vector. From 891c4ce972af97e9a0cea780a08909eaa4ee5500 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 21 Feb 2014 19:55:44 -0800 Subject: [PATCH 14/18] added invPropertyTensor and makeTensor to matUtils combined sp utils into mat utils --- SimPEG/Mesh/InnerProducts.py | 46 +------- SimPEG/Tests/test_utils.py | 50 ++++++++- SimPEG/Utils/__init__.py | 5 +- SimPEG/Utils/lomutils.py | 75 ------------- SimPEG/Utils/matutils.py | 197 ++++++++++++++++++++++++++++++++++- SimPEG/Utils/sputils.py | 41 -------- 6 files changed, 247 insertions(+), 167 deletions(-) delete mode 100644 SimPEG/Utils/sputils.py diff --git a/SimPEG/Mesh/InnerProducts.py b/SimPEG/Mesh/InnerProducts.py index 7c00967b..28befd6e 100644 --- a/SimPEG/Mesh/InnerProducts.py +++ b/SimPEG/Mesh/InnerProducts.py @@ -1,5 +1,5 @@ from scipy import sparse as sp -from SimPEG.Utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal +from SimPEG.Utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal, makePropertyTensor import numpy as np @@ -174,7 +174,7 @@ class InnerProducts(object): P011 = V3*Pxxx('fXm', 'fYp', 'fZp') P111 = V3*Pxxx('fXp', 'fYp', 'fZp') - Mu = _makeTensor(M, mu) + Mu = makePropertyTensor(M, mu) A = P000.T*Mu*P000 + P100.T*Mu*P100 P = [P000, P100] @@ -283,7 +283,7 @@ class InnerProducts(object): P011 = V*eP('eX3', 'eY2', 'eZ2') P111 = V*eP('eX3', 'eY3', 'eZ3') - Sigma = _makeTensor(M, sigma) + Sigma = makePropertyTensor(M, sigma) A = P000.T*Sigma*P000 + P100.T*Sigma*P100 + P010.T*Sigma*P010 + P110.T*Sigma*P110 P = [P000, P100, P010, P110] if M.dim == 3: @@ -313,46 +313,6 @@ class InnerProducts(object): # | |/ # node(i+1,j,k) ------ edge2(i+1,j,k) ----- node(i+1,j+1,k) -def _makeTensor(M, sigma): - if sigma is None: # default is ones - sigma = np.ones((M.nC, 1)) - - if type(sigma) is float: - sigma = np.ones((M.nC, 1))*sigma - - if M.dim == 1: - if sigma.size == M.nC: # Isotropic! - sigma = mkvc(sigma) # ensure it is a vector. - Sigma = sdiag(sigma) - else: - raise Exception('Unexpected shape of sigma') - elif M.dim == 2: - if sigma.size == M.nC: # Isotropic! - sigma = mkvc(sigma) # ensure it is a vector. - Sigma = sdiag(np.r_[sigma, sigma]) - elif sigma.shape[1] == 2: # Diagonal tensor - Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1]]) - elif sigma.shape[1] == 3: # Fully anisotropic - row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 2]))) - row2 = sp.hstack((sdiag(sigma[:, 2]), sdiag(sigma[:, 1]))) - Sigma = sp.vstack((row1, row2)) - else: - raise Exception('Unexpected shape of sigma') - elif M.dim == 3: - if sigma.size == M.nC: # Isotropic! - sigma = mkvc(sigma) # ensure it is a vector. - Sigma = sdiag(np.r_[sigma, sigma, sigma]) - elif sigma.shape[1] == 3: # Diagonal tensor - Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1], sigma[:, 2]]) - elif sigma.shape[1] == 6: # Fully anisotropic - row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 3]), sdiag(sigma[:, 4]))) - row2 = sp.hstack((sdiag(sigma[:, 3]), sdiag(sigma[:, 1]), sdiag(sigma[:, 5]))) - row3 = sp.hstack((sdiag(sigma[:, 4]), sdiag(sigma[:, 5]), sdiag(sigma[:, 2]))) - Sigma = sp.vstack((row1, row2, row3)) - else: - raise Exception('Unexpected shape of sigma') - return Sigma - def _getFacePx(M): assert M._meshType == 'TENSOR', 'Only supported for a tensor mesh' diff --git a/SimPEG/Tests/test_utils.py b/SimPEG/Tests/test_utils.py index 12d79fb4..4a68b139 100644 --- a/SimPEG/Tests/test_utils.py +++ b/SimPEG/Tests/test_utils.py @@ -1,6 +1,6 @@ -import numpy as np import unittest -from SimPEG.Utils import mkvc, ndgrid, indexCube, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal,sub2ind,ind2sub +from SimPEG.Utils import * +from SimPEG import Mesh, np, sp from SimPEG.Tests import checkDerivative @@ -96,8 +96,6 @@ class TestSequenceFunctions(unittest.TestCase): self.assertTrue(np.all(indexCube('H', nN) == np.array([10, 11, 13, 14, 19, 20, 22, 23]))) def test_invXXXBlockDiagonal(self): - import scipy.sparse as sp - a = [np.random.rand(5, 1) for i in range(4)] B = inv2X2BlockDiagonal(*a) @@ -120,6 +118,50 @@ class TestSequenceFunctions(unittest.TestCase): self.assertTrue(np.linalg.norm(Z3.todense().ravel(), 2) < 1e-12) + def test_invPropertyTensor2D(self): + M = Mesh.TensorMesh([6, 6]) + a1 = np.random.rand(M.nC) + a2 = np.random.rand(M.nC) + a3 = np.random.rand(M.nC) + prop1 = a1 + prop2 = np.c_[a1, a2] + prop3 = np.c_[a1, a2, a3] + + for prop in [4, prop1, prop2, prop3]: + b = invPropertyTensor(M, prop) + A = makePropertyTensor(M, prop) + B1 = makePropertyTensor(M, b) + B2 = invPropertyTensor(M, prop, returnMatrix=True) + + Z = B1*A - sp.identity(M.nC*2) + self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < 1e-12) + Z = B2*A - sp.identity(M.nC*2) + self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < 1e-12) + + + def test_invPropertyTensor3D(self): + M = Mesh.TensorMesh([6, 6, 6]) + a1 = np.random.rand(M.nC) + a2 = np.random.rand(M.nC) + a3 = np.random.rand(M.nC) + a4 = np.random.rand(M.nC) + a5 = np.random.rand(M.nC) + a6 = np.random.rand(M.nC) + prop1 = a1 + prop2 = np.c_[a1, a2, a3] + prop3 = np.c_[a1, a2, a3, a4, a5, a6] + + for prop in [4, prop1, prop2, prop3]: + b = invPropertyTensor(M, prop) + A = makePropertyTensor(M, prop) + B1 = makePropertyTensor(M, b) + B2 = invPropertyTensor(M, prop, returnMatrix=True) + + Z = B1*A - sp.identity(M.nC*3) + self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < 1e-12) + Z = B2*A - sp.identity(M.nC*3) + self.assertTrue(np.linalg.norm(Z.todense().ravel(), 2) < 1e-12) + if __name__ == '__main__': unittest.main() diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index b2d0d56a..f45a1039 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -1,7 +1,6 @@ -from matutils import getSubArray, mkvc, ndgrid, ind2sub, sub2ind -from sputils import spzeros, kron3, speye, sdiag, sdInv, ddx, av, avExtrap +from matutils import * from meshutils import exampleLomGird, meshTensors -from lomutils import volTetra, faceInfo, inv2X2BlockDiagonal, inv3X3BlockDiagonal, indexCube +from lomutils import volTetra, faceInfo, indexCube from interputils import interpmat from ipythonutils import easyAnimate as animate import ModelBuilder diff --git a/SimPEG/Utils/lomutils.py b/SimPEG/Utils/lomutils.py index 653cfdde..210a134d 100644 --- a/SimPEG/Utils/lomutils.py +++ b/SimPEG/Utils/lomutils.py @@ -188,78 +188,3 @@ def faceInfo(xyz, A, B, C, D, average=True, normalizeNormals=True): return N, area - -def inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33): - """ B = inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33) - - inverts a stack of 3x3 matrices - - Input: - A - a11, a12, a13, a21, a22, a23, a31, a32, a33 - - Output: - B - inverse - """ - - a11 = mkvc(a11) - a12 = mkvc(a12) - a13 = mkvc(a13) - a21 = mkvc(a21) - a22 = mkvc(a22) - a23 = mkvc(a23) - a31 = mkvc(a31) - a32 = mkvc(a32) - a33 = mkvc(a33) - - detA = a31*a12*a23 - a31*a13*a22 - a21*a12*a33 + a21*a13*a32 + a11*a22*a33 - a11*a23*a32 - - b11 = +(a22*a33 - a23*a32)/detA - b12 = -(a12*a33 - a13*a32)/detA - b13 = +(a12*a23 - a13*a22)/detA - - b21 = +(a31*a23 - a21*a33)/detA - b22 = -(a31*a13 - a11*a33)/detA - b23 = +(a21*a13 - a11*a23)/detA - - b31 = -(a31*a22 - a21*a32)/detA - b32 = +(a31*a12 - a11*a32)/detA - b33 = -(a21*a12 - a11*a22)/detA - - B = sp.vstack((sp.hstack((sdiag(b11), sdiag(b12), sdiag(b13))), - sp.hstack((sdiag(b21), sdiag(b22), sdiag(b23))), - sp.hstack((sdiag(b31), sdiag(b32), sdiag(b33))))) - - return B - - -def inv2X2BlockDiagonal(a11, a12, a21, a22): - """ B = inv2X2BlockDiagonal(a11, a12, a21, a22) - - Inverts a stack of 2x2 matrices by using the inversion formula - - inv(A) = (1/det(A)) * cof(A)^T - - Input: - A - a11, a12, a13, a21, a22, a23, a31, a32, a33 - - Output: - B - inverse - """ - - a11 = mkvc(a11) - a12 = mkvc(a12) - a21 = mkvc(a21) - a22 = mkvc(a22) - - # compute inverse of the determinant. - detAinv = 1./(a11*a22 - a21*a12) - - b11 = +detAinv*a22 - b12 = -detAinv*a12 - b21 = -detAinv*a21 - b22 = +detAinv*a11 - - B = sp.vstack((sp.hstack((sdiag(b11), sdiag(b12))), - sp.hstack((sdiag(b21), sdiag(b22))))) - - return B diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py index b2e1daf2..9e350957 100644 --- a/SimPEG/Utils/matutils.py +++ b/SimPEG/Utils/matutils.py @@ -1,5 +1,5 @@ import numpy as np - +import scipy.sparse as sp def mkvc(x, numDims=1): """Creates a vector with the number of dimension specified @@ -30,6 +30,42 @@ def mkvc(x, numDims=1): elif numDims == 3: return x.flatten(order='F')[:, np.newaxis, np.newaxis] +def sdiag(h): + """Sparse diagonal matrix""" + return sp.spdiags(mkvc(h), 0, h.size, h.size, format="csr") + +def sdInv(M): + "Inverse of a sparse diagonal matrix" + return sdiag(1/M.diagonal()) + +def speye(n): + """Sparse identity""" + return sp.identity(n, format="csr") + + +def kron3(A, B, C): + """Three kron prods""" + return sp.kron(sp.kron(A, B), C, format="csr") + + +def spzeros(n1, n2): + """spzeros""" + return sp.coo_matrix((n1, n2)).tocsr() + + +def ddx(n): + """Define 1D derivatives, inner, this means we go from n+1 to n""" + return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr") + + +def av(n): + """Define 1D averaging operator from nodes to cell-centers.""" + return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr") + +def avExtrap(n): + """Define 1D averaging operator from cell-centers to nodes.""" + Av = sp.spdiags((0.5*np.ones((n, 1))*[1, 1]).T, [-1, 0], n+1, n, format="csr") + sp.csr_matrix(([0.5,0.5],([0,n],[0,n-1])),shape=(n+1,n)) + return Av def ndgrid(*args, **kwargs): """ @@ -97,6 +133,7 @@ def ndgrid(*args, **kwargs): else: return XYZ[2], XYZ[1], XYZ[0] + def ind2sub(shape, inds): """From the given shape, returns the subscripts of the given index""" if type(inds) is not np.ndarray: @@ -104,6 +141,7 @@ def ind2sub(shape, inds): assert len(inds.shape) == 1, 'Indexing must be done as a 1D row vector, e.g. [3,6,6,...]' return np.unravel_index(inds, shape, order='F') + def sub2ind(shape, subs): """From the given shape, returns the index of the given subscript""" if type(subs) is not np.ndarray: @@ -114,6 +152,7 @@ def sub2ind(shape, subs): inds = np.ravel_multi_index(subs.T, shape, order='F') return mkvc(inds) + def getSubArray(A, ind): """subArray""" assert type(ind) == list, "ind must be a list of vectors" @@ -125,3 +164,159 @@ def getSubArray(A, ind): return A[ind[0], :, :][:, ind[1], :][:, :, ind[2]] else: raise Exception("getSubArray does not support dimension asked.") + + +def inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33, returnMatrix=True): + """ B = inv3X3BlockDiagonal(a11, a12, a13, a21, a22, a23, a31, a32, a33) + + inverts a stack of 3x3 matrices + + Input: + A - a11, a12, a13, a21, a22, a23, a31, a32, a33 + + Output: + B - inverse + """ + + a11 = mkvc(a11) + a12 = mkvc(a12) + a13 = mkvc(a13) + a21 = mkvc(a21) + a22 = mkvc(a22) + a23 = mkvc(a23) + a31 = mkvc(a31) + a32 = mkvc(a32) + a33 = mkvc(a33) + + detA = a31*a12*a23 - a31*a13*a22 - a21*a12*a33 + a21*a13*a32 + a11*a22*a33 - a11*a23*a32 + + b11 = +(a22*a33 - a23*a32)/detA + b12 = -(a12*a33 - a13*a32)/detA + b13 = +(a12*a23 - a13*a22)/detA + + b21 = +(a31*a23 - a21*a33)/detA + b22 = -(a31*a13 - a11*a33)/detA + b23 = +(a21*a13 - a11*a23)/detA + + b31 = -(a31*a22 - a21*a32)/detA + b32 = +(a31*a12 - a11*a32)/detA + b33 = -(a21*a12 - a11*a22)/detA + + if not returnMatrix: + return b11, b12, b13, b21, b22, b23, b31, b32, b33 + + return sp.vstack((sp.hstack((sdiag(b11), sdiag(b12), sdiag(b13))), + sp.hstack((sdiag(b21), sdiag(b22), sdiag(b23))), + sp.hstack((sdiag(b31), sdiag(b32), sdiag(b33))))) + + + +def inv2X2BlockDiagonal(a11, a12, a21, a22, returnMatrix=True): + """ B = inv2X2BlockDiagonal(a11, a12, a21, a22) + + Inverts a stack of 2x2 matrices by using the inversion formula + + inv(A) = (1/det(A)) * cof(A)^T + + Input: + A - a11, a12, a21, a22 + + Output: + B - inverse + """ + + a11 = mkvc(a11) + a12 = mkvc(a12) + a21 = mkvc(a21) + a22 = mkvc(a22) + + # compute inverse of the determinant. + detAinv = 1./(a11*a22 - a21*a12) + + b11 = +detAinv*a22 + b12 = -detAinv*a12 + b21 = -detAinv*a21 + b22 = +detAinv*a11 + + if not returnMatrix: + return b11, b12, b21, b22 + + return sp.vstack((sp.hstack((sdiag(b11), sdiag(b12))), + sp.hstack((sdiag(b21), sdiag(b22))))) + +def makePropertyTensor(M, sigma): + if sigma is None: # default is ones + sigma = np.ones(M.nC) + + if type(sigma) in [float, int, long]: + sigma = sigma * np.ones(M.nC) + + if M.dim == 1: + if sigma.size == M.nC: # Isotropic! + sigma = mkvc(sigma) # ensure it is a vector. + Sigma = sdiag(sigma) + else: + raise Exception('Unexpected shape of sigma') + elif M.dim == 2: + if sigma.size == M.nC: # Isotropic! + sigma = mkvc(sigma) # ensure it is a vector. + Sigma = sdiag(np.r_[sigma, sigma]) + elif sigma.shape[1] == 2: # Diagonal tensor + Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1]]) + elif sigma.shape[1] == 3: # Fully anisotropic + row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 2]))) + row2 = sp.hstack((sdiag(sigma[:, 2]), sdiag(sigma[:, 1]))) + Sigma = sp.vstack((row1, row2)) + else: + raise Exception('Unexpected shape of sigma') + elif M.dim == 3: + if sigma.size == M.nC: # Isotropic! + sigma = mkvc(sigma) # ensure it is a vector. + Sigma = sdiag(np.r_[sigma, sigma, sigma]) + elif sigma.shape[1] == 3: # Diagonal tensor + Sigma = sdiag(np.r_[sigma[:, 0], sigma[:, 1], sigma[:, 2]]) + elif sigma.shape[1] == 6: # Fully anisotropic + row1 = sp.hstack((sdiag(sigma[:, 0]), sdiag(sigma[:, 3]), sdiag(sigma[:, 4]))) + row2 = sp.hstack((sdiag(sigma[:, 3]), sdiag(sigma[:, 1]), sdiag(sigma[:, 5]))) + row3 = sp.hstack((sdiag(sigma[:, 4]), sdiag(sigma[:, 5]), sdiag(sigma[:, 2]))) + Sigma = sp.vstack((row1, row2, row3)) + else: + raise Exception('Unexpected shape of sigma') + return Sigma + + +def invPropertyTensor(M, tensor, returnMatrix=False): + + T = None + + if type(tensor) in [float, int, long]: + T = 1./tensor + + elif tensor.size == M.nC: # Isotropic! + T = 1./mkvc(tensor) # ensure it is a vector. + + elif M.dim == 2: + if tensor.shape[1] == 2: # Diagonal tensor + T = 1./tensor + elif tensor.shape[1] == 3: # Fully anisotropic + B = inv2X2BlockDiagonal(tensor[:,0], tensor[:,2], + tensor[:,2], tensor[:,1], + returnMatrix=False) + b11, b12, b21, b22 = B + T = np.c_[b11, b22, b12] + elif M.dim == 3: + if tensor.shape[1] == 3: # Diagonal tensor + T = 1./tensor + elif tensor.shape[1] == 6: # Fully anisotropic + B = inv3X3BlockDiagonal(tensor[:,0], tensor[:,3], tensor[:,4], + tensor[:,3], tensor[:,1], tensor[:,5], + tensor[:,4], tensor[:,5], tensor[:,2], + returnMatrix=False) + b11, b12, b13, b21, b22, b23, b31, b32, b33 = B + T = np.c_[b11, b22, b33, b12, b13, b23] + + if T is None: + raise Exception('Unexpected shape of tensor') + if returnMatrix: + return makePropertyTensor(M, T) + return T diff --git a/SimPEG/Utils/sputils.py b/SimPEG/Utils/sputils.py deleted file mode 100644 index f671be43..00000000 --- a/SimPEG/Utils/sputils.py +++ /dev/null @@ -1,41 +0,0 @@ -from scipy import sparse as sp -from matutils import mkvc -import numpy as np - - -def sdiag(h): - """Sparse diagonal matrix""" - return sp.spdiags(mkvc(h), 0, h.size, h.size, format="csr") - -def sdInv(M): - "Inverse of a sparse diagonal matrix" - return sdiag(1/M.diagonal()) - -def speye(n): - """Sparse identity""" - return sp.identity(n, format="csr") - - -def kron3(A, B, C): - """Three kron prods""" - return sp.kron(sp.kron(A, B), C, format="csr") - - -def spzeros(n1, n2): - """spzeros""" - return sp.coo_matrix((n1, n2)).tocsr() - - -def ddx(n): - """Define 1D derivatives, inner, this means we go from n+1 to n""" - return sp.spdiags((np.ones((n+1, 1))*[-1, 1]).T, [0, 1], n, n+1, format="csr") - - -def av(n): - """Define 1D averaging operator from nodes to cell-centers.""" - return sp.spdiags((0.5*np.ones((n+1, 1))*[1, 1]).T, [0, 1], n, n+1, format="csr") - -def avExtrap(n): - """Define 1D averaging operator from cell-centers to nodes.""" - Av = sp.spdiags((0.5*np.ones((n, 1))*[1, 1]).T, [-1, 0], n+1, n, format="csr") + sp.csr_matrix(([0.5,0.5],([0,n],[0,n-1])),shape=(n+1,n)) - return Av From ff5cd54297e35e04f007da64573a49877ec542d2 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 21 Feb 2014 21:53:44 -0800 Subject: [PATCH 15/18] clean up after sputils delete --- SimPEG/Solver.py | 3 +-- SimPEG/Utils/interputils.py | 3 +-- SimPEG/Utils/lomutils.py | 3 +-- SimPEG/Utils/meshutils.py | 3 +-- 4 files changed, 4 insertions(+), 8 deletions(-) diff --git a/SimPEG/Solver.py b/SimPEG/Solver.py index 60ca5558..6c47be27 100644 --- a/SimPEG/Solver.py +++ b/SimPEG/Solver.py @@ -1,8 +1,7 @@ import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as linalg -from Utils.matutils import mkvc -from Utils.sputils import sdiag +from Utils.matutils import mkvc, sdiag import warnings DEFAULTS = {'direct':'scipy', 'iter':'scipy', 'triangular':'fortran', 'diagonal':'python'} diff --git a/SimPEG/Utils/interputils.py b/SimPEG/Utils/interputils.py index 3f9b7e71..f1dc28b9 100644 --- a/SimPEG/Utils/interputils.py +++ b/SimPEG/Utils/interputils.py @@ -1,7 +1,6 @@ import numpy as np import scipy.sparse as sp -from sputils import spzeros -from matutils import mkvc, sub2ind +from matutils import mkvc, sub2ind, spzeros def _interp_point_1D(x, xr_i): """ diff --git a/SimPEG/Utils/lomutils.py b/SimPEG/Utils/lomutils.py index 210a134d..b641ef7f 100644 --- a/SimPEG/Utils/lomutils.py +++ b/SimPEG/Utils/lomutils.py @@ -1,7 +1,6 @@ import numpy as np from scipy import sparse as sp -from matutils import mkvc, ndgrid, sub2ind -from sputils import sdiag +from matutils import mkvc, ndgrid, sub2ind, sdiag def volTetra(xyz, A, B, C, D): diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index d7fe0836..f8add82b 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -1,7 +1,6 @@ import numpy as np from scipy import sparse as sp -from matutils import mkvc, ndgrid, sub2ind -from sputils import sdiag +from matutils import mkvc, ndgrid, sub2ind, sdiag def exampleLomGird(nC, exType): assert type(nC) == list, "nC must be a list containing the number of nodes" From aad54e64ddd70953f21911bbcb3bff9cdff9bbe5 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sat, 22 Feb 2014 13:03:05 -0800 Subject: [PATCH 16/18] Updates to plotSlice to work on irregular grid --- SimPEG/Mesh/TensorView.py | 55 ++++++++++++++++++++++++++++++++++++--- 1 file changed, 51 insertions(+), 4 deletions(-) diff --git a/SimPEG/Mesh/TensorView.py b/SimPEG/Mesh/TensorView.py index ef0957eb..e3967c4c 100644 --- a/SimPEG/Mesh/TensorView.py +++ b/SimPEG/Mesh/TensorView.py @@ -241,8 +241,27 @@ class TensorView(object): def plotSlice(self, v, vType='CC', normal='Z', ind=None, grid=False, view='real', ax=None, clim=None, showIt=False, + pcolorOpts={}, + streamOpts={'color':'k'}, gridOpts={'color':'k'} ): + + """ + Plots a slice of a 3D mesh. + + .. plot:: + + from SimPEG import * + mT = Utils.meshTensors(((2,5),(4,2),(2,5)),((2,2),(6,2),(2,2)),((2,2),(6,2),(2,2))) + M = Mesh.TensorMesh(mT) + q = np.zeros(M.vnC) + q[[4,4],[4,4],[2,6]]=[-1,1] + q = Utils.mkvc(q) + A = M.faceDiv*M.cellGrad + b = Solver(A).solve(q) + M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8}) + + """ viewOpts = ['real','imag','abs','vec'] normalOpts = ['X', 'Y', 'Z'] vTypeOpts = ['CC','F','E'] @@ -302,13 +321,31 @@ class TensorView(object): v = doSlice(v) if clim is None: clim = [v.min(),v.max()] - out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, v.T, vmin=clim[0], vmax=clim[1]),) + out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, v.T, vmin=clim[0], vmax=clim[1], **pcolorOpts),) elif view in ['vec']: U, V = doSlice(v) if clim is None: - clim = [v.min(),v.max()] - out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, 0.5*(U+V).T, vmin=clim[0], vmax=clim[1]),) - out += (plt.streamplot(tM.vectorCCx, tM.vectorCCy, U.T, V.T),) + uv = np.r_[mkvc(U), mkvc(V)] + uv = np.sqrt(uv**2) + clim = [uv.min(),uv.max()] + + # Matplotlib seems to not support irregular + # spaced vectors at the moment. So we will + # Interpolate down to a regular mesh at the + # smallest mesh size in this 2D slice. + nxi = int(tM.hx.sum()/tM.hx.min()) + nyi = int(tM.hy.sum()/tM.hy.min()) + tMi = self.__class__([np.ones(nxi)*tM.hx.sum()/nxi, + np.ones(nyi)*tM.hy.sum()/nyi]) + P = tM.getInterpolationMat(tMi.gridCC,'CC',zerosOutside=True) + Ui = P*mkvc(U) + Vi = P*mkvc(V) + Ui = tMi.r(Ui, 'CC', 'CC', 'M') + Vi = tMi.r(Vi, 'CC', 'CC', 'M') + # End Interpolation + + out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, np.sqrt(U**2+V**2).T, vmin=clim[0], vmax=clim[1], **pcolorOpts),) + out += (plt.streamplot(tMi.vectorCCx, tMi.vectorCCy, Ui.T, Vi.T, **streamOpts),) if grid: xXGrid = np.c_[tM.vectorNx,tM.vectorNx,np.nan*np.ones(tM.nNx)].flatten() @@ -512,3 +549,13 @@ class TensorView(object): return animate(fig, animateFrame, frames=len(frames)) +if __name__ == '__main__': + from SimPEG import * + mT = Utils.meshTensors(((2,5),(4,2),(2,5)),((2,2),(6,2),(2,2)),((2,2),(6,2),(2,2))) + M = Mesh.TensorMesh(mT) + q = np.zeros(M.vnC) + q[[4,4],[4,4],[2,6]]=[-1,1] + q = Utils.mkvc(q) + A = M.faceDiv*M.cellGrad + b = Solver(A).solve(q) + M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8}) From d41e6dd13df263b5f87de261a9d6b43027e9922d Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sat, 22 Feb 2014 13:26:03 -0800 Subject: [PATCH 17/18] updates to plotSlice --- SimPEG/Mesh/TensorView.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/SimPEG/Mesh/TensorView.py b/SimPEG/Mesh/TensorView.py index e3967c4c..d345f0c3 100644 --- a/SimPEG/Mesh/TensorView.py +++ b/SimPEG/Mesh/TensorView.py @@ -264,7 +264,7 @@ class TensorView(object): """ viewOpts = ['real','imag','abs','vec'] normalOpts = ['X', 'Y', 'Z'] - vTypeOpts = ['CC','F','E'] + vTypeOpts = ['CC', 'CCv','F','E'] # Some user error checking assert vType in vTypeOpts, "vType must be in ['%s']" % "','".join(vTypeOpts) @@ -296,11 +296,15 @@ class TensorView(object): def doSlice(v): if vType == 'CC': return getIndSlice(self.r(v,'CC','CC','M')) - # Now just deal with 'F' and 'E' - aveOp = 'ave' + vType + ('2CCV' if view == 'vec' else '2CC') - v = getattr(self,aveOp)*v # average to cell centers (might be a vector) - if view == 'vec': + elif vType == 'CCv': v = self.r(v.reshape((self.nC,3),order='F'),'CC','CC','M') + assert view == 'vec', 'Other types for CCv not yet supported' + else: + # Now just deal with 'F' and 'E' + aveOp = 'ave' + vType + ('2CCV' if view == 'vec' else '2CC') + v = getattr(self,aveOp)*v # average to cell centers (might be a vector) + v = self.r(v.reshape((self.nC,3),order='F'),'CC','CC','M') + if view == 'vec': outSlice = [] if 'X' not in normal: outSlice.append(getIndSlice(v[0])) if 'Y' not in normal: outSlice.append(getIndSlice(v[1])) @@ -357,6 +361,8 @@ class TensorView(object): ax.set_xlabel('y' if normal == 'X' else 'x') ax.set_ylabel('y' if normal == 'Z' else 'z') ax.set_title('Slice %d' % ind) + ax.set_xlim(*tM.vectorNx[[0,-1]]) + ax.set_ylim(*tM.vectorNy[[0,-1]]) if showIt: plt.show() return out From 33431c4e5f5e165dd773fd12a8c9e4714baba776 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sat, 22 Feb 2014 13:33:10 -0800 Subject: [PATCH 18/18] minor bug fix. --- SimPEG/Mesh/TensorView.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Mesh/TensorView.py b/SimPEG/Mesh/TensorView.py index d345f0c3..5757b00b 100644 --- a/SimPEG/Mesh/TensorView.py +++ b/SimPEG/Mesh/TensorView.py @@ -349,7 +349,7 @@ class TensorView(object): # End Interpolation out += (ax.pcolormesh(tM.vectorNx, tM.vectorNy, np.sqrt(U**2+V**2).T, vmin=clim[0], vmax=clim[1], **pcolorOpts),) - out += (plt.streamplot(tMi.vectorCCx, tMi.vectorCCy, Ui.T, Vi.T, **streamOpts),) + out += (ax.streamplot(tMi.vectorCCx, tMi.vectorCCy, Ui.T, Vi.T, **streamOpts),) if grid: xXGrid = np.c_[tM.vectorNx,tM.vectorNx,np.nan*np.ones(tM.nNx)].flatten()