From ded7a3967c9ae8420a13a8f7372f6a12bbf8e0fd Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Sun, 27 Apr 2014 20:51:03 -0700 Subject: [PATCH 1/5] Updates to survey fields class. --- SimPEG/Survey.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 7954bc21..34b803d5 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -229,7 +229,7 @@ class Fields(object): dtype = self.dtype[name] else: dtype = self.dtype - field = np.empty(self._storageShape(nP), dtype=dtype) + field = np.zeros(self._storageShape(nP), dtype=dtype) self._fields[name] = field @@ -280,17 +280,13 @@ class Fields(object): assert type(value) is dict, 'New fields must be a dictionary, if field is not specified.' newFields = value elif name in self.knownFields: - assert type(value) is np.ndarray, 'Must be set to a numpy array' newFields = {name: value} else: raise Exception('Unknown setter') for name in newFields: field = self._initStore(name) - NEWF = newFields[name] - if field.shape[1] == 1 or NEWF.ndim == 1: - NEWF = Utils.mkvc(NEWF,2) - self._setField(field, NEWF, ind) + self._setField(field, newFields[name], ind) def __getitem__(self, key): ind, name = self._indexAndNameFromKey(key) @@ -302,6 +298,8 @@ class Fields(object): return self._getField(name, ind) def _setField(self, field, val, ind): + if type(val) is np.ndarray and (field.shape[1] == 1 or val.ndim == 1): + val = Utils.mkvc(val,2) field[:,ind] = val def _getField(self, name, ind): @@ -310,6 +308,9 @@ class Fields(object): out = Utils.mkvc(out) return out + def __contains__(self, other): + return self._fields.__contains__(other) + class TimeFields(Fields): """Fancy Field Storage for time domain problems @@ -342,6 +343,8 @@ class TimeFields(Fields): def _setField(self, field, val, ind): txInd, timeInd = ind + if type(val) is np.ndarray and (val.ndim == 1): + val = Utils.mkvc(val,2) field[:,txInd,timeInd] = val def _getField(self, name, ind): From 96e45ea3b0860da2b08d5833d4ad620a125f9395 Mon Sep 17 00:00:00 2001 From: seogi Date: Mon, 28 Apr 2014 10:07:14 -0700 Subject: [PATCH 2/5] modify ubcmeshutils --- SimPEG/Utils/meshutils.py | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index 43efcc33..f893bf2d 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -1,9 +1,13 @@ import numpy as np from scipy import sparse as sp from matutils import mkvc, ndgrid, sub2ind, sdiag +<<<<<<< ded7a3967c9ae8420a13a8f7372f6a12bbf8e0fd from codeutils import asArray_N_x_Dim from codeutils import isScalar +======= +import SimPEG +>>>>>>> 6e4308fea2b6ca5756fbf45ec4ebb801125345ad def exampleLrmGrid(nC, exType): assert type(nC) == list, "nC must be a list containing the number of nodes" @@ -130,9 +134,9 @@ def readUBCTensorMesh(fileName): x0 = mesh[1][0] y0 = mesh[1][1] - z0 = -mesh[1][2] + z0 = -(hz.sum()-mesh[1][2]) - mesh3D = Mesh.TensorMesh([hx, hy, hz], np.r_[x0, y0, z0]) + mesh3D = SimPEG.Mesh.TensorMesh([hx, hy, hz], np.r_[x0, y0, z0]) return mesh3D @@ -147,7 +151,7 @@ def readUBCTensorModel(fileName, mesh): model = np.reshape(model, (mesh.nCz, mesh.nCx, mesh.nCy), order = 'F') model = model[::-1,:,:] model = np.transpose(model, (1, 2, 0)) - model = Utils.mkvc(model) + model = SimPEG.Utils.mkvc(model) return model From c814029d89af32011cadc85f446ca42257bc237a Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Mon, 28 Apr 2014 10:41:28 -0700 Subject: [PATCH 3/5] Fix commit issues in MeshUtils. --- SimPEG/Utils/meshutils.py | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index f893bf2d..1cd8bb08 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -1,14 +1,9 @@ import numpy as np from scipy import sparse as sp from matutils import mkvc, ndgrid, sub2ind, sdiag -<<<<<<< ded7a3967c9ae8420a13a8f7372f6a12bbf8e0fd from codeutils import asArray_N_x_Dim from codeutils import isScalar -======= -import SimPEG ->>>>>>> 6e4308fea2b6ca5756fbf45ec4ebb801125345ad - def exampleLrmGrid(nC, exType): assert type(nC) == list, "nC must be a list containing the number of nodes" assert len(nC) == 2 or len(nC) == 3, "nC must either two or three dimensions" @@ -136,7 +131,8 @@ def readUBCTensorMesh(fileName): y0 = mesh[1][1] z0 = -(hz.sum()-mesh[1][2]) - mesh3D = SimPEG.Mesh.TensorMesh([hx, hy, hz], np.r_[x0, y0, z0]) + from SimPEG import Mesh + mesh3D = Mesh.TensorMesh([hx, hy, hz], np.r_[x0, y0, z0]) return mesh3D @@ -151,7 +147,7 @@ def readUBCTensorModel(fileName, mesh): model = np.reshape(model, (mesh.nCz, mesh.nCx, mesh.nCy), order = 'F') model = model[::-1,:,:] model = np.transpose(model, (1, 2, 0)) - model = SimPEG.Utils.mkvc(model) + model = mkvc(model) return model From c2032c3f31f36caf18a44a647561f071a3f2f3fa Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 30 Apr 2014 10:24:06 -0700 Subject: [PATCH 4/5] ComplexMap --- SimPEG/Maps.py | 37 ++++++++++++++++++++++++++++++++++--- 1 file changed, 34 insertions(+), 3 deletions(-) diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index a32ff13b..c5f1c1a5 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -368,10 +368,41 @@ class ComboMap(IdentityMap): mi = map_i.transform(mi) return deriv +class ComplexMap(IdentityMap): + """docstring for ComplexMap""" + def __init__(self, mesh): + IdentityMap.__init__(self, mesh) + + @property + def nP(self): + return self.mesh.nC * 2 + + def transform(self, m): + nC = self.mesh.nC + return m[:nC] + m[nC:]*1j + + def transformDeriv(self, m, v=None, adjoint=False): + nC = self.mesh.nC + if v is None and adjoint is False: + return sp.hstack((sp.identity(nC), np.identity(nC,dtype=complex)*1j)) + + assert v is not None, 'Must have a vector to multiply by.' + + if adjoint is False: + return sp.hstack((sp.identity(nC), np.identity(nC,dtype=complex)*1j)) * v + elif adjoint is True: + return np.r_[v.real,v.imag] + + if __name__ == '__main__': from SimPEG import * + # mesh = Mesh.TensorMesh([10,8]) + # combo = ComboMap(mesh, [ExpMap, Vertical1DMap]) + # m = combo.example() + # print m.shape + # print combo.test(np.arange(8)) mesh = Mesh.TensorMesh([10,8]) - combo = ComboMap(mesh, [ExpMap, Vertical1DMap]) - m = combo.example() + mapping = ComplexMap(mesh) + m = mapping.example() print m.shape - print combo.test(np.arange(8)) + print mapping.test(m) From 6adcce86c8d1bb04b8a89712c82cf36ad9e80c2a Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Wed, 30 Apr 2014 11:48:10 -0700 Subject: [PATCH 5/5] Added SimPEGLinearOperator which has transpose functionality. Changed ComplexMap to use this operator. --- SimPEG/Maps.py | 32 +++++++++++++++++++------------- SimPEG/Utils/matutils.py | 11 +++++++++++ 2 files changed, 30 insertions(+), 13 deletions(-) diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index c5f1c1a5..72386695 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -369,29 +369,35 @@ class ComboMap(IdentityMap): return deriv class ComplexMap(IdentityMap): - """docstring for ComplexMap""" - def __init__(self, mesh): + """docstring for ComplexMap + + default nP is nC in the mesh times 2 [real, imag] + + """ + def __init__(self, mesh, nP=None): IdentityMap.__init__(self, mesh) + if nP is not None: + assert nP%2 == 0, 'nP must be even.' + self._nP = nP or (self.mesh.nC * 2) @property def nP(self): - return self.mesh.nC * 2 + return self._nP def transform(self, m): nC = self.mesh.nC return m[:nC] + m[nC:]*1j - def transformDeriv(self, m, v=None, adjoint=False): - nC = self.mesh.nC - if v is None and adjoint is False: - return sp.hstack((sp.identity(nC), np.identity(nC,dtype=complex)*1j)) - - assert v is not None, 'Must have a vector to multiply by.' - - if adjoint is False: - return sp.hstack((sp.identity(nC), np.identity(nC,dtype=complex)*1j)) * v - elif adjoint is True: + def transformDeriv(self, m): + nC = self.nP/2 + shp = (nC, nC*2) + def fwd(v): + return v[:nC] + v[nC:]*1j + def adj(v): return np.r_[v.real,v.imag] + return Utils.SimPEGLinearOperator(shp,fwd,adj) + + transformInverse = transformDeriv if __name__ == '__main__': diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py index 2d3b87f7..9fedfce8 100644 --- a/SimPEG/Utils/matutils.py +++ b/SimPEG/Utils/matutils.py @@ -330,3 +330,14 @@ def invPropertyTensor(M, tensor, returnMatrix=False): return makePropertyTensor(M, T) return T + + + +from scipy.sparse.linalg import LinearOperator + +class SimPEGLinearOperator(LinearOperator): + """Extends scipy.sparse.linalg.LinearOperator to have a .T function.""" + @property + def T(self): + return self.__class__((self.shape[1],self.shape[0]),self.rmatvec,rmatvec=self.matvec,matmat=self.matmat) +