Start on the global problem class.

SimPEG/EM/FDEM/SurveyFDEM.py

@@ -106,7 +106,7 @@ class Survey(SimPEG.Survey.BaseSurvey):
         SimPEG.Survey.BaseSurvey.__init__(self, **kwargs)
 
         _freqDict = {}
-        for src in srcList:
+        for src in self.srcList:
             if src.freq not in _freqDict:
                 _freqDict[src.freq] = []
             _freqDict[src.freq] += [src]

SimPEG/Exceptions.py

@@ -0,0 +1,13 @@
+
+
+class SimPEGException(Exception):
+
+    def __init__(self, reason=''):
+        self.reason = reason
+
+    def __str__(self):
+        return '%s: %s' %(self.__class__.__name__,  self.reason)
+
+
+class PairingException(SimPEGException):
+    pass

SimPEG/Maps.py

@@ -113,7 +113,7 @@ class IdentityMap(object):
             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!')
+        raise Exception('Unrecognized data type to multiply. Try a map or a numpy.ndarray! Not a %s'%type(val))
 
     def __str__(self):
         return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1])
@@ -475,7 +475,7 @@ class ActiveCells(IdentityMap):
         else:
             self.valInactive = valInactive.copy()
         self.valInactive[self.indActive] = 0
-        
+
         inds = np.nonzero(self.indActive)[0]
         self.P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(self.nC, self.nP))
 
@@ -708,7 +708,7 @@ class PolyMap(IdentityMap):
         Parameterize the model space using a polynomials in a wholespace.
 
         ..math::
-            
+
             y = \mathbf{V} c
 
         Define the model as:
@@ -752,10 +752,10 @@ class PolyMap(IdentityMap):
             else:
                 raise(Exception("Input for normal = X or Y or Z"))
         #3D
-        elif self.mesh.dim == 3: 
+        elif self.mesh.dim == 3:
             X = self.mesh.gridCC[:,0]
-            Y = self.mesh.gridCC[:,1]            
-            Z = self.mesh.gridCC[:,2]            
+            Y = self.mesh.gridCC[:,1]
+            Z = self.mesh.gridCC[:,2]
             if self.normal =='X':
                 f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
             elif self.normal =='Y':
@@ -766,43 +766,43 @@ class PolyMap(IdentityMap):
                 raise(Exception("Input for normal = X or Y or Z"))
         else:
             raise(Exception("Only supports 2D"))
-                    
+
 
         return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5)
-    
+
     def deriv(self, m):
         alpha = self.slope
         sig1,sig2, c = m[0],m[1],m[2:]
         if self.logSigma:
             sig1, sig2 = np.exp(sig1), np.exp(sig2)
         #2D
-        if self.mesh.dim == 2:            
+        if self.mesh.dim == 2:
             X = self.mesh.gridCC[:,0]
             Y = self.mesh.gridCC[:,1]
 
             if self.normal =='X':
                 f = polynomial.polyval(Y, c) - X
-                V = polynomial.polyvander(Y, len(c)-1)        
+                V = polynomial.polyvander(Y, len(c)-1)
             elif self.normal =='Y':
                 f = polynomial.polyval(X, c) - Y
-                V = polynomial.polyvander(X, len(c)-1)        
+                V = polynomial.polyvander(X, len(c)-1)
             else:
-                raise(Exception("Input for normal = X or Y or Z"))            
+                raise(Exception("Input for normal = X or Y or Z"))
         #3D
-        elif self.mesh.dim == 3: 
+        elif self.mesh.dim == 3:
             X = self.mesh.gridCC[:,0]
             Y = self.mesh.gridCC[:,1]
             Z = self.mesh.gridCC[:,2]
 
             if self.normal =='X':
                 f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
-                V = polynomial.polyvander2d(Y, Z, self.order)        
+                V = polynomial.polyvander2d(Y, Z, self.order)
             elif self.normal =='Y':
                 f = polynomial.polyval2d(X, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - Y
-                V = polynomial.polyvander2d(X, Z, self.order)        
+                V = polynomial.polyvander2d(X, Z, self.order)
             elif self.normal =='Z':
                 f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z
-                V = polynomial.polyvander2d(X, Y, self.order)        
+                V = polynomial.polyvander2d(X, Y, self.order)
             else:
                 raise(Exception("Input for normal = X or Y or Z"))
 
@@ -815,16 +815,16 @@ class PolyMap(IdentityMap):
 
         g3 = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*V
 
-        return sp.csr_matrix(np.c_[g1,g2,g3])  
+        return sp.csr_matrix(np.c_[g1,g2,g3])
 
 class SplineMap(IdentityMap):
 
     """SplineMap
 
-        Parameterize the boundary of two geological units using a spline interpolation 
+        Parameterize the boundary of two geological units using a spline interpolation
 
         ..math::
-            
+
             g = f(x)-y
 
         Define the model as:
@@ -849,7 +849,7 @@ class SplineMap(IdentityMap):
     def nP(self):
         if self.mesh.dim == 2:
             return np.size(self.pts)+2
-        elif self.mesh.dim == 3: 
+        elif self.mesh.dim == 3:
             return np.size(self.pts)*2+2
         else:
             raise(Exception("Only supports 2D and 3D"))
@@ -866,28 +866,28 @@ class SplineMap(IdentityMap):
             X = self.mesh.gridCC[:,0]
             Y = self.mesh.gridCC[:,1]
             self.spl = UnivariateSpline(self.pts, c, k=self.order, s=0)
-            if self.normal =='X':                
+            if self.normal =='X':
                 f = self.spl(Y) - X
             elif self.normal =='Y':
                 f = self.spl(X) - Y
             else:
                 raise(Exception("Input for normal = X or Y or Z"))
 
-        # 3D: 
-        # Comments: 
+        # 3D:
+        # Comments:
         # Make two spline functions and link them using linear interpolation.
         # This is not quite direct extension of 2D to 3D case
         # Using 2D interpolation  is possible
 
-        elif self.mesh.dim == 3: 
+        elif self.mesh.dim == 3:
             X = self.mesh.gridCC[:,0]
-            Y = self.mesh.gridCC[:,1]            
+            Y = self.mesh.gridCC[:,1]
             Z = self.mesh.gridCC[:,2]
 
-            npts = np.size(self.pts)            
+            npts = np.size(self.pts)
             if np.mod(c.size, 2):
                 raise(Exception("Put even points!"))
-            
+
             self.spl = {"splb":UnivariateSpline(self.pts, c[:npts], k=self.order, s=0),
                         "splt":UnivariateSpline(self.pts, c[npts:], k=self.order, s=0)}
 
@@ -902,7 +902,7 @@ class SplineMap(IdentityMap):
                 raise(Exception("Input for normal = X or Y or Z"))
         else:
             raise(Exception("Only supports 2D and 3D"))
-                    
+
 
         return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5)
 
@@ -912,7 +912,7 @@ class SplineMap(IdentityMap):
         if self.logSigma:
             sig1, sig2 = np.exp(sig1), np.exp(sig2)
         #2D
-        if self.mesh.dim == 2:            
+        if self.mesh.dim == 2:
             X = self.mesh.gridCC[:,0]
             Y = self.mesh.gridCC[:,1]
 
@@ -921,9 +921,9 @@ class SplineMap(IdentityMap):
             elif self.normal =='Y':
                 f = self.spl(X) - Y
             else:
-                raise(Exception("Input for normal = X or Y or Z"))            
+                raise(Exception("Input for normal = X or Y or Z"))
         #3D
-        elif self.mesh.dim == 3: 
+        elif self.mesh.dim == 3:
             X = self.mesh.gridCC[:,0]
             Y = self.mesh.gridCC[:,1]
             Z = self.mesh.gridCC[:,2]
@@ -931,7 +931,7 @@ class SplineMap(IdentityMap):
                 zb = self.ptsv[0]
                 zt = self.ptsv[1]
                 flines = (self.spl["splt"](Y)-self.spl["splb"](Y))*(Z-zb)/(zt-zb) + self.spl["splb"](Y)
-                f = flines - X            
+                f = flines - X
             # elif self.normal =='Y':
             # elif self.normal =='Z':
             else:
@@ -944,7 +944,7 @@ class SplineMap(IdentityMap):
             g1 = -(np.arctan(alpha*f)/np.pi + 0.5) + 1.0
             g2 = (np.arctan(alpha*f)/np.pi + 0.5)
 
-        
+
         if self.mesh.dim ==2:
             g3 = np.zeros((self.mesh.nC, self.npts))
             if self.normal =='Y':
@@ -958,7 +958,7 @@ class SplineMap(IdentityMap):
                     cb = c.copy()
                     dy = self.mesh.hy[ind]*1.5
                     ca[i] = ctemp+dy
-                    cb[i] = ctemp-dy                
+                    cb[i] = ctemp-dy
                     spla = UnivariateSpline(self.pts, ca, k=self.order, s=0)
                     splb = UnivariateSpline(self.pts, cb, k=self.order, s=0)
                     fderiv = (spla(X)-splb(X))/(2*dy)
@@ -968,7 +968,7 @@ class SplineMap(IdentityMap):
             g3 = np.zeros((self.mesh.nC, self.npts*2))
             if self.normal =='X':
                 # Here we use perturbation to compute sensitivity
-                for i in range(self.npts*2):                    
+                for i in range(self.npts*2):
                     ctemp = c[i]
                     ind = np.argmin(abs(self.mesh.vectorCCy-ctemp))
                     ca = c.copy()
@@ -982,20 +982,20 @@ class SplineMap(IdentityMap):
                         splbb = UnivariateSpline(self.pts, cb[:self.npts], k=self.order, s=0)
                         flinesa = (self.spl["splt"](Y)-splba(Y))*(Z-zb)/(zt-zb) + splba(Y) - X
                         flinesb = (self.spl["splt"](Y)-splbb(Y))*(Z-zb)/(zt-zb) + splbb(Y) - X
-                    #treat top boundary                        
+                    #treat top boundary
                     else:
                         splta = UnivariateSpline(self.pts, ca[self.npts:], k=self.order, s=0)
                         spltb = UnivariateSpline(self.pts, ca[self.npts:], k=self.order, s=0)
                         flinesa = (self.spl["splt"](Y)-splta(Y))*(Z-zb)/(zt-zb) + splta(Y) - X
-                        flinesb = (self.spl["splt"](Y)-spltb(Y))*(Z-zb)/(zt-zb) + spltb(Y) - X                        
-                    fderiv = (flinesa-flinesb)/(2*dy)                                
+                        flinesb = (self.spl["splt"](Y)-spltb(Y))*(Z-zb)/(zt-zb) + spltb(Y) - X
+                    fderiv = (flinesa-flinesb)/(2*dy)
                     g3[:,i] = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*fderiv
         else :
             raise(Exception("Not Implemented for Y and Z, your turn :)"))
-        return sp.csr_matrix(np.c_[g1,g2,g3])  
+        return sp.csr_matrix(np.c_[g1,g2,g3])
+
 
 
 
 
 
-            

SimPEG/Problem.py

@@ -1,6 +1,6 @@
 import Utils, Survey, Models, numpy as np, scipy.sparse as sp
 Solver = Utils.SolverUtils.Solver
-import Maps, Mesh
+import Maps, Mesh, Exceptions
 from Fields import Fields, TimeFields
 
 class BaseProblem(object):
@@ -18,10 +18,14 @@ class BaseProblem(object):
     Solver = Solver   #: A SimPEG Solver class.
     solverOpts = {}   #: Sovler options as a kwarg dict
 
-    mesh    = None    #: A SimPEG.Mesh instance.
-
     PropMap = None    #: A SimPEG PropertyMap class.
 
+    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)
+
     @property
     def mapping(self):
         "A SimPEG.Map instance or a property map is PropMap is not None"
@@ -32,13 +36,8 @@ class BaseProblem(object):
             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._propMapMapping = val
+            self._mapping = self.PropMap(val)
 
     @property
     def survey(self):
@@ -47,13 +46,22 @@ class BaseProblem(object):
         """
         return getattr(self, '_survey', None)
 
-    def pair(self, d):
+    def pair(self, survey):
         """Bind a survey to this problem instance using pointers."""
-        assert isinstance(d, self.surveyPair), "Data object must be an instance of a %s class."%(self.surveyPair.__name__)
-        if d.ispaired:
+        assert isinstance(survey, self.surveyPair), "Survey must be an instance of a %s class."%(self.surveyPair.__name__)
+        if survey.ispaired:
             raise Exception("The survey object is already paired to a problem. Use survey.unpair()")
-        self._survey = d
-        d._prob = self
+        try:
+            self._survey = survey
+            self._validatePairing()
+        except Exceptions.PairingException, e:
+            self._survey = None
+            raise e
+        survey._prob = self
+
+    def _validatePairing(self):
+        """Called when the pair is done, raise a SimPEG.Exceptions.PairingException if unsuccessful"""
+        pass
 
     def unpair(self):
         """Unbind a survey from this problem instance."""
@@ -214,4 +222,177 @@ class BaseTimeProblem(BaseProblem):
             del self._timeMesh
 
 
+class GlobalProblem(BaseProblem):
+    """
+
+        The GlobalProblem allows you to run a whole bunch of SubProblems,
+        potentially in parallel, potentially of different meshes.
+
+        This is handy for working with lots of sources,
+
+    """
+
+    surveyKwargs = {}
+    probKwargs   = {}
+
+    def __init__(self, SubProblem, globalMesh, mapping=None, **kwargs):
+
+        # assert isclass??(SubProblem, BaseProblem), "SubProblem must be a SimPEG.Problem.BaseProblem object."
+        self.surveyPair = SubProblem.surveyPair
+        self.PropMap = SubProblem.PropMap
+        self.mapPair = SubProblem.mapPair
+        self.SubProblem = SubProblem
+
+        Utils.setKwargs(self, **kwargs)
+        assert isinstance(globalMesh, Mesh.BaseMesh), "globalMesh must be a SimPEG.Mesh object."
+        self.globalMesh = globalMesh
+        self.mapping = mapping or Maps.IdentityMap(mesh)
+
+    @property
+    def groups(self):
+        """
+            List of lists/integers to say how the sources are grouped.
+
+            e.g.
+
+            survey.srcList = [s0,s1,s2,s3,s4]
+            groups = [ [0,4], [1,3], 2 ]
+        """
+        if getattr(self, '_groups', None) is None:
+            if not self.ispaired: return None
+            self._groups = range(self.survey.nSrc)
+        return self._groups
+    @groups.setter
+    def groups(self, val):
+        assert type(val) is list, 'This should be an list of groups'
+        if self.ispaired:
+            for g in val:
+                assert type(g) in [int, list], 'Must be an integer or a list'
+                if type(g) is int:
+                    assert g >= 0 and g < self.survey.nSrc, '%d is outside the number of sources in the surveys list'%g
+                if type(g) is list:
+                    for sg in g:
+                        assert type(g) is int, 'Must be an integer or a list'
+                        assert g >= 0 and g < self.survey.nSrc, '%d is outside the number of sources in the surveys list'%g
+            assert len(val) == len(self.survey.srcList), 'The groups must be the same length as the srcList in the survey'
+        self._groups = val
+        self._nGroups = None
+
+    @property
+    def meshes(self):
+        if getattr(self, '_meshes', None) is None:
+            if not self.ispaired: return None
+            self._meshes = [self.globalMesh]*self.nGroups
+        return self._meshes
+    @meshes.setter
+    def meshes(self, val):
+        assert type(val) is list
+        if self.ispaired:
+            assert len(val) == self.nGroups
+        self._meshes = val
+
+    @property
+    def nGroups(self):
+        if getattr(self, '_groups', None) is None:
+            return None
+        return len(self.groups)
+
+    def _validatePairing(self):
+        try:
+            self.groups = self.groups # check the assumptions for the grouping
+        except Exception, e:
+            raise Exceptions.PairingException(reason='The grouping does not match the survey')
+        if self.nGroups is not len(self.meshes):
+            raise Exceptions.PairingException(reason='The meshes are not the the same length as the number of groups')
+
+    def getSubProblem(self, ind):
+
+        assert self.ispaired, 'You must be paired to a survey'
+        assert type(ind) in [int,long] and ind >= 0 and ind < self.nGroups, 'ind must be an index into the group list'
+
+        subMesh = self.meshes[ind]
+        subMap  = Maps.IdentityMap(subMesh) # this is probably a mesh2mesh mapping?
+
+        if self.PropMap is None:
+            prob = self.SubProblem(subMesh, mapping=subMap * self.mapping, **self.probKwargs)
+        else:
+            prob = self.SubProblem(subMesh, mapping=subMap * self._propMapMapping, **self.probKwargs)
+
+        survey = self.survey.__class__(srcList=self.survey.srcList[self.groups[ind]], **self.surveyKwargs)
+        prob.pair(survey)
+
+        return prob
+
+
+if __name__ == '__main__':
+
+
+    from SimPEG import *
+    from SimPEG import EM
+    from scipy.constants import mu_0
+    from pymatsolver import MumpsSolver
+
+    cs = 10.
+    ncx, ncy, ncz = 10, 10, 10
+    npad = 4
+    freq = 1e2
+
+    hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)]
+    hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
+    hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
+    mesh = Mesh.TensorMesh([hx,hy,hz], 'CCC')
+
+    mapping = Maps.ExpMap(mesh)
+
+    x = np.linspace(-10,10,5)
+    XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
+    rxList = EM.FDEM.Rx(XYZ, 'exi')
+    Src0 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq)
+    Src1 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq)
+
+
+    prb0 = EM.FDEM.Problem_b(mesh, mapping=mapping, Solver=MumpsSolver)
+    survey = EM.FDEM.Survey([Src0])
+    prb0.pair(survey)
+    prb1 = EM.FDEM.Problem_b(mesh, mapping=mapping, Solver=MumpsSolver)
+    survey = EM.FDEM.Survey([Src1])
+    prb1.pair(survey)
+
+
+
+    sig = 1e-1
+    sigma = np.ones(mesh.nC)*sig
+    sigma[mesh.gridCC[:,2] > 0] = 1e-8
+    m = np.log(sigma)
+
+    GP = GlobalProblem(EM.FDEM.Problem_b, mesh, mapping=mapping, meshes=[mesh,mesh])
+    survey = EM.FDEM.Survey([Src0, Src1])
+    GP.pair(survey)
+
+    gp1 = GP.getSubProblem(0)
+    gp1.Solver = MumpsSolver
+
+    pu = prb0.fields(m)
+    gpu = gp1.fields(m)
+
+    bfz = mesh.r(pu[Src0, 'b'],'F','Fz','M')
+    bfz = mesh.r(gpu[Src0, 'b'],'F','Fz','M')
+    x = np.linspace(-55,55,12)
+    XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
+    P = mesh.getInterpolationMat(XYZ, 'Fz')
+
+    # an = EM.Analytics.FDEM.hzAnalyticDipoleF(x, Src0.freq, sig)
+
+    # diff = np.log10(np.abs(P*np.imag(pu[Src0, 'b']) - mu_0*np.imag(an)))
+    # diff = np.log10(np.abs(P*np.imag(gpu[Src0, 'b']) - mu_0*np.imag(an)))
+
+    import matplotlib.pyplot as plt
+    plt.plot(x,np.log10(np.abs(P*np.imag(pu[Src0, 'b']))), 'r-s')
+    plt.plot(x,np.log10(np.abs(P*np.imag(gpu[Src0, 'b']))), 'b')
+    # plt.plot(x,np.log10(np.abs(mu_0*np.imag(an))), 'r')
+    # plt.plot(x,diff,'g')
+    plt.show()
+
+
+
 

SimPEG/Survey.py

@@ -223,6 +223,8 @@ class BaseSurvey(object):
 
     @srcList.setter
     def srcList(self, value):
+        if isinstance(value, self.srcPair):
+            value = [value]
         assert type(value) is list, 'srcList must be a list'
         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'