Merge branch 'master' into parallel

broke apart smoothness and smallness terms

.bumpversion.cfg

@@ -1,4 +1,4 @@
 [bumpversion]
-current_version = 0.1.2
+current_version = 0.1.3
 files = setup.py SimPEG/__init__.py docs/conf.py
 

.travis.yml

@@ -15,7 +15,7 @@ before_install:
 
 # Install packages
 install:
-  - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython
+  - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython networkx pyzmq
   - pip install nose-cov python-coveralls
   # - pip install -r requirements.txt
   - python setup.py install

SimPEG/Directives.py

@@ -144,23 +144,69 @@ 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 UpdateReferenceModel(Parameter):

SimPEG/Dispatcher.py

@@ -0,0 +1,112 @@
+### PROTOTYPE INTERFACE FOR PARALLEL DISPATCHER ###
+
+from functools import wraps
+
+def synchronize(fn):
+    @wraps(fn)
+    def wrapper(*args, **kwargs):
+        self = args[0]
+        
+        pr = isinstance(getattr(self, '_dispatcher', None), ParallelDispatcher)
+        if pr:
+            print('Parallel stuff: (start) %(prob)s.%(fn)s'%{'prob': self.__class__.__name__, 'fn': fn.__name__})
+
+        result =  fn(*args, **kwargs)
+        
+        if pr:
+            print('Parallel stuff: ( end ) %(prob)s.%(fn)s'%{'prob': self.__class__.__name__, 'fn': fn.__name__})
+        
+        return result
+
+    return wrapper
+
+
+class BaseDispatcher(object):
+    
+    def __init__(self, *args, **kwargs):
+        print('INIT: Dispatcher!')
+
+    def pair(self, problem):
+        self._prob = problem
+        print('PAIR: Dispatcher setup...')
+
+class SerialDispatcher(BaseDispatcher):
+    
+    def __init__(self, *args, **kwargs):
+        BaseDispatcher.__init__(self, *args, **kwargs)
+        print('INIT: Serial dispatcher...')
+
+class ParallelDispatcher(BaseDispatcher):
+
+    remoteOnly = ['someotherattribute']
+    
+    def __init__(self, *args, **kwargs):
+        BaseDispatcher.__init__(self, *args, **kwargs)
+        print('INIT: Parallel dispatcher...')
+
+    def pair(self, problem):
+        BaseDispatcher.pair(self, problem)
+        print('PAIR: Parallel dispatcher setup...')
+
+    def interceptSetattr(self, prob, name, value):
+        print('SET: Parallel dispatcher set %(prob)s.%(name)s = %(value)r'
+            %{'prob': prob.__class__.__name__, 'name': name, 'value': value})
+
+        if name in self.remoteOnly:
+            print('Setting remote state...')
+        else:
+            raise AttributeError('Set local copy!')
+
+    def interceptGetattr(self, prob, name):
+        print('GET: Parallel dispatcher get %(prob)s.%(name)s'
+            %{'prob': prob.__class__.__name__, 'name': name})
+
+        if name in self.remoteOnly:
+            return '***Value from remote state***'
+        else:
+            raise AttributeError('Attribute %s not in parallel namespace!'%(name,))
+
+class StandinSurvey(object):
+    
+    def pair(self, problem):
+        self._prob = problem
+        
+class StandinProblem(object):
+    
+    def __init__(self):
+        print('INIT: Problem!')
+        self._dispatcher = SerialDispatcher()
+
+    def __setattr__(self, name, value):
+        d = getattr(self, '_dispatcher', None)
+
+        if isinstance(d, ParallelDispatcher):
+            try:
+                d.interceptSetattr(self, name, value)
+            except AttributeError:
+                super(self.__class__, self).__setattr__(name, value)
+            finally:
+                return
+
+        else:
+            super(self.__class__, self).__setattr__(name, value)
+
+    def __getattr__(self, name):
+        d = super(self.__class__, self).__getattribute__('_dispatcher')
+
+        if isinstance(d, ParallelDispatcher):
+            return d.interceptGetattr(self, name)
+    
+    def pair(self, survey, dispatcher=None):
+        
+        self._survey = survey
+        self._survey.pair(self)
+        
+        if dispatcher is not None:
+            self._dispatcher = dispatcher
+            print('PAIR: Problem setup...')
+            self._dispatcher.pair(self)
+    
+    @synchronize
+    def dosomething(self):
+        print('Doing something!')

SimPEG/Fields.py

@@ -0,0 +1,261 @@
+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.
+
+    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):
@@ -22,6 +23,8 @@ class IdentityMap(object):
             :rtype: int
             :return: number of parameters in the model
         """
+        if self.mesh is None:
+            return '*'
         return self.mesh.nC
 
     @property
@@ -32,8 +35,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 +104,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 +125,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 +161,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 +226,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 +285,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/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
 
@@ -186,6 +188,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 +213,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/Parallel.py

@@ -0,0 +1,638 @@
+from ipyparallel import Client, parallel, Reference, require, depend, interactive
+from SimPEG.Utils import CommonReducer
+import numpy as np
+import networkx
+
+DEFAULT_MPI = True
+MPI_BELLWETHERS = ['PMI_SIZE', 'OMPI_UNIVERSE_SIZE']
+
+class SuperReference(object):
+    '''
+    Object that can be called to return a reference, but
+    will only be schedulable on the correct worker(s) if
+    its 'lrank' parameter has been set.
+    '''
+
+    def __init__(self, ref, lrank=None):
+
+        if (lrank is None) or (type(lrank) is list):
+            self.rank = lrank
+        else:
+            self.rank = [lrank]
+            
+        self.ref = ref
+
+    def __call__(self, *args, **kwargs):
+
+        from ipyparallel import depend
+        from ipyparallel.error import UnmetDependency
+
+        if (self.rank is not None) and (globals().get('rank', None) not in self.rank):
+            raise UnmetDependency('Global \'rank\' does not satisfy requirements')
+
+        return self.ref(*args, **kwargs)
+
+class Endpoint(object):
+    '''
+    Object that holds the namespace of the SimPEG parallel
+    footprint on the remote workers.
+    '''
+
+    problemFactory = lambda: None   # Callable for constructing system / problem
+    surveyFactory = lambda: None    # Callable for constructing survey
+    localFields = {}                # Dictionary for storing local fields
+    globalFields = {}               # Dictionary for storing merged fields
+    localProblems = {}              # Dictionary of local subsystem / problem objects
+    localSurveys = {}               # Dictionary of local survey objects
+    functions = {}                  # Dictionary of callables to carry out modelling / etc.
+    fieldspec = None                # Dictionary of callables to setup field storage objects
+    baseSystemConfig = {}           # Base configuration for system
+
+    def setupLocalFields(self, whichfields=None):
+
+        # If no names are specified, clear all fields first
+        if whichfields is None:
+            self.localFields = {}
+
+        # If we have a 'fieldspec' object...
+        if getattr(self, 'fieldspec', None) is not None:
+            # ...either loop over the specified names, or all the fields...
+            for fn in (whichfields or self.fieldspec):
+                # ...and construct a new empty object per the 'fieldspec' constructor.
+                self.localFields[fn] = self.fieldspec[fn]()
+
+    def setupLocalSurveys(self, subConfigs):
+
+        # Loop over possible survey configurations (may differ in source terms, etc.)
+        for isub in subConfigs:
+            # For each 'isub' create a separate copy of the base configuration...
+            geom = self.baseSystemConfig['geom'].copy()
+            # ...and update with any differences...
+            geom.update(subConfigs[isub])
+            # ...then construct the Survey object and store it for later pairing.
+            self.localSurveys[isub] = self.surveyFactory(geom)
+
+    def setupLocalProblem(self, subConfig):
+
+        # Make a copy w/o the geometry information (which is used by Survey)
+        systemConfig = {key: self.baseSystemConfig[key] for key in self.baseSystemConfig if key not in ['geom']}
+
+        # Update with the configuration for this subproblem
+        systemConfig.update(subConfig)
+
+        # Create the local subproblem...
+        problem = self.problemFactory(systemConfig)
+        # ...and pair it with a corresponding survey for this 'isub' (e.g., frequency)...
+        problem.pair(self.localSurveys[subConfig['isub']])
+        # ...then store in the Endpoint for later access by the scheduler.
+        self.localProblems[subConfig['tag']] = problem
+
+
+class SystemGraph(networkx.DiGraph):
+    '''
+    NetworkX Directed Graph subclass that knows about
+    job status information, and can return a representation
+    of itself for use in interactive debugging/testing.
+    '''
+
+    @staticmethod
+    def _codeStatus(data):
+        
+        status = 0
+        
+        if 'jobs' in data:
+            status = 1 * data['jobs'][-1].ready() + 1
+            if status > 1:
+                status += 1 * (not data['jobs'][-1].successful())
+        
+        return status 
+
+    def _codeGraph(self):
+        from networkx.readwrite import json_graph
+        
+        G = networkx.DiGraph()
+        
+        for e in self.edges_iter():
+            G.add_edge(e[0], e[1])
+        
+        for n, data in self.nodes_iter(data=True):
+            G.add_node(n, status=self._codeStatus(data))
+        
+        return json_graph.node_link_data(G)
+
+    def RenderHTML(self):
+        import pkg_resources
+        from IPython.core import display
+        import time
+
+        data = str(self._codeGraph())
+        uniqueID = hash(time.time())
+
+        formatstr = {
+            'uniqueID': 'Graph%s'%uniqueID,
+            'JSONData': data,
+        }
+
+        code = pkg_resources.resource_string('SimPEG', 'Resources/Parallel/SystemGraph.html')%formatstr
+
+        return display.HTML(data=code)._repr_html_()
+
+try:
+    get_ipython().display_formatter.formatters['text/html'].for_type(SystemGraph, SystemGraph.RenderHTML)
+except NameError:
+    pass
+
+class SystemSolver(object):
+
+    def __init__(self, problem, schedule):
+
+        self.problem = problem
+        self.remote = problem.remote
+        self.schedule = schedule
+
+    def __call__(self, entry, isrcs):
+        
+        # TODO: Replace with SuperReference instances
+        fnformat = '%s.functions["%s"]'
+        fnRef = Reference(fnformat%(self.remote.endpointName, self.schedule[entry]['solve']))
+        clearRef = Reference(fnformat%(self.remote.endpointName, self.schedule[entry]['clear']))
+        reduceLabels = self.schedule[entry]['reduce']
+
+        dview = self.problem.remote.dview
+        lview = self.problem.remote.lview
+
+        chunksPerWorker = getattr(self.problem, 'chunksPerWorker', 1)
+
+        G = SystemGraph()
+
+        mainNode = 'Beginning'
+        G.add_node(mainNode)
+
+        # Parse sources
+        # TODO: Get from Survey somehow?
+        nsrc = self.problem.nsrc
+        if isrcs is None:
+            isrcs = slice(None)
+
+        elif not isinstance(isrcs, slice):
+            raise Exception('Scheduler must run over slice or None!')
+
+        # TODO: Replace w/ hook into Endpoint classes
+        systemsOnWorkers = dview['%s.localProblems.keys()'%self.remote.endpointName]
+        ids = dview['rank']
+        tags = set()
+        for ltags in systemsOnWorkers:
+            tags = tags.union(set(ltags))
+
+        clearJobs = []
+        endNodes = {}
+        tailNodes = []
+
+        for tag in tags:
+
+            tagNode = 'Head: %d, %d'%tag
+            G.add_edge(mainNode, tagNode)
+
+            relIDs = []
+            for i in xrange(len(ids)):
+
+                systems = systemsOnWorkers[i]
+                rank = ids[i]
+
+                if tag in systems:
+                    relIDs.append(rank)
+
+            systemJobs = []
+            endNodes[tag] = []
+            systemNodes = []
+
+            with lview.temp_flags(block=False):
+                iworks = 0
+                for work in self._subSlice(isrcs, int(round(chunksPerWorker*len(relIDs)))):
+                    if work:
+                        job = lview.apply(fnRef, Reference(self.remote.endpointName), tag, work)
+                        systemJobs.append(job)
+                        label = 'Compute: %d, %d, %d'%(tag[0], tag[1], iworks)
+                        systemNodes.append(label)
+                        G.add_node(label, jobs=[job], subslice=work, tag=tag)
+                        G.add_edge(tagNode, label)
+                        iworks += 1
+
+            if getattr(self.problem, 'ensembleClear', False): # True for ensemble ending, False for individual ending
+                tagNode = 'Wrap: %d, %d'%tag
+                for label in systemNodes:
+                    G.add_edge(label, tagNode)
+
+                for rank in relIDs:
+
+                    with lview.temp_flags(block=False, after=systemJobs):
+                        # TODO: Remove dependency on self._hasSystemRank, once the SuperReferences
+                        #       are able to be used. They will automatically schedule only on the
+                        #       correct (allowed) systems.
+                        job = lview.apply(clearRef, Reference(self.remote.endpointName), tag, rank)
+                        clearJobs.append(job)
+                        label = 'Wrap: %d, %d, %d'%(tag[0],tag[1], rank)
+                        G.add_node(label, jobs=[job], tag=tag, rank=rank)
+                        endNodes[tag].append(label)
+                        G.add_edge(tagNode, label)
+            else:
+
+                for i, sjob in enumerate(systemJobs):
+                    with lview.temp_flags(block=False, follow=sjob, after=sjob):
+                        job = lview.apply(clearRef, Reference(self.remote.endpointName), tag)
+                        clearJobs.append(job)
+                        label = 'Wrap: %d, %d, %d'%(tag[0],tag[1],i)
+                        G.add_node(label, jobs=[job])
+                        endNodes[tag].append(label)
+                        G.add_edge(systemNodes[i], label)
+
+            tagNode = 'Tail: %d, %d'%tag
+            for label in endNodes[tag]:
+                G.add_edge(label, tagNode)
+            tailNodes.append(tagNode)
+
+        endNode = 'End'
+        jobs = []
+        after = clearJobs
+        for label in reduceLabels:
+            job = self.problem.remote.reduceLB(Reference(self.remote.endpointName), label, after)
+            after = job
+            if job is not None:
+                jobs.append(job)
+        G.add_node(endNode, jobs=jobs)
+        for node in tailNodes:
+            G.add_edge(node, endNode)
+
+        return G
+
+    def wait(self, G):
+        self.problem.remote.lview.wait(G.node['End']['jobs'] if G.node['End']['jobs'] else (G.node[wn]['jobs'] for wn in (G.predecessors(tn)[0] for tn in G.predecessors('End'))))
+
+    # TODO: Hopefully obsoleted by SuperReference
+    @staticmethod
+    @interactive
+    def _hasSystemRank(endpoint, tag, wid):
+        global rank
+
+        return (tag in endpoint.localProblems) and (rank == wid)
+
+    @staticmethod
+    def _getChunks(problems, chunks=1):
+        nproblems = len(problems)
+        return (problems[i*nproblems // chunks: (i+1)*nproblems // chunks] for i in range(chunks))
+
+    @staticmethod
+    def _subSlice(insl,  chunks=1):
+        start = insl.start or 0
+        nproblems = insl.stop - start
+        return [slice(start + i*nproblems/chunks, start + (i+1)*nproblems/chunks) for i in xrange(chunks)]
+
+
+class RemoteInterface(object):
+
+    def __init__(self, profile=None, MPI=None, nThreads=1, bootstrap=None, endpointName='endpoint'):
+
+        # TODO: Add interface for namespace bootstrapping from
+        #       the dispatcher / problem side
+
+        if profile is not None:
+            pupdate = {'profile': profile}
+        else:
+            pupdate = {}
+
+        pclient = Client(**pupdate)
+
+        if not self._cdSame(pclient):
+            print('Could not change all workers to the same directory as the client!')
+
+        dview = pclient[:]
+        dview.block = True
+        dview.clear()
+
+        remoteSetup = '''
+        import os'''
+
+        parMPISetup = ''' 
+        from mpi4py import MPI
+        comm = MPI.COMM_WORLD
+        rank = comm.Get_rank()''' 
+
+        for command in remoteSetup.strip().split('\n'):
+            dview.execute(command.strip())
+
+        dview.scatter('rank', pclient.ids, flatten=True)
+
+        self.e0 = pclient[0]
+        self.e0.block = True
+
+        self.useMPI = False
+        MPI = DEFAULT_MPI if MPI is None else MPI
+        if MPI:
+            MPISafe = False
+
+            for var in MPI_BELLWETHERS:
+                MPISafe = MPISafe or all(dview['os.getenv("%s")'%(var,)])
+
+            if MPISafe:
+                for command in parMPISetup.strip().split('\n'):
+                    dview.execute(command.strip())
+                ranks = dview['rank']
+                reorder = [ranks.index(i) for i in xrange(len(ranks))]
+                dview = pclient[reorder]
+                dview.block = True
+                dview.activate()
+
+                # Set up necessary parts for broadcast-based communication
+                self.e0 = pclient[reorder[0]]
+                self.e0.block = True
+                self.comm = Reference('comm')
+
+            self.useMPI = MPISafe
+
+        self.pclient = pclient
+        self.dview = dview
+        self.lview = pclient.load_balanced_view()
+
+        self.nThreads = nThreads
+
+        if bootstrap is not None:
+            for command in bootstrap.strip().split('\n'):
+                dview.execute(command.strip())
+
+        self.endpointName = endpointName
+
+    @property
+    def nThreads(self):
+        return self._nThreads
+    @nThreads.setter
+    def nThreads(self, value):
+        self._nThreads = value
+        self.dview.apply(self._adjustMKLVectorization, self._nThreads)
+    
+
+    def __setitem__(self, key, item):
+
+        if self.useMPI:
+            self.e0[key] = item
+            code = 'if rank != 0: %(key)s = None\n%(key)s = comm.bcast(%(key)s, root=0)'
+            self.dview.execute(code%{'key': key})
+
+        else:
+            self.dview[key] = item
+
+    def __getitem__(self, key):
+
+        if self.useMPI:
+            code = 'temp_%(key)s = None\ntemp_%(key)s = comm.gather(%(key)s, root=%(root)d)'
+            self.dview.execute(code%{'key': key, 'root': 0})
+            item = self.e0['temp_%s'%(key,)]
+            self.e0.execute('del temp_%s'%(key,))
+
+        else:
+            item = self.dview[key]
+
+        return item
+
+    def reduceLB(self, endpoint, key, after=None):
+
+        repeat = lambda value: (value for i in xrange(len(self.pclient.ids)))
+
+        if self.useMPI:
+            with self.lview.temp_flags(block=False, after=after):
+                job = self.lview.map(self._reduceJob, xrange(len(self.pclient.ids)), repeat(0), repeat(endpoint), repeat(key))
+
+            return job
+
+    def reduce(self, key, axis=None):
+
+        if self.useMPI:
+            code = 'temp_%(key)s = comm.reduce(%(key)s, root=%(root)d)'
+            self.dview.execute(code%{'key': key, 'root': 0})
+
+            # if axis is not None:
+            #     code = 'temp_%(key)s = temp_%(key)s.sum(axis=%(axis)d)'
+            #     self.e0.execute(code%{'key': key, 'axis': axis})
+
+            item = self.e0['temp_%s'%(key,)]
+            self.dview.execute('del temp_%s'%(key,))
+
+        else:
+            item = reduce(np.add, self.dview[key])
+
+        return item
+
+    def reduceMul(self, key1, key2, axis=None):
+
+        if self.useMPI:
+            # Gather
+            code_reduce = 'temp_%(key)s = comm.reduce(%(key)s, root=%(root)d)'
+            self.dview.execute(code_reduce%{'key': key1, 'root': 0})
+            self.dview.execute(code_reduce%{'key': key2, 'root': 0})
+
+            # Multiply
+            code_mul = 'temp_%(key1)s%(key2)s = temp_%(key1)s * temp_%(key2)s'
+            self.e0.execute(code_mul%{'key1': key1, 'key2': key2})
+
+            # Potentially sum
+            if axis is not None:
+                code = 'temp_%(key1)s%(key2)s = temp_%(key1)s%(key2)s.sum(axis=%(axis)d)'
+                self.e0.execute(code%{'key1': key1, 'key2': key2, 'axis': axis})
+
+            # Pull
+            item = self.e0['temp_%(key1)s%(key2)s'%{'key1': key1, 'key2': key2}]
+
+            # Clear
+            self.dview.execute('del temp_%s'%(key1,))
+            self.dview.execute('del temp_%s'%(key2,))
+            self.e0.execute('del temp_%(key1)s%(key2)s'%{'key1': key1, 'key2': key2})
+
+        else:
+            item1 = reduce(np.add, self.dview[key1])
+            item2 = reduce(np.add, self.dview[key2])
+            item = item1 * item2
+
+        return item
+
+    def remoteMulE0(self, key1, key2, axis=None):
+
+        code_mul = 'temp_field = %(key1)s * %(key2)s'
+        self.e0.execute(code_mul%{'key1': key1, 'key2': key2})
+
+        if axis is not None:
+            code = 'temp_field = temp_field.sum(axis=%(axis)d)'
+            self.e0.execute(code%{'axis': axis})
+
+        item = self.e0['temp_field']
+        self.e0.execute('del temp_field')
+
+        return item
+
+    def remoteDifference(self, key1, key2, keyresult):
+
+        if self.useMPI:
+
+            root = 0
+
+            # Gather
+            code_reduce = 'temp_%(key)s = comm.reduce(%(key)s, root=%(root)d)'
+            self.dview.execute(code_reduce%{'key': key1, 'root': root})
+            self.dview.execute(code_reduce%{'key': key2, 'root': root})
+
+            # Difference
+            code_difference = '%(keyresult)s = temp_%(key1)s - temp_%(key2)s'
+            self.e0.execute(code_difference%{'key1': key1, 'key2': key2, 'keyresult': keyresult})
+
+            # Broadcast
+            code = 'if rank != 0: %(key)s = None\n%(key)s = comm.bcast(%(key)s, root=%(root)d)'
+            self.dview.execute(code%{'key': keyresult, 'root': root})
+
+            # Clear
+            self.e0.execute('del temp_%s'%(key1,))
+            self.e0.execute('del temp_%s'%(key2,))
+
+        else:
+            item1 = reduce(np.add, self.dview[key1])
+            item2 = reduce(np.add, self.dview[key2])
+
+            item = item1 - item2
+            self.dview[keyresult] = item
+
+    def remoteOpGatherFirst(self, op, key1, key2, keyresult):
+
+        if self.useMPI:
+
+            root = 0
+
+            # Gather
+            code_reduce = 'temp_%(key)s = comm.reduce(%(key)s, root=%(root)d)'
+            self.dview.execute(code_reduce%{'key': key1, 'root': root})
+
+            # Difference
+            code_difference = '%(keyresult)s = temp_%(key1)s %(op)s %(key2)s'
+            self.e0.execute(code_difference%{'op': op, 'key1': key1, 'key2': key2, 'keyresult': keyresult})
+
+            # Broadcast
+            code = 'if rank != 0: %(key)s = None\n%(key)s = comm.bcast(%(key)s, root=%(root)d)'
+            self.dview.execute(code%{'key': keyresult, 'root': root})
+
+            # Clear
+            self.e0.execute('del temp_%s'%(key1,))
+
+        else:
+            item1 = reduce(np.add, self.dview[key1])
+            item2 = self.e0[key2] # Assumes that any arbitrary worker has this information
+
+            item = eval('item1 %s item2'%(op,))
+            self.dview[keyresult] = item
+
+    def remoteDifferenceGatherFirst(self, *args):
+        self.remoteOpGatherFirst('-', *args)
+
+    def remoteSrcEstGatherFirst(self, keyresult, key1, key2, individual=False):
+
+        if self.useMPI:
+
+            root = 0
+
+            # # Gather
+            # code_reduce = 'temp_%(key)s = comm.reduce(%(key)s, root=%(root)d)'
+            # self.dview.execute(code_reduce%{'key': key1, 'root': root})
+
+            # SrcEst
+            if individual:
+                code_srcest = '%(keyresult)s = (%(key2)s.conj() * %(key1)s).sum(axis=1) / (%(key1)s.conj() * %(key1)s).sum(axis=1)'
+            else:
+                code_srcest = '%(keyresult)s = (%(key2)s.conj() * %(key1)s).sum() / (%(key1)s.conj() * %(key1)s).sum()'
+            self.e0.execute(code_srcest%{'key1': key1, 'key2': key2, 'keyresult': keyresult})
+
+            # Broadcast
+            code = 'if rank != %(root)d: %(key)s = None\n%(key)s = comm.bcast(%(key)s, root=%(root)d)'
+            self.dview.execute(code%{'key': keyresult, 'root': root})
+
+        else:
+
+            item1 = reduce(np.add, self.dview[key1])
+            item2 = self.e0[key2]
+
+            if individual:
+                item = (item2.conj() * item1).sum(axis=1) / (item1.conj() * item1).sum(axis=1)
+            else:
+                item = (item2.conj() * item1).sum() / (item1.conj() * item1).sum()
+
+            self.dview[keyresult] = item
+
+    def remoteApplySrc(self, keyData, keySrc):
+
+        code = '%(keyData)s = %(keySrc)s * %(keyData)s'
+        self.dview.execute(code%{'keyData': keyData, 'keySrc': keySrc})
+
+    # def normFromDifference(self, key):
+
+    #     code = 'temp_norm%(key)s = (%(key)s * %(key)s.conj()).sum(0).sum(0)'
+    #     self.e0.execute(code%{'key': key})
+    #     code = 'temp_norm%(key)s = {key: np.sqrt(temp_norm%(key)s[key]).real for key in temp_norm%(key)s.keys()}'
+    #     self.e0.execute(code%{'key': key})
+    #     result = CommonReducer(self.e0['temp_norm%s'%(key,)])
+    #     self.e0.execute('del temp_norm%s'%(key,))
+
+    #     return result
+
+    def normFromDifference(self, key):
+
+        code = 'temp_norm = (%(key)s * %(key)s.conj()).sum(0).sum(0)'
+        self.e0.execute(code%{'key': key})
+        code = 'temp_norm = {key: np.sqrt(temp_norm[key]).real for key in temp_norm}'
+        self.e0.execute(code%{'key': key})
+        result = CommonReducer(self.e0['temp_norm'])
+        self.e0.execute('del temp_norm')
+
+        return result
+
+    @staticmethod
+    @interactive
+    def _reduceJob(worker, root, endpoint, key):
+
+        from ipyparallel.error import UnmetDependency
+        if not rank == worker:
+            raise UnmetDependency
+
+        # code = '%(endpoint)s.globalFields["%(key)s"] = comm.reduce(%(endpoint)s.localFields["%(key)s"], root=%(root)d)'
+        # exec(code%{'endpoint': endpoint, 'key': key, 'root': root})
+
+        if key not in endpoint.localFields:
+            endpoint.localFields[key] = endpoint.fieldspec[key]()
+
+        endpoint.globalFields[key] = comm.reduce(endpoint.localFields[key], root=root)
+
+    @staticmethod
+    def _adjustMKLVectorization(nt=1):
+        try:
+            import mkl
+            mkl.set_num_threads(nt)
+        except ImportError:
+            pass
+
+    @staticmethod
+    def _cdSame(rc):
+        import os
+
+        dview = rc[:]
+
+        home = os.getenv('HOME')
+        cwd = os.getcwd()
+
+        @interactive
+        def cdrel(relpath):
+            import os
+            home = os.getenv('HOME')
+            fullpath = os.path.join(home, relpath)
+            try:
+                os.chdir(fullpath)
+            except OSError:
+                return False
+            else:
+                return True
+
+        if cwd.find(home) == 0:
+            relpath = cwd[len(home)+1:]
+            return all(rc[:].apply_sync(cdrel, relpath))

SimPEG/Problem.py

@@ -1,278 +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[1] == 1 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 out.shape[0] == out.size:
-            out = Utils.mkvc(out)
-        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])
-        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.ndim == 2 and fieldI.shape[1] == 1:
-                        fieldI = Utils.mkvc(fieldI)
-                    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):
     """
@@ -289,15 +18,27 @@ 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.
+
+    @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):
@@ -333,7 +74,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,261 @@
+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)
+
+    @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
+
+    @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.')
+
+
+    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

@@ -261,22 +261,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 +310,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/Resources/Parallel/SystemGraph.html

@@ -0,0 +1,75 @@
+<div id="%(uniqueID)s"></div>
+<style>
+.node {stroke: #fff; stroke-width: 1.5px;}
+.link {stroke: #999; stroke-opacity: .6;}
+</style>
+
+<script>
+var True = true;
+var False = false;
+var graph = %(JSONData)s;
+
+require.config({paths: {d3: "http://d3js.org/d3.v3.min"}});
+require(["d3"], function(d3) {
+
+    var width = 800, height = 450, radius = 8;
+
+    var color = d3.scale.category10();
+    var domain = [0, 1, 2, 3];
+    color.domain(domain);
+
+    var force = d3.layout.force()
+        .charge(-120)
+        .linkDistance(30)
+        .size([width, height]);
+
+    var svg = d3.select("#%(uniqueID)s").select("svg");
+    if (svg.empty()) {
+        svg = d3.select("#%(uniqueID)s").append("svg")
+                    .attr("width", width)
+                    .attr("height", height);
+    }
+
+    force.nodes(graph.nodes)
+        .links(graph.links)
+        .start();
+
+    var link = svg.selectAll(".link")
+        .data(graph.links)
+        .enter().append("line")
+        .attr("class", "link");
+
+    var node = svg.selectAll(".node")
+        .data(graph.nodes)
+        .enter().append("circle")
+        .attr("class", "node")
+        .attr("r", radius)
+        .style("fill", function(d) {
+            return color(d.status); 
+        })
+
+        .call(force.drag);
+
+    node.append("title")
+        .text(function(d) { return d.id; });
+
+    node.on("dblclick", function() {
+        n = d3.select(this);
+        name = n.text();
+        graph = IPython.notebook.get_selected_cell().get_text();
+        cell = IPython.notebook.insert_cell_below();
+        cell.set_text(graph + ".node['" + name + "'].get('jobs', [])");
+    })
+
+    force.on("tick", function() {
+        link.attr("x1", function(d) { return d.source.x; })
+            .attr("y1", function(d) { return d.source.y; })
+            .attr("x2", function(d) { return d.target.x; })
+            .attr("y2", function(d) { return d.target.y; });
+
+        node.attr("cx", function(d) { return d.x; })
+            .attr("cy", function(d) { return d.y; });
+    });
+
+});
+</script>

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,6 +13,7 @@ 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 = {}
@@ -119,14 +120,12 @@ class BaseSrc(object):
     rxList = None #: SimPEG Receiver List
     rxPair = BaseRx
 
-    def __init__(self, loc, srcType, rxList, **kwargs):
+    def __init__(self, rxList, **kwargs):
         assert type(rxList) is list, 'rxList must be a list'
         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.loc    = loc
-        self.srcType = srcType
+        self.uid = str(uuid.uuid4())
         self.rxList = rxList
         Utils.setKwargs(self, **kwargs)
 
@@ -146,6 +145,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:
@@ -227,6 +227,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):
@@ -360,3 +373,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])
@@ -9,14 +10,14 @@ class DataAndFieldsTest(unittest.TestCase):
         XYZ = Utils.ndgrid(x,x,np.r_[0.])
         srcLoc = np.r_[0,0,0.]
         rxList0 = Survey.BaseRx(XYZ, 'exi')
-        Src0 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0])
+        Src0 = Survey.BaseSrc([rxList0], loc=srcLoc)
         rxList1 = Survey.BaseRx(XYZ, 'bxi')
-        Src1 = Survey.BaseSrc(srcLoc, 'VMD', [rxList1])
+        Src1 = Survey.BaseSrc([rxList1], loc=srcLoc)
         rxList2 = Survey.BaseRx(XYZ, 'bxi')
-        Src2 = Survey.BaseSrc(srcLoc, 'VMD', [rxList2])
+        Src2 = Survey.BaseSrc([rxList2], loc=srcLoc)
         rxList3 = Survey.BaseRx(XYZ, 'bxi')
-        Src3 = Survey.BaseSrc(srcLoc, 'VMD', [rxList3])
-        Src4 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
+        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)
@@ -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
@@ -80,9 +62,9 @@ class DataAndFieldsTest(unittest.TestCase):
 
         b = np.random.rand(F.mesh.nF,1)
         F[self.Src0, 'b'] = b
-        self.assertTrue(np.all(F[self.Src0, 'b'] == Utils.mkvc(b)))
+        self.assertTrue(np.all(F[self.Src0, 'b'] == b))
 
-        b = np.random.rand(F.mesh.nF)
+        b = np.random.rand(F.mesh.nF,1)
         F[self.Src0, 'b'] = b
         self.assertTrue(np.all(F[self.Src0, 'b'] == b))
 
@@ -96,10 +78,10 @@ class DataAndFieldsTest(unittest.TestCase):
 
         b = np.random.rand(F.mesh.nF, 2)
         F[[self.Src0, self.Src1],'b'] = b
-        self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF,))
-        self.assertTrue(F[self.Src0,'b'].shape == (F.mesh.nF,))
-        self.assertTrue(np.all(F[self.Src0,'b'] == b[:,0]))
-        self.assertTrue(np.all(F[self.Src1,'b'] == b[:,1]))
+        self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF,1))
+        self.assertTrue(F[self.Src0,'b'].shape == (F.mesh.nF,1))
+        self.assertTrue(np.all(F[self.Src0,'b'] == Utils.mkvc(b[:,0],2)))
+        self.assertTrue(np.all(F[self.Src1,'b'] == Utils.mkvc(b[:,1],2)))
 
     def test_assertions(self):
         freq = [self.Src0, self.Src1]
@@ -122,17 +104,16 @@ class FieldsTest_Alias(unittest.TestCase):
         XYZ = Utils.ndgrid(x,x,np.r_[0.])
         srcLoc = np.r_[0,0,0.]
         rxList0 = Survey.BaseRx(XYZ, 'exi')
-        Src0 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0])
+        Src0 = Survey.BaseSrc([rxList0],loc=srcLoc)
         rxList1 = Survey.BaseRx(XYZ, 'bxi')
-        Src1 = Survey.BaseSrc(srcLoc, 'VMD', [rxList1])
+        Src1 = Survey.BaseSrc([rxList1],loc=srcLoc)
         rxList2 = Survey.BaseRx(XYZ, 'bxi')
-        Src2 = Survey.BaseSrc(srcLoc, 'VMD', [rxList2])
+        Src2 = Survey.BaseSrc([rxList2],loc=srcLoc)
         rxList3 = Survey.BaseRx(XYZ, 'bxi')
-        Src3 = Survey.BaseSrc(srcLoc, 'VMD', [rxList3])
-        Src4 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
+        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)
         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
@@ -158,7 +139,7 @@ class FieldsTest_Alias(unittest.TestCase):
 
         e = np.random.rand(F.mesh.nE,1)
         F[self.Src0, 'e'] = e
-        self.assertTrue(np.all(F[self.Src0, 'b'] == F.mesh.edgeCurl * Utils.mkvc(e)))
+        self.assertTrue(np.all(F[self.Src0, 'b'] == F.mesh.edgeCurl * e))
 
         def f():
             F[self.Src0, 'b'] = F[self.Src0, 'b']
@@ -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)
@@ -193,14 +174,14 @@ class FieldsTest_Time(unittest.TestCase):
         XYZ = Utils.ndgrid(x,x,np.r_[0.])
         srcLoc = np.r_[0,0,0.]
         rxList0 = Survey.BaseRx(XYZ, 'exi')
-        Src0 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0])
+        Src0 = Survey.BaseSrc([rxList0], loc=srcLoc)
         rxList1 = Survey.BaseRx(XYZ, 'bxi')
-        Src1 = Survey.BaseSrc(srcLoc, 'VMD', [rxList1])
+        Src1 = Survey.BaseSrc([rxList1], loc=srcLoc)
         rxList2 = Survey.BaseRx(XYZ, 'bxi')
-        Src2 = Survey.BaseSrc(srcLoc, 'VMD', [rxList2])
+        Src2 = Survey.BaseSrc([rxList2], loc=srcLoc)
         rxList3 = Survey.BaseRx(XYZ, 'bxi')
-        Src3 = Survey.BaseSrc(srcLoc, 'VMD', [rxList3])
-        Src4 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
+        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)
         prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10.,3), (20.,2)])
@@ -249,7 +230,7 @@ class FieldsTest_Time(unittest.TestCase):
 
         b = np.random.rand(F.mesh.nF,1,nT)
         F[self.Src0, 'b', 0] = b[:,:,0]
-        self.assertTrue(np.all(F[self.Src0, 'b', 0] == b[:,0,0]))
+        self.assertTrue(np.all(F[self.Src0, 'b', 0] == Utils.mkvc(b[:,0,0],2)))
 
         phi = np.random.rand(F.mesh.nC,2,nT)
         F[[self.Src0,self.Src1], 'phi'] = phi
@@ -265,10 +246,10 @@ class FieldsTest_Time(unittest.TestCase):
         self.assertTrue(F[self.Src0,'b'].shape == (F.mesh.nF,nT))
         self.assertTrue(np.all(F[self.Src0,'b'] == b[:,0,:]))
         self.assertTrue(np.all(F[self.Src1,'b'] == b[:,1,:]))
-        self.assertTrue(np.all(F[self.Src0,'b',1] == b[:,0,1]))
-        self.assertTrue(np.all(F[self.Src1,'b',1] == b[:,1,1]))
-        self.assertTrue(np.all(F[self.Src0,'b',4] == b[:,0,4]))
-        self.assertTrue(np.all(F[self.Src1,'b',4] == b[:,1,4]))
+        self.assertTrue(np.all(F[self.Src0,'b',1] == Utils.mkvc(b[:,0,1],2)))
+        self.assertTrue(np.all(F[self.Src1,'b',1] == Utils.mkvc(b[:,1,1],2)))
+        self.assertTrue(np.all(F[self.Src0,'b',4] == Utils.mkvc(b[:,0,4],2)))
+        self.assertTrue(np.all(F[self.Src1,'b',4] == Utils.mkvc(b[:,1,4],2)))
 
 
         b = np.random.rand(F.mesh.nF, 2, nT)
@@ -295,14 +276,14 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
         XYZ = Utils.ndgrid(x,x,np.r_[0.])
         srcLoc = np.r_[0,0,0.]
         rxList0 = Survey.BaseRx(XYZ, 'exi')
-        Src0 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0])
+        Src0 = Survey.BaseSrc( [rxList0],loc=srcLoc)
         rxList1 = Survey.BaseRx(XYZ, 'bxi')
-        Src1 = Survey.BaseSrc(srcLoc, 'VMD', [rxList1])
+        Src1 = Survey.BaseSrc( [rxList1],loc=srcLoc)
         rxList2 = Survey.BaseRx(XYZ, 'bxi')
-        Src2 = Survey.BaseSrc(srcLoc, 'VMD', [rxList2])
+        Src2 = Survey.BaseSrc( [rxList2],loc=srcLoc)
         rxList3 = Survey.BaseRx(XYZ, 'bxi')
-        Src3 = Survey.BaseSrc(srcLoc, 'VMD', [rxList3])
-        Src4 = Survey.BaseSrc(srcLoc, 'VMD', [rxList0, rxList1, rxList2, rxList3])
+        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)
         prob = Problem.BaseTimeProblem(mesh, timeSteps=[(10.,3), (20.,2)])
@@ -344,7 +325,7 @@ class FieldsTest_Time_Aliased(unittest.TestCase):
         self.assertTrue(np.all(F[self.Src0, 'e', :] == e[:,0,:] ))
         self.assertTrue(np.all(F[self.Src1, 'e', :] == e[:,1,:] ))
         for t in range(nT):
-            self.assertTrue(np.all(F[self.Src1, 'e', t] == e[:,1,t] ))
+            self.assertTrue(np.all(F[self.Src1, 'e', t] == Utils.mkvc(e[:,1,t],2) ))
 
         b = np.random.rand(F.mesh.nF,nT)
         F[self.Src0, 'b',:] = b
@@ -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/DataWrappers.py

@@ -0,0 +1,91 @@
+
+class CommonReducer(dict):
+    '''
+    Object based on 'dict' that implements the binary addition (obj1 + obj1) and
+    accumulation (obj += obj2). These operations pass through to the entries in
+    the commonReducer.
+    
+    Instances of commonReducer are also callable, with the syntax:
+        cr(key, value)
+    this is equivalent to cr += {key: value}.
+    ''' 
+
+    DISALLOWED = ['__getinitargs__', '__getnewargs__', '__getstate__', '__setstate__']
+
+    def __init__(self, *args, **kwargs):
+        dict.__init__(self, *args, **kwargs)
+
+    def __add__(self, other):
+        result = CommonReducer(self)
+        for key in other.keys():
+            if key in result:
+                result[key] = self[key] + other[key]
+            else:
+                result[key] = other[key]
+
+        return result
+
+    def __iadd__(self, other):
+        for key in other.keys():
+            if key in self:
+                self[key] += other[key]
+            else:
+                self[key] = other[key]
+
+        return self
+
+    def __mul__(self, other):
+        result = CommonReducer()
+        for key in other.keys():
+            if key in self:
+                result[key] = self[key] * other[key]
+
+        return result
+
+    def __sub__(self, other):
+        result = CommonReducer()
+        for key in other.keys():
+            if key in self:
+                result[key] = self[key] - other[key]
+
+        return result
+
+    def __div__(self, other):
+        result = CommonReducer()
+        for key in other.keys():
+            if key in self:
+                result[key] = self[key] / other[key]
+
+        return result
+
+    def __getattr__(self, attr):
+
+        if not attr in self.DISALLOWED and all((getattr(self[key], attr, None) is not None for key in self)):
+
+            if any((callable(getattr(self[key], attr)) for key in self)):
+
+                def wrapperFunction(*args, **kwargs):
+
+                    innerresult = CommonReducer({key: getattr(self[key], attr, None)(*args, **kwargs) for key in self})
+
+                    if not all((innerresult[key] is None for key in innerresult)):
+                        return innerresult
+
+                result = wrapperFunction
+
+            else:
+                return CommonReducer({key: getattr(self[key], attr) for key in self})
+        else:
+            raise AttributeError('\'CommonReducer\' object has no attribute \'%s\', and it could not be satisfied through cascade lookup'%attr)
+
+        return result
+
+    def copy(self):
+
+        return CommonReducer(self)
+
+    def __call__(self, key, result):
+        if key in self:
+            self[key] += result
+        else:
+            self[key] = result

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];
 

SimPEG/Utils/__init__.py

@@ -5,6 +5,7 @@ from curvutils import volTetra, faceInfo, indexCube
 from interputils import interpmat
 from ipythonutils import easyAnimate as animate
 from CounterUtils import *
+from DataWrappers import *
 import ModelBuilder
 import SolverUtils
 

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/__init__.py

@@ -13,9 +13,10 @@ import InvProblem
 import Optimization
 import Directives
 import Inversion
+import Parallel
 import Tests
 
-__version__   = '0.1.2'
+__version__   = '0.1.3'
 __author__    = 'Rowan Cockett'
 __license__   = 'MIT'
 __copyright__ = 'Copyright 2014 Rowan Cockett'

docs/conf.py

@@ -51,9 +51,9 @@ copyright = u'2013, SimPEG Developers'
 # built documents.
 #
 # The short X.Y version.
-version = '0.1.2'
+version = '0.1.3'
 # The full version, including alpha/beta/rc tags.
-release = '0.1.2'
+release = '0.1.3'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

setup.py

@@ -33,7 +33,7 @@ with open("README.rst") as f:
 
 setup(
     name = "SimPEG",
-    version = "0.1.2",
+    version = "0.1.3",
     packages = find_packages(),
     install_requires = ['numpy>=1.7',
                         'scipy>=0.13'
@@ -44,7 +44,7 @@ setup(
     long_description = LONG_DESCRIPTION,
     license = "MIT",
     keywords = "geophysics inverse problem",
-    url = "http://simpeg.3ptscience.com/",
+    url = "http://simpeg.xyz/",
     download_url = "http://github.com/simpeg/simpeg",
     classifiers=CLASSIFIERS,
     platforms = ["Windows", "Linux", "Solaris", "Mac OS-X", "Unix"],