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 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): 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 b759e50a..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" @@ -131,9 +129,10 @@ def readUBCTensorMesh(fileName): x0 = mesh[1][0] y0 = mesh[1][1] - z0 = -mesh[1][2] + 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