Fixed bugs in meshutils, writing VTR files

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

.travis.yml

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

SimPEG/DataMisfit.py

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

SimPEG/Directives.py

@@ -8,6 +8,34 @@ class InversionDirective(object):
     def __init__(self, **kwargs):
         Utils.setKwargs(self, **kwargs)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     @property
     def inversion(self):
         """This is the inversion of the InversionDirective instance."""
@@ -207,6 +235,76 @@ class SaveOutputEveryIteration(_SaveEveryIteration):
         f.close()
 
 
+class SaveOutputDictEveryIteration(_SaveEveryIteration):
+    """SaveOutputDictEveryIteration"""
+
+    def initialize(self):
+        print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName
+
+    def endIter(self):
+        # Save the data.
+        ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
+        phi_ms = 0.5*ms.dot(ms)
+        if self.reg.smoothModel == True:
+            mref = self.reg.mref
+        else:
+            mref = 0
+        mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
+        phi_mx = 0.5 * mx.dot(mx)
+        if self.prob.mesh.dim==2:
+            my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
+            phi_my = 0.5 * my.dot(my)
+        else:
+            phi_my = 'NaN'
+        if self.prob.mesh.dim==3:
+            mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
+            phi_mz = 0.5 * mz.dot(mz)
+        else:
+            phi_mz = 'NaN'
+
+
+        # Save the file as a npz
+        np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel)
+
+
+class SaveOutputDictEveryIteration(_SaveEveryIteration):
+    """SaveOutputDictEveryIteration
+
+    A directive that saves some relevant information from the inversion run to a numpy .npz dictionary file (see numpy.savez function for further info).
+
+
+
+    """
+
+    def initialize(self):
+        print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName
+
+    def endIter(self):
+        # Save the data.
+        ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
+        phi_ms = 0.5*ms.dot(ms)
+        if self.reg.smoothModel == True:
+            mref = self.reg.mref
+        else:
+            mref = 0
+        mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
+        phi_mx = 0.5 * mx.dot(mx)
+        if self.prob.mesh.dim==2:
+            my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
+            phi_my = 0.5 * my.dot(my)
+        else:
+            phi_my = 'NaN'
+        if self.prob.mesh.dim==3:
+            mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
+            phi_mz = 0.5 * mz.dot(mz)
+        else:
+            phi_mz = 'NaN'
+
+
+        # Save the file as a npz
+        np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
+
+
 
 
 # class UpdateReferenceModel(Parameter):

SimPEG/Dispatcher.py

@@ -1,112 +0,0 @@
-### 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

@@ -12,6 +12,34 @@ class Fields(object):
     aliasFields = None  #: Aliased fields, a dict with [alias, location, function], e.g. {"b":["e","F",lambda(F,e,ind)]}
     dtype = float       #: dtype is the type of the storage matrix. This can be a dictionary.
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     def __init__(self, mesh, survey, **kwargs):
         self.survey = survey
         self.mesh = mesh

SimPEG/Maps.py

@@ -17,6 +17,30 @@ class IdentityMap(object):
         Utils.setKwargs(self, **kwargs)
         self.mesh = mesh
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+
     @property
     def nP(self):
         """
@@ -289,7 +313,7 @@ class FullMap(IdentityMap):
     """
     FullMap
 
-    Given a scalar, the FullMap maps the value to the 
+    Given a scalar, the FullMap maps the value to the
     full model space.
     """
 
@@ -314,8 +338,8 @@ class FullMap(IdentityMap):
             :rtype: numpy.array
             :return: derivative of transformed model
         """
-        return np.ones([self.mesh.nC,1])     
-      
+        return np.ones([self.mesh.nC,1])
+
 
 class Vertical1DMap(IdentityMap):
     """Vertical1DMap

SimPEG/Mesh/InnerProducts.py

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

SimPEG/Optimization.py

@@ -115,6 +115,34 @@ class Minimize(object):
 
         Utils.setKwargs(self, **kwargs)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     @property
     def callback(self):
         return getattr(self, '_callback', None)

SimPEG/Parallel.py

@@ -1,638 +0,0 @@
-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

@@ -22,6 +22,34 @@ class BaseProblem(object):
 
     PropMap = None    #: A SimPEG PropertyMap class.
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     @property
     def mapping(self):
         "A SimPEG.Map instance or a property map is PropMap is not None"
@@ -32,8 +60,8 @@ class BaseProblem(object):
             val._assertMatchesPair(self.mapPair)
             self._mapping = val
         else:
-            self._mapping = self.PropMap(val) 
-    
+            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."

SimPEG/PropMaps.py

@@ -13,6 +13,34 @@ class Property(object):
         self.doc = doc
         Utils.setKwargs(self, **kwargs)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     @property
     def propertyLink(self):
         "Can be something like: ('sigma', Maps.ReciprocalMap)"
@@ -118,6 +146,36 @@ class PropModel(object):
         self.vector  = vector
         assert len(self.vector) == self.nP
 
+    # Pickleing support methods
+    # def __reduce__(self):
+    #     return (dict,{self.propMap,self.vector})
+
+    # def __getstate__(self):
+    #     '''
+    #     Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+    #     Used when doing:
+    #         pickle.dump(pickleFile,object)
+    #     '''
+    #     self.__class__ = ProbModel
+    #     odict = {}
+    #     odict['vec'] = self.__dict__['vector']
+    #     odict['pMap'] = self.__dict__['propMap']
+    #     # Return the dict
+    #     return odict
+
+    # def __setstate__(self,odict):
+    #     '''
+    #     Function that sets a pickle dictionary in to an object.
+
+    #     Used when doing:
+    #         object = pickle.load(pickleFile)
+    #     '''
+    #     # Update the dict
+    #     # Re-hook the methods to the object
+    #     self.propMap = odict['prMap']
+    #     self.vector = odict['vec']
+
     @property
     def nP(self):
         inds = []
@@ -207,6 +265,24 @@ class PropMap(object):
         else:
             raise Exception('mappings must be a dict, a mapping, or a list of tuples.')
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        pass
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        pass
 
     def setup(self, maps, slices=None):
         """
@@ -239,7 +315,7 @@ class PropMap(object):
             setattr(self, '%sMap'%name, mapping)
             setattr(self, '%sIndex'%name, slices.get(name, slice(nP, nP + mapping.nP)))
             nP += mapping.nP
-        self.nP = nP 
+        self.nP = nP
 
     @property
     def defaultInvProp(self):

SimPEG/Regularization.py

@@ -27,6 +27,34 @@ class BaseRegularization(object):
         self.mapping = mapping or Maps.IdentityMap(mesh)
         self.mapping._assertMatchesPair(self.mapPair)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
     @property
     def parent(self):
         """This is the parent of the regularization."""
@@ -311,7 +339,7 @@ class Tikhonov(BaseRegularization):
         if self.smoothModel == True:
             mD1 = self.mapping.deriv(m)
             mD2 = self.mapping.deriv(m - self.mref)
-            r1 = self.Wsmooth * ( self.mapping * (m)) 
+            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 )

SimPEG/Resources/Parallel/SystemGraph.html

@@ -1,75 +0,0 @@
-<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

@@ -19,6 +19,35 @@ class BaseRx(object):
         self._Ps = {}
         Utils.setKwargs(self, **kwargs)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
+
     @property
     def rxType(self):
         """Receiver Type"""
@@ -129,6 +158,33 @@ class BaseSrc(object):
         self.rxList = rxList
         Utils.setKwargs(self, **kwargs)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
 
     @property
     def nD(self):
@@ -153,6 +209,25 @@ class Data(object):
         if v is not None:
             self.fromvec(v)
 
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        pass
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        pass
+
     def _ensureCorrectKey(self, key):
         if type(key) is tuple:
             if len(key) is not 2:
@@ -210,11 +285,39 @@ class BaseSurvey(object):
     mtrue = None     #: True model, if data is synthetic
 
     counter = None   #: A SimPEG.Utils.Counter object
+    srcPair = BaseSrc  #: Source Pair
 
     def __init__(self, **kwargs):
         Utils.setKwargs(self, **kwargs)
 
-    srcPair = BaseSrc  #: Source Pair
+    # Pickleing support methods
+    def __getstate__(self):
+        '''
+        Method that makes the dictionary of the object pickleble, removes non-pickleble elements of the object.
+
+        Used when doing:
+            pickle.dump(pickleFile,object)
+        '''
+        odict = self.__dict__.copy()
+        # Remove fields that are not needed
+        del odict['hook']
+        del odict['setKwargs']
+        # Return the dict
+        return odict
+
+    def __setstate__(self,odict):
+        '''
+        Function that sets a pickle dictionary in to an object.
+
+        Used when doing:
+            object = pickle.load(pickleFile)
+        '''
+        # Update the dict
+        self.__dict__.update(odict)
+        # Re-hook the methods to the object
+        Utils.codeutils.hook(self,Utils.codeutils.hook)
+        Utils.codeutils.hook(self,Utils.codeutils.setKwargs)
+
 
     @property
     def srcList(self):

SimPEG/Utils/DataWrappers.py

@@ -1,91 +0,0 @@
-
-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

@@ -170,16 +170,58 @@ def scalarConductivity(ccMesh,pFunction):
 
     return mkvc(sigma)
 
+def layeredModel(ccMesh, layerTops, layerValues):
+    """
+        Define a layered model from layerTops (z-positive up)
+
+        :param numpy.array ccMesh: cell-centered mesh
+        :param numpy.array layerTops: z-locations of the tops of each layer
+        :param numpy.array layerValue: values of the property to assign for each layer (starting at the top) 
+        :rtype: numpy.array
+        :return: M, layered model on the mesh 
+    """
+
+    descending = np.linalg.norm(sorted(layerTops, reverse=True) - layerTops) < 1e-20
+
+    # TODO: put an error check to make sure that there is an ordering... needs to work with inf elts 
+    # assert ascending or descending, "Layers must be listed in either ascending or descending order"
+
+    # start from bottom up
+    if not descending:
+        zprop = np.hstack([mkvc(layerTops,2),mkvc(layerValues,2)])
+        zprop.sort(axis=0)
+        layerTops, layerValues  = zprop[::-1,0], zprop[::-1,1]
+
+    # put in vector form
+    layerTops, layerValues = mkvc(layerTops), mkvc(layerValues)
+
+    # initialize with bottom layer
+    dim = ccMesh.shape[1]
+    if dim == 3:
+        z = ccMesh[:,2]
+    elif dim == 2:
+        z = ccMesh[:,1]
+    elif dim == 1:
+        z = ccMesh[:,0]
+
+    model = np.zeros(ccMesh.shape[0])
+
+    for i, top in enumerate(layerTops):
+        zind = z <= top 
+        model[zind] = layerValues[i]
+
+    return model 
+
 
 
 def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]):
     """
-        Create a random model by convolving a kernal with a
+        Create a random model by convolving a kernel with a
         uniformly distributed model.
 
         :param int,tuple shape: shape of the model.
         :param int seed: pick which model to produce, prints the seed if you don't choose.
-        :param numpy.ndarray,list anisotropy: this is the (3 x n) blurring kernal that is used.
+        :param numpy.ndarray,list anisotropy: this is the (3 x n) blurring kernel that is used.
         :param int its: number of smoothing iterations
         :param list bounds: bounds on the model, len(list) == 2
         :rtype: numpy.ndarray

SimPEG/Utils/__init__.py

@@ -5,7 +5,6 @@ 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/meshutils.py

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

SimPEG/__init__.py

@@ -13,7 +13,6 @@ import InvProblem
 import Optimization
 import Directives
 import Inversion
-import Parallel
 import Tests
 
 __version__   = '0.1.3'