Merge branch 'dcip/ref' of https://github.com/simpeg/simpeg into dcip/ref

Conflicts:
	SimPEG/EM/Static/DC/ProblemDC.py

Confused about ...  this line

self.mesh.getFaceInnerProduct(self.curModel.rho)(u)

SimPEG/EM/Base.py

@@ -190,7 +190,13 @@ class BaseEMProblem(Problem.BaseProblem):
         """
             Derivative of :code:`MfRhoI` with respect to the model.
         """
-        return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv
+
+        dMfRhoI_dI = -self.MfRhoI**2
+        dMf_drho = self.mesh.getFaceInnerProduct(self.curModel.rho)(u)
+        drho_dm = self.curModel.rhoDeriv
+        return dMfRhoI_dI * ( dMf_drho * ( drho_dm))
+
+        # return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv
 
 class BaseEMSurvey(Survey.BaseSurvey):
 

SimPEG/EM/Static/DC/FieldsDC.py

@@ -1,5 +1,6 @@
 import SimPEG
 from SimPEG.Utils import Identity, Zero
+import numpy as np
 
 class Fields(SimPEG.Problem.Fields):
     knownFields = {}
@@ -10,8 +11,8 @@ class Fields(SimPEG.Problem.Fields):
             raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0])
 
         if adjoint:
-            return self._phiDeriv_u(src, v, adjoint), self._phiDeriv_m(src, v, adjoint)
-        return np.array(self._phiDeriv_u(src, du_dm_v, adjoint) + self._phiDeriv_m(src, v, adjoint), dtype = complex)
+            return self._phiDeriv_u(src, v, adjoint=adjoint), self._phiDeriv_m(src, v, adjoint=adjoint)
+        return np.array(self._phiDeriv_u(src, du_dm_v, adjoint) + self._phiDeriv_m(src, v, adjoint), dtype = float)
 
     def _eDeriv(self, src, du_dm_v, v, adjoint=False):
         if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None:
@@ -19,7 +20,7 @@ class Fields(SimPEG.Problem.Fields):
 
         if adjoint:
             return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint)
-        return np.array(self._eDeriv_u(src, du_dm_v, adjoint) + self._eDeriv_m(src, v, adjoint), dtype = complex)
+        return np.array(self._eDeriv_u(src, du_dm_v, adjoint) + self._eDeriv_m(src, v, adjoint), dtype = float)
 
     def _jDeriv(self, src, du_dm_v, v, adjoint=False):
         if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None:
@@ -27,7 +28,7 @@ class Fields(SimPEG.Problem.Fields):
 
         if adjoint:
             return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint)
-        return np.array(self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint), dtype = complex)
+        return np.array(self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint), dtype = float)
 
 
 class Fields_CC(Fields):
@@ -57,10 +58,10 @@ class Fields_CC(Fields):
     def _phi(self, phiSolution, srcList):
         return phiSolution
 
-    def _phiDeriv_u():
+    def _phiDeriv_u(self, src, v, adjoint = False):
         return Identity()
 
-    def _phiDeriv_m():
+    def _phiDeriv_m(self, src, v, adjoint = False):
         return Zero()
 
     def _j(self, phiSolution, srcList):
@@ -96,10 +97,10 @@ class Fields_N(Fields):
     def _phi(self, phiSolution, srcList):
         return phiSolution
 
-    def _phiDeriv_u():
+    def _phiDeriv_u(self, src, v, adjoint = False):
         return Identity()
 
-    def _phiDeriv_m():
+    def _phiDeriv_m(self, src, v, adjoint = False):
         return Zero()
 
     def _j(self, phiSolution, srcList):

SimPEG/EM/Static/DC/ProblemDC.py

@@ -4,6 +4,7 @@ from SurveyDC import Survey
 from FieldsDC import Fields, Fields_CC, Fields_N
 from SimPEG.Utils import sdiag
 import numpy as np
+from SimPEG.Utils import Zero
 
 class BaseDCProblem(BaseEMProblem):
 
@@ -22,7 +23,30 @@ class BaseDCProblem(BaseEMProblem):
         return f
 
     def Jvec(self, m, v, f=None):
-        raise NotImplementedError
+
+        if f is None:
+            f = self.fields(m)
+
+        self.curModel = m
+
+        Jv = self.dataPair(self.survey) #same size as the data
+
+        A = self.getA()
+        Ainv = self.Solver(A, **self.solverOpts)
+
+        for src in self.survey.srcList:
+            u_src = f[src, self._solutionType] # solution vector
+            dA_dm_v = self.getADeriv(u_src, v)
+            dRHS_dm_v = self.getRHSDeriv(src, v)
+            print type(dA_dm_v + dRHS_dm_v), (dA_dm_v + dRHS_dm_v).shape
+            du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v )
+
+            for rx in src.rxList:
+                df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
+                Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
+            Ainv.clean()
+        return Utils.mkvc(Jv)
 
     def Jtvec(self, m, v, f=None):
         raise NotImplementedError
@@ -126,23 +150,32 @@ class Problem3D_CC(BaseDCProblem):
 
         """
 
+        V = self.Vol
+        D = V * self.mesh.faceDiv
         # TODO: this won't work for full anisotropy
-        # V = self.Vol
-        # Div = V*self.mesh.faceDiv
         MfRhoI = self.MfRhoI
-        A = self.Div * MfRhoI * self.Div.T
+        A = D * MfRhoI * D.T
+
+        # I think we should deprecate this for DC problem.
         # if self._makeASymmetric is True:
         #     return V.T * A
         return A
 
-    def getADeriv(self, u, v, adoint=False):
-        """
+    def getADeriv(self, u, v, adjoint= False):
 
-        Product of the derivative of our system matrix with respect to the model and a vector
+        V = self.Vol
+        D = V * self.mesh.faceDiv
+        MfRhoIDeriv = self.MfRhoIDeriv
 
-        """
-        return Div*self.MfRhoIDeriv(Div.T*u)
+        if adjoint:
+            # if self._makeASymmetric is True:
+            #     v = V * v
+            return D * MfRhoIDeriv(D * v)
 
+        # I think we should deprecate this for DC problem.
+        # if self._makeASymmetric is True:
+        #     return V.T * ( D * ( MfRhoIDeriv( D.T * ( V * u ) ) * v ) )
+        return D * (MfRhoIDeriv( D.T * u ) * v)
 
     def getRHS(self):
         """
@@ -152,16 +185,21 @@ class Problem3D_CC(BaseDCProblem):
         """
 
         RHS = self.getSourceTerm()
+
+        # I think we should deprecate this for DC problem.
         # if self._makeASymmetric is True:
         #     return self.Vol.T * RHS
+
         return RHS
 
     def getRHSDeriv(self, src, v, adjoint=False):
         """
         Derivative of the right hand side with respect to the model
         """
-        qDeriv = src.evalDeriv(self, adjoint=adjoint)
-        return qDeriv
+        # TODO: add qDeriv for RHS depending on m
+        # qDeriv = src.evalDeriv(self, adjoint=adjoint)
+        # return qDeriv
+        return Zero()
 
 
 

SimPEG/EM/Static/DC/RxDC.py

@@ -36,6 +36,10 @@ class BaseRx(SimPEG.Survey.BaseRx):
         P = self.getP(mesh, self.projGLoc(f))
         return P*f[src, self.projField]
 
+    def evalDeriv(self, src, mesh, f, v, adjoint=False):
+        P = self.getP(mesh, self.projGLoc(f))
+        return P*v
+
 # DC.Rx.Dipole(locs)
 class Dipole(BaseRx):
 
@@ -50,6 +54,10 @@ class Dipole(BaseRx):
         """Number of data in the receiver."""
         return self.locs[0].shape[0]
 
+        # Not sure why ...
+        # return int(self.locs[0].size / 2)
+
+
     def getP(self, mesh, Gloc):
         if mesh in self._Ps:
             return self._Ps[mesh]
@@ -63,7 +71,6 @@ class Dipole(BaseRx):
 
         return P
 
-
 # class Pole(BaseRx):
 
 

SimPEG/EM/Static/DC/Utils.py

@@ -0,0 +1,38 @@
+import numpy as np
+
+def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
+
+    import SimPEG.EM.Static.DC as DC
+
+    elocs = np.arange(0,aSpacing*nElecs,aSpacing)
+    elocs -= (nElecs*aSpacing - aSpacing)/2
+    space = 1
+    WENNER = np.zeros((0,),dtype=int)
+    for ii in range(nElecs):
+        for jj in range(nElecs):
+            test = np.r_[jj,jj+space,jj+space*2,jj+space*3]
+            if np.any(test >= nElecs):
+                break
+            WENNER = np.r_[WENNER, test]
+        space += 1
+    WENNER = WENNER.reshape((-1,4))
+
+
+    if plotIt:
+        for i, s in enumerate('rbkg'):
+            plt.plot(elocs[WENNER[:,i]],s+'.')
+        plt.show()
+
+    # Create sources and receivers
+    i = 0
+    if in2D:
+        getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0]
+    else:
+        getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0, 0]
+    srcList = []
+    for i in range(WENNER.shape[0]):
+        rx = DC.Rx.Dipole(getLoc(i,1).reshape([1,-1]),getLoc(i,2).reshape([1,-1]))
+        src = DC.Src.Dipole([rx], getLoc(i,0),getLoc(i,3))
+        srcList += [src]
+
+    return srcList

SimPEG/EM/Static/DC/__init__.py

@@ -3,3 +3,4 @@ from SurveyDC import Survey
 import SrcDC as Src #Pole
 import RxDC as Rx
 from FieldsDC import Fields_CC
+import Utils

tests/em/static/__init__.py

@@ -0,0 +1,12 @@
+import os
+import glob
+import unittest
+
+if __name__ == '__main__':
+    test_file_strings = glob.glob('test_*.py')
+    module_strings = [str[0:len(str)-3] for str in test_file_strings]
+    suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
+              in module_strings]
+    testSuite = unittest.TestSuite(suites)
+
+    unittest.TextTestRunner(verbosity=2).run(testSuite)

tests/em/static/test_DC.py

@@ -0,0 +1,77 @@
+import unittest
+from SimPEG import *
+import SimPEG.EM.Static.DC as DC
+
+
+class DCProblemTests(unittest.TestCase):
+
+    def setUp(self):
+
+        aSpacing=2.5
+        nElecs=10
+
+        surveySize = nElecs*aSpacing - aSpacing
+        cs = surveySize/nElecs/4
+
+        mesh = Mesh.TensorMesh([
+                [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)],
+                [(cs,3, -1.3),(cs,3,1.3)],
+        #         [(cs,5, -1.3),(cs,10)]
+            ],'CN')
+
+        srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
+        survey = DC.Survey(srcList)
+        problem = DC.Problem3D_CC(mesh, mapping=[('rho', Maps.IdentityMap(mesh))])
+        problem.pair(survey)
+
+        mSynth = np.ones(mesh.nC)
+        survey.makeSyntheticData(mSynth)
+
+        # Now set up the problem to do some minimization
+        dmis = DataMisfit.l2_DataMisfit(survey)
+        reg = Regularization.Tikhonov(mesh)
+        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
+        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
+        inv = Inversion.BaseInversion(invProb)
+
+        self.inv = inv
+        self.reg = reg
+        self.p =     problem
+        self.mesh = mesh
+        self.m0 = mSynth
+        self.survey = survey
+        self.dmis = dmis
+
+    def test_misfit(self):
+        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
+        self.assertTrue(passed)
+
+    # def test_adjoint(self):
+    #     # Adjoint Test
+    #     u = np.random.rand(self.mesh.nC*self.survey.nSrc)
+    #     v = np.random.rand(self.mesh.nC)
+    #     w = np.random.rand(self.survey.dobs.shape[0])
+    #     wtJv = w.dot(self.p.Jvec(self.m0, v))
+    #     vtJtw = v.dot(self.p.Jtvec(self.m0, w))
+    #     passed = np.abs(wtJv - vtJtw) < 1e-10
+    #     print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
+    #     self.assertTrue(passed)
+
+    # def test_dataObj(self):
+    #     derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
+    #     passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
+    #     self.assertTrue(passed)
+
+
+    # def test_massMatrices(self):
+    #     Gu = np.random.rand(self.mesh.nF)
+    #     def derChk(m):
+    #         self.p.curModel = m
+    #         return [self.p.Msig * Gu, self.p.dMdsig(Gu)]
+    #     passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
+    #     self.assertTrue(passed)
+
+
+if __name__ == '__main__':
+    unittest.main()