From 9ae281f368e33897988c241d7ddacf2645767c05 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 18 Feb 2014 12:53:22 -0800 Subject: [PATCH 01/13] name space bug. --- 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 d8a215b7..b91f5fae 100644 --- a/SimPEG/Mesh/TensorView.py +++ b/SimPEG/Mesh/TensorView.py @@ -316,7 +316,7 @@ class TensorView(object): if normalizeVector: length = np.sqrt(U**2+V**2) U, V = U/length, V/length - out += (quiver( tM.gridCC[:,0], tM.gridCC[:,1], U, V),) + out += (plt.quiver( tM.gridCC[:,0], tM.gridCC[:,1], U, V),) if grid: xXGrid = np.c_[tM.vectorNx,tM.vectorNx,np.nan*np.ones(tM.nNx)].flatten() From 3f98938d3ef01241cfc2a5162c2e5d2e98af6495 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Tue, 18 Feb 2014 14:20:21 -0800 Subject: [PATCH 02/13] quiver --> streamline --- SimPEG/Mesh/TensorView.py | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/SimPEG/Mesh/TensorView.py b/SimPEG/Mesh/TensorView.py index b91f5fae..ef0957eb 100644 --- a/SimPEG/Mesh/TensorView.py +++ b/SimPEG/Mesh/TensorView.py @@ -243,7 +243,7 @@ class TensorView(object): ax=None, clim=None, showIt=False, gridOpts={'color':'k'} ): - viewOpts = ['real','imag','abs','vec','|vec|'] + viewOpts = ['real','imag','abs','vec'] normalOpts = ['X', 'Y', 'Z'] vTypeOpts = ['CC','F','E'] @@ -258,11 +258,6 @@ class TensorView(object): if ind is None: ind = int(szSliceDim/2) assert type(ind) in [int, long], 'ind must be an integer' - normalizeVector = False - if view == '|vec|': - view = 'vec' - normalizeVector = True - if ax is None: fig = plt.figure(1) fig.clf() @@ -313,10 +308,7 @@ class TensorView(object): 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]),) - if normalizeVector: - length = np.sqrt(U**2+V**2) - U, V = U/length, V/length - out += (plt.quiver( tM.gridCC[:,0], tM.gridCC[:,1], U, V),) + out += (plt.streamplot(tM.vectorCCx, tM.vectorCCy, U.T, V.T),) if grid: xXGrid = np.c_[tM.vectorNx,tM.vectorNx,np.nan*np.ones(tM.nNx)].flatten() From fdfc1c600e73cddf754f8af8e675f802a219405e Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 19 Feb 2014 18:49:30 -0800 Subject: [PATCH 03/13] 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 04/13] 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 05/13] 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 06/13] 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 07/13] 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 08/13] 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 09/13] 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 10/13] 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 11/13] 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 12/13] 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 13/13] 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 ../