Serialize the Properties.

SimPEG/Examples/EM_FDEM_SusEffects.py

@@ -1,154 +0,0 @@
-from SimPEG import *
-from SimPEG import EM
-from scipy.constants import mu_0
-
-def run(plotIt=True):
-    """
-        EM: FDEM: Effects of susceptibility
-        ===================================
-
-        When airborne freqeuncy domain EM (AFEM) survey is flown over
-        the earth including significantly susceptible bodies (magnetite-rich rocks),
-        negative data is often observed in the real part of the lowest frequency
-        (e.g. Dighem system 900 Hz). This phenomenon mostly based upon magnetization
-        occurs due to a susceptible body when the magnetic field is applied.
-
-        To clarify what is happening in the earth when we are exciting the earth with
-        a loop source in the frequency domain we run three forward modelling:
-
-            - F[:math:`\sigma`, :math:`\mu`]: Anomalous conductivity and susceptibility
-            - F[:math:`\sigma`, :math:`\mu_0`]: Anomalous conductivity
-            - F[:math:`\sigma_{air}`, :math:`\mu_0`]: primary field
-
-        We plot vector magnetic fields in the earth. For secondary fields we provide
-        F[:math:`\sigma`, :math:`\mu`]-F[:math:`\sigma`, :math:`\mu_0`]. Following
-        figure show both real and parts.
-
-    """
-    # Generate Cylindrical mesh
-    cs, ncx, ncz, npad = 5, 25, 24, 20.
-    hx = [(cs,ncx), (cs,npad,1.3)]
-    hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
-    mesh = Mesh.CylMesh([hx,1,hz], '00C')
-    sighalf = 1e-3
-    sigma = np.ones(mesh.nC)*1e-8
-    sigmahomo = sigma.copy()
-    mu = np.ones(mesh.nC)*mu_0
-    sigma[mesh.gridCC[:,-1]<0.] = sighalf
-    blkind  = np.logical_and(mesh.gridCC[:,0]<30., (mesh.gridCC[:,2]<0)&(mesh.gridCC[:,2]>-150)&(mesh.gridCC[:,2]<-50))
-    sigma[blkind] = 1e-1
-    mu[blkind] = mu_0*1.1
-    offset = 0.
-    frequency = np.r_[10., 100., 1000.]
-    rx0 = EM.FDEM.Rx(np.array([[8., 0., 30.]]), 'bzr')
-    rx1 = EM.FDEM.Rx(np.array([[8., 0., 30.]]), 'bzi')
-    srcLists = []
-    nfreq = frequency.size
-    for ifreq in range(nfreq):
-        src = EM.FDEM.Src.CircularLoop([rx0, rx1], frequency[ifreq], np.array([[0., 0., 30.]]), radius=5.)
-        srcLists.append(src)
-    survey = EM.FDEM.Survey(srcLists)
-    iMap = Maps.IdentityMap(nP=int(mesh.nC))
-    # Use PhysPropMap
-    maps = [('sigma', iMap), ('mu', iMap)]
-    prob = EM.FDEM.Problem_b(mesh, mapping=maps)
-    prob.Solver = MumpsSolver
-    survey.pair(prob)
-    m = np.r_[sigma, mu]
-    survey0 = EM.FDEM.Survey(srcLists)
-    prob0 = EM.FDEM.Problem_b(mesh, mapping=maps)
-    
-
-    try:
-        from pymatsolver import MumpsSolver
-        prb.Solver = MumpsSolver
-    except ImportError, e:
-            prb.Solver = SolverLU
-
-    survey0.pair(prob0)
-    m = np.r_[sigma, mu]
-    m0 = np.r_[sigma, np.ones(mesh.nC)*mu_0]
-    m00 = np.r_[np.ones(mesh.nC)*1e-8, np.ones(mesh.nC)*mu_0]
-    # Anomalous conductivity and susceptibility
-    F = prob.fields(m)
-    # Only anomalous conductivity
-    F0 = prob.fields(m0)
-    # Primary field
-    F00 = prob.fields(m00)
-
-    if plotIt:
-        import matplotlib.pyplot as plt
-        def vizfields(ifreq=0, primsec="secondary",realimag="real"):
-
-            titles = ["F[$\sigma$, $\mu$]", "F[$\sigma$, $\mu_0$]", "F[$\sigma$, $\mu$]-F[$\sigma$, $\mu_0$]"]
-            actind  = np.logical_and(mesh.gridCC[:,0]<200., (mesh.gridCC[:,2]>-400)&(mesh.gridCC[:,2]<200))
-
-            if primsec=="secondary":
-                bCCprim = (mesh.aveF2CCV*F00[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')
-                bCC = (mesh.aveF2CCV*F[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')-bCCprim
-                bCC0 = (mesh.aveF2CCV*F0[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')-bCCprim
-            elif primsec=="primary":
-                bCC = (mesh.aveF2CCV*F[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')
-                bCC0 = (mesh.aveF2CCV*F0[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')
-
-            XYZ = mesh.gridCC[actind,:]
-            X = XYZ[:,0].reshape((31,43), order='F')
-            Z = XYZ[:,2].reshape((31,43), order='F')
-            bx = bCC[actind,0].reshape((31,43), order='F')
-            bz = bCC[actind,1].reshape((31,43), order='F')
-            bx0 = bCC0[actind,0].reshape((31,43), order='F')
-            bz0 = bCC0[actind,1].reshape((31,43), order='F')
-
-            bxsec = (bCC[actind,0]-bCC0[actind,0]).reshape((31,43), order='F')
-            bzsec = (bCC[actind,1]-bCC0[actind,1]).reshape((31,43), order='F')
-
-            absbreal = np.sqrt(bx.real**2+bz.real**2)
-            absbimag = np.sqrt(bx.imag**2+bz.imag**2)
-            absb0real = np.sqrt(bx0.real**2+bz0.real**2)
-            absb0imag = np.sqrt(bx0.imag**2+bz0.imag**2)
-
-            absbrealsec = np.sqrt(bxsec.real**2+bzsec.real**2)
-            absbimagsec = np.sqrt(bxsec.imag**2+bzsec.imag**2)
-
-            fig = plt.figure(figsize=(15,5))
-            ax1 = plt.subplot(131)
-            ax2 = plt.subplot(132)
-            ax3 = plt.subplot(133)
-            typefield="real"
-            scale=20
-            if realimag=="real":
-                ax1.contourf(X, Z,np.log10(absbreal), 100)
-                ax1.quiver(X, Z,bx.real/absbreal,bz.real/absbreal,scale=scale,width=0.005, alpha = 0.5)
-                ax2.contourf(X, Z,np.log10(absb0real), 100)
-                ax2.quiver(X, Z,bx0.real/absb0real,bz0.real/absb0real,scale=scale,width=0.005, alpha = 0.5)
-                ax3.contourf(X, Z,np.log10(absbrealsec), 100)
-                ax3.quiver(X, Z,bxsec.real/absbrealsec,bzsec.real/absbrealsec,scale=scale,width=0.005, alpha = 0.5)
-            elif realimag=="imag":
-                ax1.contourf(X, Z,np.log10(absbimag), 100)
-                ax1.quiver(X, Z,bx.imag/absbimag,bz.imag/absbimag,scale=scale,width=0.005, alpha = 0.5)
-                ax2.contourf(X, Z,np.log10(absb0imag), 100)
-                ax2.quiver(X, Z,bx0.imag/absb0imag,bz0.imag/absb0imag,scale=scale,width=0.005, alpha = 0.5)
-                ax3.contourf(X, Z,np.log10(absbimagsec), 100)
-                ax3.quiver(X, Z,bxsec.imag/absbimagsec,bzsec.imag/absbimagsec,scale=scale,width=0.005, alpha = 0.5)
-
-            ax = [ax1, ax2, ax3]
-            ax3.text(30, 50, ("Frequency=%5.2f Hz")%(frequency[ifreq]), color="k", fontsize=18)
-            ax2.text(30, 50, primsec, color="k", fontsize=18)
-            ax1.text(30, 50, realimag, color="k", fontsize=18)
-            for i, axtemp in enumerate(ax):
-                axtemp.plot(np.r_[0, 29.75], np.r_[-50, -50], 'w', lw=3)
-
-                axtemp.plot(np.r_[29.5, 29.5], np.r_[-50, -142.5], 'w', lw=3)
-                axtemp.plot(np.r_[0, 29.5], np.r_[-142.5, -142.5], 'w', lw=3)
-                axtemp.plot(np.r_[0, 100.], np.r_[0, 0], 'w', lw=3)
-                axtemp.set_ylim(-200, 100.)
-                axtemp.set_xlim(10, 100.)
-                axtemp.set_title(titles[i])
-            plt.show()
-            return fig, ax
-        fig1, ax1 = vizfields(1, primsec="primary", realimag="real")
-        fig2, ax2 = vizfields(1, primsec="secondary", realimag="real")
-        fig4, ax4 = vizfields(1, primsec="secondary", realimag="imag")
-
-if __name__ == '__main__':
-    run()

SimPEG/Examples/__init__.py

@@ -3,7 +3,6 @@
 ##### AUTOIMPORTS #####
 import EM_FDEM_1D_Inversion
 import EM_FDEM_Analytic_MagDipoleWholespace
-import EM_FDEM_SusEffects
 import EM_TDEM_1D_Inversion
 import FLOW_Richards_1D_Celia1990
 import Forward_BasicDirectCurrent
@@ -15,9 +14,8 @@ import Mesh_QuadTree_Creation
 import Mesh_QuadTree_FaceDiv
 import Mesh_QuadTree_HangingNodes
 import Mesh_Tensor_Creation
-import MT_1D_analytic_nlayer_Earth
 
-__examples__ = ["EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_FDEM_SusEffects", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_analytic_nlayer_Earth"]
+__examples__ = ["EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation"]
 
 ##### AUTOIMPORTS #####
 

SimPEG/Optimization.py

@@ -990,18 +990,4 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
                     cgFlag = 1
                 # End CG Iterations
 
-            # Take a gradient step on the active cells if exist
-            if temp != self.xc.size:
-
-                rhs_a = (Active) * -self.g
-
-                dm_i = max( abs( delx ) )
-                dm_a = max( abs(rhs_a) )
-
-                delx = delx + rhs_a * dm_i / dm_a /10.
-
-            # Only keep gradients going in the right direction on the active set
-            indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0))
-            delx[indx] = 0.
-
             return delx

SimPEG/PropMaps.py

@@ -12,6 +12,7 @@ class Property(object):
         # Set the default after all other params are set
         self.doc = doc
         Utils.setKwargs(self, **kwargs)
+        self._kwargs = kwargs
 
     @property
     def propertyLink(self):
@@ -110,6 +111,11 @@ class Property(object):
             return getattr(self.propMap, '_%sMap'%prop.name, None)
         return property(fget=fget)
 
+    def toJSON(self):
+        out = dict(doc=self.doc)
+        for k in self._kwargs:
+            out[k] = self._kwargs[k]
+        return out
 
 
 class PropModel(object):
@@ -189,6 +195,12 @@ class _PropMapMetaClass(type):
 
         return type(name.replace('PropMap', 'PropModel'), (PropModel, ), attrs)
 
+def fromPickle(name, properties, maps, slices):
+    attrs = dict()
+    for p in properties:
+        attrs[p] = Property(**properties[p])
+    PM = type(name, (PropMap,), attrs)
+    return PM(dict(maps=maps, slices=slices))
 
 class PropMap(object):
     __metaclass__ = _PropMapMetaClass
@@ -197,6 +209,7 @@ class PropMap(object):
         """
             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'])
@@ -239,7 +252,11 @@ 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._maps   = maps
+        self._slices = slices
+
+        self.nP = nP
 
     @property
     def defaultInvProp(self):
@@ -253,9 +270,27 @@ class PropMap(object):
             setattr(self, '%sMap'%name, None)
             setattr(self, '%sIndex'%name, None)
 
+        self._maps   = None
+        self._slices = 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
+
+    def __reduce__(self):
+
+        import cPickle
+        props = dict()
+        for p in self._properties:
+            props[p] = self._properties[p].toJSON()
+        className = self.__class__.__name__
+
+
+        pickledMaps = []
+        for name, mapping in self._maps:
+            pickledMaps += [name, cPickle.dumps(mapping)]
+
+        return (fromPickle, (className, props, self._maps, self._slices))

SimPEG/Utils/codeutils.py

@@ -65,8 +65,8 @@ def setKwargs(obj, ignore=[],  **kwargs):
         else:
             raise Exception('%s attr is not recognized' % attr)
 
-    hook(obj,hook, silent=True)
-    hook(obj,setKwargs, silent=True)
+    # hook(obj,hook, silent=True)
+    # hook(obj,setKwargs, silent=True)
 
 def printTitles(obj, printers, name='Print Titles', pad=''):
     titles = ''

docs/examples/EM_FDEM_SusEffects.rst

@@ -1,41 +0,0 @@
-.. _examples_EM_FDEM_SusEffects:
-
-.. --------------------------------- ..
-..                                   ..
-..    THIS FILE IS AUTO GENEREATED   ..
-..                                   ..
-..    SimPEG/Examples/__init__.py    ..
-..                                   ..
-.. --------------------------------- ..
-
-
-EM: FDEM: Effects of susceptibility
-===================================
-
-When airborne freqeuncy domain EM (AFEM) survey is flown over
-the earth including significantly susceptible bodies (magnetite-rich rocks),
-negative data is often observed in the real part of the lowest frequency
-(e.g. Dighem system 900 Hz). This phenomenon mostly based upon magnetization
-occurs due to a susceptible body when the magnetic field is applied.
-
-To clarify what is happening in the earth when we are exciting the earth with
-a loop source in the frequency domain we run three forward modelling:
-
-    - F[:math:`\sigma`, :math:`\mu`]: Anomalous conductivity and susceptibility
-    - F[:math:`\sigma`, :math:`\mu_0`]: Anomalous conductivity
-    - F[:math:`\sigma_{air}`, :math:`\mu_0`]: primary field
-
-We plot vector magnetic fields in the earth. For secondary fields we provide
-F[:math:`\sigma`, :math:`\mu`]-F[:math:`\sigma`, :math:`\mu_0`]. Following
-figure show both real and parts.
-
-
-
-.. plot::
-
-    from SimPEG import Examples
-    Examples.EM_FDEM_SusEffects.run()
-
-.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_SusEffects.py
-    :language: python
-    :linenos:

tests/base/test_PropMaps.py

@@ -1,6 +1,7 @@
 import unittest
 from SimPEG import *
 from scipy.constants import mu_0
+import cPickle
 
 
 class MyPropMap(Maps.PropMap):
@@ -28,9 +29,10 @@ class TestPropMaps(unittest.TestCase):
         PM3 = MyPropMap({'maps':[('sigma', expMap)], 'slices':{'sigma':slice(0,3)}})
 
         for PM in [PM1,PM2,PM3]:
+            PM = cPickle.loads( cPickle.dumps(PM) )
             assert PM.defaultInvProp == 'sigma'
             assert PM.sigmaMap is not None
-            assert PM.sigmaMap is expMap
+            assert PM.sigmaMap.__class__ is expMap.__class__
             assert PM.sigmaIndex == slice(0,3)
             assert getattr(PM, 'sigma', None) is None
             assert PM.muMap is None
@@ -52,7 +54,7 @@ class TestPropMaps(unittest.TestCase):
             assert m.muDeriv is None
 
             assert np.all(m.sigmaModel == np.r_[1.,2,3])
-            assert m.sigmaMap is expMap
+            assert m.sigmaMap.__class__ is expMap.__class__
             assert np.all(m.sigma == np.exp(np.r_[1.,2,3]))
             assert m.sigmaDeriv is not None