Investigating problem with GNCG

simpegPF/Dev/Intgrl_MAG_Fwr_InterpData_Test.py

@@ -4,7 +4,7 @@ import matplotlib.pyplot as plt
 import simpegPF as PF
 import scipy.interpolate as interpolation
 import time
-from interpFFT import interpFFT
+#from interpFFT import interpFFT
 #from fwr_MAG_data import fwr_MAG_data
 
 import os

simpegPF/Dev/Intgrl_MAG_Inv_Driver.py

@@ -1,7 +1,8 @@
 import os
 
 #home_dir = 'C:\Users\dominiquef.MIRAGEOSCIENCE\Documents\GIT\SimPEG\simpegpf\simpegPF\Dev'
-home_dir = 'C:\\Users\\dominiquef.MIRAGEOSCIENCE\\ownCloud\\Research\\Modelling\\Synthetic\\Block_Gaussian_topo'
+#home_dir = 'C:\\Users\\dominiquef.MIRAGEOSCIENCE\\ownCloud\\Research\\Modelling\\Synthetic\\Block_Gaussian_topo'
+home_dir = 'C:\\Users\\dominiquef.MIRAGEOSCIENCE\\ownCloud\\Research\\Modelling\\Synthetic\\Nut_Cracker\\Induced_MAG3C'
 
 inpfile = 'PYMAG3D_inv.inp'
 
@@ -9,7 +10,7 @@ dsep = '\\'
 os.chdir(home_dir)
 
 #%%
-from SimPEG import np, Utils, mkvc
+from SimPEG import *
 import simpegPF as PF
 import pylab as plt
 
@@ -33,30 +34,36 @@ wd = dobs[:,4]
 
 ndata = rxLoc.shape[0]
 
+beta_in = 1e+2
 
 # Load in topofile or create flat surface
 if topofile == 'null':
     
-    Nx,Ny = np.meshgrid(mesh.vectorNx,mesh.vectorNy)
-    Nz = np.ones(Nx.shape) * mesh.vectorNz[-1]
-    
-    topo = np.c_[mkvc(Nx),mkvc(Ny),mkvc(Nz)]   
+    # All active
+    actv = np.ones(mesh.nC)
     
 else: 
+    
     topo = np.genfromtxt(topofile,skip_header=1)
-
-# Work with flat topogrphy for now
-actv = PF.Magnetics.getActiveTopo(mesh,topo,'N')
+    # Find the active cells
+    actv = PF.Magnetics.getActiveTopo(mesh,topo,'N')
 
 nC = int(sum(actv))
 
-# Load model file
+# Load starting model file
 if isinstance(mstart, float):
     mstart = np.ones(nC) * mstart
 else:
     mstart = Utils.meshutils.readUBCTensorModel(mstart,mesh)
     mstart = mstart[actv==1]
 
+# Load reference file
+if isinstance(mref, float):
+    mref = np.ones(nC) * mref
+else:
+    mref = Utils.meshutils.readUBCTensorModel(mref,mesh)
+    mref = mref[actv==1]
+    
 # Get magnetization vector for MOF
 if magfile=='DEFAULT':
     
@@ -73,10 +80,12 @@ midy = int(mesh.nCy/2)
 F = PF.Magnetics.Intrgl_Fwr_Op(mesh,B,M_xyz,rxLoc,actv,'tmi')
 
 # Get distance weighting function
-wr = PF.Magnetics.get_dist_wgt(mesh,rxLoc,3.,np.min(mesh.hx)/4)
+wr = PF.Magnetics.get_dist_wgt(mesh,rxLoc,actv,3.,np.min(mesh.hx)/4)
 wrMap = PF.BaseMag.WeightMap(mesh, wr)
 
-Utils.writeUBCTensorModel(home_dir+dsep+'wr.dat',mesh,wr)
+wr_out = np.zeros(mesh.nC)
+wr_out[actv==1] = wr
+Utils.writeUBCTensorModel(home_dir+dsep+'wr.dat',mesh,wr_out)
 
 # Write out the predicted
 pred = F.dot(mstart)
@@ -85,9 +94,56 @@ PF.Magnetics.writeUBCobs(home_dir + dsep + 'Pred.dat',B,M,rxLoc,pred,wd)
 #%%
 plt.figure()
 ax = plt.subplot()
-mesh.plotSlice(wr, ax = ax, normal = 'Y', ind=midx)
+mesh.plotSlice(wr_out, ax = ax, normal = 'Y', ind=midx)
 plt.title('Distance weighting')
 plt.xlabel('x');plt.ylabel('z')
 plt.gca().set_aspect('equal', adjustable='box')
 
+#%% Plot obs data
+PF.Magnetics.plot_obs_2D(rxLoc,d,wd,'Observed Data')
 
+#%% Run inversion
+prob = PF.Magnetics.MagneticIntegral(mesh, F)
+prob.solverOpts['accuracyTol'] = 1e-4
+survey = Survey.LinearSurvey()
+survey.pair(prob)
+#survey.makeSyntheticData(data, std=0.01)
+survey.dobs=d
+#survey.mtrue = model
+
+# Create pre-conditioner 
+diagA = np.sum(F**2.,axis=0) + beta_in*np.ones(nC)
+PC     = sp.spdiags(diagA**-1., 0, nC, nC);
+
+reg = Regularization.Tikhonov(mesh, mapping=wrMap)
+reg.mref = mref
+
+dmis = DataMisfit.l2_DataMisfit(survey)
+dmis.Wd = wd
+#opt = Optimization.ProjectedGNCG(maxIter=6,lower=-1.,upper=1.)
+opt.approxHinv = PC
+
+opt = Optimization.InexactGaussNewton(maxIter=6)
+invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = beta_in)
+beta = Directives.BetaSchedule(coolingFactor=8, coolingRate=2)
+#betaest = Directives.BetaEstimate_ByEig()
+target = Directives.TargetMisfit()
+inv = Inversion.BaseInversion(invProb, directiveList=[beta, target])
+reg.alpha_s =0.0025
+m0 = mstart
+
+# Run inversion
+mrec = inv.run(m0)
+
+m_out = np.ones(mesh.nC)
+m_out[actv==1] = mrec
+
+# Write result
+Utils.meshutils.writeUBCTensorModel('SimPEG_inv.sus',mesh,m_out)
+
+# Plot predicted
+pred = F.dot(mrec)
+PF.Magnetics.plot_obs_2D(rxLoc,pred,wd,'Predicted Data')
+PF.Magnetics.plot_obs_2D(rxLoc,(d-pred),wd,'Residual Data')
+
+print "Final misfit:" + str(np.sum( ((d-pred)/wd)**2. ) ) 

simpegPF/Magnetics.py

@@ -4,9 +4,22 @@ from scipy.constants import mu_0
 from MagAnalytics import spheremodel, CongruousMagBC
 
 class MagneticIntegral(Problem.BaseProblem):
-    """
-         approach using IE
-    """
+    
+    surveyPair = Survey.LinearSurvey
+
+    def __init__(self, mesh, G, **kwargs):
+        Problem.BaseProblem.__init__(self, mesh, **kwargs)
+        self.G = G
+
+    def fields(self, m):
+        return self.G.dot(m)
+
+    def Jvec(self, m, v, u=None):
+        return self.G.dot(v)
+
+    def Jtvec(self, m, v, u=None):
+        return self.G.T.dot(v)
+
 
 
 class MagneticsDiffSecondary(Problem.BaseProblem):
@@ -809,7 +822,7 @@ def dipazm_2_xyz(dip,azm_N):
 
     return M
 
-def get_dist_wgt(mesh,rxLoc,R,R0):
+def get_dist_wgt(mesh,rxLoc,actv,R,R0):
     """
     get_dist_wgt(xn,yn,zn,rxLoc,R,R0)
 
@@ -819,6 +832,7 @@ def get_dist_wgt(mesh,rxLoc,R,R0):
     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)
 
@@ -830,17 +844,34 @@ def get_dist_wgt(mesh,rxLoc,R,R0):
     @author: dominiquef
     """
 
+    # Find non-zero cells
+    inds = np.nonzero(actv)[0]
+    
+    # Create active cell projector
+    P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))),
+                      shape=(mesh.nC, len(inds)))
+            
+    # 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)
-
-    V = np.reshape(mesh.vol,hY.shape)
-    wr = np.zeros(hY.shape)
+    
+    # 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(np.sum(actv))
 
     ndata = rxLoc.shape[0]
-    count = -1;
+    count = -1
     print "Begin calculation of distance weighting for R= " + str(R)
 
     for dd in range(ndata):
@@ -926,7 +957,7 @@ def getActiveTopo(mesh,topo,flag):
 
     Created on Dec, 27th 2015
 
-    @author: dominiquef
+    @founrdo
     """
     import scipy.interpolate as interpolation
 
@@ -959,4 +990,43 @@ def getActiveTopo(mesh,topo,flag):
 
     return actv
 
+def plot_obs_2D(rxLoc,d,wd,varstr):   
+    """ Function plot_obs(rxLoc,d,wd)
+    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
+    
+    # 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')
+    
+    # Plot result
+    plt.figure()
+    plt.subplot()
+    plt.imshow(d_grid, extent=[x.min(), x.max(), y.min(), y.max()],origin = 'lower')
+    plt.colorbar(fraction=0.02)
+    plt.contour(X,Y, d_grid,10)
+    plt.scatter(rxLoc[:,0],rxLoc[:,1], c=d, s=20)
+    plt.title(varstr)
+    plt.gca().set_aspect('equal', adjustable='box')
+