From f2a8cf0a623bf2c48d44d2a69f9808f728e41629 Mon Sep 17 00:00:00 2001 From: GudniRos Date: Wed, 24 Jun 2015 11:33:14 -0700 Subject: [PATCH] Jvec working for MT1D, Jtvec getting close --- notebooks/Derivative test MT1D.ipynb | 328 ++++++++++++++------------- simpegMT/BaseMT.py | 77 ++++--- simpegMT/FieldsMT.py | 33 ++- simpegMT/ProblemMT1D/Problems.py | 13 +- simpegMT/SurveyMT.py | 39 +++- 5 files changed, 265 insertions(+), 225 deletions(-) diff --git a/notebooks/Derivative test MT1D.ipynb b/notebooks/Derivative test MT1D.ipynb index 29e5ce9c..a0ffef2c 100644 --- a/notebooks/Derivative test MT1D.ipynb +++ b/notebooks/Derivative test MT1D.ipynb @@ -159,7 +159,7 @@ "source": [ "As matrices the formulas above can be written as\n", "\\begin{align}\n", - "\\left[ \\frac{\\partial P(f(m))}{\\partial u} v \\right] = diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right] [P_e v] - diag[P_e u] diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right]^T diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right] \\left[ P_b \\frac{d f_b}{du}(v) \\frac{1}{mu_0} \\right]\n", + "\\left[ \\frac{\\partial P(f(m))}{\\partial u} v \\right] = \\left[ diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right] [P_e v] , diag[P_e u] diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right]^T diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right] \\left[ P_b \\frac{1}{mu_0} \\frac{d f_b}{du}(v) \\right] \\right]\n", "\\end{align}\n", "\n" ] @@ -170,7 +170,7 @@ "source": [ "The adjoint problem is done simliarly\n", "\\begin{align}\n", - "\\left[ \\frac{\\partial P(f(m))}{\\partial u} v \\right]^T = [P_e v]^T diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right]^T - \\left[ P_b \\frac{d f_b}{du}(v) \\frac{1}{mu_0} \\right]^T diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right] diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right]^T diag \\left[ P_e u \\right]^T\n", + "\\left[ \\frac{\\partial P(f(m))}{\\partial u} v \\right]^T = [P_e]^T diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right]^T v - \\left[ P_b \\frac{d f_b}{du} \\frac{1}{mu_0} \\right]^T diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right] diag \\left[ \\frac{1}{\\left(P_b \\frac{1}{mu_0} f_b(u)\\right)} \\right]^T diag \\left[ P_e u \\right]^T v\n", "\\end{align}\n" ] }, @@ -210,135 +210,21 @@ "src = survey.srcList[0]\n", "rx = src.rxList[0]\n", "v = np.random.randn(m1d.nN)\n", + "v0 = np.random.randn(m1d.nF+m1d.nE)\n", "u0 = np.random.randn(m1d.nN)+np.random.randn(m1d.nN)*1j\n", "f0 = problem.fieldsPair(m1d,survey)\n", "f0[src,'e_1dSolution'] = u0\n", "# f0[src,'b_1d'] = -1/(1j*simpegem.Utils.EMUtils.omega(src.freq))*m1d.nodalGrad*u0" ] }, - { - "cell_type": "code", - "execution_count": 7, - "metadata": { - "collapsed": false, - "scrolled": true - }, - "outputs": [ - { - "data": { - "text/plain": [ - "array([ -5.21576299e-01-0.28596337j, 1.33775588e+00-1.08159078j,\n", - " 1.02311300e+00-0.45033307j, -1.14140017e+00-0.53814712j,\n", - " -1.51095317e+00-0.0225078j , 2.29842631e-01-0.33435588j,\n", - " -6.13742564e-01+0.14777572j, 1.74709600e+00+0.07973247j,\n", - " -6.07949930e-01+0.53467919j, 7.44109732e-01-1.06383585j,\n", - " -6.22557065e-01-0.6604577j , 5.54530896e-01+0.49763768j,\n", - " -1.42655640e+00+1.03849457j, -8.79197413e-01-0.05384536j,\n", - " -5.91554616e-01-0.29187425j, -2.09918231e+00+0.65127732j,\n", - " -2.79577518e-01+0.08460458j, -8.49733033e-01-1.15162916j,\n", - " -5.61159802e-01+0.40323403j, -1.32856823e-01+0.83145628j,\n", - " 1.38730097e-01-0.64695336j, 4.74841673e-01+2.16259344j,\n", - " 5.20314820e-02+1.39337638j, 2.40669484e+00+1.44814886j,\n", - " 1.02271608e+00-0.64912151j, 6.52379785e-01+0.59394264j,\n", - " 2.50946314e-01-1.71467789j, -3.35061639e-01+1.80721483j,\n", - " 1.87006132e+00-0.18721108j, 1.85379814e+00+0.06414066j,\n", - " 7.66449796e-01+0.28757766j, 3.77939332e-01+0.5861157j ,\n", - " -2.01974930e+00+2.36530034j, -5.70969153e-01+0.16116942j,\n", - " 3.37395529e-01+0.76006325j, -1.12084482e-01-2.02591559j,\n", - " 1.36070180e-01-0.88380389j, 1.99378191e+00-0.50681261j,\n", - " 1.50662216e-01+1.76352214j, -4.33494912e-01-0.80584325j,\n", - " -5.30448843e-01+0.99832665j, -1.52568430e+00+0.46127907j,\n", - " -4.78999736e-01-0.52250698j, 5.07078441e-01+0.02687913j,\n", - " 1.52990246e+00+0.09587891j, 9.10105334e-01-0.22098742j,\n", - " 6.63562676e-01-0.45367061j, -1.08906319e+00-1.81627004j,\n", - " 3.41582990e-02+1.71027836j, -2.22532053e-01+1.06004075j,\n", - " 1.74917677e+00-1.41735414j, 1.06979706e+00-1.89094749j,\n", - " -3.35092655e-01+0.39943895j, 1.05206441e+00+0.80828213j,\n", - " -6.87168211e-01+1.72353937j, 1.40047301e+00-1.3727236j ,\n", - " -1.03495040e-01+0.26680256j, -1.12873662e+00-0.56805852j,\n", - " -8.57584940e-01+1.85418143j, 1.03116279e+00+0.06711151j,\n", - " 5.46075841e-01+1.07590359j, 2.67048944e-01+0.75324202j,\n", - " 9.60209324e-01-0.09519331j, -1.36504332e-01-1.35387957j,\n", - " -3.38824018e-01+0.70215028j, -4.82487190e-02-0.24027277j,\n", - " 1.55731408e+00-2.38970922j, -3.35634102e-03-1.11778716j,\n", - " 1.05996566e+00-1.54847977j, 1.37371396e+00+0.21763618j,\n", - " -2.34257640e+00+1.16351971j, -5.41859146e-01+1.0902067j ,\n", - " -1.79555485e+00+0.06162846j, 7.46946726e-01-0.87760098j,\n", - " -9.37453865e-02+0.26210336j, -6.00849545e-01+0.49005865j,\n", - " 2.47173181e-04-0.41345406j, 1.00840211e+00+1.23027562j,\n", - " 8.39952932e-01+0.43605496j, 8.65132399e-01+1.07984908j,\n", - " 2.04841885e+00-0.0851883j , 1.44345869e-01+0.4539564j ,\n", - " 9.19226689e-01-0.37858608j, -1.51448024e-01-1.2974253j ,\n", - " 5.32408479e-02-1.60301977j, -1.84229724e+00-0.15278222j,\n", - " -2.13350910e+00+1.87102448j, 9.03468928e-01+0.23498138j,\n", - " -9.26803345e-01+0.51213564j, 7.94552525e-01+0.54895852j,\n", - " -9.77798824e-01-0.49901395j, -6.05296444e-01-0.03209891j,\n", - " 1.24765034e+00-0.56443723j, -1.90295464e+00-0.21190697j,\n", - " 1.44761031e+00+0.14009266j, -9.23615842e-01+1.93414685j,\n", - " -2.89396133e-01+0.47920215j, 1.22494934e-01+1.33535626j,\n", - " -6.00157469e-01+0.61558122j, -1.40975425e-01+1.61693392j,\n", - " 7.08238379e-01+0.536495j , 5.90147090e-01+0.60041602j,\n", - " -7.09067864e-01+0.08489417j, -1.24030970e+00-0.73986538j,\n", - " -1.03488057e+00-0.09736817j, 2.89640823e-01-0.41632583j,\n", - " -4.94975749e-02-0.61519627j, 5.41641220e-01+0.77100395j,\n", - " 7.22415122e-01+1.59030908j, -5.13304114e-01-0.09690537j,\n", - " 5.41180474e-01+2.6467645j , -8.05664508e-01+1.44276035j,\n", - " -4.06768654e-01+1.58421981j, 1.33368455e+00+0.63576588j,\n", - " 4.29526056e-02+0.16497284j, 5.92169353e-01-0.76120082j,\n", - " 1.84734942e+00-0.29548982j, 3.25330682e-01+0.84903188j,\n", - " 1.08709120e-01-1.73236832j, -5.01403552e-01-0.083525j ,\n", - " 3.43892536e-01+0.3488206j , -2.42917431e+00-0.25699335j,\n", - " 5.53836789e-01-0.45887772j, -6.75365309e-01+0.0848306j ,\n", - " 4.78327738e-01-0.77164497j, -4.05543707e-01-0.04087569j,\n", - " 6.48445605e-01+1.00477682j, 4.50922906e-02-0.09451425j,\n", - " -2.97801954e-01-0.40438134j, -2.47619828e+00+0.25935539j,\n", - " 1.95520198e+00+0.93190891j, 7.52207479e-01+0.60075967j,\n", - " -8.95365798e-01+1.22241712j, -4.65516271e-01-0.10833261j,\n", - " -1.55272666e-01-1.20074828j, 5.78006258e-01+0.06705178j,\n", - " 9.35873450e-01-1.244475j , 3.08977268e-03-1.57055076j,\n", - " -1.69903140e+00-1.00183406j, -5.73059946e-01-1.96431137j,\n", - " 3.80255478e-01-0.58071467j, -3.27026898e-01-0.44312791j,\n", - " 2.66653209e-01+0.20442809j, -7.26239447e-02-0.73470193j,\n", - " -1.27190710e+00+0.6681173j , -1.66240293e+00-0.68067832j,\n", - " -1.30894453e-01+0.4142016j , 1.17940864e-01-0.29430427j,\n", - " 5.40779365e-01-0.1570116j , 4.59411844e-01-1.33449408j,\n", - " -2.86538969e-01-0.46958326j, -2.91593793e-02-1.25227084j,\n", - " 2.45440216e-02+0.07172228j, -8.17237335e-02-0.71911812j,\n", - " 1.74142448e+00-0.26649839j, -4.31150193e-01-0.79560216j,\n", - " -7.04116352e-01-1.53969088j, -4.88139638e-01-0.59075677j,\n", - " -1.61059835e-01-0.62180822j, -1.75918449e-01-0.31350392j,\n", - " -2.28918438e-01-1.25884104j, -3.94756005e-01-0.08464047j,\n", - " 3.23000023e-01+0.87737735j, 1.05597723e-02-0.74774027j,\n", - " -2.25509504e+00+2.81870339j, -9.95775953e-01-1.21962899j,\n", - " -8.59196834e-01+1.4407182j , -6.38861601e-01+0.67564119j,\n", - " 1.99028382e-01+2.32353483j, -1.15758435e+00-0.70250997j,\n", - " -9.18642125e-01-0.72232602j, 1.00415977e+00+0.40547797j,\n", - " -2.61807248e-01-1.27673902j, -1.39773535e+00+1.04765733j,\n", - " 8.93091651e-01-0.67650682j, -1.26058855e+00+1.24753141j,\n", - " -4.41649966e-01-0.64575466j, -1.35393159e+00+0.42326246j,\n", - " 1.35067221e-01+1.44186272j, 3.95784662e-01-0.99344302j,\n", - " 3.22900419e-01+1.13931166j])" - ] - }, - "execution_count": 7, - "metadata": {}, - "output_type": "execute_result" - } - ], - "source": [ - "u0" - ] - }, { "cell_type": "raw", "metadata": {}, - "source": [ - "#rx.projectFieldsDeriv(src,m1d,fields,v)\n" - ] + "source": [] }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 7, "metadata": { "collapsed": false, "scrolled": true @@ -351,12 +237,13 @@ "==================== checkDerivative ====================\n", "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", "---------------------------------------------------------\n", - " 0 1.00e-01 5.703e-05 1.929e-06 nan\n", - " 1 1.00e-02 5.877e-06 1.954e-08 1.994\n", - " 2 1.00e-03 5.894e-07 1.956e-10 1.999\n", - " 3 1.00e-04 5.896e-08 1.957e-12 2.000\n", + " 0 1.00e-01 3.893e-05 5.627e-08 nan\n", + " 1 1.00e-02 3.898e-06 5.632e-10 2.000\n", + " 2 1.00e-03 3.898e-07 5.633e-12 2.000\n", + " 3 1.00e-04 3.898e-08 5.633e-14 2.000\n", + " 4 1.00e-05 3.898e-09 5.637e-16 2.000\n", "========================= PASS! =========================\n", - "Not just a pretty face Gudni\n", + "Yay passed!\n", "\n" ] }, @@ -366,7 +253,7 @@ "True" ] }, - "execution_count": 8, + "execution_count": 7, "metadata": {}, "output_type": "execute_result" } @@ -376,9 +263,30 @@ "def fun(u):\n", " f = problem.fieldsPair(m1d,survey)\n", " f[src,'e_1dSolution'] = u\n", - "# f[src,'b_1d'] = -(m1d.nodalGrad*u)/(1j*simpegem.Utils.EMUtils.omega(src.freq))\n", " return rx.projectFields(src,m1d,f), lambda t: rx.projectFieldsDeriv(src,m1d,f0,t)\n", - "simpeg.Tests.checkDerivative(fun,u0,num=4,plotIt=False)" + "simpeg.Tests.checkDerivative(fun,u0,num=5,plotIt=False)" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "array(-0.005683351098616839)" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "rx.projectFieldsDeriv(src,m1d,f0,u0)" ] }, { @@ -387,6 +295,28 @@ "metadata": { "collapsed": false }, + "outputs": [ + { + "data": { + "text/plain": [ + "array([[ 0.00374478]])" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "rx.projectFields(src,m1d,f0)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, "outputs": [], "source": [ "# Test the Jvec derivative." @@ -394,7 +324,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 11, "metadata": { "collapsed": false }, @@ -407,25 +337,26 @@ "iter h |ft-f0| |ft-f0-h*J0*dx| Order\n", "---------------------------------------------------------\n", "Project at freq: 1.000e+02\n", - "Project at freq: 1.000e+02\n" + "Project at freq: 1.000e+02\n", + " 0 1.00e-01 7.283e-06 9.326e-09 nan\n", + "Project at freq: 1.000e+02\n", + " 1 1.00e-02 7.290e-07 9.338e-11 1.999\n", + "Project at freq: 1.000e+02\n", + " 2 1.00e-03 7.290e-08 9.339e-13 2.000\n", + "========================= PASS! =========================\n", + "Testing is important.\n", + "\n" ] }, { - "ename": "ValueError", - "evalue": "dimension mismatch", - "output_type": "error", - "traceback": [ - "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", - "\u001b[1;31mValueError\u001b[0m Traceback (most recent call last)", - "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[0;32m 13\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 14\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0msurvey\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdpred\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;32mlambda\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m:\u001b[0m \u001b[0mproblem\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mJvec\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 15\u001b[1;33m \u001b[0msimpeg\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mTests\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcheckDerivative\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mu0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnum\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mplotIt\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", - "\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Tests/TestUtils.pyc\u001b[0m in \u001b[0;36mcheckDerivative\u001b[1;34m(fctn, x0, num, plotIt, dx, expectedOrder, tolerance, eps, ax)\u001b[0m\n\u001b[0;32m 256\u001b[0m \u001b[1;31m# 1st order Taylor\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 257\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0minspect\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0misfunction\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mJ0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 258\u001b[1;33m \u001b[0mE1\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0ml2norm\u001b[0m\u001b[1;33m(\u001b[0m \u001b[0mft\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mf0\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mh\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mi\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mJ0\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdx\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 259\u001b[0m \u001b[1;32melse\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 260\u001b[0m \u001b[1;31m# We assume it is a numpy.ndarray\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m(x)\u001b[0m\n\u001b[0;32m 12\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0mfun\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 13\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 14\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0msurvey\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdpred\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;32mlambda\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m:\u001b[0m \u001b[0mproblem\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mJvec\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 15\u001b[0m \u001b[0msimpeg\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mTests\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcheckDerivative\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mu0\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnum\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;36m3\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mplotIt\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/BaseMT.pyc\u001b[0m in \u001b[0;36mJvec\u001b[1;34m(self, m, v, f)\u001b[0m\n\u001b[0;32m 64\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 65\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 66\u001b[1;33m \u001b[0mJv\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mrx\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mP\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mdu_dm\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 67\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 68\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mJv\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/BaseMT.pyc\u001b[0m in \u001b[0;36m\u001b[1;34m(v)\u001b[0m\n\u001b[0;32m 61\u001b[0m \u001b[0mdu_dm\u001b[0m \u001b[1;33m+=\u001b[0m \u001b[0mdf_dm\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 62\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 63\u001b[1;33m \u001b[0mP\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;32mlambda\u001b[0m \u001b[0mv\u001b[0m\u001b[1;33m:\u001b[0m \u001b[0mrx\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mprojectFieldsDeriv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmesh\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mf\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mv\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;31m# wrt u, also have wrt m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 64\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 65\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/SurveyMT.py\u001b[0m in \u001b[0;36mprojectFieldsDeriv\u001b[1;34m(self, src, mesh, f, v, adjoint)\u001b[0m\n\u001b[0;32m 141\u001b[0m \u001b[1;31m# ex = Pex*mkvc(f[src,'e_1d'],2)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 142\u001b[0m \u001b[1;31m# bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 143\u001b[1;33m \u001b[0mderiv_complex\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1.\u001b[0m\u001b[1;33m/\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPbx\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'b_1d'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mmu_0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPex\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPex\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'e_1d'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1.\u001b[0m\u001b[1;33m/\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPbx\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'b_1d'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mmu_0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1.\u001b[0m\u001b[1;33m/\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPbx\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'b_1d'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mmu_0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPbx\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_bDeriv_u\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mmu_0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 144\u001b[0m \u001b[1;31m# elif self.projType is 'Z2D\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 145\u001b[0m \u001b[1;32melif\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mprojType\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;34m'Z3D'\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;32m/home/gudni/anaconda/lib/python2.7/site-packages/scipy/sparse/base.pyc\u001b[0m in \u001b[0;36m__mul__\u001b[1;34m(self, other)\u001b[0m\n\u001b[0;32m 325\u001b[0m \u001b[1;31m# dense row or column vector\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 326\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mother\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m \u001b[1;33m!=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mN\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mand\u001b[0m \u001b[0mother\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m \u001b[1;33m!=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mN\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 327\u001b[1;33m \u001b[1;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'dimension mismatch'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 328\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 329\u001b[0m \u001b[0mresult\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_mul_vector\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mravel\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mother\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[1;31mValueError\u001b[0m: dimension mismatch" - ] + "data": { + "text/plain": [ + "True" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" } ], "source": [ @@ -442,24 +373,64 @@ "# survey = prb.survey\n", "def fun(x):\n", " \n", - " return survey.dpred(m0), lambda x: problem.Jvec(m0, x)\n", - "simpeg.Tests.checkDerivative(fun, u0, num=3, plotIt=False)" + " return survey.dpred(x), lambda x: problem.Jvec(m0, x)\n", + "simpeg.Tests.checkDerivative(fun, m0, num=3, plotIt=False)" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 12, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "### Adjoint test" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We have \n", + "\\begin{align}\n", + "Jvec =&" + ] + }, + { + "cell_type": "code", + "execution_count": 13, "metadata": { "collapsed": false }, "outputs": [], "source": [ - "problem.getADeriv_m(freq,fields[src,'e_1dSolution'],v)" + "TOL = 1e-4\n", + "FLR = 1e-20\n", + "\n", + "def adjointTest(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", + "\n", + " v = np.random.rand(survey.nD)\n", + " # print prb.PropMap.PropModel.nP\n", + " w = np.random.rand(problem.mesh.nC)\n", + "\n", + " vJw = v.dot(problem.Jvec(m, w, u))\n", + " wJtv = w.dot(problem.Jtvec(m, v, u))\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" ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 14, "metadata": { "collapsed": false }, @@ -468,21 +439,53 @@ "name": "stdout", "output_type": "stream", "text": [ - "> \u001b[1;32m/home/gudni/anaconda/lib/python2.7/site-packages/scipy/sparse/base.py\u001b[0m(327)\u001b[0;36m__mul__\u001b[1;34m()\u001b[0m\n", - "\u001b[1;32m 326 \u001b[1;33m \u001b[1;32mif\u001b[0m \u001b[0mother\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m \u001b[1;33m!=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mN\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mand\u001b[0m \u001b[0mother\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m \u001b[1;33m!=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mN\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[0m\u001b[1;32m--> 327 \u001b[1;33m \u001b[1;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'dimension mismatch'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[0m\u001b[1;32m 328 \u001b[1;33m\u001b[1;33m\u001b[0m\u001b[0m\n", + "Adjoint E formulation - e\n" + ] + }, + { + "ename": "IndexError", + "evalue": "tuple index out of range", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mIndexError\u001b[0m Traceback (most recent call last)", + "\u001b[1;32m\u001b[0m in \u001b[0;36m\u001b[1;34m()\u001b[0m\n\u001b[1;32m----> 1\u001b[1;33m \u001b[0madjointTest\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'E'\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'e'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m", + "\u001b[1;32m\u001b[0m in \u001b[0;36madjointTest\u001b[1;34m(fdemType, comp)\u001b[0m\n\u001b[0;32m 16\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 17\u001b[0m \u001b[0mvJw\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mv\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mproblem\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mJvec\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mw\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mu\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 18\u001b[1;33m \u001b[0mwJtv\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mw\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mproblem\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mJtvec\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mm\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mv\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mu\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 19\u001b[0m \u001b[0mtol\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmax\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mTOL\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m10\u001b[0m\u001b[1;33m**\u001b[0m\u001b[0mint\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlog10\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvJw\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mFLR\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 20\u001b[0m \u001b[1;32mprint\u001b[0m \u001b[0mvJw\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mwJtv\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvJw\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mwJtv\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mtol\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvJw\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mwJtv\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m<\u001b[0m \u001b[0mtol\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/BaseMT.pyc\u001b[0m in \u001b[0;36mJtvec\u001b[1;34m(self, m, v, f)\u001b[0m\n\u001b[0;32m 102\u001b[0m \u001b[0mPTv\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mrx\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mprojectFieldsDeriv\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmesh\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mf\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mv\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mrx\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0madjoint\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mTrue\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;31m# wrt u, need possibility wrt m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 103\u001b[0m \u001b[1;31m# Get the\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 104\u001b[1;33m \u001b[0mdA_duIT\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mATinv\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mPTv\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 105\u001b[0m \u001b[0mdA_dmT\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetADeriv_m\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mu_src\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdA_duIT\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0madjoint\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 106\u001b[0m \u001b[0mdRHS_dmT\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetRHSDeriv_m\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdA_duIT\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0madjoint\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Utils/SolverUtils.pyc\u001b[0m in \u001b[0;36m__mul__\u001b[1;34m(self, b)\u001b[0m\n\u001b[0;32m 35\u001b[0m \u001b[1;32mraise\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m'Can only multiply by a numpy array.'\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 36\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 37\u001b[1;33m \u001b[1;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m1\u001b[0m \u001b[1;32mor\u001b[0m \u001b[0mb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 38\u001b[0m \u001b[0mb\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mflatten\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 39\u001b[0m \u001b[1;31m# Just one RHS\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[1;31mIndexError\u001b[0m: tuple index out of range" + ] + } + ], + "source": [ + "adjointTest('E','e')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": false, + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "> \u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpeg/SimPEG/Utils/SolverUtils.py\u001b[0m(37)\u001b[0;36m__mul__\u001b[1;34m()\u001b[0m\n", + "\u001b[1;32m 36 \u001b[1;33m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[0m\u001b[1;32m---> 37 \u001b[1;33m \u001b[1;32mif\u001b[0m \u001b[0mlen\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m1\u001b[0m \u001b[1;32mor\u001b[0m \u001b[0mb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;36m1\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;36m1\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[0m\u001b[1;32m 38 \u001b[1;33m \u001b[0mb\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mb\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mflatten\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[0m\n", "ipdb> u\n", - "> \u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/SurveyMT.py\u001b[0m(143)\u001b[0;36mprojectFieldsDeriv\u001b[1;34m()\u001b[0m\n", - "\u001b[1;32m 142 \u001b[1;33m \u001b[1;31m# bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[0m\u001b[1;32m--> 143 \u001b[1;33m \u001b[0mderiv_complex\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1.\u001b[0m\u001b[1;33m/\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPbx\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'b_1d'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mmu_0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPex\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;33m-\u001b[0m \u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPex\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'e_1d'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1.\u001b[0m\u001b[1;33m/\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPbx\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'b_1d'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mmu_0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mT\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mUtils\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0msdiag\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;36m1.\u001b[0m\u001b[1;33m/\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPbx\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mmkvc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;34m'b_1d'\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;36m2\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mmu_0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mPbx\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0mf\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_bDeriv_u\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[1;33m,\u001b[0m\u001b[0mv\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m/\u001b[0m\u001b[0mmu_0\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", - "\u001b[0m\u001b[1;32m 144 \u001b[1;33m \u001b[1;31m# elif self.projType is 'Z2D\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "> \u001b[1;32m/media/gudni/ExtraDrive1/Codes/python/simpegmt/simpegMT/BaseMT.py\u001b[0m(104)\u001b[0;36mJtvec\u001b[1;34m()\u001b[0m\n", + "\u001b[1;32m 103 \u001b[1;33m \u001b[1;31m# Get the\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[0m\u001b[1;32m--> 104 \u001b[1;33m \u001b[0mdA_duIT\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mATinv\u001b[0m \u001b[1;33m*\u001b[0m \u001b[0mPTv\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", + "\u001b[0m\u001b[1;32m 105 \u001b[1;33m \u001b[0mdA_dmT\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mgetADeriv_m\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mfreq\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mu_src\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mdA_duIT\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0madjoint\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mTrue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n", "\u001b[0m\n", - "ipdb> Pex*v\n", - "*** ValueError: dimension mismatch\n", - "ipdb> v.shape\n", - "(180,)\n" + "ipdb> PTv\n", + "array(-6.9926641874700175e-06)\n" ] } ], @@ -498,7 +501,8 @@ }, "outputs": [], "source": [ - "problem.getA" + "v = np.random.randn(1,1) + 1j* np.random.randn(1,1)\n", + "rx.projectFieldsDeriv(src, problem.mesh, fields, v, adjoint=True)" ] }, { diff --git a/simpegMT/BaseMT.py b/simpegMT/BaseMT.py index 31b28300..370f6476 100644 --- a/simpegMT/BaseMT.py +++ b/simpegMT/BaseMT.py @@ -2,7 +2,8 @@ from simpegEM.FDEM import BaseFDEMProblem from SurveyMT import SurveyMT from DataMT import DataMT from FieldsMT import FieldsMT -from SimPEG import SolverLU as SimpegSolver +from SimPEG import SolverLU as SimpegSolver, mkvc +import numpy as np class BaseMTProblem(BaseFDEMProblem): @@ -24,18 +25,32 @@ class BaseMTProblem(BaseFDEMProblem): # Might need to add more stuff here. def Jvec(self, m, v, f=None): + """ + Function to calculate the data sensitivities dD/dm times a vector. + + :param numpy.ndarray (nC, 1) - conductive model + :param numpy.ndarray (nC, 1) - random vector + :param MTfields object (optional) - MT fields object, if not given it is calculated + :rtype: MTdata object + :return: Data sensitivities wrt m + """ + + # Calculate the fields if f is None: f = self.fields(m) - + # Set current model self.curModel = m - + # Initiate the Jv object Jv = self.dataPair(self.survey) + # Loop all the frequenies for freq in self.survey.freqs: dA_du = self.getA(freq) # dA_duI = self.Solver(dA_du, **self.solverOpts) for src in self.survey.getSrcByFreq(freq): + # We need fDeriv_m = df/du*du/dm + df/dm + # Construct du/dm, it requires a solve ftype = self._fieldType + 'Solution' u_src = f[src, ftype] dA_dm = self.getADeriv_m(freq, u_src, v) @@ -44,28 +59,23 @@ class BaseMTProblem(BaseFDEMProblem): du_dm = dA_duI * ( - dA_dm ) else: du_dm = dA_duI * ( - dA_dm + dRHS_dm ) + # Calculate the projection derivatives for rx in src.rxList: - # df_duFun = u.deriv_u(rx.fieldsUsed, m) - if 'e' in self._fieldType: - projField = 'b' - elif 'b' in self._fieldType: - projField = 'e' - df_duFun = getattr(f, '_%sDeriv_u'%projField, None) - df_du = df_duFun(src, du_dm, adjoint=False) - if df_du is not None: - du_dm = df_du + # Get the stacked derivative + # df_duFun = getattr(f, '_fDeriv_u', None) + # df_dmFun = getattr(f, '_fDeriv_m', None) + # df_dm = df_dmFun(src,v,adjoint=False) + # if df_dm is None: + # fDeriv_m = df_duFun(src, du_dm, adjoint=False) + # else: + # fDeriv_m = df_duFun(src, du_dm, adjoint=False) + df_dm + # Not needed for now. Since PDeriv does this currently. - df_dmFun = getattr(f, '_%sDeriv_m'%projField, None) - df_dm = df_dmFun(src, v, adjoint=False) - if df_dm is not None: - du_dm += df_dm - - P = lambda v: rx.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m - - - Jv[src, rx] = P(du_dm) - - return Utils.mkvc(Jv) + # Get the projection derivative + PDeriv = lambda v: rx.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m + Jv[src, rx] = PDeriv(du_dm) + # Return the vectorized sensitivities + return mkvc(Jv) def Jtvec(self, m, v, f=None): if f is None: @@ -88,29 +98,18 @@ class BaseMTProblem(BaseFDEMProblem): u_src = f[src, ftype] for rx in src.rxList: + # Get the adjoint projectFieldsDeriv PTv = rx.projectFieldsDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt u, need possibility wrt m - - df_duTFun = getattr(f, '_%sDeriv_u'%rx.projField, None) - df_duT = df_duTFun(src, PTv, adjoint=True) - if df_duT is not None: - dA_duIT = ATinv * df_duT - else: - dA_duIT = ATinv * PTv - + # Get the + dA_duIT = ATinv * PTv dA_dmT = self.getADeriv_m(freq, u_src, dA_duIT, adjoint=True) - dRHS_dmT = self.getRHSDeriv_m(src, dA_duIT, adjoint=True) - + # Make du_dmT if dRHS_dmT is None: du_dmT = - dA_dmT else: du_dmT = -dA_dmT + dRHS_dmT - - df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None) - dfT_dm = df_dmFun(src, PTv, adjoint=True) - if dfT_dm is not None: - du_dmT += dfT_dm - + # Select the correct component real_or_imag = rx.projComp if real_or_imag == 'real': Jtv += du_dmT.real diff --git a/simpegMT/FieldsMT.py b/simpegMT/FieldsMT.py index 0fe0f3be..b9832309 100644 --- a/simpegMT/FieldsMT.py +++ b/simpegMT/FieldsMT.py @@ -80,10 +80,11 @@ class FieldsMT_1D(FieldsMT): return - 1./(1j*omega(src.freq)) * (C * v) def _bSecondaryDeriv_m(self, src, v, adjoint = False): - S_mDeriv, _ = src.evalDeriv(self.survey.prob, adjoint) - S_mDeriv = S_mDeriv(v) - if S_mDeriv is not None: - return 1./(1j * omega(src.freq)) * S_mDeriv + # Doesn't depend on m + # _, S_eDeriv = src.evalDeriv(self.survey.prob, adjoint) + # S_eDeriv = S_eDeriv(v) + # if S_eDeriv is not None: + # return 1./(1j * omega(src.freq)) * S_eDeriv return None def _bDeriv_u(self, src, v, adjoint=False): @@ -94,3 +95,27 @@ class FieldsMT_1D(FieldsMT): # Assuming the primary does not depend on the model return self._bSecondaryDeriv_m(src, v, adjoint) + def _fDeriv_u(self, src, v, adjoint=False): + """ + Derivative of the fields object wrt u. + + :param MTsrc src: MT source + :param numpy.ndarray v: random vector of f_sol.size + This function stacks the fields derivatives appropriately + + return a vector of size (nreEle+nrbEle) + """ + + de_du = v #Utils.spdiag(np.ones((self.nF,))) + db_du = self._bDeriv_u(src, v, adjoint) + # Return the stack + # This doesn't work... + return np.vstack((de_du,db_du)) + + def _fDeriv_m(self, src, v, adjoint=False): + """ + Derivative of the fields object wrt m. + + This function stacks the fields derivatives appropriately + """ + return None \ No newline at end of file diff --git a/simpegMT/ProblemMT1D/Problems.py b/simpegMT/ProblemMT1D/Problems.py index 7ffe46fa..a667cc7c 100644 --- a/simpegMT/ProblemMT1D/Problems.py +++ b/simpegMT/ProblemMT1D/Problems.py @@ -60,7 +60,8 @@ class eForm_psField(BaseMTProblem): if adjoint: return 1j * omega(freq) * ( dMf_dsig.T * v ) # Note: output has to be nN/nF, not nC/nE. - return 1j * omega(freq) * ( (dMf_dsig * dMf_dsig.T)**(1/2) * v) + # v should be nC + return 1j * omega(freq) * ( dMf_dsig * v ) def getRHS(self, freq): """ @@ -106,17 +107,11 @@ class eForm_psField(BaseMTProblem): # Store the fields Src = self.survey.getSrcByFreq(freq)[0] - # Calculate total e + # NOTE: only store the e_solution(secondary), all other components calculated in the fields object + F[Src, 'e_1dSolution'] = e_s[:,1] # Only storing the yx polarization as 1d - e = Src.ePrimary(self) + e_s - - # Store the fields - # NOTE: only store - F[Src, 'e_1dSolution'] = e[:,1] # Only storing the yx polarization as 1d - # F[Src, 'e_py'] = 0*e[:,0] # Note curl e = -iwb so b = -curl e /iw # b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) ) - # F[Src, 'b_px'] = 0*b[:,0] # F[Src, 'b_1d'] = b[:,1] if self.verbose: print 'Ran for {:f} seconds'.format(time.time()-startTime) diff --git a/simpegMT/SurveyMT.py b/simpegMT/SurveyMT.py index 491966dd..96f6a902 100644 --- a/simpegMT/SurveyMT.py +++ b/simpegMT/SurveyMT.py @@ -85,7 +85,7 @@ class RxMT(Survey.BaseRx): def projectFields(self, src, mesh, f): ''' - Project the fields and return the + Project the fields and return the correct data. ''' if self.projType is 'Z1D': @@ -131,32 +131,48 @@ class RxMT(Survey.BaseRx): def projectFieldsDeriv(self, src, mesh, f, v, adjoint=False): """ The derivative of the projection wrt u + + :param MTsrc src: MT source + :param TensorMesh mesh: Mesh defining the topology of the problem + :param MTfields f: MT fields object of the source + :param numpy.ndarray v: Random vector of size """ real_or_imag = self.projComp + if not adjoint: if self.projType is 'Z1D': Pex = mesh.getInterpolationMat(self.locs,'Fx') Pbx = mesh.getInterpolationMat(self.locs,'Ex') # ex = Pex*mkvc(f[src,'e_1d'],2) # bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0 - deriv_complex = Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Pex*v) - Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))*(Pbx*f._bDeriv_u(src,v)/mu_0) + dP_de = Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Pex*v) + dP_db = - Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))*(Pbx*f._bDeriv_u(src,v)/mu_0) + PDeriv_complex = np.sum((dP_de,dP_db)) # elif self.projType is 'Z2D elif self.projType is 'Z3D': - pass - Pv = getattr(deriv_complex, real_or_imag) + raise NotImplementedError('Has not be implement for full impedance tensor') + Pv = np.array(getattr(PDeriv_complex, real_or_imag)) elif adjoint: - raise NotImplementedError('must be real or imag') - + if self.projType is 'Z1D': + Pex = mesh.getInterpolationMat(self.locs,'Fx') + Pbx = mesh.getInterpolationMat(self.locs,'Ex') + # 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 + dP_dbTv = -mkvc(f._bDeriv_u(src,db_duv,adjoint=True)*v,2) + PDeriv_complex = np.sum((dP_deTv,dP_dbTv)) + elif self.projType is 'Z3D': + raise NotImplementedError('must be real or imag') + Pv = np.array(getattr(PDeriv_complex, real_or_imag)) return Pv -# Note: Might need to add tests to make sure that both polarization have the same rxList. - ############### ### Sources ### ############### -# Note: Should like inheret from FDEM + class srcMT(SrcFDEM): # Survey.BaseSrc): ''' Sources for the MT problem. @@ -238,15 +254,16 @@ class srcMT_polxy_1Dprimary(srcMT): if problem.mesh.dim == 1: # Need to use the faceInnerProduct MsigmaDeriv = problem.mesh.getFaceInnerProductDeriv(problem.curModel.sigma)(self.ePrimary(problem)[:,-1]) * problem.curModel.sigmaDeriv - MsigmaDeriv = ( MsigmaDeriv * MsigmaDeriv.T)**2 + # MsigmaDeriv = ( MsigmaDeriv * MsigmaDeriv.T)**2 if problem.mesh.dim == 2: pass if problem.mesh.dim == 3: MsigmaDeriv = problem.MeSigmaDeriv(self.ePrimary(problem)) if adjoint: + # return MsigmaDeriv.T * v else: - # Moved the v in front to make the multi work + # v should be nC size return MsigmaDeriv * v