From d233e40a95e9c9f7e75b9d7880d02825ac72e1a2 Mon Sep 17 00:00:00 2001 From: GudniRos Date: Wed, 24 Jun 2015 19:55:10 -0700 Subject: [PATCH] Jvec adjoint test not working, but all other derivatives adjoint tests are working. --- notebooks/Derivative test MT1D.ipynb | 302 ++++++++++++++++++++++++--- simpegMT/BaseMT.py | 6 +- simpegMT/SurveyMT.py | 4 +- 3 files changed, 275 insertions(+), 37 deletions(-) diff --git a/notebooks/Derivative test MT1D.ipynb b/notebooks/Derivative test MT1D.ipynb index 9e8e36c3..cde75b75 100644 --- a/notebooks/Derivative test MT1D.ipynb +++ b/notebooks/Derivative test MT1D.ipynb @@ -237,13 +237,13 @@ "==================== checkDerivative ====================\n", "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", "---------------------------------------------------------\n", - " 0 1.00e-01 9.779e-05 3.219e-06 nan\n", - " 1 1.00e-02 1.007e-05 3.145e-08 2.010\n", - " 2 1.00e-03 1.010e-06 3.137e-10 2.001\n", - " 3 1.00e-04 1.010e-07 3.136e-12 2.000\n", - " 4 1.00e-05 1.010e-08 3.136e-14 2.000\n", + " 0 1.00e-01 1.076e-05 3.966e-07 nan\n", + " 1 1.00e-02 1.040e-06 3.875e-09 2.010\n", + " 2 1.00e-03 1.037e-07 3.866e-11 2.001\n", + " 3 1.00e-04 1.037e-08 3.865e-13 2.000\n", + " 4 1.00e-05 1.037e-09 3.865e-15 2.000\n", "========================= PASS! =========================\n", - "Yay passed!\n", + "The test be workin!\n", "\n" ] }, @@ -277,7 +277,7 @@ { "data": { "text/plain": [ - "array([ 0.00017028])" + "array([ -2.07544500e-05])" ] }, "execution_count": 8, @@ -299,7 +299,7 @@ { "data": { "text/plain": [ - "array([[ 0.0014539]])" + "array([[-0.00014517]])" ] }, "execution_count": 9, @@ -338,13 +338,15 @@ "---------------------------------------------------------\n", "Project at freq: 1.000e+02\n", "Project at freq: 1.000e+02\n", - " 0 1.00e-01 8.466e-06 5.557e-08 nan\n", + " 0 1.00e-01 1.142e-06 2.072e-08 nan\n", "Project at freq: 1.000e+02\n", - " 1 1.00e-02 8.416e-07 5.519e-10 2.003\n", + " 1 1.00e-02 1.161e-07 2.039e-10 2.007\n", "Project at freq: 1.000e+02\n", - " 2 1.00e-03 8.411e-08 5.515e-12 2.000\n", + " 2 1.00e-03 1.162e-08 2.036e-12 2.001\n", + "Project at freq: 1.000e+02\n", + " 3 1.00e-04 1.163e-09 2.036e-14 2.000\n", "========================= PASS! =========================\n", - "Not just a pretty face Gudni\n", + "Testing is important.\n", "\n" ] }, @@ -374,7 +376,7 @@ "def fun(x):\n", " \n", " return survey.dpred(x), lambda x: problem.Jvec(m0, x)\n", - "simpeg.Tests.checkDerivative(fun, m0, num=3, plotIt=False)" + "simpeg.Tests.checkDerivative(fun, m0, num=4, plotIt=False)" ] }, { @@ -403,12 +405,195 @@ "metadata": { "collapsed": false }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Adjoint e formulation - projectFieldsDeriv\n", + "-9.0191451349e-06 -9.0191451349e-06 3.38813178902e-21 1e-09 True\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Run a test\n", + "TOL = 1e-4\n", + "FLR = 1e-20\n", + "\n", + "def projectFieldsAdjointTest(fdemType, comp):\n", + " print 'Adjoint %s formulation - %s' % (fdemType, comp)\n", + "\n", + " m = np.log(np.ones(problem.mesh.nC)*0.01)\n", + " if True:\n", + " m = m + np.random.randn(problem.mesh.nC)*0.01*1e-1 \n", + "\n", + " u = problem.fields(m)\n", + " v = np.random.randn(1)#+np.random.randn(1)*1j\n", + " # print prb.PropMap.PropModel.nP\n", + " w = np.random.randn(m1d.nN)#+np.random.randn(m1d.nN)*1j\n", + "\n", + " vJw = v.dot(rx.projectFieldsDeriv(src,m1d,f0,w))\n", + " wJtv = w.dot(rx.projectFieldsDeriv(src,m1d,f0,v,adjoint=True))\n", + " tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) \n", + " print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol\n", + " return np.abs(vJw - wJtv) < tol\n", + "projectFieldsAdjointTest('e','projectFieldsDeriv')" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Adjoint test e formulation - getADeriv_m\n", + "(508851.474801-63591.6147417j) (508851.474801-63591.6147417j) (-1.7462298274e-10-5.82076609135e-11j) 10.0 True\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Run a test\n", + "TOL = 1e-4\n", + "FLR = 1e-20\n", + "\n", + "def getADeriv_mAdjointTest():\n", + " print 'Adjoint test e formulation - getADeriv_m' \n", + "\n", + " m = np.log(np.ones(problem.mesh.nC)*0.01)\n", + " if True:\n", + " m = m + np.random.randn(problem.mesh.nC)*0.01*1e-1 \n", + "\n", + " u = problem.fields(m)\n", + " v = np.random.randn(m1d.nN)#+np.random.randn(1)*1j\n", + " # print prb.PropMap.PropModel.nP\n", + " w = np.random.randn(m1d.nC)#+np.random.randn(m1d.nN)*1j\n", + "\n", + " vJw = v.dot(problem.getADeriv_m(freq,u0,w))\n", + " wJtv = w.dot(problem.getADeriv_m(freq,u0,v,adjoint=True))\n", + " tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) \n", + " print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol\n", + " return np.abs(vJw - wJtv) < tol\n", + "getADeriv_mAdjointTest()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Adjoint test e formulation - getRHSDeriv_m\n", + "(-29427.4754714-17701.4445852j) (-29427.4754714-17701.4445852j) (2.91038304567e-11+0j) 1.0 True\n" + ] + }, + { + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Run a test\n", + "TOL = 1e-4\n", + "FLR = 1e-20\n", + "\n", + "def getRHSDeriv_mAdjointTest():\n", + " print 'Adjoint test e formulation - getRHSDeriv_m'\n", + "\n", + " m = np.log(np.ones(problem.mesh.nC)*0.01)\n", + " if True:\n", + " m = m + np.random.randn(problem.mesh.nC)*0.01*1e-1 \n", + "\n", + " u = problem.fields(m)\n", + " v = np.random.randn(m1d.nN)#+np.random.randn(1)*1j\n", + " # print prb.PropMap.PropModel.nP\n", + " w = np.random.randn(m1d.nC)#+np.random.randn(m1d.nN)*1j\n", + "\n", + " vJw = v.dot(problem.getRHSDeriv_m(freq,w))\n", + " wJtv = w.dot(problem.getRHSDeriv_m(freq,v,adjoint=True))\n", + " tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) \n", + " print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol\n", + " return np.abs(vJw - wJtv) < tol\n", + "getRHSDeriv_mAdjointTest( )" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2\n" + ] + } + ], + "source": [ + "simpeg.mkvc(np.random.randn(survey.nD)+np.random.randn(survey.nD)*1j,2)\n", + "\n", + "print survey.nD" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": { + "collapsed": false + }, "outputs": [], "source": [ "TOL = 1e-4\n", "FLR = 1e-20\n", "\n", - "def adjointTest(fdemType, comp):\n", + "def JvecAdjointTest(fdemType, comp):\n", " print 'Adjoint %s formulation - %s' % (fdemType, comp)\n", "\n", " m = np.log(np.ones(problem.mesh.nC)*0.01)\n", @@ -430,7 +615,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 18, "metadata": { "collapsed": false }, @@ -440,7 +625,7 @@ "output_type": "stream", "text": [ "Adjoint E formulation - e\n", - "-4.65700174631e-05 3.27055690471e-05 -7.92755865102e-05 1e-08 False\n" + "-1.7867035849e-05 1.42750293876e-05 -3.21420652366e-05 1e-08 False\n" ] }, { @@ -449,25 +634,63 @@ "False" ] }, - "execution_count": 14, + "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ - "adjointTest('E','e')" + "JvecAdjointTest('E','e')" ] }, { "cell_type": "code", - "execution_count": 17, + "execution_count": 19, "metadata": { "collapsed": false, "scrolled": true }, - "outputs": [], + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "ERROR: No traceback has been produced, nothing to debug.\n" + ] + } + ], "source": [ - "# %debug" + "%debug" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "collapsed": false, + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Project at freq: 1.000e+02\n" + ] + }, + { + "data": { + "text/plain": [ + "array([ 9.80523303e-06, -1.98372645e-03])" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "survey.dpred(sigma)" ] }, { @@ -481,20 +704,23 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 21, "metadata": { "collapsed": false }, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], + "outputs": [ + { + "data": { + "text/plain": [ + "<180x180 sparse matrix of type ''\n", + "\twith 180 stored elements in Compressed Sparse Row format>" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "problem.mesh.getEdgeInnerProductDeriv(problem.curModel.sigma)(u0[1::])" ] @@ -514,6 +740,18 @@ "display_name": "Python 2", "language": "python", "name": "python2" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.9" } }, "nbformat": 4, diff --git a/simpegMT/BaseMT.py b/simpegMT/BaseMT.py index bc630a5d..f1092f03 100644 --- a/simpegMT/BaseMT.py +++ b/simpegMT/BaseMT.py @@ -106,15 +106,15 @@ class BaseMTProblem(BaseFDEMProblem): dRHS_dmT = self.getRHSDeriv_m(freq, dA_duIT, adjoint=True) # Make du_dmT if dRHS_dmT is None: - du_dmT = - dA_dmT + du_dmT = -dA_dmT else: du_dmT = -dA_dmT + dRHS_dmT # Select the correct component real_or_imag = rx.projComp if real_or_imag == 'real': - Jtv += du_dmT.real + Jtv += du_dmT.real elif real_or_imag == 'imag': - Jtv += - du_dmT.real + Jtv += -du_dmT.real else: raise Exception('Must be real or imag') diff --git a/simpegMT/SurveyMT.py b/simpegMT/SurveyMT.py index 3256015b..c440a812 100644 --- a/simpegMT/SurveyMT.py +++ b/simpegMT/SurveyMT.py @@ -161,8 +161,8 @@ class RxMT(Survey.BaseRx): # ex = Pex*mkvc(f[src,'e_1d'],2) # bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0 dP_deTv = mkvc(Pex.T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*v,2) - db_duv = Pbx.T/mu_0*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2)).T*v - dP_dbTv = -mkvc(f._bDeriv_u(src,db_duv,adjoint=True),2) + db_duv = -Pbx.T/mu_0*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2)).T*v + dP_dbTv = mkvc(f._bDeriv_u(src,db_duv,adjoint=True),2) PDeriv_complex = np.sum(np.hstack((dP_deTv,dP_dbTv)),1) # raise Exception('Debug error') elif self.projType is 'Z3D':