diff --git a/simpegEM/TDEM/SurveyTDEM.py b/simpegEM/TDEM/SurveyTDEM.py index 722846f5..0e04f10a 100644 --- a/simpegEM/TDEM/SurveyTDEM.py +++ b/simpegEM/TDEM/SurveyTDEM.py @@ -60,8 +60,8 @@ class FieldsTDEM(Survey.TimeFields): e = self[:,'e',i+1] else: e = np.zeros(nE if nTx == 1 else (nE, nTx)) - u = np.r_[u, b, e] - return u + u = np.concatenate((u, b, e)) + return Utils.mkvc(u) class TxTDEM(Survey.BaseTx): rxPair = RxTDEM diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py index 01c6a50b..9622af14 100644 --- a/simpegEM/TDEM/TDEM_b.py +++ b/simpegEM/TDEM/TDEM_b.py @@ -64,7 +64,7 @@ class ProblemTDEM_b(BaseTDEMProblem): p = self.survey.projectFieldsDeriv(u, v=v, adjoint=True) y = self.solveAht(m, p) w = self.Gtvec(m, y, u) - return - w + return - mkvc(w) def Gvec(self, m, vec, u=None): """ @@ -107,10 +107,13 @@ class ProblemTDEM_b(BaseTDEMProblem): if u is None: u = self.fields(m) nTx, nE = self.survey.nTx, self.mesh.nE - tmp = np.zeros(nE if nTx == 1 else (nE,nTx)) + tmp = np.zeros(nE) # Here we can do internal multiplications of Gt*v and then multiply by MsigDeriv.T in one go. for i in range(1,self.nT+1): - tmp += vec[:,'e',i]*u[:,'e',i] + vu = vec[:,'e',i]*u[:,'e',i] + if nTx > 1: + vu = vu.sum(axis=1) + tmp += vu curModel = self.mapping.transform(m) p = -mkvc(self.mapping.transformDeriv(m).T*self.mesh.getEdgeInnerProductDeriv(curModel).T*tmp) @@ -137,6 +140,12 @@ class ProblemTDEM_b(BaseTDEMProblem): def solveAht(self, m, p): + # Mini Example: + # + # nT = 3, len(times) == 4, fields stored in F[:,:,1:4] + # + # 0 is held for initial conditions (this shifts the storage by +1) + # ^ # fLoc 0 1 2 3 # |-----|-----|-----| # tInd 0 1 2 / / diff --git a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py index 841cef25..36176944 100644 --- a/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py +++ b/simpegEM/Tests/test_TDEM_b_DerivAdjoint.py @@ -197,7 +197,7 @@ class TDEM_bDerivTests(unittest.TestCase): for i in range(prb.nT): f[:,'b',i] = np.random.rand(mesh.nF, 1) f[:,'e',i] = np.random.rand(mesh.nE, 1) - d_vec = np.random.rand(survey.nD, survey.nTx).flatten() + d_vec = np.random.rand(survey.nD) d = Survey.Data(survey,v=d_vec) # Check that d.T*Q*f = f.T*Q.T*d diff --git a/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py b/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py new file mode 100644 index 00000000..477952e1 --- /dev/null +++ b/simpegEM/Tests/test_TDEM_b_MultiTx_DerivAdjoint.py @@ -0,0 +1,150 @@ +import unittest +from SimPEG import * +import simpegEM as EM + +plotIt = False + +class TDEM_bDerivTests(unittest.TestCase): + + def setUp(self): + + cs = 5. + ncx = 20 + ncy = 6 + npad = 20 + hx = [(cs,ncx), (cs,npad,1.3)] + hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)] + mesh = Mesh.CylMesh([hx,1,hy], '00C') + + active = mesh.vectorCCz<0. + activeMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz) + mapping = Maps.ComboMap(mesh, + [Maps.ExpMap, Maps.Vertical1DMap, activeMap]) + + rxOffset = 40. + rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz') + tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx]) + rx2 = EM.TDEM.RxTDEM(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), 'bz') + tx2 = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx2]) + + survey = EM.TDEM.SurveyTDEM([tx,tx2]) + + self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping) + # self.prb.timeSteps = [1e-5] + self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)] + # self.prb.timeSteps = [(1e-05, 100)] + + self.sigma = np.ones(mesh.nCz)*1e-8 + self.sigma[mesh.vectorCCz<0] = 1e-1 + self.sigma = np.log(self.sigma[active]) + + self.prb.pair(survey) + self.mesh = mesh + + def test_DerivG(self): + """ + Test the derivative of c with respect to sigma + """ + + # Random model and perturbation + sigma = np.random.rand(self.prb.mapping.nP) + + f = self.prb.fields(sigma) + dm = 1000*np.random.rand(self.prb.mapping.nP) + h = 0.01 + + derChk = lambda m: [self.prb.AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()] + print '\ntest_DerivG' + passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) + self.assertTrue(passed) + + def test_Deriv_dUdM(self): + + prb = self.prb + prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] + mesh = self.mesh + sigma = self.sigma + + dm = 10*np.random.rand(prb.mapping.nP) + f = prb.fields(sigma) + + derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()] + print '\n' + print 'test_Deriv_dUdM' + passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20) + self.assertTrue(passed) + + def test_Deriv_J(self): + + prb = self.prb + prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)] + mesh = self.mesh + sigma = self.sigma + + # d_sig = 0.8*sigma #np.random.rand(mesh.nCz) + d_sig = 10*np.random.rand(prb.mapping.nP) + + + derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)] + print '\n' + print 'test_Deriv_J' + passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20) + self.assertTrue(passed) + + def test_projectAdjoint(self): + prb = self.prb + survey = prb.survey + mesh = self.mesh + + # Generate random fields and data + f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + for i in range(prb.nT): + f[:,'b',i] = np.random.rand(mesh.nF, 1) + f[:,'e',i] = np.random.rand(mesh.nE, 1) + d_vec = np.random.rand(survey.nD) + d = Survey.Data(survey,v=d_vec) + + # Check that d.T*Q*f = f.T*Q.T*d + V1 = d_vec.dot(survey.projectFieldsDeriv(None, v=f).tovec()) + V2 = f.tovec().dot(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec()) + + self.assertLess((V1-V2)/np.abs(V1), 1e-6) + + def test_adjointGvecVsGtvec(self): + mesh = self.mesh + prb = self.prb + + m = np.random.rand(prb.mapping.nP) + sigma = np.random.rand(prb.mapping.nP) + + u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + for i in range(1,prb.nT+1): + u[:,'b',i] = np.random.rand(mesh.nF, 2) + u[:,'e',i] = np.random.rand(mesh.nE, 2) + + v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey) + for i in range(1,prb.nT+1): + v[:,'b',i] = np.random.rand(mesh.nF, 2) + v[:,'e',i] = np.random.rand(mesh.nE, 2) + + V1 = m.dot(prb.Gtvec(sigma, v, u)) + V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec()) + self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) + + def test_adjointJvecVsJtvec(self): + mesh = self.mesh + prb = self.prb + sigma = self.sigma + + m = np.random.rand(prb.mapping.nP) + d = np.random.rand(prb.survey.nD) + + V1 = d.dot(prb.Jvec(sigma, m)) + V2 = m.dot(prb.Jtvec(sigma, d)) + print 'AdjointTest', V1, V2 + self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6) + + + +if __name__ == '__main__': + unittest.main()