From 18546d038b2f4d25bfb2c452ea5e81dc1e017df6 Mon Sep 17 00:00:00 2001 From: seogi Date: Fri, 28 Feb 2014 16:44:12 -0800 Subject: [PATCH] resolved conflicts --- docs/api_PF.rst | 25 +- simpegPF/BaseMag.py | 63 +++- simpegPF/Magnetics.py | 291 +++++++++++------ simpegPF/Tests/test_forward_PFproblem.py | 17 +- simpegPF/Tests/test_sensitivity_PFproblem.py | 316 +++++++++++++++++++ simpegPF/notebooks/Jacobian.ipynb | 200 ++++++------ 6 files changed, 695 insertions(+), 217 deletions(-) create mode 100644 simpegPF/Tests/test_sensitivity_PFproblem.py diff --git a/docs/api_PF.rst b/docs/api_PF.rst index d045dd08..9c8fe7e6 100644 --- a/docs/api_PF.rst +++ b/docs/api_PF.rst @@ -111,13 +111,13 @@ In applied geophysics, which means in practice, it is common to refer to measure .. math :: - \triangle\vec{B} = |\vec{B}_0-\vec{B}_s|-|\vec{B}_0| \approx |\vec{B}_s|cos \theta + \triangle\vec{B} = |\hat{B}_o-\vec{B}_s|-|\hat{B}_o| \approx |\vec{B}_s|cos \theta -where \\(\\theta\\) is the angle between total and anomalous fields. Equivalently, we can use the vector dot product to show that the anomalous field is aproximately equal to the projection of that field onto the direction of the inducing field. Using this approach we would write +where \\(\\theta\\) is the angle between total and anomalous fields, \\(\\hat{B}_o\\) is the unit vector for \\(\\vec{B}_o\\). Equivalently, we can use the vector dot product to show that the anomalous field is aproximately equal to the projection of that field onto the direction of the inducing field. Using this approach we would write .. math :: - \triangle\vec{B} = |\vec{B}_s|cos \theta = |\vec{B}_o||\vec{B}_s|cos \theta = \vec{B}_o \cdot \vec{B}_s + \triangle\vec{B} = |\vec{B}_s|cos \theta = |\hat{B}_o||\vec{B}_s|cos \theta = \hat{B}_o \cdot \vec{B}_s This is important because, in practice we usually use a total field magnetometer (like a proton precession or optically pumped sensor), which can measure only that part of the anomalous field which is in the direction of the earth's main field. @@ -128,6 +128,19 @@ Sphere in a whole space Forward problem =============== +Differential equation approach +------------------------------ + + .. math :: + + \mathbf{A}\mathbf{u} = \mathbf{rhs} + + \mathbf{A} = \Div(\MfMui)^{-1}\Div^{T} + + \mathbf{rhs} = \Div(\MfMui)^{-1}\mathbf{M}^f_{\mu_0^{-1}}\mathbf{B}_0 - \Div\mathbf{B}_0+\diag(v)\mathbf{D} \mathbf{P}_{out}^T \mathbf{B}_{sBC} + + \mathbf{B}_s = (\MfMui)^{-1}\mathbf{M}^f_{\mu_0^{-1}}\mathbf{B}_0-\mathbf{B}_0 -(\MfMui)^{-1}\Div^T \mathbf{u} + Inverse problem @@ -143,6 +156,12 @@ Mag Differential eq. approach :undoc-members: :inherited-members: +.. autoclass:: simpegPF.BaseMag.BaseMagData + :show-inheritance: + :members: + :undoc-members: + :inherited-members: + Mag Integral eq. approach ************************* diff --git a/simpegPF/BaseMag.py b/simpegPF/BaseMag.py index 1e53f12e..c7b3ea28 100644 --- a/simpegPF/BaseMag.py +++ b/simpegPF/BaseMag.py @@ -11,10 +11,14 @@ class BaseMagData(Data.BaseData): def __init__(self, **kwargs): Data.BaseData.__init__(self, **kwargs) - #TODO: change to inc, dec, intensity - def setBackgroundField(self, x=1., y=0., z=0.): - # Primary field in x-direction (background) - self.B0 = np.r_[x,y,z] + + def setBackgroundField(self, Inc, Dec, Btot): + + Bx = Btot*np.cos(Inc/180.*np.pi)*np.sin(Dec/180.*np.pi) + By = Btot*np.cos(Inc/180.*np.pi)*np.cos(Dec/180.*np.pi) + Bz = -Btot*np.sin(Inc/180.*np.pi) + + self.B0 = np.r_[Bx,By,Bz] @property def Qfx(self): @@ -34,36 +38,63 @@ class BaseMagData(Data.BaseData): self._Qfz = self.prob.mesh.getInterpolationMat(self.rxLoc,'Fz') return self._Qfz - def projectFields(self, B): + def projectFields(self, u): """ This function projects the fields onto the data space. + Esepcially, here for we use total magnetic intensity (TMI) data, + which is common in practice. + + First we project our B on to data location .. math:: - d_\\text{pred} = \mathbf{P} u(m) + + \mathbf{B}_{rec} = \mathbf{P} \mathbf{B} + + then we take the dot product between B and b_0 + + .. math :: + + \\text{TMI} = \\vec{B}_s \cdot \hat{B}_0 + """ + #TODO: There can be some different tyes of data like |B| or B - # bfx = self.Qfx*B - # bfy = self.Qfy*B - bfz = self.Qfz*B - return bfz + bfx = self.Qfx*u['B'] + bfy = self.Qfy*u['B'] + bfz = self.Qfz*u['B'] - # return np.sqrt(bfx**2 + bfy**2 + bfz**2) - + # Generate unit vector + B0 = self.prob.data.B0 + Bot = np.sqrt(B0[0]**2+B0[1]**2+B0[2]**2) + box = B0[0]/Bot + boy = B0[1]/Bot + boz = B0[2]/Bot + + # return bfx*box + bfx*boy + bfx*boz + return bfx*box + bfy*boy + bfz*boz - # return np.sqrt(bfx**2 + bfy**2 + bfz**2) @Utils.count def projectFieldsDeriv(self, B): """ This function projects the fields onto the data space. - .. math:: - \\frac{\partial d_\\text{pred}}{\partial u} = \mathbf{P} + \\frac{\partial d_\\text{pred}}{\partial \mathbf{B}} = \mathbf{P} + + Esepcially, this function is for TMI data type + """ - return self.Qfz + # Generate unit vector + B0 = self.prob.data.B0 + Bot = np.sqrt(B0[0]**2+B0[1]**2+B0[2]**2) + box = B0[0]/Bot + boy = B0[1]/Bot + boz = B0[2]/Bot + + return self.Qfx*box+self.Qfy*boy+self.Qfz*boz def projectFieldsAsVector(self, B): diff --git a/simpegPF/Magnetics.py b/simpegPF/Magnetics.py index eccf6786..38de9b6d 100644 --- a/simpegPF/Magnetics.py +++ b/simpegPF/Magnetics.py @@ -41,6 +41,7 @@ class MagneticsDiffSecondary(Problem.BaseProblem): self._MfMu0 = self.mesh.getFaceMass(1/mu_0) # self._MfMu0 = self.mesh.getFaceInnerProduct(1/mu_0) + @Utils.requires('data') def getB0(self): b0 = self.data.B0 B0 = np.r_[b0[0]*np.ones(self.mesh.nFx), @@ -49,17 +50,21 @@ class MagneticsDiffSecondary(Problem.BaseProblem): return B0 def getRHS(self, m): + """ + + .. math :: + \mathbf{rhs} = \Div(\MfMui)^{-1}\mathbf{M}^f_{\mu_0^{-1}}\mathbf{B}_0 - \Div\mathbf{B}_0+\diag(v)\mathbf{D} \mathbf{P}_{out}^T \mathbf{B}_{sBC} + + """ B0 = self.getB0() Dface = self.mesh.faceDiv Mc = Utils.sdiag(self.mesh.vol) chi = self.model.transform(m, asMu=False) - Bbc, const = CongruousMagBC(self.mesh, self.data.B0, chi) - self.Bbc_const = const + Bbc, Bbc_const = CongruousMagBC(self.mesh, self.data.B0, chi) self.Bbc = Bbc - #TODO: put congrous BC back in - # return self._Div*self.MfMuI*self.MfMu0*B0 - self._Div*B0 #+ Mc*Dface*self._Pout.T*Bbc + self.Bbc_const = Bbc_const return self._Div*self.MfMuI*self.MfMu0*B0 - self._Div*B0 + Mc*Dface*self._Pout.T*Bbc def getA(self, m): @@ -68,41 +73,39 @@ class MagneticsDiffSecondary(Problem.BaseProblem): The A matrix has the form: - .. math :: - - \mathbf{A}\mathbf{u} = \mathbf{rhs} - - \mathbf{A} = - \Div(\MfMui)^{-1}\Div^{T} - - \mathbf{rhs} = - \Div(\MfMui)^{-1}\mathbf{M}^f_{\\frac{1}{\mu_0}}\mathbf{B}_0 + \Div\mathbf{B}_0-\diag(v)\mathbf{D} \mathbf{P}_{out}^T \mathbf{B}_{sBC} + .. math :: + \mathbf{A} = \Div(\MfMui)^{-1}\Div^{T} """ return self._Div*self.MfMuI*self._Div.T def fields(self, m): + """ + Return magnetic potential (u) and flux (B) + u: defined on the cell center [nC x 1] + B: defined on the cell center [nF x 1] + + After we compute u, then we update B. + + .. math :: + + \mathbf{B}_s = (\MfMui)^{-1}\mathbf{M}^f_{\mu_0^{-1}}\mathbf{B}_0-\mathbf{B}_0 -(\MfMui)^{-1}\Div^T \mathbf{u} + + """ self.makeMassMatrices(m) A = self.getA(m) rhs = self.getRHS(m) - m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(1/A.diagonal())) u, info = sp.linalg.bicgstab(A, rhs, tol=1e-6, maxiter=1000, M=m1) - B0 = self.getB0() - B = self.MfMuI*self.MfMu0*B0-B0-self.MfMuI*self._Div.T*u - #TODO: Create a mag fields object class. - # F = self.getInitialFields() - # e.g. {'B': B, 'u': u} - return {'B': B, 'u': u} - # return self.forward(m, self.getRHS, self.calcFields, F=F) - @Utils.timeIt - def Jvec(self, m, v, u=None): + def Jvec(self, m, v, fields=None): """ Computing Jacobian multiplied by vector @@ -122,7 +125,7 @@ class MagneticsDiffSecondary(Problem.BaseProblem): \\frac{\delta \mathbf{u}}{\delta \mathbf{m}} = - [\\nabla_u \mathbf{C}(\mathbf{u})]^{-1}\\nabla_m \mathbf{C}(\mathbf{m}) With some linear algebra we can have - + .. math :: \\nabla_u \mathbf{C}(\mathbf{u}) = \mathbf{A} @@ -138,66 +141,177 @@ class MagneticsDiffSecondary(Problem.BaseProblem): \dMfMuI = \diag(\MfMui)^{-1}_{vec} \mathbf{Av}_{F2CC}^T\diag(\mathbf{v})\diag(\\frac{1}{\mu^2}) \\frac{\partial \mathbf{rhs}(\mathbf{m})}{\partial \mathbf{m}} = \\frac{\partial \mathbf{\mu}}{\partial \mathbf{m}} \left[ - \Div \diag(\M^f_{\mu_{0}^{-1} \mathbf{B}_0}) \dMfMuI \\right] - \diag(\mathbf{v})\mathbf{D} \mathbf{P}_{out}^T\\frac{\partial B_{sBC}}{\partial \mathbf{m}} + \Div \diag(\M^f_{\mu_{0}^{-1}}\mathbf{B}_0) \dMfMuI \\right] - \diag(\mathbf{v})\mathbf{D} \mathbf{P}_{out}^T\\frac{\partial B_{sBC}}{\partial \mathbf{m}} + + In the end, + + .. math :: + + \\frac{\delta \mathbf{u}}{\delta \mathbf{m}} = + - [ \mathbf{A} ]^{-1}\left[ \\frac{\partial \mathbf{A}}{\partial \mathbf{m}}(\mathbf{m})\mathbf{u} + - \\frac{\partial \mathbf{rhs}(\mathbf{m})}{\partial \mathbf{m}} \\right] + + A little tricky point here is we are not interested in potential (u), but interested in magnetic flux (B). + Thus, we need sensitivity for B. Now we take derivative of B w.r.t m and have + + .. math :: + + \\frac{\delta \mathbf{B}} {\delta \mathbf{m}} = \\frac{\partial \mathbf{\mu} } {\partial \mathbf{m} } + \left[ + \diag(\M^f_{\mu_{0}^{-1} } \mathbf{B}_0) \dMfMuI \\ + - \diag (\Div^T\mathbf{u})\dMfMuI + \\right ] + + - (\MfMui)^{-1}\Div^T\\frac{\delta\mathbf{u}}{\delta \mathbf{m}} + + Finally we evaluate the above, but we should remember that + + .. note :: + + We only want to evalute + + .. math :: + + \mathbf{J}\mathbf{v} = \\frac{\delta \mathbf{P}\mathbf{B}} {\delta \mathbf{m}}\mathbf{v} + + Since forming sensitivity matrix is very expensive in that this monster is "big" and "dense" matrix!! """ - if u is None: - u = self.fields(m) - - #TODO: B, u = u['B'], u['u'] - - B, u = u['B'], u['u'] + if fields is None: + fields = self.fields(m) + B, u = fields['B'], fields['u'] mu = self.model.transform(m, asMu=True) - dmudm = self.model.transform(m, asMu=True) - - P = self.data.projectFieldsDeriv(u) - - A = self.getA(m) - dCdu = A - # (Av_m)^-1 - # -(Av_m)^-2 * MfMu_dm * d/dm(1/mu(m)) - # (Av_m)^-2 * MfMu_dm * diag(mu(m)^-2) * mT_dm - - #TODO: only works for diagonal MfMui - # Some chain rule! - - - # harm_dm = Utils.sdiag(self.MfMui.diagonal()**(-2)) - # MfMu_dm = self.mesh.getFaceMassDeriv() - # dmuI_dm = Utils.sdiag(mu**(-2)) - # mT_dm = self.model.transformDeriv(m, asMu=True) - - getFIPconst = 1./3 - MfMuIvec = 1/self.MfMui.diagonal()*getFIPconst - dMfMuI = Utils.sdiag(MfMuIvec**2)*self.mesh.aveF2CC.T*Utils.sdiag(self.mesh.vol*1./mu**2) + dmudm = self.model.transformDeriv(m, asMu=True) + dmdmu = Utils.sdiag(1/(dmudm.diagonal())) + vol = self.mesh.vol Div = self._Div - # lots-o-bracket for vector multiplication first! - # MfMu_dmXv = harm_dm * ( MfMu_dm * ( dmuI_dm * ( mT_dm * v ) ) ) - #dCdm_A = D * ( Utils.sdiag( D.T * u ) * MfMu_dmXv ) - - dCdm_A = dmudm*Div * ( Utils.sdiag( Div.T * u * dMfMuI ) ) - - # rhs = D * MfMuI * MfMu0 * B0 - + Dface = self.mesh.faceDiv + P = self.data.projectFieldsDeriv(B) # Projection matrix B0 = self.getB0() - #TODO: add congrous stuff - dCdm_RHS = dmudm* Div * Utils.sdiag( self.MfMu0*B0 ) * dMfMuI - Utils.sdiag(self.mesh.vol)*self.mesh.faceDiv*self.Pout.T*self.Bbc*self.Bbc_const + MfMuIvec = 1/self.MfMui.diagonal() + dMfMuI = Utils.sdiag(MfMuIvec**2)*self.mesh.aveF2CC.T*Utils.sdiag(vol*1./mu**2) + # A = self._Div*self.MfMuI*self._Div.T + # RHS = Div*MfMuI*MfMu0*B0 - Div*B0 + Mc*Dface*Pout.T*Bbc + # C(m,u) = A*m-rhs + # dudm = -(dCdu)^(-1)dCdm - # c(m,u) = A(m)u - rhs(m) - dCdm = dCdm_A - dCdm_RHS + dCdu = self.getA(m) + dCdm_A = Div * ( Utils.sdiag( Div.T * u )* dMfMuI *dmudm ) + dCdm_RHS1 = Div * (Utils.sdiag( self.MfMu0*B0 ) * dMfMuI) + temp1 = (Dface*(self._Pout.T*self.Bbc_const*self.Bbc)) + dCdm_RHS2v = (Utils.sdiag(vol)*temp1)*np.inner(vol, v) + dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v + dCdm_v = dCdm_A*v - dCdm_RHSv - solve = Solver(dCdu) + m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(1/dCdu.diagonal())) + sol, info = sp.linalg.bicgstab(dCdu, dCdm_v, tol=1e-8, maxiter=1000, M=m1) - #TODO: Multiply by the dP(u(m))/du - # We transformed u in our fields object. - # ( dBdu * + dBdm(u) ) - Jv = - P * solve.solve(dCdm) - return Utils.mkvc(Jv) + if info > 0: + raise Exception ("Iterative solver did not work well") + + # B = self.MfMuI*self.MfMu0*B0-B0-self.MfMuI*self._Div.T*u + # dBdm = d\mudm*dBd\mu + + dudm = -sol + dBdmv = ( Utils.sdiag(self.MfMu0*B0)*(dMfMuI * (dmudm*v)) \ + - Utils.sdiag(Div.T*u)*(dMfMuI* (dmudm*v)) \ + - self.MfMuI*(Div.T* (dudm)) ) + + return Utils.mkvc(P*dBdmv) + + @Utils.timeIt + def Jtvec(self, m, v, fields=None): + """ + Computing Jacobian^T multiplied by vector. + + .. math :: + + (\\frac{\delta \mathbf{P}\mathbf{B}} {\delta \mathbf{m}})^{T} = \left[ \mathbf{P}_{deriv}\\frac{\partial \mathbf{\mu} } {\partial \mathbf{m} } + \left[ + \diag(\M^f_{\mu_{0}^{-1} } \mathbf{B}_0) \dMfMuI \\ + - \diag (\Div^T\mathbf{u})\dMfMuI + \\right ]\\right]^{T} + + - \left[\mathbf{P}_{deriv}(\MfMui)^{-1}\Div^T\\frac{\delta\mathbf{u}}{\delta \mathbf{m}} \\right]^{T} + + where + + .. math :: + + \mathbf{P}_{derv} = \\frac{\partial \mathbf{P}}{\partial\mathbf{B}} + + .. note :: + + Here we only want to compute + + .. math :: + + \mathbf{J}^{T}\mathbf{v} = (\\frac{\delta \mathbf{P}\mathbf{B}} {\delta \mathbf{m}})^{T} \mathbf{v} + + """ + if fields is None: + fields = self.fields(m) + + B, u = fields['B'], fields['u'] + mu = self.model.transform(m, asMu=True) + dmudm = self.model.transformDeriv(m, asMu=True) + dmdmu = Utils.sdiag(1/(dmudm.diagonal())) + + vol = self.mesh.vol + Div = self._Div + Dface = self.mesh.faceDiv + P = self.data.projectFieldsDeriv(B) # Projection matrix + B0 = self.getB0() + + MfMuIvec = 1/self.MfMui.diagonal() + dMfMuI = Utils.sdiag(MfMuIvec**2)*self.mesh.aveF2CC.T*Utils.sdiag(vol*1./mu**2) + + # A = self._Div*self.MfMuI*self._Div.T + # RHS = Div*MfMuI*MfMu0*B0 - Div*B0 + Mc*Dface*Pout.T*Bbc + # C(m,u) = A*m-rhs + # dudm = -(dCdu)^(-1)dCdm + + dCdu = self.getA(m) + s = Div * ( self.MfMuI.T * ( P.T*v ) ) + + m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(1/(dCdu.T).diagonal())) + sol, info = sp.linalg.bicgstab(dCdu.T, s, tol=1e-8, maxiter=1000, M=m1) + + if info > 0: + raise Exception ("Iterative solver did not work well") + + # dCdm_A = Div * ( Utils.sdiag( Div.T * u )* dMfMuI *dmudm ) + dCdm_Atsol = ( dMfMuI.T*( Utils.sdiag( Div.T * u ) * (Div.T * dmudm)) ) * sol + + # dCdm_RHS1 = Div * (Utils.sdiag( self.MfMu0*B0 ) * dMfMuI) + dCdm_RHS1tsol = (dMfMuI.T*( Utils.sdiag( self.MfMu0*B0 ) ) * Div.T * dmudm) * sol + + # temp1 = (Dface*(self._Pout.T*self.Bbc_const*self.Bbc)) + # dCdm_RHS2v = (Utils.sdiag(vol)*temp1)*np.inner(vol, v) + temp1sol = ( Dface.T*( Utils.sdiag(vol)*sol ) ) + temp2 = self.Bbc_const*(self._Pout.T*self.Bbc).T + dCdm_RHS2tsol = vol*np.inner(temp2, temp1sol) + dCdm_RHStsol = dCdm_RHS1tsol - dCdm_RHS2tsol + + # dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v + # dCdm_v = dCdm_A*v - dCdm_RHSv + + Ctv = dCdm_Atsol - dCdm_RHStsol + + # B = self.MfMuI*self.MfMu0*B0-B0-self.MfMuI*self._Div.T*u + # dBdm = d\mudm*dBd\mu + # dPBdm^T*v = Atemp^T*P^T*v - Btemp^T*P^T*v - Ctv + + Atemp = Utils.sdiag(self.MfMu0*B0)*(dMfMuI * (dmudm)) + Btemp = Utils.sdiag(Div.T*u)*(dMfMuI* (dmudm)) + Jtv = Atemp.T*(P.T*v) - Btemp.T*(P.T*v) - Ctv + + return Utils.mkvc(Jtv) @@ -218,7 +332,12 @@ if __name__ == '__main__': # mu = (1.+chi)*mu_0 data = BaseMag.BaseMagData() - data.setBackgroundField(x=1., y=1., z=0.) + Inc = 90. + Dec = 0. + Btot = 51000 + + data.setBackgroundField(Inc, Dec, Btot) + xr = np.linspace(-300, 300, 41) yr = np.linspace(-300, 300, 41) X, Y = np.meshgrid(xr, yr) @@ -230,30 +349,14 @@ if __name__ == '__main__': prob.pair(data) - B = prob.fields(chi) - # mesh.plotSlice(B, 'F', view='vec', showIt=True) - - dpred = data.dpred(chi, u=B) - - - # ################## - # # Test J - # ################## - - # d_chi = 0.8*chi #np.random.rand(mesh.nCz) - # d_sph_ind = spheremodel(mesh, 0., 0., -100., 50) - # d_chi[d_sph_ind] = 0.02 - - # from SimPEG.Tests import checkDerivative - - # derChk = lambda m: [prob.data.dpred(m), lambda mx: -prob.Jvec(chi, mx)] - # print '\n' - # passed = checkDerivative(derChk, chi, plotIt=False, dx=d_chi, num=2) - - - # # plt.pcolor(X, Y, dpred.reshape(X.shape, order='F')) - - # # plt.show() + dpred = data.dpred(chi) + fig = plt.figure( figsize = (8,5) ) + ax = plt.subplot(111) + dat = plt.imshow(np.reshape(dpred, (xr.size, yr.size), order='F'), extent=[min(xr), max(xr), min(yr), max(yr)]) + plt.colorbar(dat, ax = ax) + plt.show() + + diff --git a/simpegPF/Tests/test_forward_PFproblem.py b/simpegPF/Tests/test_forward_PFproblem.py index e786a34d..1948c28b 100644 --- a/simpegPF/Tests/test_forward_PFproblem.py +++ b/simpegPF/Tests/test_forward_PFproblem.py @@ -4,7 +4,7 @@ import matplotlib.pyplot as plt import simpegPF as PF -class MagProblemTests(unittest.TestCase): +class MagFwdProblemTests(unittest.TestCase): def setUp(self): @@ -29,7 +29,13 @@ class MagProblemTests(unittest.TestCase): def test_anal_forward(self): data = PF.BaseMag.BaseMagData() - data.setBackgroundField(x=1., y=1., z=0.) + + Inc = 90. + Dec = 0. + Btot = 51000 + + b0 = PF.MagAnalytics.IDTtoxyz(Inc, Dec, Btot) + data.setBackgroundField(Inc, Dec, Btot) xr = np.linspace(-300, 300, 41) yr = np.linspace(-300, 300, 41) X, Y = np.meshgrid(xr, yr) @@ -38,9 +44,10 @@ class MagProblemTests(unittest.TestCase): data.rxLoc = rxLoc self.prob.pair(data) - B = self.prob.fields(self.chi) + u = self.prob.fields(self.chi) + B = u['B'] - bxa,bya,bza = PF.MagAnalytics.MagSphereAnalFunA(rxLoc[:,0],rxLoc[:,1],rxLoc[:,2],100.,0.,0.,0.,0.01,np.array([1.,1.,0.]),'secondary') + bxa,bya,bza = PF.MagAnalytics.MagSphereAnalFunA(rxLoc[:,0],rxLoc[:,1],rxLoc[:,2],100.,0.,0.,0.,0.01, b0,'secondary') dpred = data.projectFieldsAsVector(B) err = np.linalg.norm(dpred-np.r_[bxa, bya, bza])/np.linalg.norm(np.r_[bxa, bya, bza]) @@ -50,6 +57,6 @@ class MagProblemTests(unittest.TestCase): else: print "Anaytic test is passed" pass - + if __name__ == '__main__': unittest.main() diff --git a/simpegPF/Tests/test_sensitivity_PFproblem.py b/simpegPF/Tests/test_sensitivity_PFproblem.py new file mode 100644 index 00000000..de1a94cf --- /dev/null +++ b/simpegPF/Tests/test_sensitivity_PFproblem.py @@ -0,0 +1,316 @@ +import unittest +from SimPEG import * +from simpegPF import BaseMag +import matplotlib.pyplot as plt +import simpegPF as PF +from scipy.constants import mu_0 + + +class MagSensProblemTests(unittest.TestCase): + + def setUp(self): + + hxind = ((5,25,1.3),(21, 25.),(5,25,1.3)) + hyind = ((5,25,1.3),(21, 25.),(5,25,1.3)) + hzind = ((5,25,1.3),(20, 25.),(5,25,1.3)) + hx, hy, hz = Utils.meshTensors(hxind, hyind, hzind) + M = Mesh.TensorMesh([hx, hy, hz], [-hx.sum()/2,-hy.sum()/2,-hz.sum()/2]) + chibkg = 0.001 + chiblk = 0.01 + chi = np.ones(M.nC)*chibkg + + Inc = 90. + Dec = 0. + Btot = 51000 + + b0 = PF.MagAnalytics.IDTtoxyz(Inc, Dec, Btot) + sph_ind = PF.MagAnalytics.spheremodel(M, 0., 0., 0., 100) + chi[sph_ind] = chiblk + model = PF.BaseMag.BaseMagModel(M) + + data = BaseMag.BaseMagData() + + data.setBackgroundField(Inc, Dec, Btot) + xr = np.linspace(-300, 300, 41) + yr = np.linspace(-300, 300, 41) + X, Y = np.meshgrid(xr, yr) + Z = np.ones((xr.size, yr.size))*150 + rxLoc = np.c_[Utils.mkvc(X), Utils.mkvc(Y), Utils.mkvc(Z)] + data.rxLoc = rxLoc + + + + prob = PF.Magnetics.MagneticsDiffSecondary(M, model) + prob.pair(data) + dpre = data.dpred(chi) + + fields = prob.fields(chi) + self.u = fields['u'] + self.B = fields['B'] + + self.data = data + self.model = model + self.prob = prob + self.M = M + self.chi = chi + + + + def test_mass(self): + print '\n >>Derivative for MfMuI works.' + mu = self.model.transform(self.chi, asMu=True) + def MfmuI(mu): + + chi = mu/mu_0-1 + self.prob.makeMassMatrices(chi) + vol = self.prob.mesh.vol + aveF2CC = self.prob.mesh.aveF2CC + MfMuI = self.prob.MfMuI.diagonal() + + return MfMuI + + def dMfmuI(mu, v): + + chi = mu/mu_0-1 + self.prob.makeMassMatrices(chi) + vol = self.prob.mesh.vol + aveF2CC = self.prob.mesh.aveF2CC + MfMuI = self.prob.MfMuI.diagonal() + dMfMuI = Utils.sdiag(MfMuI**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) + + return dMfMuI*v + + d_mu = mu*0.8 + derChk = lambda m: [MfmuI(m), lambda mx: dMfmuI(self.chi, mx)] + passed = Tests.checkDerivative(derChk, mu, num=4, dx = d_mu, plotIt=False) + + + self.assertTrue(passed) + + + def test_dCdm_Av(self): + print '\n >>Derivative for Cm_A.' + Div = self.prob._Div + vol = self.prob.mesh.vol + aveF2CC = self.prob.mesh.aveF2CC + + def Cm_A(chi): + dmudm = self.model.transformDeriv(chi, asMu=True) + u = self.u + # chi = mu/mu_0-1 + self.prob.makeMassMatrices(chi) + mu = self.model.transform(self.chi, asMu=True) + A = self.prob.getA(self.chi) + MfMuIvec = 1/self.prob.MfMui.diagonal() + dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) + + Cm_A = A*u + + return Cm_A + + def dCdm_A(chi, v): + + dmudm = self.model.transformDeriv(chi, asMu=True) + u = self.u + self.prob.makeMassMatrices(chi) + mu = self.model.transform(self.chi, asMu=True) + A = self.prob.getA(self.chi) + MfMuIvec = 1/self.prob.MfMui.diagonal() + dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) + + Cm_A = A*u + dCdm_A = Div * ( Utils.sdiag( Div.T * u )* dMfMuI *dmudm ) + + return dCdm_A*v + + d_chi = self.chi*0.8 + derChk = lambda m: [Cm_A(m), lambda mx: dCdm_A(self.chi, mx)] + passed = Tests.checkDerivative(derChk, self.chi, num=4, dx = d_chi, plotIt=False) + self.assertTrue(passed) + + + def test_dCdmu_RHS(self): + print '\n >>Derivative for Cm_RHS.' + u = self.u + Div = self.prob._Div + mu = self.model.transform(self.chi, asMu=True) + vol = self.prob.mesh.vol + Mc = Utils.sdiag(vol) + aveF2CC = self.prob.mesh.aveF2CC + B0 = self.prob.getB0() + Dface = self.prob.mesh.faceDiv + + def Cm_RHS(chi): + + self.prob.makeMassMatrices(chi) + dmudm = self.model.transformDeriv(chi, asMu=True) + dchidmu = Utils.sdiag(1/(dmudm.diagonal())) + Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.data.B0, chi) + MfMuIvec = 1/self.prob.MfMui.diagonal() + dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) + RHS1 = Div*self.prob.MfMuI*self.prob.MfMu0*B0 + RHS2 = Mc*Dface*self.prob._Pout.T*Bbc + RHS = RHS1 + RHS2 + Div*B0 + + return RHS + + + def dCdm_RHS(chi, v): + + + self.prob.makeMassMatrices(chi) + dmudm = self.model.transformDeriv(chi, asMu=True) + dmdmu = Utils.sdiag(1/(dmudm.diagonal())) + Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.data.B0, chi) + MfMuIvec = 1/self.prob.MfMui.diagonal() + dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) + dCdm_RHS1 = Div * (Utils.sdiag( self.prob.MfMu0*B0 ) * dMfMuI) + temp1 = (Dface*(self.prob._Pout.T*Bbc_const*Bbc)) + dCdm_RHS2v = (Utils.sdiag(vol)*temp1)*np.inner(vol, v) + dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v + + return dCdm_RHSv + + d_chi = self.chi*0.8 + derChk = lambda m: [Cm_RHS(m), lambda mx: dCdm_RHS(self.chi, mx)] + passed = Tests.checkDerivative(derChk, self.chi, num=4, dx = d_chi, plotIt=False) + self.assertTrue(passed) + + + def test_dudm(self): + print ">> Derivative test for dudm" + u = self.u + Div = self.prob._Div + mu = self.model.transform(self.chi, asMu=True) + vol = self.prob.mesh.vol + Mc = Utils.sdiag(vol) + aveF2CC = self.prob.mesh.aveF2CC + B0 = self.prob.getB0() + Dface = self.prob.mesh.faceDiv + + def ufun(chi): + u = self.prob.fields(chi)['u'] + return u + + def dudm(chi, v): + + chi = mu/mu_0-1 + self.prob.makeMassMatrices(chi) + u = self.u + dmudm = self.model.transformDeriv(chi, asMu=True) + dmdmu = Utils.sdiag(1/(dmudm.diagonal())) + Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.data.B0, chi) + MfMuIvec = 1/self.prob.MfMui.diagonal() + dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) + dCdu = self.prob.getA(chi) + dCdm_A = Div * ( Utils.sdiag( Div.T * u )* dMfMuI *dmudm ) + dCdm_RHS1 = Div * (Utils.sdiag( self.prob.MfMu0*B0 ) * dMfMuI) + temp1 = (Dface*(self.prob._Pout.T*Bbc_const*Bbc)) + dCdm_RHS2v = (Utils.sdiag(vol)*temp1)*np.inner(vol, v) + dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v + dCdm_v = dCdm_A*v - dCdm_RHSv + m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(1/dCdu.diagonal())) + sol, info = sp.linalg.bicgstab(dCdu, dCdm_v, tol=1e-8, maxiter=1000, M=m1) + + dudm = -sol + + return dudm + + d_chi = 10.0*self.chi #np.random.rand(mesh.nCz) + d_sph_ind = PF.MagAnalytics.spheremodel(self.prob.mesh, 0., 0., -50., 50) + d_chi[d_sph_ind] = 0.1 + + derChk = lambda m: [ufun(m), lambda mx: dudm(self.chi, mx)] + # TODO: I am not sure why the order get worse as step decreases .. --; + passed = Tests.checkDerivative(derChk, self.chi, num=2, dx = d_chi, plotIt=False) + self.assertTrue(passed) + + + def test_dBdm(self): + print ">> Derivative test for dBdm" + u = self.u + Div = self.prob._Div + mu = self.model.transform(self.chi, asMu=True) + vol = self.prob.mesh.vol + Mc = Utils.sdiag(vol) + aveF2CC = self.prob.mesh.aveF2CC + B0 = self.prob.getB0() + Dface = self.prob.mesh.faceDiv + + def Bfun(chi): + B = self.prob.fields(chi)['B'] + return B + + def dBdm(chi, v): + + chi = mu/mu_0-1 + self.prob.makeMassMatrices(chi) + u = self.u + dmudm = self.model.transformDeriv(chi, asMu=True) + dmdmu = Utils.sdiag(1/(dmudm.diagonal())) + Bbc, Bbc_const = PF.MagAnalytics.CongruousMagBC(self.prob.mesh, self.data.B0, chi) + MfMuIvec = 1/self.prob.MfMui.diagonal() + dMfMuI = Utils.sdiag(MfMuIvec**2)*aveF2CC.T*Utils.sdiag(vol*1./mu**2) + dCdu = self.prob.getA(chi) + dCdm_A = Div * ( Utils.sdiag( Div.T * u )* dMfMuI *dmudm ) + dCdm_RHS1 = Div * (Utils.sdiag( self.prob.MfMu0*B0 ) * dMfMuI) + temp1 = (Dface*(self.prob._Pout.T*Bbc_const*Bbc)) + dCdm_RHS2v = (Utils.sdiag(vol)*temp1)*np.inner(vol, v) + dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v + dCdm_v = dCdm_A*v - dCdm_RHSv + m1 = sp.linalg.interface.aslinearoperator(Utils.sdiag(1/dCdu.diagonal())) + sol, info = sp.linalg.bicgstab(dCdu, dCdm_v, tol=1e-8, maxiter=1000, M=m1) + + dudm = -sol + dBdmv = ( Utils.sdiag(self.prob.MfMu0*B0)*(dMfMuI * (dmudm*v)) \ + - Utils.sdiag(Div.T*u)*(dMfMuI* (dmudm*v)) \ + - self.prob.MfMuI*(Div.T* (dudm)) ) + + return dBdmv + + d_chi = 10.0*self.chi #np.random.rand(mesh.nCz) + d_sph_ind = PF.MagAnalytics.spheremodel(self.prob.mesh, 0., 0., -50., 50) + d_chi[d_sph_ind] = 0.1 + + derChk = lambda m: [Bfun(m), lambda mx: dBdm(self.chi, mx)] + # TODO: I am not sure why the order get worse as step decreases .. --; + passed = Tests.checkDerivative(derChk, self.chi, num=3, dx = d_chi, plotIt=False) + self.assertTrue(passed) + + + + def test_Jvec(self): + print ">> Derivative test for Jvec" + mu = self.model.transform(self.chi, asMu=True) + + d_chi = 10.0*self.chi #np.random.rand(mesh.nCz) + d_sph_ind = PF.MagAnalytics.spheremodel(self.prob.mesh, 0., 0., -50., 50) + d_chi[d_sph_ind] = 0.1 + + a = self.prob.Jvec(self.chi, d_chi) + + derChk = lambda m: [self.data.dpred(m), lambda mx: self.prob.Jvec(self.chi, mx)] + # TODO: I am not sure why the order get worse as step decreases .. --; + passed = Tests.checkDerivative(derChk, self.chi, num=2, dx = d_chi, plotIt=False) + self.assertTrue(passed) + + def test_Jtvec(self): + print ">> Derivative test for Jtvec" + mu = self.model.transform(self.chi, asMu=True) + dobs = self.data.dpred(self.chi) + + def misfit (m, dobs): + dpre = self.data.dpred(m) + misfit = 0.5*np.linalg.norm(dpre-dobs)**2 + residual = dpre-dobs + dmisfit = self.prob.Jtvec(self.chi, residual) + + return misfit, dmisfit + + # TODO: I am not sure why the order get worse as step decreases .. --; + derChk = lambda m: misfit(m, dobs) + passed = Tests.checkDerivative(derChk, self.chi, num=4, plotIt=False) + self.assertTrue(passed) + +if __name__ == '__main__': + unittest.main() diff --git a/simpegPF/notebooks/Jacobian.ipynb b/simpegPF/notebooks/Jacobian.ipynb index 8b8994c4..3c435223 100644 --- a/simpegPF/notebooks/Jacobian.ipynb +++ b/simpegPF/notebooks/Jacobian.ipynb @@ -28,31 +28,30 @@ ] } ], - "prompt_number": 3 + "prompt_number": 15 }, { "cell_type": "code", "collapsed": false, - "input": [ - "hxind = ((5,25,1.3),(21, 25.),(5,25,1.3))\n", - "hyind = ((5,25,1.3),(21, 25.),(5,25,1.3))\n", - "hzind = ((5,25,1.3),(20, 25.),(5,25,1.3))\n", - "hx, hy, hz = Utils.meshTensors(hxind, hyind, hzind)\n", - "M3 = Mesh.TensorMesh([hx, hy, hz], [-sum(hx)/2,-sum(hy)/2,-sum(hz)/2])" - ], + "input": [], "language": "python", "metadata": {}, "outputs": [], - "prompt_number": 4 + "prompt_number": 16 }, { "cell_type": "code", "collapsed": false, "input": [ "import matplotlib.pyplot as plt\n", - "hxind = ((5,25,1.3),(41, 12.5),(5,25,1.3))\n", - "hyind = ((5,25,1.3),(41, 12.5),(5,25,1.3))\n", - "hzind = ((5,25,1.3),(40, 12.5),(5,25,1.3))\n", + "# hxind = ((5,25,1.3),(41, 12.5),(5,25,1.3))\n", + "# hyind = ((5,25,1.3),(41, 12.5),(5,25,1.3))\n", + "# hzind = ((5,25,1.3),(40, 12.5),(5,25,1.3))\n", + "\n", + "hxind = ((5,25,1.3),(21, 25.),(5,25,1.3))\n", + "hyind = ((5,25,1.3),(21, 25.),(5,25,1.3))\n", + "hzind = ((5,25,1.3),(20, 25.),(5,25,1.3))\n", + "\n", "hx, hy, hz = Utils.meshTensors(hxind, hyind, hzind)\n", "mesh = Mesh.TensorMesh([hx, hy, hz], [-hx.sum()/2,-hy.sum()/2,-hz.sum()/2])\n", "\n", @@ -80,50 +79,22 @@ "language": "python", "metadata": {}, "outputs": [], - "prompt_number": 5 + "prompt_number": 50 }, { "cell_type": "code", "collapsed": false, "input": [ - "phi, B = prob.fields(chi)\n", + "fields = prob.fields(chi)\n", + "u = fields['u']\n", + "B = fields['B']\n", "# mesh.plotSlice(B, 'F', view='vec', showIt=True)\n", - "dpred = data.dpred(chi, u=B)" + "dpred = data.dpred(chi, fields)" ], "language": "python", "metadata": {}, "outputs": [], - "prompt_number": 4 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "test = {'B': B, 'u': phi}" - ], - "language": "python", - "metadata": {}, - "outputs": [], - "prompt_number": 5 - }, - { - "cell_type": "code", - "collapsed": false, - "input": [ - "print test['B']" - ], - "language": "python", - "metadata": {}, - "outputs": [ - { - "output_type": "stream", - "stream": "stdout", - "text": [ - "[ 0.001 0.00032989 0.00039923 ..., 0. 0. 0. ]\n" - ] - } - ], - "prompt_number": 7 + "prompt_number": 51 }, { "cell_type": "markdown", @@ -147,11 +118,11 @@ "output_type": "stream", "stream": "stdout", "text": [ - "6.99232804893\n" + "9.55511263546\n" ] } ], - "prompt_number": 6 + "prompt_number": 52 }, { "cell_type": "code", @@ -165,21 +136,21 @@ { "metadata": {}, "output_type": "pyout", - "prompt_number": 7, + "prompt_number": 53, "text": [ - "[]" + "[]" ] }, { "metadata": {}, "output_type": "display_data", - "png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAD9CAYAAACm2+DgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHNxJREFUeJzt3X9QVNf9//HXIsw4hIQEq4jgBBQCrMAuCm7TYLIWDYmm\nJpq0NTGMyRCHMkkTO36chkkzwXaqsY3toOYPazU2TUttZ1L101GiTNxJNFHGFGorxjFGwvLLCBEQ\nfwTZPZ8//LrfEK6IDQtr8nzM7MzuPWf3vu9R7ou95yxrM8YYAQDwJWEjXQAAIDQREAAASwQEAMAS\nAQEAsERAAAAsERAAAEsDBsTFixflcrnkdDplt9tVWlrar8+f/vQnORwOZWVl6a677tLhw4cDbZWV\nlUpLS1NKSopWr1499NUDAILGdq3PQZw/f16RkZHq7e1VXl6eXnnlFeXl5QXa33//fdntdkVHR6uy\nslJlZWU6cOCAfD6fUlNTVVVVpfj4eOXm5qqiokLp6elBPygAwFd3zUtMkZGRkqSenh75fD7FxMT0\nab/zzjsVHR0tSXK5XGpsbJQkVVdXKzk5WYmJiYqIiNDChQu1ffv2oa4fABAk1wwIv98vp9Op2NhY\nzZw5U3a7/ap9N23apDlz5kiSmpqaNHHixEBbQkKCmpqahqBkAMBwCL9Wh7CwMNXW1qqzs1MFBQXy\neDxyu939+u3du1ebN2/W/v37JUk2m21QBQy2HwCgr2D/paRBr2KKjo7W3LlzdejQoX5thw8f1pIl\nS7Rjxw7ddtttkqT4+Hh5vd5AH6/Xq4SEBMvXNsaE/O2ll14a8Rqokxqpkzqv3IbDgAHR1tamjo4O\nSdKFCxe0Z88eZWdn9+nT0NCgBQsW6I033lBycnJge05Ojo4fP676+nr19PRo69atmjdvXhAOAQAQ\nDANeYmppadHixYvl9/vl9/tVWFio/Px8bdiwQZJUXFysn//85zpz5oxKSkokSREREaqurlZ4eLjW\nr1+vgoIC+Xw+FRUVsYIJAG4g11zmGvQCbLZhe7v0VVxt7iXUUOfQuRFqlKhzqN0odQ7HuZOAAIAb\n0HCcO/lTGwAASwQEAMASAQEAsERAAAAsERAAAEsEBADAEgEBALBEQAAALBEQAABLBAQAwBIBAQCw\nREAAACwREAAASwQEAMASAQEAsERAAAAsERAAAEsEBADA0oABcfHiRblcLjmdTtntdpWWlvbr8+GH\nH+rOO+/U6NGjtWbNmj5tiYmJysrKUnZ2tqZPnz60lQMAgip8oMbRo0dr7969ioyMVG9vr/Ly8rRv\n3z7l5eUF+owZM0br1q3Ttm3b+j3fZrPJ4/EoJiZm6CsHAATVNS8xRUZGSpJ6enrk8/n6nezHjh2r\nnJwcRUREWD4/2F+qDQAIjgHfQUiS3+/X1KlTdeLECZWUlMhutw/6xW02m2bNmqVRo0apuLhYS5Ys\nsexXVlYWuO92u+V2uwe9DwD4JvB4PPJ4PMO6T5sZ5K/4nZ2dKigo0Msvv2x5Al+xYoWioqK0bNmy\nwLaWlhbFxcXp9OnTmj17ttatW6cZM2b0LcBm410GAFyn4Th3DnoVU3R0tObOnatDhw4N+sXj4uIk\nXb4MNX/+fFVXV19/hQCAETFgQLS1tamjo0OSdOHCBe3Zs0fZ2dmWfb+cZOfPn9fZs2clSefOndPu\n3buVmZk5FDUDAIbBgHMQLS0tWrx4sfx+v/x+vwoLC5Wfn68NGzZIkoqLi9Xa2qrc3Fx1dXUpLCxM\n5eXlqqur06effqoFCxZIknp7e7Vo0SLde++9wT8iAMCQGPQcRNAKYA4CAK5bSM1BAAC+WQgIAIAl\nAgIAYImAAABYIiAAAJYICACAJQICAGCJgAAAWCIgAACWCAgAgCUCAgBgiYAAAFgiIAAAlggIAIAl\nAgIAYImAAABYIiAAAJYICACAJQICAGBpwIC4ePGiXC6XnE6n7Ha7SktL+/X58MMPdeedd2r06NFa\ns2ZNn7bKykqlpaUpJSVFq1evHtrKAQBBZTPX+Nbr8+fPKzIyUr29vcrLy9Mrr7yivLy8QPvp06f1\nySefaNu2bbrtttu0bNkySZLP51NqaqqqqqoUHx+v3NxcVVRUKD09vW8Bw/DF2wDwdTMc585rXmKK\njIyUJPX09Mjn8ykmJqZP+9ixY5WTk6OIiIg+26urq5WcnKzExERFRERo4cKF2r59+xCWDgAIpvBr\ndfD7/Zo6dapOnDihkpIS2e32Qb1wU1OTJk6cGHickJCggwcPWvYtKysL3He73XK73YPaBwB8U3g8\nHnk8nmHd5zUDIiwsTLW1ters7FRBQYE8Hs+gTuA2m23QRXwxIAAA/X35l+cVK1YEfZ+DXsUUHR2t\nuXPn6tChQ4PqHx8fL6/XG3js9XqVkJBw/RUCAEbEgAHR1tamjo4OSdKFCxe0Z88eZWdnW/b98mRJ\nTk6Ojh8/rvr6evX09Gjr1q2aN2/eEJUNAAi2AS8xtbS0aPHixfL7/fL7/SosLFR+fr42bNggSSou\nLlZra6tyc3PV1dWlsLAwlZeXq66uTlFRUVq/fr0KCgrk8/lUVFTUbwUTACB0XXOZa9ALYJkrAFy3\nkFjmCgD4ZiIgAACWCAgAgCUCAgBgiYAAAFgiIAAAlggIAIAlAgIAYImAAABYIiAAAJYICACAJQIC\nAGCJgAAAWCIgAACWCAgAgCUCAgBgiYAAAFgiIAAAlggIAIClAQPi4sWLcrlccjqdstvtKi0ttez3\n7LPPKiUlRQ6HQzU1NYHtiYmJysrKUnZ2tqZPnz60lQMAgip8oMbRo0dr7969ioyMVG9vr/Ly8rRv\n3z7l5eUF+uzcuVMfffSRjh8/roMHD6qkpEQHDhyQdPlLtT0ej2JiYoJ7FACAIXfNS0yRkZGSpJ6e\nHvl8vn4n+x07dmjx4sWSJJfLpY6ODp06dSrQbowZynoBAMNkwHcQkuT3+zV16lSdOHFCJSUlstvt\nfdqbmpo0ceLEwOOEhAQ1NTUpNjZWNptNs2bN0qhRo1RcXKwlS5ZY7qOsrCxw3+12y+12/3dHAwBf\nUx6PRx6PZ1j3ec2ACAsLU21trTo7O1VQUCCPx9PvBH61dwn79u3ThAkTdPr0ac2ePVtpaWmaMWNG\nv35fDAgAQH9f/uV5xYoVQd/noFcxRUdHa+7cuTp06FCf7fHx8fJ6vYHHjY2Nio+PlyRNmDBBkjR2\n7FjNnz9f1dXVQ1EzAGAYDBgQbW1t6ujokCRduHBBe/bsUXZ2dp8+8+bN0+uvvy5JOnDggG699VbF\nxsbq/PnzOnv2rCTp3Llz2r17tzIzM4NxDACAIBjwElNLS4sWL14sv98vv9+vwsJC5efna8OGDZKk\n4uJizZkzRzt37lRycrJuuukmvfbaa5Kk1tZWLViwQJLU29urRYsW6d577w3y4QAAhorNjPAyI5vN\nxkonALhOw3Hu5JPUAABLBAQAwBIBAQCwREAAACwREAAASwQEAMASAQEAsERAAAAsERAAAEsEBADA\nEgEBALBEQAAALBEQAABLBAQAwBIBAQCwREAAACwREAAASwQEAMASAQEAsDRgQFy8eFEul0tOp1N2\nu12lpaWW/Z599lmlpKTI4XCopqYmsL2yslJpaWlKSUnR6tWrh7ZyAEBQDRgQo0eP1t69e1VbW6vD\nhw9r79692rdvX58+O3fu1EcffaTjx4/rd7/7nUpKSiRJPp9PzzzzjCorK1VXV6eKigodPXo0eEcC\nABhS4dfqEBkZKUnq6emRz+dTTExMn/YdO3Zo8eLFkiSXy6WOjg61trbq5MmTSk5OVmJioiRp4cKF\n2r59u9LT0/vt43//96seBgB8M9x+u5SVNTz7umZA+P1+TZ06VSdOnFBJSYnsdnuf9qamJk2cODHw\nOCEhQU1NTWpubu63/eDBg5b7+J//KQvcHzPGrTFj3Nd5GADw9dbe7lF7u0fJyVJu7vDs85oBERYW\nptraWnV2dqqgoEAej0dut7tPH2PMVyri2LGyr/R8APj6c/+/22UrVqwI+h4HvYopOjpac+fO1aFD\nh/psj4+Pl9frDTxubGxUQkJCv+1er1cJCQlDUDIAYDgMGBBtbW3q6OiQJF24cEF79uxRdnZ2nz7z\n5s3T66+/Lkk6cOCAbr31VsXGxionJ0fHjx9XfX29enp6tHXrVs2bNy9IhwEAGGoDXmJqaWnR4sWL\n5ff75ff7VVhYqPz8fG3YsEGSVFxcrDlz5mjnzp1KTk7WTTfdpNdee+3yC4eHa/369SooKJDP51NR\nUZHlBDUAIDTZzFedQPiqBdhsX3kOAwC+aYbj3MknqQEAlggIAIAlAgIAYImAAABYIiAAAJYICACA\nJQICAGCJgAAAWCIgAACWCAgAgCUCAgBgiYAAAFgiIAAAlggIAIAlAgIAYImAAABYIiAAAJYICACA\nJQICAGBpwIDwer2aOXOmpkyZooyMDK1du7ZfnzNnzmj+/PlyOBxyuVw6cuRIoC0xMVFZWVnKzs7W\n9OnTh756AEDQ2MwA33rd2tqq1tZWOZ1OdXd3a9q0adq2bZvS09MDfZYvX65bbrlFL774oo4dO6an\nn35aVVVVkqSkpCR98MEHiomJuXoBw/DF2wDwdTMc584B30GMHz9eTqdTkhQVFaX09HQ1Nzf36XP0\n6FHNnDlTkpSamqr6+nqdPn060M7JHwBuTOGD7VhfX6+amhq5XK4+2x0Oh958803l5eWpurpan3zy\niRobGzV27FjZbDbNmjVLo0aNUnFxsZYsWWL52mVlZYH7brdbbrf7vzoYAPi68ng88ng8w7rPAS8x\nXdHd3S23262f/exneuihh/q0nT17Vs8995xqamqUmZmpDz/8UL///e+VlZWl5uZmTZgwQadPn9bs\n2bO1bt06zZgxo28BXGICgOs2HOfOawbEpUuX9MADD+j+++/X0qVLr/mCSUlJ+ve//62oqKg+21es\nWKGoqCgtW7asbwEEBABctxGfgzDGqKioSHa7/arh0NnZqZ6eHknSxo0bdc899ygqKkrnz5/X2bNn\nJUnnzp3T7t27lZmZOcTlAwCCZcA5iP379+uNN94ILFWVpJUrV6qhoUGSVFxcrLq6Oj3xxBOy2WzK\nyMjQpk2bJEmnTp3S/PnzJUm9vb1atGiR7r333mAeCwBgCA1qDiKoBXCJCQCu24hfYgIAfHMREAAA\nSwQEAMASAQEAsERAAAAsERAAAEsEBADAEgEBALBEQAAALBEQAABLBAQAwBIBAQCwREAAACwREAAA\nSwQEAMASAQEAsERAAAAsERAAAEsEBADA0oAB4fV6NXPmTE2ZMkUZGRlau3Ztvz5nzpzR/Pnz5XA4\n5HK5dOTIkUBbZWWl0tLSlJKSotWrVw999QCAoLGZAb71urW1Va2trXI6neru7ta0adO0bds2paen\nB/osX75ct9xyi1588UUdO3ZMTz/9tKqqquTz+ZSamqqqqirFx8crNzdXFRUVfZ4rDc8XbwPA181w\nnDsHfAcxfvx4OZ1OSVJUVJTS09PV3Nzcp8/Ro0c1c+ZMSVJqaqrq6+v16aefqrq6WsnJyUpMTFRE\nRIQWLlyo7du3B+kwAABDLXywHevr61VTUyOXy9Vnu8Ph0Jtvvqm8vDxVV1frk08+UWNjo5qamjRx\n4sRAv4SEBB08eNDytcvKygL33W633G739R0FAHzNeTweeTyeYd3noAKiu7tbjzzyiMrLyxUVFdWn\n7fnnn9dzzz2n7OxsZWZmKjs7W6NGjZLNZht0EV8MCABAf1/+5XnFihVB3+c1A+LSpUt6+OGH9fjj\nj+uhhx7q137zzTdr8+bNgcdJSUmaPHmyLly4IK/XG9ju9XqVkJAwRGUDAIJtwDkIY4yKiopkt9u1\ndOlSyz6dnZ3q6emRJG3cuFH33HOPoqKilJOTo+PHj6u+vl49PT3aunWr5s2bN/RHAAAIigHfQezf\nv19vvPGGsrKylJ2dLUlauXKlGhoaJEnFxcWqq6vTE088IZvNpoyMDG3atOnyC4eHa/369SooKJDP\n51NRUVG/FUwAgNA14DLXYSmAZa4AcN1GfJkrAOCbi4AAAFgiIAAAlggIAIAlAgIAYImAAABYIiAA\nAJYICACAJQICAGCJgAAAWCIgAACWCAgAgCUCAgBgiYAAAFgiIAAAlggIAIAlAgIAYImAAABYGjAg\nvF6vZs6cqSlTpigjI0Nr167t16etrU333XefnE6nMjIytGXLlkBbYmJi4Pusp0+fPuTFAwCCZ8Dv\npG5tbVVra6ucTqe6u7s1bdo0bdu2Tenp6YE+ZWVl+vzzz7Vq1Sq1tbUpNTVVp06dUnh4uJKSkvTB\nBx8oJibm6gXwndQAcN1G/Dupx48fL6fTKUmKiopSenq6mpub+/SJi4tTV1eXJKmrq0tjxoxReHh4\noJ2TPwDcmAY9B1FfX6+amhq5XK4+25csWaIjR45owoQJcjgcKi8vD7TZbDbNmjVLOTk52rhx49BV\nDQAIuvBrd5G6u7v1yCOPqLy8XFFRUX3aVq5cKafTKY/HoxMnTmj27Nn617/+pZtvvln79+9XXFyc\nTp8+rdmzZystLU0zZszo9/plZWWB+263W263+ysdFAB83Xg8Hnk8nmHd54BzEJJ06dIlPfDAA7r/\n/vu1dOnSfu1z5szRCy+8oLvuukuSlJ+fr9WrVysnJ6dPvxUrVigqKkrLli3rWwBzEABw3UZ8DsIY\no6KiItntdstwkKS0tDRVVVVJkk6dOqVjx45p0qRJOn/+vM6ePStJOnfunHbv3q3MzMwhLh8AECwD\nvoPYt2+f7r77bmVlZclms0m6fEmpoaFBklRcXKy2tjY9+eSTamhokN/vV2lpqR577DF9/PHHWrBg\ngSSpt7dXixYtUmlpaf8CeAcBANdtOM6d17zEFGwEBABcvxG/xAQA+OYiIAAAlggIAIAlAgIAYImA\nAABYIiAAAJYICACAJQICAGCJgAAAWCIgAACWCAgAgCUCAgBgiYAAAFgiIAAAlggIAIAlAgIAYImA\nAABYIiAAAJYICACApQEDwuv1aubMmZoyZYoyMjK0du3afn3a2tp03333yel0KiMjQ1u2bAm0VVZW\nKi0tTSkpKVq9evWQFz+cPB7PSJcwKNQ5dG6EGiXqHGo3Sp3DYcCAiIiI0G9/+1sdOXJEBw4c0Kuv\nvqqjR4/26bN+/XplZ2ertrZWHo9Hy5YtU29vr3w+n5555hlVVlaqrq5OFRUV/Z57I7lR/tNQ59C5\nEWqUqHOo3Sh1DocBA2L8+PFyOp2SpKioKKWnp6u5ublPn7i4OHV1dUmSurq6NGbMGIWHh6u6ulrJ\nyclKTExURESEFi5cqO3btwfpMAAAQy18sB3r6+tVU1Mjl8vVZ/uSJUv03e9+VxMmTNDZs2f117/+\nVZLU1NSkiRMnBvolJCTo4MGDQ1Q2ACDozCCcPXvWTJs2zfz973/v1/aLX/zCPPfcc8YYYz766COT\nlJRkurq6zN/+9jfz1FNPBfr98Y9/NM8880y/50vixo0bN27/xS3YrvkO4tKlS3r44Yf1+OOP66GH\nHurX/t577+mFF16QJE2ePFlJSUk6duyYEhIS5PV6A/28Xq8SEhL6Pf9yRgAAQs2AcxDGGBUVFclu\nt2vp0qWWfdLS0lRVVSVJOnXqlI4dO6ZJkyYpJydHx48fV319vXp6erR161bNmzdv6I8AABAUNjPA\nr/D79u3T3XffraysLNlsNknSypUr1dDQIEkqLi5WW1ubnnzySTU0NMjv96u0tFSPPfaYJGnXrl1a\nunSpfD6fioqKVFpaOgyHBAAYEkG/iDWAXbt2mdTUVJOcnGxefvnlYdnn7bffbjIzM43T6TS5ubnG\nGGPa29vNrFmzTEpKipk9e7Y5c+ZMoP/KlStNcnKySU1NNW+99VZg+6FDh0xGRoZJTk42zz77bGD7\nxYsXzQ9+8AOTnJxsXC6Xqa+vH1RdTz75pBk3bpzJyMgIbBuuurZs2WJSUlJMSkqK+cMf/nBdNb70\n0ksmPj7eOJ1O43Q6zc6dO0e0RmOMaWhoMG6329jtdjNlyhRTXl4ekuN5tTpDbUwvXLhgpk+fbhwO\nh0lPTzfPP/98SI7n1eoMtfE0xpje3l7jdDrNAw88EJJjecWIBURvb6+ZPHmyOXnypOnp6TEOh8PU\n1dUFfb+JiYmmvb29z7bly5eb1atXG2OMefnll81Pf/pTY4wxR44cMQ6Hw/T09JiTJ0+ayZMnG7/f\nb4wxJjc31xw8eNAYY8z9999vdu3aZYwx5tVXXzUlJSXGGGP+8pe/mB/+8IeDquudd94x//znP/uc\nfIejrvb2djNp0iRz5swZc+bMmcD9wdZYVlZm1qxZ06/vSNVojDEtLS2mpqbGGHN5gcUdd9xh6urq\nQm48r1ZnKI7puXPnjDHGXLp0ybhcLvPuu++G3Hherc5QHM81a9aYxx57zHzve98zxoTez/oVIxYQ\n7733nikoKAg8XrVqlVm1alXQ95uYmGja2tr6bEtNTTWtra3GmMs/tKmpqcaYy8n9xXc2BQUF5v33\n3zfNzc0mLS0tsL2iosIUFxcH+hw4cMAYc/k/6be+9a1B13by5Mk+J9/hqOvPf/6z+dGPfhR4TnFx\nsamoqBh0jWVlZeaVV17p128ka/yyBx980OzZsyckx9OqzlAe03PnzpmcnBzzn//8J6TH84t1htp4\ner1ek5+fb95+++3AO4hQHcsR+1tMVp+TaGpqCvp+bTabZs2apZycHG3cuFHS5cn12NhYSVJsbKxO\nnTolSWpubu6z8upKjV/eHh8fH6j9i8cVHh6u6OhoffbZZ/9VrcGuq729/aqvdT3WrVsnh8OhoqIi\ndXR0hFSNX/z8TiiP55U6v/3tb0sKvTH1+/1yOp2KjY0N/PmdUBxPqzpDbTx/8pOf6Ne//rXCwv7/\n6TcUx1IawT/Wd2XSe7jt379fNTU12rVrl1599VW9++67/eoaqdoGEqp1lZSU6OTJk6qtrVVcXJyW\nLVs20iUFdHd36+GHH1Z5ebluvvnmPm2hNJ7d3d165JFHVF5erqioqJAc07CwMNXW1qqxsVHvvPOO\n9u7d26c9VMbzy3V6PJ6QGs9//OMfGjdunLKzs6+6xD9UxlIawYCIj48f1OckhlpcXJwkaezYsZo/\nf76qq6sVGxur1tZWSVJLS4vGjRtnWWNjY6MSEhIUHx+vxsbGftuvPOfKKq/e3l51dnYqJibmv6o1\n2HWNGTPmK/87jBs3LvAf+qmnnlJ1dXVI1Hjl8zuFhYWBz++E4nhafc4oVMdUkqKjozV37lx98MEH\nITmeX67z0KFDITWe7733nnbs2KGkpCQ9+uijevvtt1VYWBi6YzngBaggunTpkpk0aZI5efKk+fzz\nz4dlkvrcuXOmq6vLGGNMd3e3+c53vmPeeusts3z58sB1vlWrVvWbIPr888/Nxx9/bCZNmhSYIJo+\nfbo5cOCA8fv9/SaIrlznq6ioGPQktTH9r+8PR13t7e0mKSnJnDlzxnz22WeB+4Otsbm5OXD/N7/5\njXn00UdHvEa/328KCwvN0qVL+2wPtfG8Wp2hNqanT58OtJ0/f97MmDHDVFVVhdx4Xq3OlpaWkBrP\nKzweT2AOItTG8ooRXea6c+dOc8cdd5jJkyeblStXBn1/H3/8sXE4HMbhcJgpU6YE9tne3m7y8/Mt\nl5j98pe/NJMnTzapqammsrIysP3KErPJkyebH//4x4HtFy9eNN///vcDS8xOnjw5qNoWLlxo4uLi\nTEREhElISDCbN28etro2b95skpOTTXJystmyZcuga9y0aZMpLCw0mZmZJisryzz44IOBibaRqtEY\nY959911js9mMw+EILG3ctWtXyI2nVZ07d+4MuTE9fPiwyc7ONg6Hw2RmZppf/epXxpjh+7n5qnWG\n2nhe4fF4AquYQm0srxjwg3IAgG8uvlEOAGCJgAAAWCIgAACWCAgAgCUCAgBgiYAAAFj6P0f5Arjx\nzIiGAAAAAElFTkSuQmCC\n", + "png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD9CAYAAABTJWtQAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHaxJREFUeJzt3X9Qm/XhB/B3KNz1MC1Kh0AJGwgUCIEkLTTnpBoEiv0h\nSquutiJ2jMsxnXbXr6ec8wbzhmXTm9T6R9fVn93Rzp1r667FwklmWwUOC6srbUcrSPhpqeU3LE3y\n+f7Rbz9fIxioJj6svl93uSPP50Oe95PS502e5wlRCSEEiIiIAPgpHYCIiOYOlgIREUksBSIiklgK\nREQksRSIiEhiKRARkeSxFCYnJ2EymWAwGKDValFSUjJlzp///Gfo9XqkpKTgtttuw8mTJ+VYdXU1\nEhISEBcXh4qKCu+nJyIir1LN9D6F8fFxBAYGwuFwID09HS+88ALS09Pl+EcffQStVougoCBUV1ej\ntLQU9fX1cDqdiI+PR21tLSIiIpCWloaqqiokJib6fKOIiOibmfHwUWBgIADAbrfD6XQiODjYbfzW\nW29FUFAQAMBkMqGrqwsA0NjYiNjYWERFRSEgIAAbNmzAgQMHvJ2fiIi8aMZScLlcMBgMCA0NRUZG\nBrRa7dfO3b17N1avXg0A6O7uRmRkpBzTaDTo7u72QmQiIvIV/5km+Pn5oaWlBUNDQ8jJyYHVaoXZ\nbJ4yr66uDq+++iqOHz8OAFCpVLMKMNt5RETkzhd/pWjWVx8FBQVhzZo1aGpqmjJ28uRJFBUV4eDB\ng7jpppsAABEREbDZbHKOzWaDRqOZ9rGFEHPq9utf/1rxDP8NmeZqLmZipu9DLl/xWAoDAwMYHBwE\nAExMTKCmpgZGo9FtTmdnJ9atW4c9e/YgNjZWLk9NTUVbWxs6Ojpgt9uxb98+5Obm+mATiIjIWzwe\nPurt7UVBQQFcLhdcLhfy8/ORmZmJnTt3AgAsFgt+85vf4NKlSyguLgYABAQEoLGxEf7+/tixYwdy\ncnLgdDpRWFjIK4+IiOY6obA5EGGKuro6pSNMMRczCTE3czHT7DDT7M3FXL7ad874PgVfU6lUPj0+\nRkR0PfLVvpN/5oKIiCSWAhERSSwFIiKSWApERCSxFIiISGIpEBGRxFIgIiKJpUBERBJLgYiIJJYC\nERFJLAUiIpJYCkREJLEUiIhIYikQEZHEUiAiIomlQEREEkuBiIgklgIREUkeS2FychImkwkGgwFa\nrRYlJSVT5pw5cwa33nor5s+fjxdffNFtLCoqCikpKTAajVi+fLl3kxMRkdf5exqcP38+6urqEBgY\nCIfDgfT0dBw7dgzp6elyzqJFi/Dyyy9j//79U75fpVLBarUiODjY+8mJiMjrZjx8FBgYCACw2+1w\nOp1TdvAhISFITU1FQEDAtN/viw+WJiIi3/D4SgEAXC4Xli5divPnz6O4uBharXbWD65SqZCVlYV5\n8+bBYrGgqKho2nmlpaXya7PZDLPZPOt1EBF9H1itVlitVp+vRyVm+av80NAQcnJysG3btml32mVl\nZVCr1di6datc1tvbi/DwcFy4cAHZ2dl4+eWXsWLFCvcAKhVfTRARXSNf7TtnffVRUFAQ1qxZg6am\nplk/eHh4OIArh5jy8vLQ2Nh47QmJiOg747EUBgYGMDg4CACYmJhATU0NjEbjtHO/2ljj4+MYGRkB\nAIyNjeHIkSNITk72RmYiIvIRj+cUent7UVBQAJfLBZfLhfz8fGRmZmLnzp0AAIvFgr6+PqSlpWF4\neBh+fn6orKxEa2srPv/8c6xbtw4A4HA4sGnTJqxcudL3W0RERN/YrM8p+CwAzykQEV0zxc8pEBHR\n9Y+lQEREEkuBiIgklgIREUksBSIiklgKREQksRSIiEhiKRARkcRSICIiiaVAREQSS4GIiCSWAhER\nSSwFIiKSWApERCSxFIiISGIpEBGRxFIgIiKJpUBERBJLgYiIJI+lMDk5CZPJBIPBAK1Wi5KSkilz\nzpw5g1tvvRXz58/Hiy++6DZWXV2NhIQExMXFoaKiwrvJiYjI61Rihk9+Hh8fR2BgIBwOB9LT0/HC\nCy8gPT1djl+4cAGfffYZ9u/fj5tuuglbt24FADidTsTHx6O2thYRERFIS0tDVVUVEhMT3QP46MOn\niYiuZ77ad854+CgwMBAAYLfb4XQ6ERwc7DYeEhKC1NRUBAQEuC1vbGxEbGwsoqKiEBAQgA0bNuDA\ngQNejE5ERN7mP9MEl8uFpUuX4vz58yguLoZWq53VA3d3dyMyMlLe12g0aGhomHZuaWmp/NpsNsNs\nNs9qHURE3xdWqxVWq9Xn65mxFPz8/NDS0oKhoSHk5OTAarXOaqetUqlmHeLLpUBERFN99RfmsrIy\nn6xn1lcfBQUFYc2aNWhqaprV/IiICNhsNnnfZrNBo9Fce0IiIvrOeCyFgYEBDA4OAgAmJiZQU1MD\no9E47dyvnvBITU1FW1sbOjo6YLfbsW/fPuTm5nopNhER+YLHw0e9vb0oKCiAy+WCy+VCfn4+MjMz\nsXPnTgCAxWJBX18f0tLSMDw8DD8/P1RWVqK1tRVqtRo7duxATk4OnE4nCgsLp1x5REREc8uMl6T6\nPAAvSSUiumaKXZJKRETfHywFIiKSWApERCSxFIiISGIpEBGRxFIgIiKJpUBERBJLgYiIJJYCERFJ\nLAUiIpJYCkREJLEUiIhIYikQEZHEUiAiIomlQEREEkuBiIgklgIREUksBSIiklgKREQkeSyFyclJ\nmEwmGAwGaLValJSUTDvv8ccfR1xcHPR6PZqbm+XyqKgopKSkwGg0Yvny5d5NTkREXufvaXD+/Pmo\nq6tDYGAgHA4H0tPTcezYMaSnp8s5hw4dwrlz59DW1oaGhgYUFxejvr4ewJUPlrZarQgODvbtVhAR\nkVfMePgoMDAQAGC32+F0Oqfs4A8ePIiCggIAgMlkwuDgIPr7++W4EMKbeYmIyIc8vlIAAJfLhaVL\nl+L8+fMoLi6GVqt1G+/u7kZkZKS8r9Fo0N3djdDQUKhUKmRlZWHevHmwWCwoKiqadh2lpaXya7PZ\nDLPZ/M22hojoOmW1WmG1Wn2+nhlLwc/PDy0tLRgaGkJOTg6sVuuUnfbXvRo4duwYFi9ejAsXLiA7\nOxsJCQlYsWLFlHlfLgUiIprqq78wl5WV+WQ9s776KCgoCGvWrEFTU5Pb8oiICNhsNnm/q6sLERER\nAIDFixcDAEJCQpCXl4fGxkZvZCYiIh/xWAoDAwMYHBwEAExMTKCmpgZGo9FtTm5uLt58800AQH19\nPW688UaEhoZifHwcIyMjAICxsTEcOXIEycnJvtgGIiLyEo+Hj3p7e1FQUACXywWXy4X8/HxkZmZi\n586dAACLxYLVq1fj0KFDiI2NxQ033IDXXnsNANDX14d169YBABwOBzZt2oSVK1f6eHOIiOjbUAmF\nLw9SqVS8QomI6Br5at/JdzQTEZHEUiAiIomlQEREEkuBiIgklgIREUksBSIiklgKREQksRSIiEhi\nKRARkcRSICIiiaVAREQSS4GIiCSWAhERSSwFIiKSWApERCSxFIiISGIpEBGRxFIgIiKJpUBERJLH\nUpicnITJZILBYIBWq0VJScm08x5//HHExcVBr9ejublZLq+urkZCQgLi4uJQUVHh3eREROR1Hkth\n/vz5qKurQ0tLC06ePIm6ujocO3bMbc6hQ4dw7tw5tLW14Y9//COKi4sBAE6nE4899hiqq6vR2tqK\nqqoqnD592ndbQkRE35r/TBMCAwMBAHa7HU6nE8HBwW7jBw8eREFBAQDAZDJhcHAQfX19aG9vR2xs\nLKKiogAAGzZswIEDB5CYmDhlHe+++203g4jo+rdgAWA2+3YdM5aCy+XC0qVLcf78eRQXF0Or1bqN\nd3d3IzIyUt7XaDTo7u5GT0/PlOUNDQ3TruN//qdUfr1okRmLFpmvcTOIiK5vFy9a4XBYsXq1b9cz\nYyn4+fmhpaUFQ0NDyMnJgdVqhfkrVSWE+FYhzp4t/VbfT0R0/TP/3+2KsrIyn6xl1lcfBQUFYc2a\nNWhqanJbHhERAZvNJu93dXVBo9FMWW6z2aDRaLwQmYiIfMVjKQwMDGBwcBAAMDExgZqaGhiNRrc5\nubm5ePPNNwEA9fX1uPHGGxEaGorU1FS0tbWho6MDdrsd+/btQ25uro82g4iIvMHj4aPe3l4UFBTA\n5XLB5XIhPz8fmZmZ2LlzJwDAYrFg9erVOHToEGJjY3HDDTfgtddeu/LA/v7YsWMHcnJy4HQ6UVhY\nOO1JZiIimjtU4tueEPi2AVSqb31Ogojo+8ZX+06+o5mIiCSWAhERSSwFIiKSWApERCSxFIiISGIp\nEBGRxFIgIiKJpUBERBJLgYiIJJYCERFJLAUiIpJYCkREJLEUiIhIYikQEZHEUiAiIomlQEREEkuB\niIgklgIREUksBSIikjyWgs1mQ0ZGBpKSkqDT6bB9+/Ypcy5duoS8vDzo9XqYTCacOnVKjkVFRSEl\nJQVGoxHLly/3fnoiIvIqlfDwyc99fX3o6+uDwWDA6Ogoli1bhv379yMxMVHOefLJJ7Fw4UI8++yz\nOHv2LB599FHU1tYCAKKjo/Hxxx8jODj46wP46MOniYiuZ77ad3p8pRAWFgaDwQAAUKvVSExMRE9P\nj9uc06dPIyMjAwAQHx+Pjo4OXLhwQY5zh09E9N/Df7YTOzo60NzcDJPJ5LZcr9fjnXfeQXp6Ohob\nG/HZZ5+hq6sLISEhUKlUyMrKwrx582CxWFBUVDTtY5eWlsqvzWYzzGbzN9oYIqLrldVqhdVq9fl6\nPB4+ump0dBRmsxm/+tWvcO+997qNjYyM4IknnkBzczOSk5Nx5swZ/OlPf0JKSgp6enqwePFiXLhw\nAdnZ2Xj55ZexYsUK9wA8fEREdM18te+csRQuX76MtWvXYtWqVdiyZcuMDxgdHY1PPvkEarXabXlZ\nWRnUajW2bt3qHoClQER0zRQ5pyCEQGFhIbRa7dcWwtDQEOx2OwBg165duOOOO6BWqzE+Po6RkREA\nwNjYGI4cOYLk5GQvxyciIm/yeE7h+PHj2LNnj7ysFADKy8vR2dkJALBYLGhtbcUjjzwClUoFnU6H\n3bt3AwD6+/uRl5cHAHA4HNi0aRNWrlzpy20hIqJvaVbnFHwagIePiIiumSKHj4iI6PuFpUBERBJL\ngYiIJJYCERFJLAUiIpJYCkREJLEUiIhIYikQEZHEUiAiIomlQEREEkuBiIgklgIREUksBSIiklgK\nREQksRSIiEhiKRARkcRSICIiiaVAREQSS4GIiCSPpWCz2ZCRkYGkpCTodDps3759ypxLly4hLy8P\ner0eJpMJp06dkmPV1dVISEhAXFwcKioqvJ+eiIi8SiU8fPJzX18f+vr6YDAYMDo6imXLlmH//v1I\nTEyUc5588kksXLgQzz77LM6ePYtHH30UtbW1cDqdiI+PR21tLSIiIpCWloaqqiq37wV89+HTRETX\nM1/tOz2+UggLC4PBYAAAqNVqJCYmoqenx23O6dOnkZGRAQCIj49HR0cHPv/8czQ2NiI2NhZRUVEI\nCAjAhg0bcODAAa9vABEReY//bCd2dHSgubkZJpPJbbler8c777yD9PR0NDY24rPPPkNXVxe6u7sR\nGRkp52k0GjQ0NEz72KWlpfJrs9kMs9l8bVtBRHSds1qtsFqtPl/PrEphdHQU9913HyorK6FWq93G\nnn76aTzxxBMwGo1ITk6G0WjEvHnzoFKpZh3iy6VARERTffUX5rKyMp+sZ8ZSuHz5MtavX4+HHnoI\n995775TxBQsW4NVXX5X3o6OjERMTg4mJCdhsNrncZrNBo9F4KTYREfmCx3MKQggUFhZCq9Viy5Yt\n084ZGhqC3W4HAOzatQt33HEH1Go1UlNT0dbWho6ODtjtduzbtw+5ubne3wIiIvIaj68Ujh8/jj17\n9iAlJQVGoxEAUF5ejs7OTgCAxWJBa2srHnnkEahUKuh0OuzevfvKA/v7Y8eOHcjJyYHT6URhYeGU\nK4+IiGhu8XhJ6ncSgJekEhFdM0UuSSUiou8XlgIREUksBSIiklgKREQksRSIiEhiKRARkcRSICIi\niaVAREQSS4GIiCSWAhERSSwFIiKSWApERCSxFIiISGIpEBGRxFIgIiKJpUBERBJLgYiIJJYCERFJ\nHkvBZrMhIyMDSUlJ0Ol02L59+5Q5AwMDuOuuu2AwGKDT6fD666/LsaioKPn5zsuXL/d6eCIi8i6P\nn9Hc19eHvr4+GAwGjI6OYtmyZdi/fz8SExPlnNLSUvznP//B888/j4GBAcTHx6O/vx/+/v6Ijo7G\nxx9/jODg4K8PwM9oJiK6Zop8RnNYWBgMBgMAQK1WIzExET09PW5zwsPDMTw8DAAYHh7GokWL4O/v\nL8e5wyci+u8x63MKHR0daG5uhslkclteVFSEU6dOYfHixdDr9aisrJRjKpUKWVlZSE1Nxa5du7yX\nmoiIfMJ/5inA6Ogo7rvvPlRWVkKtVruNlZeXw2AwwGq14vz588jOzsY///lPLFiwAMePH0d4eDgu\nXLiA7OxsJCQkYMWKFVMev7S0VH5tNpthNpu/1UYREV1vrFYrrFarz9fj8ZwCAFy+fBlr167FqlWr\nsGXLlinjq1evxjPPPIPbbrsNAJCZmYmKigqkpqa6zSsrK4NarcbWrVvdA/CcAhHRNVPknIIQAoWF\nhdBqtdMWAgAkJCSgtrYWANDf34+zZ8/illtuwfj4OEZGRgAAY2NjOHLkCJKTk70cn4iIvMnjK4Vj\nx47h9ttvR0pKClQqFYArh4s6OzsBABaLBQMDA9i8eTM6OzvhcrlQUlKCjRs34tNPP8W6desAAA6H\nA5s2bUJJScnUAHylQER0zXy175zx8JGvsRSIiK6dIoePiIjo+4WlQEREEkuBiIgklgIREUksBSIi\nklgKREQksRSIiEhiKRARkcRSICIiiaVAREQSS4GIiCSWAhERSSwFIiKSWApERCSxFIiISGIpEBGR\nxFIgIiKJpUBERBJLgYiIJI+lYLPZkJGRgaSkJOh0Omzfvn3KnIGBAdx1110wGAzQ6XR4/fXX5Vh1\ndTUSEhIQFxeHiooKr4f3FavVqnSEKeZiJmBu5mKm2WGm2ZuruXzBYykEBATgD3/4A06dOoX6+nq8\n8sorOH36tNucHTt2wGg0oqWlBVarFVu3boXD4YDT6cRjjz2G6upqtLa2oqqqasr3zlVz8QdgLmYC\n5mYuZpodZpq9uZrLFzyWQlhYGAwGAwBArVYjMTERPT09bnPCw8MxPDwMABgeHsaiRYvg7++PxsZG\nxMbGIioqCgEBAdiwYQMOHDjgo80gIiJv8J/txI6ODjQ3N8NkMrktLyoqwp133onFixdjZGQEf/nL\nXwAA3d3diIyMlPM0Gg0aGhq8FJuIiHxCzMLIyIhYtmyZ+Nvf/jZl7LnnnhNPPPGEEEKIc+fOiejo\naDE8PCzefvtt8bOf/UzOe+utt8Rjjz025fsB8MYbb7zx9g1uvjDjK4XLly9j/fr1eOihh3DvvfdO\nGf/www/xzDPPAABiYmIQHR2Ns2fPQqPRwGazyXk2mw0ajWbK91/pBSIimgs8nlMQQqCwsBBarRZb\ntmyZdk5CQgJqa2sBAP39/Th79ixuueUWpKamoq2tDR0dHbDb7di3bx9yc3O9vwVEROQ1KuHhV/Vj\nx47h9ttvR0pKClQqFQCgvLwcnZ2dAACLxYKBgQFs3rwZnZ2dcLlcKCkpwcaNGwEAhw8fxpYtW+B0\nOlFYWIiSkpLvYJOIiOgb88lBqVk6fPiwiI+PF7GxsWLbtm1efezNmzeLm2++Weh0Orns4sWLIisr\nS8TFxYns7Gxx6dIlOVZeXi5iY2NFfHy8eO+99+TypqYmodPpRGxsrHj88cfl8snJSfHAAw+I2NhY\nYTKZREdHx4yZOjs7hdlsFlqtViQlJYnKykrFc01MTIjly5cLvV4vEhMTxdNPP614pi9zOBzCYDCI\ntWvXzolcP/rRj0RycrIwGAwiLS1tTmS6dOmSWL9+vUhISBCJiYmivr5e0UxnzpwRBoNB3hYuXCgq\nKysVf57Ky8uFVqsVOp1OPPjgg2JyclLxTEII8dJLLwmdTieSkpLESy+9JIRQ9mdKsVJwOBwiJiZG\ntLe3C7vdLvR6vWhtbfXa43/wwQfixIkTbqXw5JNPioqKCiGEENu2bRNPPfWUEEKIU6dOCb1eL+x2\nu2hvbxcxMTHC5XIJIYRIS0sTDQ0NQgghVq1aJQ4fPiyEEOKVV14RxcXFQggh9u7dK37yk5/MmKm3\nt1c0NzcLIa6cvF+yZIlobW1VPNfY2JgQQojLly8Lk8kkjh49qnimq1588UWxceNGcffddwshlP83\njIqKEhcvXnRbpnSmhx9+WOzevVsIceXfcHBwUPFMVzmdThEWFiY6OzsVzdTe3i6io6PF5OSkEEKI\nBx54QLz++uuKP0+ffPKJ0Ol0YmJiQjgcDpGVlSXOnTunaC7FSuHDDz8UOTk58v7zzz8vnn/+ea+u\no7293a0U4uPjRV9fnxDiyg46Pj5eCHGleb/8SiUnJ0d89NFHoqenRyQkJMjlVVVVwmKxyDn19fVC\niCv/EX/wgx9cc7577rlH1NTUzJlcY2NjIjU1VfzrX/+aE5lsNpvIzMwU77//vnyloHSuqKgoMTAw\n4LZMyUyDg4MiOjp6ynKln6er3nvvPZGenq54posXL4olS5aIL774Qly+fFmsXbtWHDlyRPHn6e23\n3xaFhYXy/nPPPScqKioUzaXY3z6a7n0M3d3dPl1nf38/QkNDAQChoaHo7+8HAPT09LhdGXU1y1eX\nR0REyIxfzu/v74+goCB88cUXs87y5fd9KJ3L5XLBYDAgNDRU/lkTpTMBwC9/+Uv8/ve/h5/f//+Y\nKp1LpVIhKysLqamp2LVrl+KZ2tvbERISgs2bN2Pp0qUoKirC2NiY4s/TVXv37sWDDz6o+PMUHByM\nrVu34oc//CEWL16MG2+8EdnZ2Yo/TzqdDkePHsUXX3yB8fFxHDp0CF1dXYrmUqwUrp64VnL9SmUY\nHR3F+vXrUVlZiQULFiiey8/PDy0tLejq6sIHH3yAuro6xTP9/e9/x8033wyj0fi1ly0rkev48eNo\nbm7G4cOH8corr+Do0aOKZnI4HDhx4gR+/vOf48SJE7jhhhuwbds2RTNdZbfb8e677+L++++fMvZd\nZzp//jxeeukldHR0oKenB6Ojo9izZ4+imYArV28+9dRTWLlyJVatWgWDwYB58+YpmkuxUoiIiJjV\n+xi8KTQ0FH19fQCA3t5e3HzzzdNm6erqgkajQUREBLq6uqYsv/o9V6/CcjgcGBoaQnBw8IwZrr7v\nIz8/X77vYy7kAoCgoCCsWbMGH3/8seKZPvzwQxw8eBDR0dF48MEH8f777yM/P1/xXOHh4QCAkJAQ\n5OXlobGxUdFMGo0GGo0GaWlpAID77rsPJ06cQFhYmOI/U4cPH8ayZcsQEhICQNmf86amJvz4xz+W\nf4Zn3bp1+Oijj+bE8/TTn/4UTU1N+Mc//oGbbroJS5YsUfS5UqwUlHgfQ25uLt544w0AwBtvvCF3\nyrm5udi7dy/sdjva29vR1taG5cuXIywsDAsXLkRDQwOEEHjrrbdwzz33THmsv/71r8jMzJxx/eJr\n3vehZK6BgQEMDg4CACYmJlBTUwOj0aj4c1VeXg6bzYb29nbs3bsXd955J9566y1Fc42Pj2NkZAQA\nMDY2hiNHjiA5OVnRTGFhYYiMjMS///1vAEBtbS2SkpJw9913K/rvBwBVVVXy0NFXH+e7zpSQkID6\n+npMTExACIHa2lpotdo58Tx9/vnnAIDOzk6888472Lhxo7L//2Y8E+JDhw4dEkuWLBExMTGivLzc\nq4+9YcMGER4eLgICAoRGoxGvvvqquHjxosjMzJz2Mq/f/va3IiYmRsTHx4vq6mq5/OplXjExMeIX\nv/iFXD45OSnuv/9+eZlXe3v7jJmOHj0qVCqV0Ov18nK9w4cPK5rr5MmTwmg0Cr1eL5KTk8Xvfvc7\nIYRQ/Ln6MqvVKq8+UjLXp59+KvR6vdDr9SIpKUn+zCr9XLW0tIjU1FSRkpIi8vLyxODgoOKZRkdH\nxaJFi8Tw8LBcpnSmiooKeUnqww8/LOx2u+KZhBBixYoVQqvVCr1eL95//33FnyuPb14jIqLvF37y\nGhERSSwFIiKSWApERCSxFIiISGIpEBGRxFIgIiLpfwGexzRJjCWLBQAAAABJRU5ErkJggg==\n", "text": [ - "" + "" ] } ], - "prompt_number": 7 + "prompt_number": 53 }, { "cell_type": "code", @@ -201,7 +172,7 @@ "language": "python", "metadata": {}, "outputs": [], - "prompt_number": 8 + "prompt_number": 54 }, { "cell_type": "code", @@ -221,7 +192,7 @@ "==================== checkDerivative ====================\n", "iter\th\t\t|J0-Jt|\t\t|J0+h*dJ'*dx-Jt|\tOrder\n", "---------------------------------------------------------\n", - "0\t1.00e-01\t7.240e-08\t\t1.904e-22\t\tnan" + "0\t1.00e-01\t4.299e-09\t\t1.152e-23\t\tnan" ] }, { @@ -229,7 +200,7 @@ "stream": "stdout", "text": [ "\n", - "1\t1.00e-02\t7.240e-09\t\t1.009e-22\t\t0.276" + "1\t1.00e-02\t4.299e-10\t\t7.723e-24\t\t0.174" ] }, { @@ -237,7 +208,7 @@ "stream": "stdout", "text": [ "\n", - "2\t1.00e-03\t7.240e-10\t\t2.166e-22\t\t-0.332" + "2\t1.00e-03\t4.299e-11\t\t1.170e-23\t\t-0.180" ] }, { @@ -245,7 +216,7 @@ "stream": "stdout", "text": [ "\n", - "3\t1.00e-04\t7.240e-11\t\t9.316e-23\t\t0.366" + "3\t1.00e-04\t4.299e-12\t\t8.749e-24\t\t0.126" ] }, { @@ -253,7 +224,7 @@ "stream": "stdout", "text": [ "\n", - "4\t1.00e-05\t7.240e-12\t\t2.209e-22\t\t-0.375" + "4\t1.00e-05\t4.299e-13\t\t1.107e-23\t\t-0.102" ] }, { @@ -262,35 +233,72 @@ "text": [ "\n", "========================= PASS! =========================\n", - "You are awesome.\n", + "Once upon a time, a happy little test passed.\n", "\n" ] }, { "metadata": {}, "output_type": "pyout", - "prompt_number": 9, + "prompt_number": 55, "text": [ "True" ] } ], - "prompt_number": 9 + "prompt_number": 55 }, { "cell_type": "code", "collapsed": false, "input": [ "def derBbc(chi):\n", - " Bbc, const = CongruousMagBC(mesh, np.array([1., 0., 0.]), chi)\n", - " bc = Utils.sdiag(prob.mesh.vol)*prob.mesh.faceDiv*prob._Pout.T*Bbc \n", - " dBbc = lambda x: Utils.sdiag(prob.mesh.vol)*prob.mesh.faceDiv*prob._Pout.T*const*Bbc*(np.inner(prob.mesh.vol, x))\n", + " Bbc, const = CongruousMagBC(mesh, np.array([1., 1., 0.]), chi)\n", + " bc = Utils.sdiag(prob.mesh.vol)*prob.mesh.faceDiv*prob._Pout.T*Bbc \n", + " print const, np.linalg.norm(Bbc)\n", + " dBbc = lambda x: (Utils.sdiag(prob.mesh.vol)*(prob.mesh.faceDiv*prob._Pout.T*const*Bbc))*(np.inner(prob.mesh.vol, x)) \n", " return bc, dBbc" ], "language": "python", "metadata": {}, "outputs": [], - "prompt_number": 8 + "prompt_number": 56 + }, + { + "cell_type": "code", + "collapsed": false, + "input": [ + "a, b = derBbc(chi)\n", + "print np.linalg.norm(a), np.linalg.norm(b(chi))\n", + "print np.linalg.norm(chi)\n", + "print np.array([1., 1., 0.])\n", + "print mesh" + ], + "language": "python", + "metadata": {}, + "outputs": [ + { + "output_type": "stream", + "stream": "stdout", + "text": [ + "7.65830184557e-07 0.0659506783728\n", + "55.718105908 55.6914561973\n", + "1.6089045963\n", + "[ 1. 1. 0.]\n", + " ---- 3-D TensorMesh ---- \n", + " x0: -488.58\n", + " y0: -488.58\n", + " z0: -476.08\n", + " nCx: 31\n", + " nCy: 31\n", + " nCz: 30\n", + " hx: 71.40, 54.93, 42.25, 32.50, 23*25.00, 32.50, 42.25, 54.93, 71.40\n", + " hy: 71.40, 54.93, 42.25, 32.50, 23*25.00, 32.50, 42.25, 54.93, 71.40\n", + " hz: 71.40, 54.93, 42.25, 32.50, 22*25.00, 32.50, 42.25, 54.93, 71.40\n" + ] + } + ], + "prompt_number": 57 }, { "cell_type": "code", @@ -310,85 +318,79 @@ "==================== checkDerivative ====================\n", "iter\th\t\t|J0-Jt|\t\t|J0+h*dJ'*dx-Jt|\tOrder\n", "---------------------------------------------------------\n", - "0\t1.00e-01\t2.212e+00\t\t8.604e-05\t\tnan" + "7.65830184557e-07" ] }, { "output_type": "stream", "stream": "stdout", "text": [ - "\n", - "1\t1.00e-02\t2.212e-01\t\t8.604e-07\t\t2.000" + " 0.0659506783728\n", + "7.0907489096e-07" ] }, { "output_type": "stream", "stream": "stdout", "text": [ - "\n", - "2\t1.00e-03\t2.212e-02\t\t8.574e-09\t\t2.001" + " 0.0712240073538\n", + "0\t1.00e-01\t4.455e+00\t\t1.705e-04\t\tnan\n", + "7.59749260143e-07" ] }, { "output_type": "stream", "stream": "stdout", "text": [ - "\n", - "3\t1.00e-04\t2.212e-03\t\t9.857e-11\t\t1.939" + " 0.0664780294308\n", + "1\t1.00e-02\t4.455e-01\t\t1.705e-06\t\t2.000\n", + "7.65217717349e-07" ] }, { "output_type": "stream", "stream": "stdout", "text": [ - "\n", - "4\t1.00e-05\t2.212e-04\t\t1.424e-11\t\t0.840" + " 0.0660034136602\n", + "2\t1.00e-03\t4.455e-02\t\t1.706e-08\t\t2.000\n", + "7.65768893742e-07" ] }, { "output_type": "stream", "stream": "stdout", "text": [ - "\n", + " 0.0659559519034\n", + "3\t1.00e-04\t4.455e-03\t\t1.749e-10\t\t1.989\n", + "7.65824055034e-07" + ] + }, + { + "output_type": "stream", + "stream": "stdout", + "text": [ + " 0.0659512057259\n", + "4\t1.00e-05\t4.455e-04\t\t2.506e-12\t\t1.844\n", "========================= PASS! =========================\n", - "Awesome, Seogi, just awesome.\n", + "The test be workin!\n", "\n" ] }, { "metadata": {}, "output_type": "pyout", - "prompt_number": 9, + "prompt_number": 58, "text": [ "True" ] } ], - "prompt_number": 9 + "prompt_number": 58 }, { "cell_type": "code", "collapsed": false, - "input": [ - "##################\n", - "# Test J\n", - "##################\n", - "\n", - "d_chi = 0.8*chi #np.random.rand(mesh.nCz)\n", - "d_sph_ind = spheremodel(mesh, 0., 0., -100., 50)\n", - "d_chi[d_sph_ind] = 0.02\n", - "\n", - "from SimPEG.Tests import checkDerivative\n", - "\n", - "derChk = lambda m: [prob.data.dpred(m), lambda mx: -prob.Jvec(chi, mx)]\n", - "print '\\n'\n", - "passed = checkDerivative(derChk, chi, plotIt=False, dx=d_chi, num=2)\n", - "\n", - "\n", - "# plt.pcolor(X, Y, dpred.reshape(X.shape, order='F'))\n", - "\n", - "# plt.show()" - ], + "input": [], "language": "python", "metadata": {}, "outputs": []