From 5239e57c69d175bde44f1f6d0b3a6835450429d4 Mon Sep 17 00:00:00 2001
From: rowanc1 <rowan@visiblegeology.com>
Date: Wed, 12 Feb 2014 16:59:56 -0800
Subject: [PATCH] testing the sensitivities

---
 simpegEM/TDEM/BaseTDEM.py                |  2 +
 simpegEM/TDEM/TDEM_b.py                  | 17 --------
 simpegEM/Tests/test_forward_EMproblem.py | 53 +++++++++++++++++++++++-
 3 files changed, 54 insertions(+), 18 deletions(-)

diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py
index 180f7b27..4e9f905a 100644
--- a/simpegEM/TDEM/BaseTDEM.py
+++ b/simpegEM/TDEM/BaseTDEM.py
@@ -135,6 +135,8 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
         self.makeMassMatrices(m)
 
         F = self.getInitialFields()
+        #TODO: Split next code to forward and adjoint.
+        # fields would call forward
         dtFact = None
         for tInd, t in enumerate(self.times):
             dt = self.getDt(tInd)
diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py
index 2b4f9832..fa007377 100644
--- a/simpegEM/TDEM/TDEM_b.py
+++ b/simpegEM/TDEM/TDEM_b.py
@@ -131,23 +131,6 @@ if __name__ == '__main__':
     prb.pair(dat)
     sigma = np.random.rand(mesh.nCz)
 
-    f = FieldsTDEM(prb.mesh, 1, prb.times.size, 'b')
-    for i in range(f.nTimes):
-        f.set_b(np.zeros((mesh.nF, 1)), i)
-        f.set_e(np.random.rand(mesh.nE, 1), i)
-
-    Ahf = prb.AhVec(sigma, f)
-    f_test = prb.solveAh(sigma, Ahf)
-
-    print np.linalg.norm(f.fieldVec() - f_test.fieldVec())
-
-    e0 = f.get_e(0)
-    e1 = f_test.get_e(0)
-    b0 = f.get_b(0)
-    b1 = f_test.get_b(0)
-    plt.semilogy(np.abs(e0))
-    plt.semilogy(np.abs(e1),'r')
-    plt.show()
 
 
 
diff --git a/simpegEM/Tests/test_forward_EMproblem.py b/simpegEM/Tests/test_forward_EMproblem.py
index 2eadab2f..7ef0dadd 100644
--- a/simpegEM/Tests/test_forward_EMproblem.py
+++ b/simpegEM/Tests/test_forward_EMproblem.py
@@ -3,6 +3,7 @@ from SimPEG import *
 import simpegEM as EM
 from scipy.constants import mu_0
 from simpegEM.Utils.Ana import hzAnalyticDipoleT
+import matplotlib.pyplot as plt
 
 
 class TDEM_bTests(unittest.TestCase):
@@ -66,6 +67,8 @@ class TDEM_bDerivTests(unittest.TestCase):
        self.sigma = np.ones(mesh.nCz)*1e-8
        self.sigma[mesh.vectorCCz<0] = 0.1
        self.prb.pair(self.dat)
+       self.mesh = mesh
+
 
     def test_AhVec(self):
         """
@@ -83,7 +86,7 @@ class TDEM_bDerivTests(unittest.TestCase):
             self.assertTrue(np.linalg.norm(Ahu.get_e(i))/np.linalg.norm(u.get_e(i)) < 1.e-2)
 
     def test_AhVecVSMat_OneTS(self):
-        
+
         self.prb.setTimes([1e-5], [1])
 
         sigma = np.ones(self.prb.mesh.nCz)*1e-8
@@ -165,6 +168,54 @@ class TDEM_bDerivTests(unittest.TestCase):
         # Assuming that the gradient is exact to machine precision
         self.assertTrue(b<1e-16)
 
+    def test_Deriv_dUdM(self):
+
+        prb = self.prb
+        prb.setTimes([1e-5, 1e-4, 1e-3], [10, 10, 10])
+        mesh = self.mesh
+        sigma = self.sigma
+
+        d_sig = sigma.copy() #np.random.rand(mesh.nCz)
+        d_sig[d_sig==1e-8] = 0
+
+        num = 10
+        error = np.zeros(num)
+        order = 0
+        hv = np.logspace(-1.2,-3, num)
+        for i, h in enumerate(hv):
+            f = prb.fields(sigma)
+            fstep = prb.fields(sigma + h*d_sig)
+            dcdm = prb.G(sigma, h*d_sig, u=f) # TODO: make negative!?!?
+            dudm = prb.solveAh(sigma, dcdm)
+
+            linear = np.linalg.norm(f.fieldVec() - fstep.fieldVec())
+            quad = np.linalg.norm(f.fieldVec() - fstep.fieldVec() - dudm.fieldVec())
+            error[i] = quad
+            if i > 0:
+                order = np.log(error[i]/error[i-1])/np.log(hv[i]/hv[i-1])
+
+            # print np.log(linearB/quadB)/np.log(h)
+            print h, linear, quad, order
+
+        self.assertTrue(order > 1.8)
+
+    def test_Deriv_J(self):
+
+        prb = self.prb
+        prb.setTimes([1e-5, 1e-4, 1e-3], [10, 10, 10])
+        mesh = self.mesh
+        sigma = self.sigma
+
+        d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
+        d_sig[d_sig==1e-8] = 0
+
+
+        derChk = lambda m: [prb.data.dpred(m), lambda mx: -prb.J(sigma, mx)]
+        passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=2, eps=1e-20)
+        self.assertTrue(passed)
+
+
+
 
 if __name__ == '__main__':
     unittest.main()