Merge branch 'master' of https://github.com/simpeg/simpeg

Conflicts:
	SimPEG/Utils/meshutils.py

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)

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):

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)
+

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