Fixed bugs in meshutils, writing VTR files

Allow fileName to be None, which will output the VTKobject without save a file.

SimPEG/DataMisfit.py

@@ -14,6 +14,34 @@ class BaseDataMisfit(object):
     debug   = False  #: Print debugging information
     counter = None   #: Set this to a SimPEG.Utils.Counter() if you want to count things
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     def __init__(self, survey, **kwargs):
         assert survey.ispaired, 'The survey must be paired to a problem.'
         if isinstance(survey, Survey.BaseSurvey):

SimPEG/Directives.py

@@ -8,6 +8,34 @@ class InversionDirective(object):
     def __init__(self, **kwargs):
         Utils.setKwargs(self, **kwargs)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     @property
     def inversion(self):
         """This is the inversion of the InversionDirective instance."""
@@ -144,23 +172,139 @@ class BetaSchedule(InversionDirective):
             if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter
             self.invProb.beta /= self.coolingFactor
 
-
-
-class SaveModelEveryIteration(InversionDirective):
-    """SaveModelEveryIteration"""
+class TargetMisfit(InversionDirective):
 
     @property
-    def modelName(self):
-        if getattr(self, '_modelName', None) is None:
-            from datetime import datetime
-            self._modelName = 'inversionModel-%s'%datetime.now().strftime('%Y-%m-%d')
-        return self._modelName
-    @modelName.setter
-    def modelName(self, value):
-        self._modelName = value
+    def target(self):
+        if getattr(self, '_target', None) is None:
+            self._target = self.survey.nD
+        return self._target
+    @target.setter
+    def target(self, val):
+        self._target = val
 
     def endIter(self):
-        np.save('%03d-%s' % (self.opt.iter, self.modelName), self.opt.xc)
+        if self.invProb.phi_d < self.target:
+            self.opt.stopNextIteration = True
+
+
+
+class _SaveEveryIteration(InversionDirective):
+    @property
+    def name(self):
+        if getattr(self, '_name', None) is None:
+            self._name = 'InversionModel'
+        return self._name
+    @name.setter
+    def name(self, value):
+        self._name = value
+
+    @property
+    def fileName(self):
+        if getattr(self, '_fileName', None) is None:
+            from datetime import datetime
+            self._fileName = '%s-%s'%(self.name, datetime.now().strftime('%Y-%m-%d-%H-%M'))
+        return self._fileName
+    @fileName.setter
+    def fileName(self, value):
+        self._fileName = value
+
+
+class SaveModelEveryIteration(_SaveEveryIteration):
+    """SaveModelEveryIteration"""
+
+    def initialize(self):
+        print "SimPEG.SaveModelEveryIteration will save your models as: '###-%s.npy'"%self.fileName
+
+    def endIter(self):
+        np.save('%03d-%s' % (self.opt.iter, self.fileName), self.opt.xc)
+
+
+class SaveOutputEveryIteration(_SaveEveryIteration):
+    """SaveModelEveryIteration"""
+
+    def initialize(self):
+        print "SimPEG.SaveOutputEveryIteration will save your inversion progress as: '###-%s.txt'"%self.fileName
+        f = open(self.fileName+'.txt', 'w')
+        f.write("  #     beta     phi_d     phi_m       f\n")
+        f.close()
+
+    def endIter(self):
+        f = open(self.fileName+'.txt', 'a')
+        f.write(' %3d %1.4e %1.4e %1.4e %1.4e\n'%(self.opt.iter, self.invProb.beta, self.invProb.phi_d, self.invProb.phi_m, self.opt.f))
+        f.close()
+
+
+class SaveOutputDictEveryIteration(_SaveEveryIteration):
+    """SaveOutputDictEveryIteration"""
+
+    def initialize(self):
+        print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName
+
+    def endIter(self):
+        # Save the data.
+        ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
+        phi_ms = 0.5*ms.dot(ms)
+        if self.reg.smoothModel == True:
+            mref = self.reg.mref
+        else:
+            mref = 0
+        mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
+        phi_mx = 0.5 * mx.dot(mx)
+        if self.prob.mesh.dim==2:
+            my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
+            phi_my = 0.5 * my.dot(my)
+        else:
+            phi_my = 'NaN'
+        if self.prob.mesh.dim==3:
+            mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
+            phi_mz = 0.5 * mz.dot(mz)
+        else:
+            phi_mz = 'NaN'
+
+
+        # Save the file as a npz
+        np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel)
+
+
+class SaveOutputDictEveryIteration(_SaveEveryIteration):
+    """SaveOutputDictEveryIteration
+
+    A directive that saves some relevant information from the inversion run to a numpy .npz dictionary file (see numpy.savez function for further info).
+
+
+
+    """
+
+    def initialize(self):
+        print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName
+
+    def endIter(self):
+        # Save the data.
+        ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
+        phi_ms = 0.5*ms.dot(ms)
+        if self.reg.smoothModel == True:
+            mref = self.reg.mref
+        else:
+            mref = 0
+        mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
+        phi_mx = 0.5 * mx.dot(mx)
+        if self.prob.mesh.dim==2:
+            my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
+            phi_my = 0.5 * my.dot(my)
+        else:
+            phi_my = 'NaN'
+        if self.prob.mesh.dim==3:
+            mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
+            phi_mz = 0.5 * mz.dot(mz)
+        else:
+            phi_mz = 'NaN'
+
+
+        # Save the file as a npz
+        np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
+
+
 
 
 # class UpdateReferenceModel(Parameter):

SimPEG/Fields.py

@@ -0,0 +1,289 @@
+import Utils, numpy as np, scipy.sparse as sp
+
+class Fields(object):
+    """Fancy Field Storage
+
+        u[:,'phi'] = phi
+        print u[src0,'phi']
+
+    """
+
+    knownFields = None  #: Known fields,   a dict with locations,         e.g. {"e": "E", "phi": "CC"}
+    aliasFields = None  #: Aliased fields, a dict with [alias, location, function], e.g. {"b":["e","F",lambda(F,e,ind)]}
+    dtype = float       #: dtype is the type of the storage matrix. This can be a dictionary.
+
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
+    def __init__(self, mesh, survey, **kwargs):
+        self.survey = survey
+        self.mesh = mesh
+        Utils.setKwargs(self, **kwargs)
+        self._fields = {}
+
+        if self.knownFields is None:
+            raise Exception('knownFields cannot be set to None')
+        if self.aliasFields is None:
+            self.aliasFields = {}
+
+        allFields = [k for k in self.knownFields] + [a for a in self.aliasFields]
+        assert len(allFields) == len(set(allFields)), 'Aliased fields and Known Fields have overlapping definitions.'
+        self.startup()
+
+    def startup(self):
+        pass
+
+    @property
+    def approxSize(self):
+        """The approximate cost to storing all of the known fields."""
+        sz = 0.0
+        for f in self.knownFields:
+            loc =self.knownFields[f]
+            sz += np.array(self._storageShape(loc)).prod()*8.0/(1024**2)
+        return "%e MB"%sz
+
+    def _storageShape(self, loc):
+        nSrc = self.survey.nSrc
+
+        nP = {'CC': self.mesh.nC,
+              'N':  self.mesh.nN,
+              'F':  self.mesh.nF,
+              'E':  self.mesh.nE}[loc]
+
+        return (nP, nSrc)
+
+    def _initStore(self, name):
+        if name in self._fields:
+            return self._fields[name]
+
+        assert name in self.knownFields, 'field name is not known.'
+
+        loc = self.knownFields[name]
+
+        if type(self.dtype) is dict:
+            dtype = self.dtype[name]
+        else:
+            dtype = self.dtype
+        field = np.zeros(self._storageShape(loc), dtype=dtype)
+
+        self._fields[name] = field
+
+        return field
+
+    def _srcIndex(self, srcTestList):
+        if type(srcTestList) is slice:
+            ind = srcTestList
+        else:
+            ind = self.survey.getSourceIndex(srcTestList)
+        return ind
+
+    def _nameIndex(self, name, accessType):
+
+        if type(name) is slice:
+            assert name == slice(None,None,None), 'Fancy field name slicing is not supported... yet.'
+            name = None
+
+        if name is None:
+            return
+        if accessType=='set' and name not in self.knownFields:
+            if name in self.aliasFields:
+                raise KeyError("Invalid field name (%s) for setter, you can't set an aliased property"%name)
+            else:
+                raise KeyError('Invalid field name (%s) for setter'%name)
+
+        elif accessType=='get' and (name not in self.knownFields and name not in self.aliasFields):
+            raise KeyError('Invalid field name (%s) for getter'%name)
+        return name
+
+    def _indexAndNameFromKey(self, key, accessType):
+        if type(key) is not tuple:
+            key = (key,)
+        if len(key) == 1:
+            key += (None,)
+
+        assert len(key) == 2, 'must be [Src, fieldName]'
+
+        srcTestList, name = key
+        name = self._nameIndex(name, accessType)
+        ind = self._srcIndex(srcTestList)
+        return ind, name
+
+    def __setitem__(self, key, value):
+        ind, name = self._indexAndNameFromKey(key, 'set')
+        if name is None:
+            freq = key
+            assert type(value) is dict, 'New fields must be a dictionary, if field is not specified.'
+            newFields = value
+        elif name in self.knownFields:
+            newFields = {name: value}
+        else:
+            raise Exception('Unknown setter')
+
+        for name in newFields:
+            field = self._initStore(name)
+            self._setField(field, newFields[name], name, ind)
+
+    def __getitem__(self, key):
+        ind, name = self._indexAndNameFromKey(key, 'get')
+        if name is None:
+            out = {}
+            for name in self._fields:
+                out[name] = self._getField(name, ind)
+            return out
+        return self._getField(name, ind)
+
+    def _setField(self, field, val, name, ind):
+        if isinstance(val, np.ndarray) and (field.shape[0] == field.size or val.ndim == 1):
+            val = Utils.mkvc(val,2)
+        field[:,ind] = val
+
+    def _getField(self, name, ind):
+        if name in self._fields:
+            out = self._fields[name][:,ind]
+        else:
+            # Aliased fields
+            alias, loc, func = self.aliasFields[name]
+
+            srcII = np.array(self.survey.srcList)[ind]
+            srcII = srcII.tolist()
+
+            if type(func) is str:
+                assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.'
+                func = getattr(self, func)
+            out = func(self._fields[alias][:,ind], srcII)
+        if out.shape[0] == out.size or out.ndim == 1:
+            out = Utils.mkvc(out,2)
+        return out
+
+    def __contains__(self, other):
+        if other in self.aliasFields:
+            other = self.aliasFields[other][0]
+        return self._fields.__contains__(other)
+
+
+class TimeFields(Fields):
+    """Fancy Field Storage for time domain problems
+
+        u[:,'phi', timeInd] = phi
+        print u[src0,'phi']
+
+    """
+
+    def _storageShape(self, loc):
+        nP = {'CC': self.mesh.nC,
+              'N':  self.mesh.nN,
+              'F':  self.mesh.nF,
+              'E':  self.mesh.nE}[loc]
+        nSrc = self.survey.nSrc
+        nT = self.survey.prob.nT + 1
+        return (nP, nSrc, nT)
+
+    def _indexAndNameFromKey(self, key, accessType):
+        if type(key) is not tuple:
+            key = (key,)
+        if len(key) == 1:
+            key += (None,)
+        if len(key) == 2:
+            key += (slice(None,None,None),)
+
+        assert len(key) == 3, 'must be [Src, fieldName, times]'
+
+        srcTestList, name, timeInd = key
+
+        name = self._nameIndex(name, accessType)
+        srcInd = self._srcIndex(srcTestList)
+
+        return (srcInd, timeInd), name
+
+    def _correctShape(self, name, ind, deflate=False):
+        srcInd, timeInd = ind
+        if name in self.knownFields:
+            loc = self.knownFields[name]
+        else:
+            loc = self.aliasFields[name][1]
+        nP, total_nSrc, total_nT = self._storageShape(loc)
+        nSrc = np.ones(total_nSrc, dtype=bool)[srcInd].sum()
+        nT  = np.ones(total_nT, dtype=bool)[timeInd].sum()
+        shape = nP, nSrc, nT
+        if deflate:
+             shape = tuple([s for s in shape if s > 1])
+        if len(shape) == 1:
+            shape = shape + (1,)
+        return shape
+
+    def _setField(self, field, val, name, ind):
+        srcInd, timeInd = ind
+        shape = self._correctShape(name, ind)
+        if Utils.isScalar(val):
+            field[:,srcInd,timeInd] = val
+            return
+        if val.size != np.array(shape).prod():
+            raise ValueError('Incorrect size for data.')
+        correctShape = field[:,srcInd,timeInd].shape
+        field[:,srcInd,timeInd] = val.reshape(correctShape, order='F')
+
+    def _getField(self, name, ind):
+        srcInd, timeInd = ind
+
+        if name in self._fields:
+            out = self._fields[name][:,srcInd,timeInd]
+        else:
+            # Aliased fields
+            alias, loc, func = self.aliasFields[name]
+            if type(func) is str:
+                assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.'
+                func = getattr(self, func)
+            pointerFields = self._fields[alias][:,srcInd,timeInd]
+            pointerShape = self._correctShape(alias, ind)
+            pointerFields = pointerFields.reshape(pointerShape, order='F')
+
+            timeII = np.arange(self.survey.prob.nT + 1)[timeInd]
+            srcII  = np.array(self.survey.srcList)[srcInd]
+            srcII  = srcII.tolist()
+
+            if timeII.size == 1:
+                pointerShapeDeflated = self._correctShape(alias, ind, deflate=True)
+                pointerFields = pointerFields.reshape(pointerShapeDeflated, order='F')
+                out = func(pointerFields, srcII, timeII)
+            else: #loop over the time steps
+                nT = pointerShape[2]
+                out = range(nT)
+                for i, TIND_i in enumerate(timeII):
+                    fieldI = pointerFields[:,:,i]
+                    if fieldI.shape[0] == fieldI.size:
+                        fieldI = Utils.mkvc(fieldI, 2)
+                    out[i] = func(fieldI, srcII, TIND_i)
+                    if out[i].ndim == 1:
+                        out[i] = out[i][:,np.newaxis,np.newaxis]
+                    elif out[i].ndim == 2:
+                        out[i] = out[i][:,:,np.newaxis]
+                out = np.concatenate(out, axis=2)
+
+        shape = self._correctShape(name, ind, deflate=True)
+        return out.reshape(shape, order='F')
+

SimPEG/Maps.py

@@ -1,5 +1,6 @@
 import Utils, numpy as np, scipy.sparse as sp
 from Tests import checkDerivative
+from PropMaps import PropMap, Property
 
 
 class IdentityMap(object):
@@ -16,12 +17,38 @@ class IdentityMap(object):
         Utils.setKwargs(self, **kwargs)
         self.mesh = mesh
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+
     @property
     def nP(self):
         """
             :rtype: int
             :return: number of parameters in the model
         """
+        if self.mesh is None:
+            return '*'
         return self.mesh.nC
 
     @property
@@ -32,8 +59,11 @@ class IdentityMap(object):
             :rtype: (int,int)
             :return: shape of the operator as a tuple
         """
+        if self.mesh is None:
+            return ('*', self.nP)
         return (self.mesh.nC, self.nP)
 
+
     def _transform(self, m):
         """
             Changes the model into the physical property.
@@ -98,17 +128,17 @@ class IdentityMap(object):
 
     def __mul__(self, val):
         if isinstance(val, IdentityMap):
-            if not self.shape[1] == val.shape[0]:
+            if not (self.shape[1] == '*' or val.shape[0] == '*') and not self.shape[1] == val.shape[0]:
                 raise ValueError('Dimension mismatch in %s and %s.' % (str(self), str(val)))
             return ComboMap([self, val])
         elif isinstance(val, np.ndarray):
-            if not self.shape[1] == val.shape[0]:
+            if not self.shape[1] == '*' and not self.shape[1] == val.shape[0]:
                 raise ValueError('Dimension mismatch in %s and np.ndarray%s.' % (str(self), str(val.shape)))
             return self._transform(val)
         raise Exception('Unrecognized data type to multiply. Try a map or a numpy.ndarray!')
 
     def __str__(self):
-        return "%s(%d,%d)" % (self.__class__.__name__, self.shape[0], self.shape[1])
+        return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1])
 
 class ComboMap(IdentityMap):
     """Combination of various maps."""
@@ -119,10 +149,10 @@ class ComboMap(IdentityMap):
         self.maps = []
         for ii, m in enumerate(maps):
             assert isinstance(m, IdentityMap), 'Unrecognized data type, inherit from an IdentityMap or ComboMap!'
-            if ii > 0 and not self.shape[1] == m.shape[0]:
+            if ii > 0 and not (self.shape[1] == '*' or m.shape[0] == '*') and not self.shape[1] == m.shape[0]:
                 prev = self.maps[-1]
-                errArgs = (prev.__name__, prev.shape[0], prev.shape[1], m.__name__, m.shape[0], m.shape[1])
-                raise ValueError('Dimension mismatch in map[%s] (%i, %i) and map[%s] (%i, %i).' % errArgs)
+                errArgs = (prev.__class__.__name__, prev.shape[0], prev.shape[1], m.__class__.__name__, m.shape[0], m.shape[1])
+                raise ValueError('Dimension mismatch in map[%s] (%s, %s) and map[%s] (%s, %s).' % errArgs)
 
             if isinstance(m, ComboMap):
                 self.maps += m.maps
@@ -155,7 +185,7 @@ class ComboMap(IdentityMap):
         return deriv
 
     def __str__(self):
-        return 'ComboMap[%s]%s' % (' * '.join([m.__str__() for m in self.maps]), str(self.shape))
+        return 'ComboMap[%s](%s,%s)' % (' * '.join([m.__str__() for m in self.maps]), self.shape[0], self.shape[1])
 
 
 class ExpMap(IdentityMap):
@@ -220,6 +250,26 @@ class ExpMap(IdentityMap):
         """
         return Utils.sdiag(np.exp(Utils.mkvc(m)))
 
+class ReciprocalMap(IdentityMap):
+    """
+        Reciprocal mapping. For example, electrical resistivity and conductivity.
+
+        .. math::
+
+            \\rho = \\frac{1}{\sigma}
+
+    """
+    def _transform(self, m):
+        return 1.0 / Utils.mkvc(m)
+
+    def inverse(self, D):
+        return 1.0 / Utils.mkvc(m)
+
+    def deriv(self, m):
+        # TODO: if this is a tensor, you might have a problem.
+        return Utils.sdiag( - Utils.mkvc(m)**(-2) )
+
+
 
 class LogMap(IdentityMap):
     """
@@ -259,6 +309,37 @@ class LogMap(IdentityMap):
     def inverse(self, m):
         return np.exp(Utils.mkvc(m))
 
+class FullMap(IdentityMap):
+    """
+    FullMap
+
+    Given a scalar, the FullMap maps the value to the
+    full model space.
+    """
+
+    def __init__(self,mesh,**kwargs):
+        IdentityMap.__init__(self, mesh,**kwargs)
+
+    @property
+    def nP(self):
+        return 1
+
+    def _transform(self, m):
+        """
+            :param m: model (scalar)
+            :rtype: numpy.array
+            :return: transformed model
+        """
+        return np.ones(self.mesh.nC)*m
+
+    def deriv(self, m):
+        """
+            :param numpy.array m: model
+            :rtype: numpy.array
+            :return: derivative of transformed model
+        """
+        return np.ones([self.mesh.nC,1])
+
 
 class Vertical1DMap(IdentityMap):
     """Vertical1DMap

SimPEG/Mesh/InnerProducts.py

@@ -33,7 +33,7 @@ class InnerProducts(object):
         return self._getInnerProduct('E', prop=prop, invProp=invProp, invMat=invMat, doFast=doFast)
 
     def _getInnerProduct(self, projType, prop=None, invProp=False, invMat=False, doFast=True):
-        """
+        """r
             :param str projType: 'F' for faces 'E' for edges
             :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6))
             :param bool invProp: inverts the material property

SimPEG/Mesh/View.py

@@ -328,6 +328,7 @@ class TensorView(object):
             v = getattr(np,view)(v) # e.g. np.real(v)
             if clim is None:
                 clim = [v.min(),v.max()]
+            v = np.ma.masked_where(np.isnan(v), v)
             out += (ax.pcolormesh(self.vectorNx, self.vectorNy, v.T, vmin=clim[0], vmax=clim[1], **pcolorOpts),)
         elif view in ['vec']:
             U, V = self.r(v.reshape((self.nC,-1), order='F'), 'CC', 'CC', 'M')

SimPEG/Models.py

@@ -31,5 +31,5 @@ class Model(np.ndarray):
     @property
     def transformDeriv(self):
         if getattr(self, '_transformDeriv', None) is None:
-            self.deriv = self.mapping.deriv(self.view(np.ndarray))
-        return self.deriv
+            self._transformDeriv = self.mapping.deriv(self.view(np.ndarray))
+        return self._transformDeriv

SimPEG/Optimization.py

@@ -97,6 +97,8 @@ class Minimize(object):
     tolG = 1e-1        #: Tolerance on gradient norm
     eps = 1e-5         #: Small value
 
+    stopNextIteration = False #: Stops the optimization program nicely.
+
     debug   = False    #: Print debugging information
     debugLS = False    #: Print debugging information for the line-search
 
@@ -113,6 +115,34 @@ class Minimize(object):
 
         Utils.setKwargs(self, **kwargs)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     @property
     def callback(self):
         return getattr(self, '_callback', None)
@@ -186,6 +216,7 @@ class Minimize(object):
                 xt, caught = self.modifySearchDirectionBreak(p)
                 if not caught: return self.xc
             self.doEndIteration(xt)
+            if self.stopNextIteration: break
 
         self.printDone()
         self.finish()
@@ -210,6 +241,7 @@ class Minimize(object):
 
         self.iter = 0
         self.iterLS = 0
+        self.stopNextIteration = False
 
         x0 = self.projection(x0)  # ensure that we start of feasible.
         self.x0 = x0

SimPEG/Problem.py

@@ -1,280 +1,7 @@
 import Utils, Survey, Models, numpy as np, scipy.sparse as sp
 Solver = Utils.SolverUtils.Solver
 import Maps, Mesh
-
-
-class Fields(object):
-    """Fancy Field Storage
-
-        u[:,'phi'] = phi
-        print u[src0,'phi']
-
-    """
-
-    knownFields = None  #: Known fields,   a dict with locations,         e.g. {"e": "E", "phi": "CC"}
-    aliasFields = None  #: Aliased fields, a dict with [alias, location, function], e.g. {"b":["e","F",lambda(F,e,ind)]}
-    dtype = float       #: dtype is the type of the storage matrix. This can be a dictionary.
-
-    def __init__(self, mesh, survey, **kwargs):
-        self.survey = survey
-        self.mesh = mesh
-        Utils.setKwargs(self, **kwargs)
-        self._fields = {}
-
-        if self.knownFields is None:
-            raise Exception('knownFields cannot be set to None')
-        if self.aliasFields is None:
-            self.aliasFields = {}
-
-        allFields = [k for k in self.knownFields] + [a for a in self.aliasFields]
-        assert len(allFields) == len(set(allFields)), 'Aliased fields and Known Fields have overlapping definitions.'
-        self.startup()
-
-    def startup(self):
-        pass
-
-    @property
-    def approxSize(self):
-        """The approximate cost to storing all of the known fields."""
-        sz = 0.0
-        for f in self.knownFields:
-            loc =self.knownFields[f]
-            sz += np.array(self._storageShape(loc)).prod()*8.0/(1024**2)
-        return "%e MB"%sz
-
-    def _storageShape(self, loc):
-        nSrc = self.survey.nSrc
-
-        nP = {'CC': self.mesh.nC,
-              'N':  self.mesh.nN,
-              'F':  self.mesh.nF,
-              'E':  self.mesh.nE}[loc]
-
-        return (nP, nSrc)
-
-    def _initStore(self, name):
-        if name in self._fields:
-            return self._fields[name]
-
-        assert name in self.knownFields, 'field name is not known.'
-
-        loc = self.knownFields[name]
-
-        if type(self.dtype) is dict:
-            dtype = self.dtype[name]
-        else:
-            dtype = self.dtype
-        field = np.zeros(self._storageShape(loc), dtype=dtype)
-
-        self._fields[name] = field
-
-        return field
-
-    def _srcIndex(self, srcTestList):
-        if type(srcTestList) is slice:
-            ind = srcTestList
-        else:
-            if type(srcTestList) is not list:
-                srcTestList = [srcTestList]
-            for srcTest in srcTestList:
-                if srcTest not in self.survey.srcList:
-                    raise KeyError('Invalid Source, not in survey list.')
-
-            ind = np.in1d(self.survey.srcList, srcTestList)
-        return ind
-
-    def _nameIndex(self, name, accessType):
-
-        if type(name) is slice:
-            assert name == slice(None,None,None), 'Fancy field name slicing is not supported... yet.'
-            name = None
-
-        if name is None:
-            return
-        if accessType=='set' and name not in self.knownFields:
-            if name in self.aliasFields:
-                raise KeyError("Invalid field name (%s) for setter, you can't set an aliased property"%name)
-            else:
-                raise KeyError('Invalid field name (%s) for setter'%name)
-
-        elif accessType=='get' and (name not in self.knownFields and name not in self.aliasFields):
-            raise KeyError('Invalid field name (%s) for getter'%name)
-        return name
-
-    def _indexAndNameFromKey(self, key, accessType):
-        if type(key) is not tuple:
-            key = (key,)
-        if len(key) == 1:
-            key += (None,)
-
-        assert len(key) == 2, 'must be [Src, fieldName]'
-
-        srcTestList, name = key
-        name = self._nameIndex(name, accessType)
-        ind = self._srcIndex(srcTestList)
-        return ind, name
-
-    def __setitem__(self, key, value):
-        ind, name = self._indexAndNameFromKey(key, 'set')
-        if name is None:
-            freq = key
-            assert type(value) is dict, 'New fields must be a dictionary, if field is not specified.'
-            newFields = value
-        elif name in self.knownFields:
-            newFields = {name: value}
-        else:
-            raise Exception('Unknown setter')
-
-        for name in newFields:
-            field = self._initStore(name)
-            self._setField(field, newFields[name], name, ind)
-
-    def __getitem__(self, key):
-        ind, name = self._indexAndNameFromKey(key, 'get')
-        if name is None:
-            out = {}
-            for name in self._fields:
-                out[name] = self._getField(name, ind)
-            return out
-        return self._getField(name, ind)
-
-    def _setField(self, field, val, name, ind):
-        if isinstance(val, np.ndarray) and (field.shape[0] == field.size or val.ndim == 1):
-            val = Utils.mkvc(val,2)
-        field[:,ind] = val
-
-    def _getField(self, name, ind):
-        if name in self._fields:
-            out = self._fields[name][:,ind]
-        else:
-            # Aliased fields
-            alias, loc, func = self.aliasFields[name]
-
-            srcII   = np.array(self.survey.srcList)[ind]
-            if isinstance(srcII, np.ndarray):
-                srcII = srcII.tolist()
-            if len(srcII) == 1:
-                srcII = srcII[0]
-
-            if type(func) is str:
-                assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.'
-                func = getattr(self, func)
-            out = func(self._fields[alias][:,ind], srcII)
-        if isinstance(out, np.ndarray) and (out.shape[0] == out.size or out.ndim == 1):
-            out = Utils.mkvc(out,2)
-        return out
-
-    def __contains__(self, other):
-        if other in self.aliasFields:
-            other = self.aliasFields[other][0]
-        return self._fields.__contains__(other)
-
-
-class TimeFields(Fields):
-    """Fancy Field Storage for time domain problems
-
-        u[:,'phi', timeInd] = phi
-        print u[src0,'phi']
-
-    """
-
-    def _storageShape(self, loc):
-        nP = {'CC': self.mesh.nC,
-              'N':  self.mesh.nN,
-              'F':  self.mesh.nF,
-              'E':  self.mesh.nE}[loc]
-        nSrc = self.survey.nSrc
-        nT = self.survey.prob.nT + 1
-        return (nP, nSrc, nT)
-
-    def _indexAndNameFromKey(self, key, accessType):
-        if type(key) is not tuple:
-            key = (key,)
-        if len(key) == 1:
-            key += (None,)
-        if len(key) == 2:
-            key += (slice(None,None,None),)
-
-        assert len(key) == 3, 'must be [Src, fieldName, times]'
-
-        srcTestList, name, timeInd = key
-
-        name = self._nameIndex(name, accessType)
-        srcInd = self._srcIndex(srcTestList)
-
-        return (srcInd, timeInd), name
-
-    def _correctShape(self, name, ind, deflate=False):
-        srcInd, timeInd = ind
-        if name in self.knownFields:
-            loc = self.knownFields[name]
-        else:
-            loc = self.aliasFields[name][1]
-        nP, total_nSrc, total_nT = self._storageShape(loc)
-        nSrc = np.ones(total_nSrc, dtype=bool)[srcInd].sum()
-        nT  = np.ones(total_nT, dtype=bool)[timeInd].sum()
-        shape = nP, nSrc, nT
-        if deflate:
-             shape = tuple([s for s in shape if s > 1])
-        if len(shape) == 1:
-            shape = shape + (1,)
-        return shape
-
-    def _setField(self, field, val, name, ind):
-        srcInd, timeInd = ind
-        shape = self._correctShape(name, ind)
-        if Utils.isScalar(val):
-            field[:,srcInd,timeInd] = val
-            return
-        if val.size != np.array(shape).prod():
-            raise ValueError('Incorrect size for data.')
-        correctShape = field[:,srcInd,timeInd].shape
-        field[:,srcInd,timeInd] = val.reshape(correctShape, order='F')
-
-    def _getField(self, name, ind):
-        srcInd, timeInd = ind
-
-        if name in self._fields:
-            out = self._fields[name][:,srcInd,timeInd]
-        else:
-            # Aliased fields
-            alias, loc, func = self.aliasFields[name]
-            if type(func) is str:
-                assert hasattr(self, func), 'The alias field function is a string, but it does not exist in the Fields class.'
-                func = getattr(self, func)
-            pointerFields = self._fields[alias][:,srcInd,timeInd]
-            pointerShape = self._correctShape(alias, ind)
-            pointerFields = pointerFields.reshape(pointerShape, order='F')
-
-            timeII = np.arange(self.survey.prob.nT + 1)[timeInd]
-            srcII   = np.array(self.survey.srcList)[srcInd]
-            if isinstance(srcII, np.ndarray):
-                srcII = srcII.tolist()
-            if len(srcII) == 1:
-                srcII = srcII[0]
-
-            if timeII.size == 1:
-                pointerShapeDeflated = self._correctShape(alias, ind, deflate=True)
-                pointerFields = pointerFields.reshape(pointerShapeDeflated, order='F')
-                out = func(pointerFields, srcII, timeII)
-            else: #loop over the time steps
-                nT = pointerShape[2]
-                out = range(nT)
-                for i, TIND_i in enumerate(timeII):
-                    fieldI = pointerFields[:,:,i]
-                    if fieldI.shape[0] == fieldI.size:
-                        fieldI = Utils.mkvc(fieldI,2)
-                    out[i] = func(fieldI, srcII, TIND_i)
-                    if out[i].ndim == 1:
-                        out[i] = out[i][:,np.newaxis,np.newaxis]
-                    elif out[i].ndim == 2:
-                        out[i] = out[i][:,:,np.newaxis]
-                out = np.concatenate(out, axis=2)
-
-        shape = self._correctShape(name, ind, deflate=True)
-        return out.reshape(shape, order='F')
-
-
+from Fields import Fields, TimeFields
 
 class BaseProblem(object):
     """
@@ -291,15 +18,55 @@ class BaseProblem(object):
     Solver = Solver   #: A SimPEG Solver class.
     solverOpts = {}   #: Sovler options as a kwarg dict
 
-    mapping = None    #: A SimPEG.Map instance.
     mesh    = None    #: A SimPEG.Mesh instance.
 
+    PropMap = None    #: A SimPEG PropertyMap class.
+
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
+    @property
+    def mapping(self):
+        "A SimPEG.Map instance or a property map is PropMap is not None"
+        return getattr(self, '_mapping', None)
+    @mapping.setter
+    def mapping(self, val):
+        if self.PropMap is None:
+            val._assertMatchesPair(self.mapPair)
+            self._mapping = val
+        else:
+            self._mapping = self.PropMap(val)
+
     def __init__(self, mesh, mapping=None, **kwargs):
         Utils.setKwargs(self, **kwargs)
         assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object."
         self.mesh = mesh
         self.mapping = mapping or Maps.IdentityMap(mesh)
-        self.mapping._assertMatchesPair(self.mapPair)
 
     @property
     def survey(self):
@@ -335,7 +102,10 @@ class BaseProblem(object):
     def curModel(self, value):
         if value is self.curModel:
             return # it is the same!
-        self._curModel = Models.Model(value, self.mapping)
+        if self.PropMap is not None:
+            self._curModel = self.mapping(value)
+        else:
+            self._curModel = Models.Model(value, self.mapping)
         for prop in self.deleteTheseOnModelUpdate:
             if hasattr(self, prop):
                 delattr(self, prop)

SimPEG/PropMaps.py

@@ -0,0 +1,337 @@
+import Utils, Maps, numpy as np, scipy.sparse as sp
+
+class Property(object):
+
+    name           = ''
+    doc            = ''
+
+    defaultVal     = None
+    defaultInvProp = False
+
+    def __init__(self, doc, **kwargs):
+        # Set the default after all other params are set
+        self.doc = doc
+        Utils.setKwargs(self, **kwargs)
+
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
+    @property
+    def propertyLink(self):
+        "Can be something like: ('sigma', Maps.ReciprocalMap)"
+        return getattr(self, '_propertyLink', None)
+    @propertyLink.setter
+    def propertyLink(self, value):
+        assert type(value) is tuple and len(value) == 2 and type(value[0]) is str and issubclass(value[1], Maps.IdentityMap), 'Use format: ("%s", Maps.ReciprocalMap)'%self.name
+        self._propertyLink = value
+
+    def _getMapProperty(self):
+        prop = self
+        def fget(self):
+            return getattr(self, '_%sMap'%prop.name, None)
+        def fset(self, val):
+            if prop.propertyLink is not None:
+                linkName, linkMap = prop.propertyLink
+                assert getattr(self, '%sMap'%linkName, None) is None, 'Cannot set both sides of a linked property.'
+            # TODO: Check if the mapping can be correct
+            setattr(self, '_%sMap'%prop.name, val)
+        return property(fget=fget, fset=fset, doc=prop.doc)
+
+    def _getIndexProperty(self):
+        prop = self
+        def fget(self):
+            return getattr(self, '_%sIndex'%prop.name, slice(None))
+        def fset(self, val):
+            setattr(self, '_%sIndex'%prop.name, val)
+        return property(fget=fget, fset=fset, doc=prop.doc)
+
+    def _getProperty(self):
+        prop = self
+        def fget(self):
+            mapping = getattr(self, '%sMap'%prop.name)
+            if mapping is None and prop.propertyLink is None:
+                return prop.defaultVal
+
+            if mapping is None and prop.propertyLink is not None:
+                linkName, linkMapClass = prop.propertyLink
+                linkMap = linkMapClass(None)
+                if getattr(self, '%sMap'%linkName, None) is None:
+                    return prop.defaultVal
+                m = getattr(self, '%s'%linkName)
+                return linkMap * m
+
+            m = getattr(self, '%sModel'%prop.name)
+            return mapping * m
+        return property(fget=fget)
+
+    def _getModelDerivProperty(self):
+        prop = self
+        def fget(self):
+            mapping = getattr(self, '%sMap'%prop.name)
+            if mapping is None and prop.propertyLink is None:
+                return None
+
+            if mapping is None and prop.propertyLink is not None:
+                linkName, linkMapClass = prop.propertyLink
+                linkedMap = getattr(self, '%sMap'%linkName)
+                if linkedMap is None:
+                    return None
+                linkMap = linkMapClass(None) * linkedMap
+                m = getattr(self, '%s'%linkName)
+                return linkMap.deriv( m )
+
+            m = getattr(self, '%sModel'%prop.name)
+            return mapping.deriv( m )
+        return property(fget=fget)
+
+    def _getModelProperty(self):
+        prop = self
+        def fget(self):
+            mapping = getattr(self, '%sMap'%prop.name)
+            if mapping is None:
+                return None
+            index = getattr(self.propMap, '%sIndex'%prop.name)
+            return self.vector[index]
+        return property(fget=fget)
+
+    def _getModelProjProperty(self):
+        prop = self
+        def fget(self):
+            mapping = getattr(self, '%sMap'%prop.name)
+            if mapping is None:
+                return None
+            inds = getattr(self.propMap, '%sIndex'%prop.name)
+            if type(inds) is slice:
+                inds = range(*inds.indices(self.nP))
+            nI, nP = len(inds),self.nP
+            return sp.csr_matrix((np.ones(nI), (range(nI), inds) ), shape=(nI, nP))
+        return property(fget=fget)
+
+    def _getModelMapProperty(self):
+        prop = self
+        def fget(self):
+            return getattr(self.propMap, '_%sMap'%prop.name, None)
+        return property(fget=fget)
+
+
+
+class PropModel(object):
+    def __init__(self, propMap, vector):
+        self.propMap = propMap
+        self.vector  = vector
+        assert len(self.vector) == self.nP
+
+    # Pickleing support methods
+    # def __reduce__(self):
+    #     return (dict,{self.propMap,self.vector})
+
+    # def __getstate__(self):
+    #     '''
+    #     Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+    #     Used when doing:
+    #         pickle.dump(pickleFile,object)
+    #     '''
+    #     self.__class__ = ProbModel
+    #     odict = {}
+    #     odict['vec'] = self.__dict__['vector']
+    #     odict['pMap'] = self.__dict__['propMap']
+    #     # Return the dict
+    #     return odict
+
+    # def __setstate__(self,odict):
+    #     '''
+    #     Function that sets a pickle dictionary in to an object.
+
+    #     Used when doing:
+    #         object = pickle.load(pickleFile)
+    #     '''
+    #     # Update the dict
+    #     # Re-hook the methods to the object
+    #     self.propMap = odict['prMap']
+    #     self.vector = odict['vec']
+
+    @property
+    def nP(self):
+        inds = []
+        if getattr(self, '_nP', None) is None:
+            for name in self.propMap._properties:
+                index = getattr(self.propMap, '%sIndex'%name, None)
+                if index is not None:
+                    if type(index) is slice:
+                        inds += range(*index.indices(len(self.vector)))
+                    else:
+                        inds += list(index)
+            self._nP = len(set(inds))
+        return self._nP
+
+    def __contains__(self, val):
+        return val in self.propMap
+
+
+
+_PROPMAPCLASSREGISTRY = {}
+
+class _PropMapMetaClass(type):
+    def __new__(cls, name, bases, attrs):
+        assert name.endswith('PropMap'), 'Please use convention: ___PropMap, e.g. ElectromagneticPropMap'
+        _properties = {}
+        for base in bases:
+            for baseProp in getattr(base, '_properties', {}):
+                _properties[baseProp] = base._properties[baseProp]
+        keys = [key for key in attrs]
+        for attr in keys:
+            if isinstance(attrs[attr], Property):
+                attrs[attr].name = attr
+                attrs[attr + 'Map'  ] = attrs[attr]._getMapProperty()
+                attrs[attr + 'Index'] = attrs[attr]._getIndexProperty()
+                _properties[attr] = attrs[attr]
+                attrs.pop(attr)
+
+        attrs['_properties'] = _properties
+
+        defaultInvProps = []
+        for p in _properties:
+            prop = _properties[p]
+            if prop.defaultInvProp:
+                defaultInvProps += [p]
+            if prop.propertyLink is not None:
+                assert prop.propertyLink[0] in _properties, "You can only link to things that exist: '%s' is trying to link to '%s'"%(prop.name, prop.propertyLink[0])
+        if len(defaultInvProps) > 1:
+            raise Exception('You have more than one default inversion property: %s' % defaultInvProps)
+
+        newClass = super(_PropMapMetaClass, cls).__new__(cls, name, bases, attrs)
+
+        newClass.PropModel = cls.createPropModelClass(newClass, name, _properties)
+
+        _PROPMAPCLASSREGISTRY[name] = newClass
+        return newClass
+
+    def createPropModelClass(self, name, _properties):
+
+        attrs = dict()
+
+        for attr in _properties:
+            prop = _properties[attr]
+
+            attrs[attr          ] = prop._getProperty()
+            attrs[attr + 'Map'  ] = prop._getModelMapProperty()
+            attrs[attr + 'Proj' ] = prop._getModelProjProperty()
+            attrs[attr + 'Model'] = prop._getModelProperty()
+            attrs[attr + 'Deriv'] = prop._getModelDerivProperty()
+
+        return type(name.replace('PropMap', 'PropModel'), (PropModel, ), attrs)
+
+
+class PropMap(object):
+    __metaclass__ = _PropMapMetaClass
+
+    def __init__(self, mappings):
+        """
+            PropMap takes a multi parameter model and maps it to the equivalent PropModel
+        """
+        if type(mappings) is dict:
+            assert np.all([k in ['maps', 'slices'] for k in mappings]), 'Dict must only have properties "maps" and "slices"'
+            self.setup(mappings['maps'], slices=mappings['slices'])
+        elif type(mappings) is list:
+            self.setup(mappings)
+        elif isinstance(mappings, Maps.IdentityMap):
+            self.setup([(self.defaultInvProp, mappings)])
+        else:
+            raise Exception('mappings must be a dict, a mapping, or a list of tuples.')
+
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        pass
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        pass
+
+    def setup(self, maps, slices=None):
+        """
+            Sets up the maps and slices for the PropertyMap
+
+
+            :param list maps: [('sigma', sigmaMap), ('mu', muMap), ...]
+            :param list slices: [('sigma', slice(0,nP)), ('mu', [1,2,5,6]), ...]
+
+        """
+        assert np.all([
+                type(m) is tuple and
+                len(m)==2 and
+                type(m[0]) is str and
+                m[0] in self._properties and
+                isinstance(m[1], Maps.IdentityMap)
+                for m in maps]), "Use signature: [%s]" % (', '.join(["('%s', %sMap)"%(p,p) for p in self._properties]))
+        if slices is None:
+            slices = dict()
+        else:
+            assert np.all([
+                s in self._properties and
+                (type(slices[s]) in [slice, list] or isinstance(slices[s], np.ndarray))
+                for s in slices]), 'Slices must be for each property'
+
+        self.clearMaps()
+
+        nP = 0
+        for name, mapping in maps:
+            setattr(self, '%sMap'%name, mapping)
+            setattr(self, '%sIndex'%name, slices.get(name, slice(nP, nP + mapping.nP)))
+            nP += mapping.nP
+        self.nP = nP
+
+    @property
+    def defaultInvProp(self):
+        for name in self._properties:
+            p = self._properties[name]
+            if p.defaultInvProp:
+                return p.name
+
+    def clearMaps(self):
+        for name in self._properties:
+            setattr(self, '%sMap'%name, None)
+            setattr(self, '%sIndex'%name, None)
+
+    def __call__(self, vec):
+        return self.PropModel(self, vec)
+
+    def __contains__(self, val):
+        activeMaps = [name for name in self._properties if getattr(self, '%sMap'%name) is not None]
+        return val in activeMaps

SimPEG/Regularization.py

@@ -27,6 +27,34 @@ class BaseRegularization(object):
         self.mapping = mapping or Maps.IdentityMap(mesh)
         self.mapping._assertMatchesPair(self.mapPair)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     @property
     def parent(self):
         """This is the parent of the regularization."""
@@ -261,22 +289,30 @@ class Tikhonov(BaseRegularization):
         return self._Wzz
 
     @property
-    def W(self):
-        """Full regularization matrix W"""
-        if getattr(self, '_W', None) is None:
-            wlist = (self.Ws, self.Wx, self.Wxx)
+    def Wsmooth(self):
+        """Full smoothness regularization matrix W"""
+        if getattr(self, '_Wsmooth', None) is None:
+            wlist = (self.Wx, self.Wxx)
             if self.mesh.dim > 1:
                 wlist += (self.Wy, self.Wyy)
             if self.mesh.dim > 2:
                 wlist += (self.Wz, self.Wzz)
+            self._Wsmooth = sp.vstack(wlist)
+        return self._Wsmooth
+
+    @property
+    def W(self):
+        """Full regularization matrix W"""
+        if getattr(self, '_W', None) is None:
+            wlist = (self.Ws, self.Wsmooth)
             self._W = sp.vstack(wlist)
         return self._W
 
     @Utils.timeIt
     def eval(self, m):
         if self.smoothModel == True:
-            r1 = self.W * ( self.mapping * (m) )
-            r2 = self.Ws * ( self.mapping * (self.mref) )
+            r1 = self.Wsmooth * ( self.mapping * (m) )
+            r2 = self.Ws * ( self.mapping * (m - self.mref) )
             return 0.5*(r1.dot(r1)+r2.dot(r2))
         elif self.smoothModel == False:
             r = self.W * ( self.mapping * (m - self.mref) )
@@ -302,12 +338,12 @@ class Tikhonov(BaseRegularization):
         """
         if self.smoothModel == True:
             mD1 = self.mapping.deriv(m)
-            mD2 = self.mapping.deriv(self.mref)
-            r1 = self.W * ( self.mapping * (m) )
-            r2 = self.Ws * ( self.mapping * (self.mref) )
-            out1 = mD1.T * ( self.W.T * r1 )
+            mD2 = self.mapping.deriv(m - self.mref)
+            r1 = self.Wsmooth * ( self.mapping * (m))
+            r2 = self.Ws * ( self.mapping * (m - self.mref) )
+            out1 = mD1.T * ( self.Wsmooth.T * r1 )
             out2 = mD2.T * ( self.Ws.T * r2 )
-            out = out1-out2
+            out = out1+out2
         elif self.smoothModel == False:
             mD = self.mapping.deriv(m - self.mref)
             r = self.W * ( self.mapping * (m - self.mref) )

SimPEG/Survey.py

@@ -1,4 +1,4 @@
-import Utils, numpy as np, scipy.sparse as sp
+import Utils, numpy as np, scipy.sparse as sp, uuid
 
 
 class BaseRx(object):
@@ -13,11 +13,41 @@ class BaseRx(object):
     storeProjections = True #: Store calls to getP (organized by mesh)
 
     def __init__(self, locs, rxType, **kwargs):
+        self.uid = str(uuid.uuid4())
         self.locs = locs
         self.rxType = rxType
         self._Ps = {}
         Utils.setKwargs(self, **kwargs)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
+
     @property
     def rxType(self):
         """Receiver Type"""
@@ -124,10 +154,37 @@ class BaseSrc(object):
         for rx in rxList:
             assert isinstance(rx, self.rxPair), 'rxList must be a %s'%self.rxPair.__name__
         assert len(set(rxList)) == len(rxList), 'The rxList must be unique'
-
+        self.uid = str(uuid.uuid4())
         self.rxList = rxList
         Utils.setKwargs(self, **kwargs)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
 
     @property
     def nD(self):
@@ -144,6 +201,7 @@ class Data(object):
     """Fancy data storage by Src and Rx"""
 
     def __init__(self, survey, v=None):
+        self.uid = str(uuid.uuid4())
         self.survey = survey
         self._dataDict = {}
         for src in self.survey.srcList:
@@ -151,6 +209,25 @@ class Data(object):
         if v is not None:
             self.fromvec(v)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        pass
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        pass
+
     def _ensureCorrectKey(self, key):
         if type(key) is tuple:
             if len(key) is not 2:
@@ -208,11 +285,39 @@ class BaseSurvey(object):
     mtrue = None     #: True model, if data is synthetic
 
     counter = None   #: A SimPEG.Utils.Counter object
+    srcPair = BaseSrc  #: Source Pair
 
     def __init__(self, **kwargs):
         Utils.setKwargs(self, **kwargs)
 
-    srcPair = BaseSrc  #: Source Pair
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
 
     @property
     def srcList(self):
@@ -225,6 +330,19 @@ class BaseSurvey(object):
         assert np.all([isinstance(src, self.srcPair) for src in value]), 'All sources must be instances of %s' % self.srcPair.__name__
         assert len(set(value)) == len(value), 'The srcList must be unique'
         self._srcList = value
+        self._sourceOrder = dict()
+        [self._sourceOrder.setdefault(src.uid, ii) for ii, src in enumerate(self._srcList)]
+
+    def getSourceIndex(self, sources):
+        if type(sources) is not list:
+            sources = [sources]
+        for src in sources:
+            if getattr(src,'uid',None) is None:
+                raise KeyError('Source does not have a uid: %s'%str(src))
+        inds = map(lambda src: self._sourceOrder.get(src.uid, None), sources)
+        if None in inds:
+            raise KeyError('Some of the sources specified are not in this survey. %s'%str(inds))
+        return inds
 
     @property
     def prob(self):
@@ -358,3 +476,4 @@ class BaseSurvey(object):
         noise = std*abs(self.dtrue)*np.random.randn(*self.dtrue.shape)
         self.dobs = self.dtrue+noise
         self.std = self.dobs*0 + std
+        return self.dobs

SimPEG/Tests/test_Fields.py

@@ -1,7 +1,8 @@
 import unittest
 from SimPEG import *
 
-class DataAndFieldsTest(unittest.TestCase):
+
+class FieldsTest(unittest.TestCase):
 
     def setUp(self):
         mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
@@ -26,25 +27,6 @@ class DataAndFieldsTest(unittest.TestCase):
         self.mesh = mesh
         self.XYZ = XYZ
 
-    def test_overlappingFields(self):
-        self.assertRaises(AssertionError, Problem.Fields, self.F.mesh, self.F.survey,
-                            knownFields={'b':'F'},
-                            aliasFields={'b':['b',(lambda F, b, ind: b)]})
-
-    def test_data(self):
-        V = []
-        for src in self.D.survey.srcList:
-            for rx in src.rxList:
-                v = np.random.rand(rx.nD)
-                V += [v]
-                self.D[src, rx] = v
-                self.assertTrue(np.all(v == self.D[src, rx]))
-        V = np.concatenate(V)
-        self.assertTrue(np.all(V == Utils.mkvc(self.D)))
-
-        D2 = Survey.Data(self.D.survey, V)
-        self.assertTrue(np.all(Utils.mkvc(D2) == Utils.mkvc(self.D)))
-
     def test_contains(self):
         F = self.F
         nSrc = F.survey.nSrc
@@ -55,10 +37,10 @@ class DataAndFieldsTest(unittest.TestCase):
         self.assertTrue('b' not in F)
         self.assertTrue('e' in F)
 
-    def test_uniqueSrcs(self):
-        srcs = self.D.survey.srcList
-        srcs += [srcs[0]]
-        self.assertRaises(AssertionError, Survey.BaseSurvey, srcList=srcs)
+    def test_overlappingFields(self):
+        self.assertRaises(AssertionError, Problem.Fields, self.F.mesh, self.F.survey,
+                            knownFields={'b':'F'},
+                            aliasFields={'b':['b',(lambda F, b, ind: b)]})
 
     def test_SetGet(self):
         F = self.F
@@ -132,7 +114,6 @@ class FieldsTest_Alias(unittest.TestCase):
         Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3],loc=srcLoc)
         srcList = [Src0,Src1,Src2,Src3,Src4]
         survey = Survey.BaseSurvey(srcList=srcList)
-        self.D = Survey.Data(survey)
         self.F = Problem.Fields(mesh, survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',(lambda e, ind: self.F.mesh.edgeCurl * e)]})
         self.Src0 = Src0
         self.Src1 = Src1
@@ -166,7 +147,7 @@ class FieldsTest_Alias(unittest.TestCase):
 
     def test_aliasFunction(self):
         def alias(e, ind):
-            self.assertTrue(ind is self.Src0)
+            self.assertTrue(ind[0] is self.Src0)
             return self.F.mesh.edgeCurl * e
         F = Problem.Fields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
         e = np.random.rand(F.mesh.nE,1)
@@ -362,7 +343,7 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
         count = [0]
         def alias(e, srcInd, timeInd):
             count[0] += 1
-            self.assertTrue(srcInd is self.Src0)
+            self.assertTrue(srcInd[0] is self.Src0)
             return self.F.mesh.edgeCurl * e
         F = Problem.TimeFields(self.F.mesh, self.F.survey, knownFields={'e':'E'}, aliasFields={'b':['e','F',alias]})
         e = np.random.rand(F.mesh.nE,1,nT)

SimPEG/Tests/test_PropMaps.py

@@ -0,0 +1,193 @@
+import unittest
+from SimPEG import *
+from scipy.constants import mu_0
+
+
+class MyPropMap(Maps.PropMap):
+    sigma = Maps.Property("Electrical Conductivity", defaultInvProp=True)
+    mu    = Maps.Property("Mu", defaultVal=mu_0)
+
+class MyReciprocalPropMap(Maps.PropMap):
+    sigma = Maps.Property("Electrical Conductivity", defaultInvProp=True, propertyLink=('rho',   Maps.ReciprocalMap))
+    rho   = Maps.Property("Electrical Resistivity",                       propertyLink=('sigma', Maps.ReciprocalMap))
+    mu    = Maps.Property("Mu", defaultVal=mu_0, propertyLink=('mui',   Maps.ReciprocalMap))
+    mui    = Maps.Property("Mu", defaultVal=1./mu_0, propertyLink=('mu',   Maps.ReciprocalMap))
+
+
+class TestPropMaps(unittest.TestCase):
+
+    def setUp(self):
+        pass
+
+    def test_setup(self):
+        expMap = Maps.ExpMap(Mesh.TensorMesh((3,)))
+        assert expMap.nP == 3
+
+        PM1 = MyPropMap(expMap)
+        PM2 = MyPropMap([('sigma', expMap)])
+        PM3 = MyPropMap({'maps':[('sigma', expMap)], 'slices':{'sigma':slice(0,3)}})
+
+        for PM in [PM1,PM2,PM3]:
+            assert PM.defaultInvProp == 'sigma'
+            assert PM.sigmaMap is not None
+            assert PM.sigmaMap is expMap
+            assert PM.sigmaIndex == slice(0,3)
+            assert getattr(PM, 'sigma', None) is None
+            assert PM.muMap is None
+            assert PM.muIndex is None
+
+            assert 'sigma' in PM
+            assert 'mu' not in PM
+            assert 'mui' not in PM
+
+            m = PM(np.r_[1.,2,3])
+
+            assert 'sigma' in m
+            assert 'mu' not in m
+            assert 'mui' not in m
+
+            assert m.mu == mu_0
+            assert m.muModel is None
+            assert m.muMap is None
+            assert m.muDeriv is None
+
+            assert np.all(m.sigmaModel == np.r_[1.,2,3])
+            assert m.sigmaMap is expMap
+            assert np.all(m.sigma == np.exp(np.r_[1.,2,3]))
+            assert m.sigmaDeriv is not None
+
+            assert m.nP == 3
+
+    def test_slices(self):
+        expMap = Maps.ExpMap(Mesh.TensorMesh((3,)))
+        PM = MyPropMap({'maps':[('sigma', expMap)], 'slices':{'sigma':[2,1,0]}})
+        assert PM.sigmaIndex == [2,1,0]
+        m = PM(np.r_[1.,2,3])
+        assert np.all(m.sigmaModel == np.r_[3,2,1])
+        assert np.all(m.sigma == np.exp(np.r_[3,2,1]))
+
+    def test_multiMap(self):
+        m = Mesh.TensorMesh((3,))
+        expMap = Maps.ExpMap(m)
+        iMap = Maps.IdentityMap(m)
+        PM = MyPropMap([('sigma', expMap), ('mu', iMap)])
+
+        pm = PM(np.r_[1.,2,3,4,5,6])
+
+        assert pm.nP == 6
+
+        assert 'sigma' in PM
+        assert 'mu' in PM
+        assert 'mui' not in PM
+
+        assert 'sigma' in pm
+        assert 'mu' in pm
+        assert 'mui' not in pm
+
+        assert np.all(pm.sigmaModel == [1.,2,3])
+        assert np.all(pm.sigma == np.exp([1.,2,3]))
+        assert np.all(pm.muModel == [4.,5,6])
+        assert np.all(pm.mu == [4.,5,6])
+
+
+    def test_multiMapCompressed(self):
+        m = Mesh.TensorMesh((3,))
+        expMap = Maps.ExpMap(m)
+        iMap = Maps.IdentityMap(m)
+        PM = MyPropMap({'maps':[('sigma', expMap), ('mu', iMap)],'slices':{'mu':[0,1,2]}})
+
+        pm = PM(np.r_[1,2.,3])
+
+        assert pm.nP == 3
+
+        assert 'sigma' in PM
+        assert 'mu' in PM
+        assert 'mui' not in PM
+
+        assert 'sigma' in pm
+        assert 'mu' in pm
+        assert 'mui' not in pm
+
+        assert np.all(pm.sigmaModel == [1,2,3])
+        assert np.all(pm.sigma == np.exp([1,2,3]))
+        assert np.all(pm.muModel == [1,2,3])
+        assert np.all(pm.mu == [1,2,3])
+
+    def test_Projections(self):
+        m = Mesh.TensorMesh((3,))
+        iMap = Maps.IdentityMap(m)
+        PM = MyReciprocalPropMap([('sigma', iMap)])
+        v = np.r_[1,2.,3]
+        pm = PM(v)
+
+        assert pm.sigmaProj is not None
+        assert pm.rhoProj   is None
+        assert pm.muProj    is None
+        assert pm.muiProj   is None
+
+        assert np.all(pm.sigmaProj * v == pm.sigmaModel)
+
+    def test_Links(self):
+        m = Mesh.TensorMesh((3,))
+        expMap = Maps.ExpMap(m)
+        iMap = Maps.IdentityMap(m)
+        PM = MyReciprocalPropMap([('sigma', iMap)])
+        pm = PM(np.r_[1,2.,3])
+        # print pm.sigma
+        # print pm.sigmaMap
+        assert np.all(pm.sigma == [1,2,3])
+        assert np.all(pm.rho == 1./np.r_[1,2,3])
+        assert pm.sigmaMap is iMap
+        assert pm.rhoMap is None
+        assert pm.sigmaDeriv is not None
+        assert pm.rhoDeriv is not None
+
+        assert 'sigma' in PM
+        assert 'rho' not in PM
+        assert 'mu' not in PM
+        assert 'mui' not in PM
+
+
+        assert 'sigma' in pm
+        assert 'rho' not in pm
+        assert 'mu' not in pm
+        assert 'mui' not in pm
+
+        assert pm.mu == mu_0
+        assert pm.mui == 1.0/mu_0
+        assert pm.muMap is None
+        assert pm.muDeriv is None
+        assert pm.muiMap is None
+        assert pm.muiDeriv is None
+
+        PM = MyReciprocalPropMap([('rho', iMap)])
+        pm = PM(np.r_[1,2.,3])
+        # print pm.sigma
+        # print pm.sigmaMap
+        assert np.all(pm.sigma == 1./np.r_[1,2,3])
+        assert np.all(pm.rho == [1,2,3])
+        assert pm.sigmaMap is None
+        assert pm.rhoMap is iMap
+        assert pm.sigmaDeriv is not None
+        assert pm.rhoDeriv is not None
+
+        assert 'sigma' not in PM
+        assert 'rho' in PM
+        assert 'mu' not in PM
+        assert 'mui' not in PM
+
+
+        assert 'sigma' not in pm
+        assert 'rho' in pm
+        assert 'mu' not in pm
+        assert 'mui' not in pm
+
+        self.assertRaises(AssertionError, MyReciprocalPropMap, [('rho', iMap), ('sigma', iMap)])
+        self.assertRaises(AssertionError, MyReciprocalPropMap, [('sigma', iMap), ('rho', iMap)])
+
+        MyReciprocalPropMap([('sigma', iMap), ('mu', iMap)]) # This should be fine
+
+if __name__ == '__main__':
+    unittest.main()
+
+

SimPEG/Tests/test_SurveyAndData.py

@@ -0,0 +1,53 @@
+import unittest
+from SimPEG import *
+
+class TestData(unittest.TestCase):
+
+    def setUp(self):
+        mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
+        x = np.linspace(5,10,3)
+        XYZ = Utils.ndgrid(x,x,np.r_[0.])
+        srcLoc = np.r_[0,0,0.]
+        rxList0 = Survey.BaseRx(XYZ, 'exi')
+        Src0 = Survey.BaseSrc([rxList0], loc=srcLoc)
+        rxList1 = Survey.BaseRx(XYZ, 'bxi')
+        Src1 = Survey.BaseSrc([rxList1], loc=srcLoc)
+        rxList2 = Survey.BaseRx(XYZ, 'bxi')
+        Src2 = Survey.BaseSrc([rxList2], loc=srcLoc)
+        rxList3 = Survey.BaseRx(XYZ, 'bxi')
+        Src3 = Survey.BaseSrc([rxList3], loc=srcLoc)
+        Src4 = Survey.BaseSrc([rxList0, rxList1, rxList2, rxList3], loc=srcLoc)
+        srcList = [Src0,Src1,Src2,Src3,Src4]
+        survey = Survey.BaseSurvey(srcList=srcList)
+        self.D = Survey.Data(survey)
+
+    def test_data(self):
+        V = []
+        for src in self.D.survey.srcList:
+            for rx in src.rxList:
+                v = np.random.rand(rx.nD)
+                V += [v]
+                self.D[src, rx] = v
+                self.assertTrue(np.all(v == self.D[src, rx]))
+        V = np.concatenate(V)
+        self.assertTrue(np.all(V == Utils.mkvc(self.D)))
+
+        D2 = Survey.Data(self.D.survey, V)
+        self.assertTrue(np.all(Utils.mkvc(D2) == Utils.mkvc(self.D)))
+
+    def test_uniqueSrcs(self):
+        srcs = self.D.survey.srcList
+        srcs += [srcs[0]]
+        self.assertRaises(AssertionError, Survey.BaseSurvey, srcList=srcs)
+
+    def test_sourceIndex(self):
+        survey = self.D.survey
+        srcs = survey.srcList
+        assert survey.getSourceIndex([srcs[1],srcs[0]]) == [1,0]
+        assert survey.getSourceIndex([srcs[1],srcs[2],srcs[2]]) == [1,2,2]
+        SrcNotThere = Survey.BaseSrc(srcs[0].rxList, loc=np.r_[0,0,0])
+        self.assertRaises(KeyError, survey.getSourceIndex, [SrcNotThere])
+        self.assertRaises(KeyError, survey.getSourceIndex, [srcs[1],srcs[2],SrcNotThere])
+
+if __name__ == '__main__':
+    unittest.main()

SimPEG/Tests/test_maps.py

@@ -6,8 +6,8 @@ from scipy.sparse.linalg import dsolve
 
 TOL = 1e-14
 
-MAPS_TO_TEST_2D = ["CircleMap", "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting"]
-MAPS_TO_TEST_3D = [             "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting"]
+MAPS_TO_TEST_2D = ["CircleMap", "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting", "FullMap"]
+MAPS_TO_TEST_3D = [             "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting", "FullMap"]
 
 class MapTests(unittest.TestCase):
 
@@ -30,8 +30,8 @@ class MapTests(unittest.TestCase):
             self.assertTrue(maps.test())
 
 
-    def test_transforms_logMap(self):
-        # Note that log maps can be kinda finicky, so we are being explicit about the random seed.
+    def test_transforms_logMap_reciprocalMap(self):
+        # Note that log/reciprocal maps can be kinda finicky, so we are being explicit about the random seed.
         v2 = np.r_[ 0.40077291, 0.14410044, 0.58452314, 0.96323738, 0.01198519, 0.79754415]
         dv2 = np.r_[ 0.80653921, 0.13132446, 0.4901117, 0.03358737, 0.65473762, 0.44252488]
         v3 = np.r_[ 0.96084865, 0.34385186, 0.39430044, 0.81671285, 0.65929109, 0.2235217, 0.87897526, 0.5784033, 0.96876393, 0.63535864, 0.84130763, 0.22123854]
@@ -41,6 +41,11 @@ class MapTests(unittest.TestCase):
         maps = Maps.LogMap(self.mesh3)
         self.assertTrue(maps.test(v3, dx=dv3))
 
+        maps = Maps.ReciprocalMap(self.mesh2)
+        self.assertTrue(maps.test(v2, dx=dv2))
+        maps = Maps.ReciprocalMap(self.mesh3)
+        self.assertTrue(maps.test(v3, dx=dv3))
+
     def test_Mesh2MeshMap(self):
         maps = Maps.Mesh2Mesh([self.mesh22, self.mesh2])
         self.assertTrue(maps.test())

SimPEG/Tests/test_regularization.py

@@ -22,7 +22,7 @@ class RegularizationTests(unittest.TestCase):
             mapping = r.mapPair(self.mesh2)
             reg = r(self.mesh2, mapping=mapping)
             m = np.random.rand(mapping.nP)
-            reg.mref = m[:]*0
+            reg.mref = m[:]*np.mean(m)
             passed = checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False)
             self.assertTrue(passed)
 

SimPEG/Utils/ModelBuilder.py

@@ -3,10 +3,27 @@ import scipy.ndimage as ndi
 import scipy.sparse as sp
 from matutils import mkvc
 
-
-def getIndecesBlock(p0,p1,ccMesh):
+def addBlock(gridCC, modelCC, p0, p1, blockProp):
     """
-        Creates a vector containing the block indexes in the cell centerd mesh.
+        Add a block to an exsisting cell centered model, modelCC
+
+        :param numpy.array, gridCC: mesh.gridCC is the cell centered grid
+        :param numpy.array, modelCC: cell centered model
+        :param numpy.array, p0: bottom, southwest corner of block
+        :param numpy.array, p1: top, northeast corner of block
+        :blockProp float, blockProp: property to assign to the model 
+
+        :return numpy.array, modelBlock: model with block
+    """
+    ind = getIndicesBlock(p0, p1, gridCC)
+    modelBlock = modelCC.copy()
+    modelBlock[ind] = blockProp
+    return modelBlock
+
+
+def getIndicesBlock(p0,p1,ccMesh):
+    """
+        Creates a vector containing the block indices in the cell centers mesh.
         Returns a tuple
 
         The block is defined by the points
@@ -78,7 +95,7 @@ def defineBlock(ccMesh,p0,p1,vals=[0,1]):
         vals[1]  conductivity of the ground
     """
     sigma = np.zeros(ccMesh.shape[0]) + vals[1]
-    ind   = getIndecesBlock(p0,p1,ccMesh)
+    ind   = getIndicesBlock(p0,p1,ccMesh)
 
     sigma[ind] = vals[0]
 
@@ -132,7 +149,7 @@ def defineTwoLayers(ccMesh,depth,vals=[0,1]):
     # The depth is always defined on the last one.
     p1[len(p1)-1] -= depth
 
-    ind   = getIndecesBlock(p0,p1,ccMesh)
+    ind   = getIndicesBlock(p0,p1,ccMesh)
 
     sigma[ind] = vals[0];
 
@@ -153,16 +170,58 @@ def scalarConductivity(ccMesh,pFunction):
 
     return mkvc(sigma)
 
+def layeredModel(ccMesh, layerTops, layerValues):
+    """
+        Define a layered model from layerTops (z-positive up)
+
+        :param numpy.array ccMesh: cell-centered mesh
+        :param numpy.array layerTops: z-locations of the tops of each layer
+        :param numpy.array layerValue: values of the property to assign for each layer (starting at the top) 
+        :rtype: numpy.array
+        :return: M, layered model on the mesh 
+    """
+
+    descending = np.linalg.norm(sorted(layerTops, reverse=True) - layerTops) < 1e-20
+
+    # TODO: put an error check to make sure that there is an ordering... needs to work with inf elts 
+    # assert ascending or descending, "Layers must be listed in either ascending or descending order"
+
+    # start from bottom up
+    if not descending:
+        zprop = np.hstack([mkvc(layerTops,2),mkvc(layerValues,2)])
+        zprop.sort(axis=0)
+        layerTops, layerValues  = zprop[::-1,0], zprop[::-1,1]
+
+    # put in vector form
+    layerTops, layerValues = mkvc(layerTops), mkvc(layerValues)
+
+    # initialize with bottom layer
+    dim = ccMesh.shape[1]
+    if dim == 3:
+        z = ccMesh[:,2]
+    elif dim == 2:
+        z = ccMesh[:,1]
+    elif dim == 1:
+        z = ccMesh[:,0]
+
+    model = np.zeros(ccMesh.shape[0])
+
+    for i, top in enumerate(layerTops):
+        zind = z <= top 
+        model[zind] = layerValues[i]
+
+    return model 
+
 
 
 def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]):
     """
-        Create a random model by convolving a kernal with a
+        Create a random model by convolving a kernel with a
         uniformly distributed model.
 
         :param int,tuple shape: shape of the model.
         :param int seed: pick which model to produce, prints the seed if you don't choose.
-        :param numpy.ndarray,list anisotropy: this is the (3 x n) blurring kernal that is used.
+        :param numpy.ndarray,list anisotropy: this is the (3 x n) blurring kernel that is used.
         :param int its: number of smoothing iterations
         :param list bounds: bounds on the model, len(list) == 2
         :rtype: numpy.ndarray

SimPEG/Utils/codeutils.py

@@ -58,9 +58,11 @@ def hook(obj, method, name=None, overwrite=False, silent=False):
         print 'Method '+name+' was not overwritten.'
 
 
-def setKwargs(obj, **kwargs):
+def setKwargs(obj, ignore=[],  **kwargs):
     """Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist."""
     for attr in kwargs:
+        if attr in ignore:
+            continue
         if hasattr(obj, attr):
             setattr(obj, attr, kwargs[attr])
         else:

SimPEG/Utils/meshutils.py

@@ -149,7 +149,7 @@ def readUBCTensorModel(fileName, mesh):
 
         Input:
         :param fileName, path to the UBC GIF mesh file to read
-        :param mesh, TensorMesh object, mesh that coresponds to the model 
+        :param mesh, TensorMesh object, mesh that coresponds to the model
 
         Output:
         :return numpy array, model with TensorMesh ordered
@@ -170,7 +170,7 @@ def writeUBCTensorMesh(fileName, mesh):
 
         :param str fileName: File to write to
         :param simpeg.Mesh.TensorMesh mesh: The mesh
-        
+
     """
     assert mesh.dim == 3
     s = ''
@@ -205,7 +205,6 @@ def writeUBCTensorModel(fileName, mesh, model):
 
     np.savetxt(fileName, modelMatTR.ravel())
 
-
 def readVTRFile(fileName):
     """
         Read VTK Rectilinear (vtr xml file) and return SimPEG Tensor mesh and model
@@ -216,7 +215,7 @@ def readVTRFile(fileName):
         Output:
         :return SimPEG TensorMesh object
         :return SimPEG model dictionary
-        
+
     """
     # Import
     from vtk import vtkXMLRectilinearGridReader as vtrFileReader
@@ -296,84 +295,87 @@ def writeVTRFile(fileName,mesh,model=None):
     vtkObj.SetZCoordinates(numpy_to_vtk(vZ,deep=1))
 
     # Assign the model('s) to the object
-    for item in model.iteritems():
-        # Convert numpy array
-        vtkDoubleArr = numpy_to_vtk(item[1],deep=1)
-        vtkDoubleArr.SetName(item[0])
-        vtkObj.GetCellData().AddArray(vtkDoubleArr)
-    # Set the active scalar
-    vtkObj.GetCellData().SetActiveScalars(model.keys()[0])
-    vtkObj.Update()
+    if model is not None:
+        for item in model.iteritems():
+            # Convert numpy array
+            vtkDoubleArr = numpy_to_vtk(item[1],deep=1)
+            vtkDoubleArr.SetName(item[0])
+            vtkObj.GetCellData().AddArray(vtkDoubleArr)
+        # Set the active scalar
+        vtkObj.GetCellData().SetActiveScalars(model.keys()[0])
 
 
     # Check the extension of the fileName
-    ext = os.path.splitext(fileName)[1]
-    if ext is '':
-        fileName = fileName + '.vtr'
-    elif ext not in '.vtr':
-        raise IOError('{:s} is an incorrect extension, has to be .vtr')
-    # Write the file.
-    vtrWriteFilter = rectWriter()
-    vtrWriteFilter.SetInput(vtkObj)
-    vtrWriteFilter.SetFileName(fileName)
-    vtrWriteFilter.Update()
+    if fileName is not None:
+        ext = os.path.splitext(fileName)[1]
+        if ext is '':
+            fileName = fileName + '.vtr'
+        elif ext not in '.vtr':
+            raise IOError('{:s} is an incorrect extension, has to be .vtr')
+        # Write the file.
 
+        vtrWriteFilter = rectWriter()
+        vtrWriteFilter.SetInputData(vtkObj)
+        vtrWriteFilter.SetFileName(fileName)
+        vtrWriteFilter.Update()
+    else:
+        return vtkObj
 
 def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
     """
         Extracts Core Mesh from Global mesh
         xyzlim: 2D array [ndim x 2]
         mesh: SimPEG mesh
-        This function ouputs: 
+        This function ouputs:
             - actind: corresponding boolean index from global to core
-            - meshcore: core SimPEG mesh  
+            - meshcore: core SimPEG mesh
         Warning: 1D and 2D has not been tested
     """
     from SimPEG import Mesh
     if mesh.dim ==1:
         xyzlim = xyzlim.flatten()
         xmin, xmax = xyzlim[0], xyzlim[1]
-        
-        xind = np.logical_and(mesh.vectorCCx>xmin, mesh.vectorCCx<xmax)        
-        
+
+        xind = np.logical_and(mesh.vectorCCx>xmin, mesh.vectorCCx<xmax)
+
         xc = mesh.vectorCCx[xind]
 
         hx = mesh.hx[xind]
-        
+
         x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5]
-        
+
         meshCore = Mesh.TensorMesh([hx, hy] ,x0=x0)
-        
+
         actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax)
-        
+
     elif mesh.dim ==2:
         xmin, xmax = xyzlim[0,0], xyzlim[0,1]
         ymin, ymax = xyzlim[1,0], xyzlim[1,1]
 
         yind = np.logical_and(mesh.vectorCCy>ymin, mesh.vectorCCy<ymax)
-        zind = np.logical_and(mesh.vectorCCz>zmin, mesh.vectorCCz<zmax)        
+        zind = np.logical_and(mesh.vectorCCz>zmin, mesh.vectorCCz<zmax)
 
         xc = mesh.vectorCCx[xind]
         yc = mesh.vectorCCy[yind]
 
         hx = mesh.hx[xind]
         hy = mesh.hy[yind]
-        
+
         x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5]
-        
+
         meshCore = Mesh.TensorMesh([hx, hy] ,x0=x0)
-        
+
         actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax) \
                & (mesh.gridCC[:,1]>ymin) & (mesh.gridCC[:,1]<ymax) \
-        
+
     elif mesh.dim==3:
         xmin, xmax = xyzlim[0,0], xyzlim[0,1]
         ymin, ymax = xyzlim[1,0], xyzlim[1,1]
         zmin, zmax = xyzlim[2,0], xyzlim[2,1]
-        
+
         xind = np.logical_and(mesh.vectorCCx>xmin, mesh.vectorCCx<xmax)
         yind = np.logical_and(mesh.vectorCCy>ymin, mesh.vectorCCy<ymax)
-        zind = np.logical_and(mesh.vectorCCz>zmin, mesh.vectorCCz<zmax)        
+        zind = np.logical_and(mesh.vectorCCz>zmin, mesh.vectorCCz<zmax)
 
         xc = mesh.vectorCCx[xind]
         yc = mesh.vectorCCy[yind]
@@ -382,19 +384,19 @@ def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'):
         hx = mesh.hx[xind]
         hy = mesh.hy[yind]
         hz = mesh.hz[zind]
-        
+
         x0 = [xc[0]-hx[0]*0.5, yc[0]-hy[0]*0.5, zc[0]-hz[0]*0.5]
-        
+
         meshCore = Mesh.TensorMesh([hx, hy, hz] ,x0=x0)
-        
+
         actind = (mesh.gridCC[:,0]>xmin) & (mesh.gridCC[:,0]<xmax) \
                & (mesh.gridCC[:,1]>ymin) & (mesh.gridCC[:,1]<ymax) \
                & (mesh.gridCC[:,2]>zmin) & (mesh.gridCC[:,2]<zmax)
-                
+
     else:
         raise(Exception("Not implemented!"))
-    
-    
+
+
     return actind, meshCore