From d8863c7c4c577fe4f0f0449378c04b72d5104641 Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 2 May 2014 12:02:46 -0700 Subject: [PATCH] Put changes back in which seogi reversed. Changed dependent property. --- SimPEG/Maps.py | 43 ++++++++++++++++++++++++++++++++++++--- SimPEG/Survey.py | 27 +++++++++++++++++++++--- SimPEG/Utils/codeutils.py | 2 +- SimPEG/Utils/matutils.py | 11 ++++++++++ SimPEG/Utils/meshutils.py | 7 +++---- 5 files changed, 79 insertions(+), 11 deletions(-) diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index a32ff13b..72386695 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -368,10 +368,47 @@ class ComboMap(IdentityMap): mi = map_i.transform(mi) return deriv +class ComplexMap(IdentityMap): + """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._nP + + def transform(self, m): + nC = self.mesh.nC + return m[:nC] + m[nC:]*1j + + 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__': 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) diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 34b803d5..450f7464 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -57,6 +57,7 @@ class BaseTimeRx(BaseRx): """SimPEG Receiver Object""" times = None #: Times when the receivers were active. + projTLoc = 'N' def __init__(self, locs, times, rxType, **kwargs): self.times = times @@ -67,6 +68,26 @@ class BaseTimeRx(BaseRx): """Number of data in the receiver.""" return self.locs.shape[0] * len(self.times) + def getSpatialP(self, mesh): + """ + Returns the spatial projection matrix. + + .. note:: + + This is not stored in memory, but is created on demand. + """ + return mesh.getInterpolationMat(self.locs, self.projGLoc) + + def getTimeP(self, timeMesh): + """ + Returns the time projection matrix. + + .. note:: + + This is not stored in memory, but is created on demand. + """ + return timeMesh.getInterpolationMat(self.times, self.projTLoc) + def getP(self, mesh, timeMesh): """ Returns the projection matrices as a @@ -75,13 +96,13 @@ class BaseTimeRx(BaseRx): .. note:: - Projection matrices are stored as a dictionary (mesh, timeMesh) + Projection matrices are stored as a dictionary (mesh, timeMesh) if storeProjections is True """ if (mesh, timeMesh) in self._Ps: return self._Ps[(mesh, timeMesh)] - Ps = mesh.getInterpolationMat(self.locs, self.projGLoc) - Pt = timeMesh.getInterpolationMat(self.times, 'N') + Ps = self.getSpatialP(mesh) + Pt = self.getTimeP(timeMesh) P = sp.kron(Pt, Ps) if self.storeProjections: diff --git a/SimPEG/Utils/codeutils.py b/SimPEG/Utils/codeutils.py index 5762117d..c1d898a6 100644 --- a/SimPEG/Utils/codeutils.py +++ b/SimPEG/Utils/codeutils.py @@ -128,7 +128,7 @@ def callHooks(match, mainFirst=False): def dependentProperty(name, value, children, doc): def fget(self): return getattr(self,name,value) def fset(self, val): - if isScalar(val) and getattr(self,name,value) == val: + if (isScalar(val) and getattr(self,name,value) == val) or val is getattr(self,name,value): return # it is the same! for child in children: if hasattr(self, child): 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) + diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index dafbe0bb..1cd8bb08 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -3,8 +3,6 @@ from scipy import sparse as sp from matutils import mkvc, ndgrid, sub2ind, sdiag from codeutils import asArray_N_x_Dim from codeutils import isScalar -import SimPEG - def exampleLrmGrid(nC, exType): assert type(nC) == list, "nC must be a list containing the number of nodes" @@ -133,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 @@ -148,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