from SimPEG import * import BaseMag as MAG from scipy.constants import mu_0 from MagAnalytics import spheremodel, CongruousMagBC import re class Problem3D_Integral(Problem.BaseProblem): #surveyPair = Survey.LinearSurvey forwardOnly = False #: Store the forward matrix by default, otherwise just compute d actInd = None #: Active cell indices provided M = None #: Magnetization matrix provided, otherwise all induced rtype = 'tmi' #: Receiver type either "tmi" | "xyz" def __init__(self, mesh, mapping=None, **kwargs): Problem.BaseProblem.__init__(self, mesh, mapping=mapping, **kwargs) def fwr_ind(self): # Add forward function # kappa = self.curModel.kappa TODO kappa = self.mapping*self.curModel if self.forwardOnly: if getattr(self, 'actInd', None) is not None: if self.actInd.dtype=='bool': inds = np.asarray([inds for inds, elem in enumerate(self.actInd, 1) if elem], dtype = int) - 1 else: inds = self.actInd else: inds = np.asarray(range(self.mesh.nC)) nC = len(inds) # Create active cell projector P = sp.csr_matrix((np.ones(nC),(inds, range(nC))), shape=(self.mesh.nC, nC)) # Create vectors of nodal location (lower and upper coners for each cell) xn = self.mesh.vectorNx; yn = self.mesh.vectorNy; zn = self.mesh.vectorNz; yn2,xn2,zn2 = np.meshgrid(yn[1:], xn[1:], zn[1:]) yn1,xn1,zn1 = np.meshgrid(yn[0:-1], xn[0:-1], zn[0:-1]) Yn = P.T*np.c_[mkvc(yn1), mkvc(yn2)] Xn = P.T*np.c_[mkvc(xn1), mkvc(xn2)] Zn = P.T*np.c_[mkvc(zn1), mkvc(zn2)] survey = self.survey rxLoc = survey.srcField.rxList[0].locs ndata = rxLoc.shape[0] # Loop through all observations and create forward operator (ndata-by-nC) print "Begin calculation forward calculations... G not stored: " # If assumes uniform magnetization direction if getattr(self, 'M', None) is None: M = dipazm_2_xyz(np.ones(nC) * survey.srcField.param[1],np.ones(nC) * survey.srcField.param[2]) Mx = Utils.sdiag(M[:,0]*survey.srcField.param[0]) My = Utils.sdiag(M[:,1]*survey.srcField.param[0]) Mz = Utils.sdiag(M[:,2]*survey.srcField.param[0]) Mxyz = sp.vstack((Mx,My,Mz)) if self.rtype == 'tmi': # Convert Bdecination from north to cartesian D = (450.-float(survey.srcField.param[2]))%360. I = survey.srcField.param[1] # Projection matrix Ptmi = mkvc(np.r_[np.cos(np.deg2rad(I))*np.cos(np.deg2rad(D)), np.cos(np.deg2rad(I))*np.sin(np.deg2rad(D)), np.sin(np.deg2rad(I))],2).T fwr_d = np.zeros(self.survey.nRx) else: fwr_d = np.zeros(3*self.survey.nRx) # Add counter to dsiplay progress. Good for large problems count = -1; for ii in range(ndata): tx, ty, tz = get_T_mat(Xn,Yn,Zn,rxLoc[ii,:]) if self.rtype =='tmi': fwr_d[ii] = (Ptmi.dot(np.vstack((tx,ty,tz)))*Mxyz).dot(kappa) elif self.rtype =='xyz': fwr_d[ii] = (tx*Mxyz).dot(kappa) fwr_d[ii+ndata] = (ty*Mxyz).dot(kappa) fwr_d[ii+2*ndata] = (tz*Mxyz).dot(kappa) # Display progress count = progress(ii,count,ndata) print "Done 100% ...forward operator completed!!\n" return fwr_d else: return self.G.dot(kappa) def fwr_rem(self): #TODO check if we are inverting for M return self.G.dot(self.mapping(m)) def fields(self, m, **kwargs): self.curModel = m if self.rtype == 'tmi': total = np.zeros(self.survey.nRx) else: total = np.zeros(3*self.survey.nRx) induced = self.fwr_ind() # rem = self.rem if induced is not None: total += induced return total # return self.G.dot(self.mapping*(m)) def Jvec(self, m, v, f=None): dmudm = self.mapping.deriv(m) return self.G.dot(dmudm*v) def Jtvec(self, m, v, f=None): dmudm = self.mapping.deriv(m) return dmudm.T * (self.G.T.dot(v)) @property def G(self): if not self.ispaired: raise Exception('Need to pair!') if getattr(self,'_G', None) is None: self._G = self.Intrgl_Fwr_Op() return self._G # @property # def Grem(self): # if not self.ispaired: # raise Exception('Need to pair!') # if getattr(self,'_Grem', None) is None: # self._Grem = Intrgl_Fwr_Op('full') # return self._Grem def Intrgl_Fwr_Op(self, Magnetization="ind"): """ Magnetic forward operator in integral form flag = 'ind' | 'full' 1- ind : Magnetization fixed by user 3- full: Full tensor matrix stored with shape([3*ndata, 3*nc]) Return _G = Linear forward modeling operation Created on March, 13th 2016 @author: dominiquef """ # Find non-zero cells #inds = np.nonzero(actv)[0] if getattr(self, 'actInd', None) is not None: if self.actInd.dtype=='bool': inds = np.asarray([inds for inds, elem in enumerate(self.actInd, 1) if elem], dtype = int) - 1 else: inds = self.actInd else: inds = np.asarray(range(self.mesh.nC)) nC = len(inds) # Create active cell projector P = sp.csr_matrix((np.ones(nC),(inds, range(nC))), shape=(self.mesh.nC, nC)) # Create vectors of nodal location (lower and upper coners for each cell) xn = self.mesh.vectorNx; yn = self.mesh.vectorNy; zn = self.mesh.vectorNz; yn2,xn2,zn2 = np.meshgrid(yn[1:], xn[1:], zn[1:]) yn1,xn1,zn1 = np.meshgrid(yn[0:-1], xn[0:-1], zn[0:-1]) Yn = P.T*np.c_[mkvc(yn1), mkvc(yn2)] Xn = P.T*np.c_[mkvc(xn1), mkvc(xn2)] Zn = P.T*np.c_[mkvc(zn1), mkvc(zn2)] survey = self.survey rxLoc = survey.srcField.rxList[0].locs ndata = rxLoc.shape[0] # Pre-allocate space and create magnetization matrix if required if (Magnetization=='ind'): # # If assumes uniform magnetization direction # if M.shape != (nC,3): # print 'Magnetization vector must be Nc x 3' # return if getattr(self, 'M', None) is None: M = dipazm_2_xyz(np.ones(nC) * survey.srcField.param[1],np.ones(nC) * survey.srcField.param[2]) Mx = Utils.sdiag(M[:,0]*survey.srcField.param[0]) My = Utils.sdiag(M[:,1]*survey.srcField.param[0]) Mz = Utils.sdiag(M[:,2]*survey.srcField.param[0]) Mxyz = sp.vstack((Mx,My,Mz)) if survey.srcField.rxList[0].rxType == 'tmi': G = np.zeros((ndata, nC)) # Convert Bdecination from north to cartesian D = (450.-float(survey.srcField.param[2]))%360. I = survey.srcField.param[1] # Projection matrix Ptmi = mkvc(np.r_[np.cos(np.deg2rad(I))*np.cos(np.deg2rad(D)), np.cos(np.deg2rad(I))*np.sin(np.deg2rad(D)), np.sin(np.deg2rad(I))],2).T; elif survey.srcField.rxList[0].rxType == 'xyz': G = np.zeros((int(3*ndata), nC)) elif Magnetization == 'full': G = np.zeros((int(3*ndata), int(3*nC))) else: print """Flag must be either 'ind' | 'full', please revised""" return # Loop through all observations and create forward operator (ndata-by-nC) print "Begin calculation of forward operator: " + Magnetization # Add counter to dsiplay progress. Good for large problems count = -1; for ii in range(ndata): tx, ty, tz = get_T_mat(Xn,Yn,Zn,rxLoc[ii,:]) if Magnetization == 'ind': if survey.srcField.rxList[0].rxType =='tmi': G[ii,:] = Ptmi.dot(np.vstack((tx,ty,tz)))*Mxyz elif survey.srcField.rxList[0].rxType =='xyz': G[ii,:] = tx*Mxyz G[ii+ndata,:] = ty*Mxyz G[ii+2*ndata,:] = tz*Mxyz elif Magnetization == 'full': G[ii,:] = tx G[ii+ndata,:] = ty G[ii+2*ndata,:] = tz # Display progress count = progress(ii,count,ndata) print "Done 100% ...forward operator completed!!\n" return G class Problem3D_DiffSecondary(Problem.BaseProblem): """ Secondary field approach using differential equations! """ surveyPair = MAG.BaseMagSurvey modelPair = MAG.BaseMagMap def __init__(self, model, mapping=None, **kwargs): Problem.BaseProblem.__init__(self, model, mapping=mapping, **kwargs) Pbc, Pin, self._Pout = \ self.mesh.getBCProjWF('neumann', discretization='CC') Dface = self.mesh.faceDiv Mc = Utils.sdiag(self.mesh.vol) self._Div = Mc*Dface*Pin.T*Pin @property def MfMuI(self): return self._MfMuI @property def MfMui(self): return self._MfMui @property def MfMu0(self): return self._MfMu0 def makeMassMatrices(self, m): mu = self.mapping*m self._MfMui = self.mesh.getFaceInnerProduct(1./mu)/self.mesh.dim # self._MfMui = self.mesh.getFaceInnerProduct(1./mu) #TODO: this will break if tensor mu self._MfMuI = Utils.sdiag(1./self._MfMui.diagonal()) self._MfMu0 = self.mesh.getFaceInnerProduct(1./mu_0)/self.mesh.dim # self._MfMu0 = self.mesh.getFaceInnerProduct(1/mu_0) @Utils.requires('survey') def getB0(self): b0 = self.survey.B0 B0 = np.r_[b0[0]*np.ones(self.mesh.nFx), b0[1]*np.ones(self.mesh.nFy), b0[2]*np.ones(self.mesh.nFz)] 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) mu = self.mapping*m chi = mu/mu_0-1 #temporary fix Bbc, Bbc_const = CongruousMagBC(self.mesh, self.survey.B0, chi) self.Bbc = Bbc self.Bbc_const = Bbc_const # return self._Div*self.MfMuI*self.MfMu0*B0 - self._Div*B0 + Mc*Dface*self._Pout.T*Bbc return self._Div*self.MfMuI*self.MfMu0*B0 - self._Div*B0 def getA(self, m): """ GetA creates and returns the A matrix for the Magnetics problem The A matrix has the form: .. 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 return {'B': B, 'u': u} @Utils.timeIt def Jvec(self, m, v, u=None): """ Computing Jacobian multiplied by vector By setting our problem as .. math :: \mathbf{C}(\mathbf{m}, \mathbf{u}) = \mathbf{A}\mathbf{u} - \mathbf{rhs} = 0 And taking derivative w.r.t m .. math :: \\nabla \mathbf{C}(\mathbf{m}, \mathbf{u}) = \\nabla_m \mathbf{C}(\mathbf{m}) \delta \mathbf{m} + \\nabla_u \mathbf{C}(\mathbf{u}) \delta \mathbf{u} = 0 \\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} \\nabla_m \mathbf{C}(\mathbf{m}) = \\frac{\partial \mathbf{A}}{\partial \mathbf{m}}(\mathbf{m})\mathbf{u} - \\frac{\partial \mathbf{rhs}(\mathbf{m})}{\partial \mathbf{m}} .. math :: \\frac{\partial \mathbf{A}}{\partial \mathbf{m}}(\mathbf{m})\mathbf{u} = \\frac{\partial \mathbf{\mu}}{\partial \mathbf{m}} \left[\Div \diag (\Div^T \mathbf{u}) \dMfMuI \\right] \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}} 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) B, u = u['B'], u['u'] mu = self.mapping*(m) dmudm = self.mapping.deriv(m) dchidmu = Utils.sdiag(1/mu_0*np.ones(self.mesh.nC)) vol = self.mesh.vol Div = self._Div Dface = self.mesh.faceDiv P = self.survey.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) 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, dchidmu*dmudm*v) #temporary fix # dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v dCdm_RHSv = dCdm_RHS1*(dmudm*v) 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-6, maxiter=1000, M=m1) if info > 0: print "Iterative solver did not work well (Jvec)" # 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, u=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 u is None: u = self.fields(m) B, u = u['B'], u['u'] mu = self.mapping*(m) dmudm = self.mapping.deriv(m) dchidmu = Utils.sdiag(1/mu_0*np.ones(self.mesh.nC)) vol = self.mesh.vol Div = self._Div Dface = self.mesh.faceDiv P = self.survey.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-6, maxiter=1000, M=m1) if info > 0: print "Iterative solver did not work well (Jtvec)" # 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_Atsol = ( dmudm.T * dMfMuI.T*( Utils.sdiag( Div.T * u ) * Div.T ) ) * sol # dCdm_RHS1 = Div * (Utils.sdiag( self.MfMu0*B0 ) * dMfMuI) # dCdm_RHS1tsol = (dMfMuI.T*( Utils.sdiag( self.MfMu0*B0 ) ) * Div.T * dmudm) * sol dCdm_RHS1tsol = ( dmudm.T * dMfMuI.T*( Utils.sdiag( self.MfMu0*B0 ) ) * Div.T ) * sol # temp1 = (Dface*(self._Pout.T*self.Bbc_const*self.Bbc)) temp1sol = ( Dface.T*( Utils.sdiag(vol)*sol ) ) temp2 = self.Bbc_const*(self._Pout.T*self.Bbc).T # dCdm_RHS2v = (Utils.sdiag(vol)*temp1)*np.inner(vol, dchidmu*dmudm*v) dCdm_RHS2tsol = (dmudm.T*dchidmu.T*vol)*np.inner(temp2, temp1sol) # dCdm_RHSv = dCdm_RHS1*(dmudm*v) + dCdm_RHS2v #temporary fix # dCdm_RHStsol = dCdm_RHS1tsol - dCdm_RHS2tsol dCdm_RHStsol = dCdm_RHS1tsol # 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) def MagneticsDiffSecondaryInv(mesh, model, data, **kwargs): """ Inversion module for MagneticsDiffSecondary """ from SimPEG import Optimization, Regularization, Parameters, ObjFunction, Inversion prob = MagneticsDiffSecondary(mesh, model) miter = kwargs.get('maxIter', 10) if prob.ispaired: prob.unpair() if data.ispaired: data.unpair() prob.pair(data) # Create an optimization program opt = Optimization.InexactGaussNewton(maxIter=miter) opt.bfgsH0 = Solver(sp.identity(model.nP),flag='D') # Create a regularization program reg = Regularization.Tikhonov(model) # Create an objective function beta = Parameters.BetaSchedule(beta0=1e0) obj = ObjFunction.BaseObjFunction(data, reg, beta=beta) # Create an inversion object inv = Inversion.BaseInversion(obj, opt) return inv, reg def Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,flag): """ Forward model magnetic data using integral equation INPUT: mesh = SimPEG.TensorMesh B = Inducing field parameter [Binc, Bdecl, B0] M = Magnetization matrix [Minc, Mdecl] -90:90, 0:360 rxLox = Observation location informat [obsx, obsy, obsz] model = Model associated with mesh actv = Active cells from topo (from getActiveTopo) flag = Data type "tmi" | "xyz" OUTPUT: dobs =Observation array in format [obsx, obsy, obsz, data] Created on Oct 7, 2015 @author: dominiquef """ if actv.dtype=='bool': inds = np.asarray([inds for inds, elem in enumerate(actv, 1) if elem], dtype = int) - 1 else: inds = actv nC = len(inds) P = sp.csr_matrix((np.ones(nC),(inds, range(nC))), shape=(mesh.nC, nC)) xn = mesh.vectorNx; yn = mesh.vectorNy; zn = mesh.vectorNz; yn2,xn2,zn2 = np.meshgrid(yn[1:], xn[1:], zn[1:]) yn1,xn1,zn1 = np.meshgrid(yn[0:-1], xn[0:-1], zn[0:-1]) Yn = P.T*np.c_[mkvc(yn1), mkvc(yn2)] Xn = P.T*np.c_[mkvc(xn1), mkvc(xn2)] Zn = P.T*np.c_[mkvc(zn1), mkvc(zn2)] nC = len(inds) ndata = rxLoc.shape[0] # Convert declination from north to cartesian Md = (450.-float(M[1]))%360. # Create magnetization matrix mx = np.cos(np.deg2rad(M[0])) * np.cos(np.deg2rad(Md)) my = np.cos(np.deg2rad(M[0])) * np.sin(np.deg2rad(Md)) mz = np.sin(np.deg2rad(M[0])) Mx = Utils.sdiag(np.ones([nC])*mx*B[2]) My = Utils.sdiag(np.ones([nC])*my*B[2]) Mz = Utils.sdiag(np.ones([nC])*mz*B[2]) #matplotlib.pyplot.spy(scipy.sparse.csr_matrix(Mx)) #plt.show() Mxyz = sp.vstack((Mx,My,Mz)); #%% Create TMI projector # Convert Bdecination from north to cartesian D = (450.-float(B[1]))%360. if flag=='tmi': Ptmi = mkvc(np.r_[np.cos(np.deg2rad(B[0]))*np.cos(np.deg2rad(D)),np.cos(np.deg2rad(B[0]))*np.sin(np.deg2rad(D)),np.sin(np.deg2rad(B[0]))],2).T; d = np.zeros(ndata) elif flag=='xyz': d = np.zeros(int(3*ndata)) # Loop through all observations and create forward operator (ndata-by-nC) print "Begin forward modeling " +str(int(ndata)) + " data points..." # Add counter to dsiplay progress. count = -1 for ii in range(ndata): tx, ty, tz = get_T_mat(Xn,Yn,Zn,rxLoc[ii,:]) Gxyz = np.vstack((tx,ty,tz))*Mxyz # Remove non-active cells if flag=='xyz': d[ii::ndata] = mkvc(Gxyz.dot(P.T*model)) elif flag=='tmi': d[ii] = Ptmi.dot(Gxyz.dot(P.T*model)) # Display progress count = progress(ii,count,ndata) print "Done 100% ...forward modeling completed!!\n" return d #def Intrgl_Fwr_Op(mesh,B,M,rxLoc,actv,flag): # """ # # Magnetic forward operator in integral form # # INPUT: # mesh = Mesh in SimPEG format # B = Inducing field parameter [Binc, Bdecl, B0] # M = Magnetization information # [OPTIONS] # 1- [Minc, Mdecl] : Assumes uniform magnetization orientation # 2- [mx1,mx2,..., my1,...,mz1] : cell-based defined magnetization direction # 3- diag(M): Block diagonal matrix with [Mx, My, Mz] along the diagonal # # rxLox = Observation location informat [obsx, obsy, obsz] # # flag = 'tmi' | 'xyz' | 'full' # [OPTIONS] # 1- tmi : Magnetization direction used and data are projected onto the # inducing field direction F.shape([ndata, nc]) # # 2- xyz : Magnetization direction used and data are given in 3-components # F.shape([3*ndata, nc]) # # 3- full: Full tensor matrix stored with shape([3*ndata, 3*nc]) # # OUTPUT: # F = Linear forward modeling operation # # Created on Dec, 20th 2015 # # @author: dominiquef # # """ # # Find non-zero cells # #inds = np.nonzero(actv)[0] # if actv.dtype=='bool': # inds = np.asarray([inds for inds, elem in enumerate(actv, 1) if elem], dtype = int) - 1 # else: # inds = actv # # nC = len(inds) # # # Create active cell projector # P = sp.csr_matrix((np.ones(nC),(inds, range(nC))), # shape=(mesh.nC, nC)) # # # Create vectors of nodal location (lower and upper coners for each cell) # xn = mesh.vectorNx; # yn = mesh.vectorNy; # zn = mesh.vectorNz; # # yn2,xn2,zn2 = np.meshgrid(yn[1:], xn[1:], zn[1:]) # yn1,xn1,zn1 = np.meshgrid(yn[0:-1], xn[0:-1], zn[0:-1]) # # Yn = P.T*np.c_[mkvc(yn1), mkvc(yn2)] # Xn = P.T*np.c_[mkvc(xn1), mkvc(xn2)] # Zn = P.T*np.c_[mkvc(zn1), mkvc(zn2)] # # ndata = rxLoc.shape[0] # # # Convert Bdecination from north to cartesian # D = (450.-float(B[1]))%360. # # # # Pre-allocate space and create magnetization matrix if required # if (flag=='tmi') | (flag == 'xyz'): # # If assumes uniform magnetization direction # if M.shape != (nC,3): # # print 'Magnetization vector must be Nc x 3' # return # # # Mx = Utils.sdiag(M[:,0]*B[2]) # My = Utils.sdiag(M[:,1]*B[2]) # Mz = Utils.sdiag(M[:,2]*B[2]) # # Mxyz = sp.vstack((Mx,My,Mz)) # # # # if flag == 'tmi': # F = np.zeros((ndata, nC)) # # # Projection matrix # Ptmi = mkvc(np.r_[np.cos(np.deg2rad(B[0]))*np.cos(np.deg2rad(D)), # np.cos(np.deg2rad(B[0]))*np.sin(np.deg2rad(D)), # np.sin(np.deg2rad(B[0]))],2).T; # # elif flag == 'xyz': # # F = np.zeros((int(3*ndata), nC)) # # elif flag == 'full': # F = np.zeros((int(3*ndata), int(3*nC))) # # # else: # print """Flag must be either 'tmi' | 'xyz' | 'full', please revised""" # return # # # # Loop through all observations and create forward operator (ndata-by-nC) # print "Begin calculation of forward operator: " + flag # # # Add counter to dsiplay progress. Good for large problems # count = -1; # for ii in range(ndata): # # # tx, ty, tz = get_T_mat(Xn,Yn,Zn,rxLoc[ii,:]) # # if flag=='tmi': # F[ii,:] = Ptmi.dot(np.vstack((tx,ty,tz)))*Mxyz # # elif flag == 'xyz': # F[ii,:] = tx*Mxyz # F[ii+ndata,:] = ty*Mxyz # F[ii+2*ndata,:] = tz*Mxyz # # elif flag == 'full': # F[ii,:] = tx # F[ii+ndata,:] = ty # F[ii+2*ndata,:] = tz # # # # Display progress # count = progress(ii,count,ndata) # # print "Done 100% ...forward operator completed!!\n" # # return F def get_T_mat(Xn, Yn, Zn, rxLoc): """ Load in the active nodes of a tensor mesh and computes the magnetic tensor for a given observation location rxLoc[obsx, obsy, obsz] INPUT: Xn, Yn, Zn: Node location matrix for the lower and upper most corners of all cells in the mesh shape[nC,2] M OUTPUT: Tx = [Txx Txy Txz] Ty = [Tyx Tyy Tyz] Tz = [Tzx Tzy Tzz] where each elements have dimension 1-by-nC. Only the upper half 5 elements have to be computed since symetric. Currently done as for-loops but will eventually be changed to vector indexing, once the topography has been figured out. Created on Oct, 20th 2015 @author: dominiquef """ eps = 1e-10 # add a small value to the locations to avoid /0 nC = Xn.shape[0] # Pre-allocate space for 1D array Tx = np.zeros((1,3*nC)) Ty = np.zeros((1,3*nC)) Tz = np.zeros((1,3*nC)) dz2 = rxLoc[2] - Zn[:,0] + eps dz1 = rxLoc[2] - Zn[:,1] + eps dy2 = Yn[:,1] - rxLoc[1] + eps dy1 = Yn[:,0] - rxLoc[1] + eps dx2 = Xn[:,1] - rxLoc[0] + eps dx1 = Xn[:,0] - rxLoc[0] + eps R1 = ( dy2**2 + dx2**2 ) R2 = ( dy2**2 + dx1**2 ) R3 = ( dy1**2 + dx2**2 ) R4 = ( dy1**2 + dx1**2 ) arg1 = np.sqrt( dz2**2 + R2 ) arg2 = np.sqrt( dz2**2 + R1 ) arg3 = np.sqrt( dz1**2 + R1 ) arg4 = np.sqrt( dz1**2 + R2 ) arg5 = np.sqrt( dz2**2 + R3 ) arg6 = np.sqrt( dz2**2 + R4 ) arg7 = np.sqrt( dz1**2 + R4 ) arg8 = np.sqrt( dz1**2 + R3 ) Tx[0,0:nC] = np.arctan2( dy1 * dz2 , ( dx2 * arg5 ) ) +\ - np.arctan2( dy2 * dz2 , ( dx2 * arg2 ) ) +\ np.arctan2( dy2 * dz1 , ( dx2 * arg3 ) ) +\ - np.arctan2( dy1 * dz1 , ( dx2 * arg8 ) ) +\ np.arctan2( dy2 * dz2 , ( dx1 * arg1 ) ) +\ - np.arctan2( dy1 * dz2 , ( dx1 * arg6 ) ) +\ np.arctan2( dy1 * dz1 , ( dx1 * arg7 ) ) +\ - np.arctan2( dy2 * dz1 , ( dx1 * arg4 ) ); Ty[0,0:nC] = np.log( ( dz2 + arg2 ) / (dz1 + arg3 ) ) +\ -np.log( ( dz2 + arg1 ) / (dz1 + arg4 ) ) +\ np.log( ( dz2 + arg6 ) / (dz1 + arg7 ) ) +\ -np.log( ( dz2 + arg5 ) / (dz1 + arg8 ) ); Ty[0,nC:2*nC] = np.arctan2( dx1 * dz2 , ( dy2 * arg1 ) ) +\ - np.arctan2( dx2 * dz2 , ( dy2 * arg2 ) ) +\ np.arctan2( dx2 * dz1 , ( dy2 * arg3 ) ) +\ - np.arctan2( dx1 * dz1 , ( dy2 * arg4 ) ) +\ np.arctan2( dx2 * dz2 , ( dy1 * arg5 ) ) +\ - np.arctan2( dx1 * dz2 , ( dy1 * arg6 ) ) +\ np.arctan2( dx1 * dz1 , ( dy1 * arg7 ) ) +\ - np.arctan2( dx2 * dz1 , ( dy1 * arg8 ) ); R1 = (dy2**2 + dz1**2); R2 = (dy2**2 + dz2**2); R3 = (dy1**2 + dz1**2); R4 = (dy1**2 + dz2**2); Ty[0,2*nC:] = np.log( ( dx1 + np.sqrt( dx1**2 + R1 ) ) / (dx2 + np.sqrt( dx2**2 + R1 ) ) ) +\ -np.log( ( dx1 + np.sqrt( dx1**2 + R2 ) ) / (dx2 + np.sqrt( dx2**2 + R2 ) ) ) +\ np.log( ( dx1 + np.sqrt( dx1**2 + R4 ) ) / (dx2 + np.sqrt( dx2**2 + R4 ) ) ) +\ -np.log( ( dx1 + np.sqrt( dx1**2 + R3 ) ) / (dx2 + np.sqrt( dx2**2 + R3 ) ) ); R1 = (dx2**2 + dz1**2); R2 = (dx2**2 + dz2**2); R3 = (dx1**2 + dz1**2); R4 = (dx1**2 + dz2**2); Tx[0,2*nC:] = np.log( ( dy1 + np.sqrt( dy1**2 + R1 ) ) / (dy2 + np.sqrt( dy2**2 + R1 ) ) ) +\ -np.log( ( dy1 + np.sqrt( dy1**2 + R2 ) ) / (dy2 + np.sqrt( dy2**2 + R2 ) ) ) +\ np.log( ( dy1 + np.sqrt( dy1**2 + R4 ) ) / (dy2 + np.sqrt( dy2**2 + R4 ) ) ) +\ -np.log( ( dy1 + np.sqrt( dy1**2 + R3 ) ) / (dy2 + np.sqrt( dy2**2 + R3 ) ) ); Tz[0,2*nC:] = -( Ty[0,nC:2*nC] + Tx[0,0:nC] ); Tz[0,nC:2*nC] = Ty[0,2*nC:]; Tx[0,nC:2*nC] = Ty[0,0:nC]; Tz[0,0:nC] = Tx[0,2*nC:]; Tx = Tx/(4*np.pi); Ty = Ty/(4*np.pi); Tz = Tz/(4*np.pi); return Tx,Ty,Tz def progress(iter, prog, final): """ progress(iter,prog,final) Function measuring the progress of a process and print to screen the %. Useful to estimate the remaining runtime of a large problem. Created on Dec, 20th 2015 @author: dominiquef """ arg = np.floor(float(iter)/float(final)*10.); if arg > prog: strg = "Done " + str(arg*10) + " %" print strg prog = arg; return prog def dipazm_2_xyz(dip, azm_N): """ dipazm_2_xyz(dip,azm_N) Function converting degree angles for dip and azimuth from north to a 3-components in cartesian coordinates. INPUT dip : Value or vector of dip from horizontal in DEGREE azm_N : Value or vector of azimuth from north in DEGREE OUTPUT M : [n-by-3] Array of xyz components of a unit vector in cartesian Created on Dec, 20th 2015 @author: dominiquef """ nC = len(azm_N) M = np.zeros((nC,3)) # Modify azimuth from North to Cartesian-X azm_X = (450.- np.asarray(azm_N)) % 360. D = np.deg2rad(np.asarray(dip)) I = np.deg2rad(azm_X) M[:,0] = np.cos(D) * np.cos(I) ; M[:,1] = np.cos(D) * np.sin(I) ; M[:,2] = np.sin(D) ; return M def get_dist_wgt(mesh, rxLoc, actv, R, R0): """ get_dist_wgt(xn,yn,zn,rxLoc,R,R0) Function creating a distance weighting function required for the magnetic inverse problem. INPUT xn, yn, zn : Node location rxLoc : Observation locations [obsx, obsy, obsz] actv : Active cell vector [0:air , 1: ground] R : Decay factor (mag=3, grav =2) R0 : Small factor added (default=dx/4) OUTPUT wr : [nC] Vector of distance weighting Created on Dec, 20th 2015 @author: dominiquef """ # Find non-zero cells if actv.dtype=='bool': inds = np.asarray([inds for inds, elem in enumerate(actv, 1) if elem], dtype=int) - 1 else: inds = actv nC = len(inds) # Create active cell projector P = sp.csr_matrix((np.ones(nC),(inds, range(nC))), shape=(mesh.nC, nC)) # Geometrical constant p = 1/np.sqrt(3); # Create cell center location Ym,Xm,Zm = np.meshgrid(mesh.vectorCCy, mesh.vectorCCx, mesh.vectorCCz) hY,hX,hZ = np.meshgrid(mesh.hy, mesh.hx, mesh.hz) # Rmove air cells Xm = P.T*mkvc(Xm) Ym = P.T*mkvc(Ym) Zm = P.T*mkvc(Zm) hX = P.T*mkvc(hX) hY = P.T*mkvc(hY) hZ = P.T*mkvc(hZ) V = P.T * mkvc(mesh.vol) wr = np.zeros(nC) ndata = rxLoc.shape[0] count = -1 print "Begin calculation of distance weighting for R= " + str(R) for dd in range(ndata): nx1 = (Xm - hX * p - rxLoc[dd,0])**2 nx2 = (Xm + hX * p - rxLoc[dd,0])**2 ny1 = (Ym - hY * p - rxLoc[dd,1])**2 ny2 = (Ym + hY * p - rxLoc[dd,1])**2 nz1 = (Zm - hZ * p - rxLoc[dd,2])**2 nz2 = (Zm + hZ * p - rxLoc[dd,2])**2 R1 = np.sqrt(nx1 + ny1 + nz1) R2 = np.sqrt(nx1 + ny1 + nz2) R3 = np.sqrt(nx2 + ny1 + nz1) R4 = np.sqrt(nx2 + ny1 + nz2) R5 = np.sqrt(nx1 + ny2 + nz1) R6 = np.sqrt(nx1 + ny2 + nz2) R7 = np.sqrt(nx2 + ny2 + nz1) R8 = np.sqrt(nx2 + ny2 + nz2) temp = (R1 + R0)**-R + (R2 + R0)**-R + (R3 + R0)**-R + (R4 + R0)**-R + (R5 + R0)**-R + (R6 + R0)**-R + (R7 + R0)**-R + (R8 + R0)**-R wr = wr + (V*temp/8.)**2. count = progress(dd,count,ndata) wr = np.sqrt(wr)/V wr = mkvc(wr) wr = np.sqrt(wr/(np.max(wr))) print "Done 100% ...distance weighting completed!!\n" return wr def writeUBCobs(filename, survey, d): """ writeUBCobs(filename,B,M,rxLoc,d,wd) Function writing an observation file in UBC-MAG3D format. INPUT filename : Name of out file including directory survey flag : dobs | dpred OUTPUT Obsfile Created on Dec, 27th 2015 @author: dominiquef """ B = survey.srcField.param rxLoc = survey.srcField.rxList[0].locs wd = survey.std data = np.c_[rxLoc , d , wd] with file(filename,'w') as fid: fid.write('%6.2f %6.2f %6.2f\n' %(B[1], B[2], B[0]) ) fid.write('%6.2f %6.2f %6.2f\n' %(B[1], B[2], 1) ) fid.write('%i\n' %len(d) ) np.savetxt(fid, data, fmt='%e',delimiter=' ',newline='\n') print "Observation file saved to: " + filename def getActiveTopo(mesh, topo, flag): """ getActiveTopo(mesh,topo) Function creates an active cell model from topography INPUT mesh : Mesh in SimPEG format topo : Scatter points defining topography [x,y,z] OUTPUT actv : Active cell model Created on Dec, 27th 2015 @founrdo """ import scipy.interpolate as interpolation print "Please remove this function! use SimPEG.Utils.surface2ind_topo(mesh, topo, gridLoc='CC')" if (flag=='N'): Zn = np.zeros((mesh.nNx,mesh.nNy)) # wght = np.zeros((mesh.nNx,mesh.nNy)) cx = mesh.vectorNx cy = mesh.vectorNy F = interpolation.NearestNDInterpolator(topo[:,0:2],topo[:,2]) [Y,X] = np.meshgrid(cy,cx) Zn = F(X,Y) actv = np.zeros((mesh.nCx, mesh.nCy, mesh.nCz)) if (flag=='N'): Nz = mesh.vectorNz[1:] for jj in range(mesh.nCy): for ii in range(mesh.nCx): temp = [kk for kk in range(len(Nz)) if np.all(Zn[ii:(ii+2),jj:(jj+2)] > Nz[kk]) ] actv[ii,jj,temp] = 1 actv = mkvc(actv==1) inds = np.asarray([inds for inds, elem in enumerate(actv, 1) if elem], dtype = int) - 1 return inds def plot_obs_2D(rxLoc,d=None ,varstr='Mag Obs', vmin=None, vmax=None, levels=None, fig=None): """ Function plot_obs(rxLoc,d) Generate a 2d interpolated plot from scatter points of data INPUT rxLoc : Observation locations [x,y,z] d : Data vector wd : Uncertainty vector OUTPUT figure() Created on Dec, 27th 2015 @author: dominiquef """ from scipy.interpolate import griddata import pylab as plt # Plot result if fig is None: fig = plt.figure() ax = plt.subplot() plt.scatter(rxLoc[:,0],rxLoc[:,1], c='k', s=10) if d is not None: if (vmin is None): vmin = d.min() if (vmax is None): vmax = d.max() # Create grid of points x = np.linspace(rxLoc[:,0].min(), rxLoc[:,0].max(), 100) y = np.linspace(rxLoc[:,1].min(), rxLoc[:,1].max(), 100) X, Y = np.meshgrid(x,y) # Interpolate d_grid = griddata(rxLoc[:,0:2],d,(X,Y), method ='linear') plt.imshow(d_grid, extent=[x.min(), x.max(), y.min(), y.max()], origin='lower', vmin=vmin, vmax=vmax, cmap="plasma") plt.colorbar(fraction=0.02) if levels is None: plt.contour(X,Y, d_grid, 10, vmin=vmin, vmax=vmax, cmap="plasma") else: plt.contour(X,Y, d_grid, levels=levels, colors='r', vmin=vmin, vmax=vmax, cmap="plasma") plt.title(varstr) plt.gca().set_aspect('equal', adjustable='box') return fig def read_MAGfwr_inp(input_file): """Read input files for forward modeling MAG data with integral form INPUT: input_file: File name containing the forward parameter OUTPUT: mshfile obsfile modfile magfile topofile # All files should be in the working directory, otherwise the path must # be specified. Created on Jul 17, 2013 @author: dominiquef """ fid = open(input_file,'r') line = fid.readline() l_input = line.split('!') mshfile = l_input[0].rstrip() line = fid.readline() l_input = line.split('!') obsfile = l_input[0].rstrip() line = fid.readline() l_input = line.split('!') modfile = l_input[0].rstrip() line = fid.readline() l_input = line.split('!') if l_input=='null': magfile = [] else: magfile = l_input[0].rstrip() line = fid.readline() l_input = line.split('!') if l_input=='null': topofile = [] else: topofile = l_input[0].rstrip() return mshfile, obsfile, modfile, magfile, topofile if __name__ == '__main__': import matplotlib.pyplot as plt hxind = ((5,25,1.3),(41, 12.5),(5,25,1.3)) hyind = ((5,25,1.3),(41, 12.5),(5,25,1.3)) hzind = ((5,25,1.3),(40, 12.5),(5,25,1.3)) hx, hy, hz = Utils.meshTensors(hxind, hyind, hzind) mesh = Mesh.TensorMesh([hx, hy, hz], [-hx.sum()/2,-hy.sum()/2,-hz.sum()/2]) chibkg = 0. chiblk = 0.01 chi = np.ones(mesh.nC)*chibkg sph_ind = spheremodel(mesh, 0., 0., 0., 100) chi[sph_ind] = chiblk model = BaseMag.BaseMagModel(mesh) # mu = (1.+chi)*mu_0 data = BaseMag.BaseMagData() 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) Z = np.ones((xr.size, yr.size))*150 rxLoc = np.c_[Utils.mkvc(X), Utils.mkvc(Y), Utils.mkvc(Z)] data.rxLoc = rxLoc prob = MagneticsDiffSecondary(mesh, model) prob.pair(data) 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()