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 ../ diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 9717affa..9a68d12a 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -325,7 +325,7 @@ class TensorMesh(BaseTensorMesh, TensorView, DiffOperators, InnerProducts): # --------------- Methods --------------------- - def isInside(self, pts): + def isInside(self, pts, locType='N'): """ Determines if a set of points are inside a mesh. @@ -334,15 +334,23 @@ class TensorMesh(BaseTensorMesh, 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 @@ -362,8 +370,21 @@ class TensorMesh(BaseTensorMesh, 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: @@ -375,7 +396,9 @@ class TensorMesh(BaseTensorMesh, 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/Mesh/View.py b/SimPEG/Mesh/View.py index ef670d83..b53cf074 100644 --- a/SimPEG/Mesh/View.py +++ b/SimPEG/Mesh/View.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 += (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() diff --git a/SimPEG/Model.py b/SimPEG/Model.py index 5ea87720..f6d4d6d4 100644 --- a/SimPEG/Model.py +++ b/SimPEG/Model.py @@ -174,12 +174,109 @@ 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 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(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(self.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.""" 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): diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index ed3d03d1..4fbf5d4b 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -59,7 +59,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 +79,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): 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/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: diff --git a/SimPEG/Tests/test_utils.py b/SimPEG/Tests/test_utils.py index fea231f2..2732c311 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]))) @@ -93,7 +106,7 @@ class TestSequenceFunctions(unittest.TestCase): sp.hstack((sdiag(a[2]), sdiag(a[3]))))) Z2 = B*A - sp.eye(10, 10) - self.assertTrue(np.linalg.norm(Z2.todense().ravel(), 2) < 1e-12) + self.assertTrue(np.linalg.norm(Z2.todense().ravel(), 2) < 1e-10) a = [np.random.rand(5, 1) for i in range(9)] B = inv3X3BlockDiagonal(*a) @@ -104,7 +117,7 @@ class TestSequenceFunctions(unittest.TestCase): Z3 = B*A - sp.eye(15, 15) - self.assertTrue(np.linalg.norm(Z3.todense().ravel(), 2) < 1e-12) + self.assertTrue(np.linalg.norm(Z3.todense().ravel(), 2) < 1e-10) diff --git a/SimPEG/Utils/interputils.py b/SimPEG/Utils/interputils.py index 38308492..3f9b7e71 100644 --- a/SimPEG/Utils/interputils.py +++ b/SimPEG/Utils/interputils.py @@ -20,9 +20,17 @@ 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 + + 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): @@ -70,14 +78,17 @@ 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 - return Q.tocsr() + + for I, V in zip(inds, vals): + Q[i, I] += V + # Q[i, mkvc(inds)] = vals + + return Q @@ -91,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 @@ -112,9 +120,13 @@ def _interpmat2D(locs, x, y): (1-dx2/Dx)*(1-dy1/Dy), (1-dx2/Dx)*(1-dy2/Dy)] - Q[i, mkvc(inds)] = vals - return Q.tocsr() + for I, V in zip(mkvc(inds), vals): + Q[i, I] += V + # Q[i, mkvc(inds)] = vals + + + return Q @@ -129,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 @@ -160,20 +168,21 @@ 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.tocsr() + 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) + Q = Utils.interpmat(locs, x) plt.plot(x, fun(x), 'bs-') plt.plot(dense, fun(dense), 'y:') plt.plot(locs, Q*fun(x), 'mo') 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"""