Update grav code and example

simpegPF/Dev/GRAV/Intgrl_GRAV_Inv_Driver.py

@@ -2,42 +2,41 @@
 from SimPEG import *
 import simpegPF as PF
 import pylab as plt
-
 import os
-
-#home_dir = 'C:\Users\dominiquef.MIRAGEOSCIENCE\Documents\GIT\SimPEG\simpegpf\simpegPF\Dev'
-#home_dir = 'C:\\Users\\dominiquef.MIRAGEOSCIENCE\\ownCloud\\Research\\Modelling\\Synthetic\\Nut_Cracker\\Induced_MAG3C'
-#home_dir = 'C:\\Users\\dominiquef.MIRAGEOSCIENCE\\ownCloud\\Research\\Modelling\\Synthetic\\Block_Gaussian_topo\\GRAV'
-#home_dir ='C:\\Users\\dominiquef.MIRAGEOSCIENCE\\ownCloud\\Research\\CraigModel'
 home_dir = '.\\'
-
 inpfile = 'PYGRAV3D_inv.inp'
-
 dsep = '\\'
 os.chdir(home_dir)
-## New scripts to be added to basecode
+plt.close('all')
-#from fwr_MAG_data import fwr_MAG_data
-#from read_MAGfwr_inp import read_MAGfwr_inp
 
+#%% User input
+# Initial beta
+beta_in = 1e+2
+
+# Treshold values for compact norm
+eps_p = 1e-2 # Small model values
+eps_q = 1e-2 # Small model gradient
+
+# Plotting parameter
+vmin = -0.1
+vmax = 0.2
 #%%
 # Read input file
 [mshfile, obsfile, topofile, mstart, mref, wgtfile, chi, alphas, bounds, lpnorms] = PF.Gravity.read_GRAVinv_inp(home_dir + dsep + inpfile)
 
 # Load mesh file
 mesh = Mesh.TensorMesh.readUBC(mshfile)
-#mesh = Utils.meshutils.readUBCTensorMesh(mshfile) 
 
 # Load in observation file
 survey = PF.Gravity.readUBCgravObs(obsfile)
 
+# Get obs location and data
 rxLoc = survey.srcField.rxList[0].locs
 d = survey.dobs
 wd = survey.std
 
 ndata = survey.srcField.rxList[0].locs.shape[0]
 
-beta_in = 1e+1
-
 # Load in topofile or create flat surface
 if topofile == 'null':
     
@@ -53,7 +52,7 @@ else:
 nC = len(actv)
 
 # Create active map to go from reduce set to full
-actvMap = Maps.ActiveCells(mesh, actv, -100)
+actvMap = Maps.InjectActiveCells(mesh, actv, -100)
 
 # Creat reduced identity map
 idenMap = Maps.IdentityMap(nP = nC)
@@ -73,7 +72,7 @@ else:
     mref = Mesh.TensorMesh.readModelUBC(mesh,mref)
     mref = mref[actv]
     
-# Get index of the center
+# Get index of the center for plotting
 midx = int(mesh.nCx/2)
 midy = int(mesh.nCy/2)
 
@@ -85,46 +84,47 @@ PF.Gravity.plot_obs_2D(survey,'Observed Data')
 prob = PF.Gravity.GravityIntegral(mesh, mapping = idenMap, actInd = actv)
 prob.solverOpts['accuracyTol'] = 1e-4
 
-#survey = Survey.LinearSurvey()
 survey.pair(prob)
-#survey.makeSyntheticData(data, std=0.01)
+
-#survey.dobs=d
+# Write out the predicted file and generate the forward operator
-#survey.mtrue = model
-# Write out the predicted
 pred = prob.fields(mstart)
 
 PF.Gravity.writeUBCobs(home_dir + dsep + 'Pred0.dat',survey,pred)
 
-wr = np.sum(prob.G**2.,axis=0)**0.5
+# Make depth weighting
+wr = np.sum(prob.G**2.,axis=0)**0.5 / mesh.vol[actv]
 wr = ( wr/np.max(wr) )
 wr_out = actvMap * wr
 
-plt.figure()
+#%% Plot depth weighting
-ax = plt.subplot()
+#plt.figure()
-mesh.plotSlice(actvMap*mstart, ax = ax, normal = 'Y', ind=midx+1 ,clim = (-1e-1, mstart.max()))
+#ax = plt.subplot()
-plt.title('Distance weighting')
+#mesh.plotSlice(actvMap*mstart, ax = ax, normal = 'Y', ind=midx+1 ,clim = (-1e-1, mstart.max()))
-plt.xlabel('x');plt.ylabel('z')
+#plt.title('Distance weighting')
-plt.gca().set_aspect('equal', adjustable='box')
+#plt.xlabel('x');plt.ylabel('z')
+#plt.gca().set_aspect('equal', adjustable='box')
 
+#%% Create inversion objects
+
+# First start with an l2 regularization
 reg = Regularization.Simple(mesh, indActive = actv, mapping = idenMap)
 reg.mref = mref
 reg.wght = wr
-#reg.alpha_s = 1.
+
 
 # Create pre-conditioner 
 diagA = np.sum(prob.G**2.,axis=0) + beta_in*(reg.W.T*reg.W).diagonal()*wr
 PC     = Utils.sdiag(diagA**-1.)
 
-
+# Data misfit function
 dmis = DataMisfit.l2_DataMisfit(survey)
-dmis.Wd = 1/wd
+dmis.Wd = 1./wd
-opt = Optimization.ProjectedGNCG(maxIter=10,lower=0.,upper=1., maxIterCG= 20, tolCG = 1e-3)
+opt = Optimization.ProjectedGNCG(maxIter=20,lower=bounds[0],upper=bounds[1], maxIterCG= 20, tolCG = 1e-3)
 opt.approxHinv = PC
 
-# opt = Optimization.InexactGaussNewton(maxIter=6)
+
 invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = beta_in)
 beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
-#betaest = Directives.BetaEstimate_ByEig()
 target = Directives.TargetMisfit()
 
 inv = Inversion.BaseInversion(invProb, directiveList=[beta,target])
@@ -147,32 +147,40 @@ pred = prob.fields(mrec)
 
 print "Final misfit:" + str(np.sum( ((d-pred)/wd)**2. ) ) 
 
-#%% Plot out a section of the model
+#%% Plot out sections of the smooth model
 
 yslice = midx+1
 
 plt.figure()
 ax = plt.subplot(221)
-mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-5, clim = (-1e-3, mrec.max()))
+mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-5, clim = (vmin,vmax))
 plt.plot(np.array([mesh.vectorCCx[0],mesh.vectorCCx[-1]]), np.array([mesh.vectorCCy[yslice],mesh.vectorCCy[yslice]]),c='w',linestyle = '--')
 plt.title('Z: ' + str(mesh.vectorCCz[-5]) + ' m')
 plt.xlabel('x');plt.ylabel('z')
 plt.gca().set_aspect('equal', adjustable='box')
 
 ax = plt.subplot(222)
-mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-8, clim = (-1e-3, mrec.max()))
+mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-8, clim = (vmin,vmax))
 plt.plot(np.array([mesh.vectorCCx[0],mesh.vectorCCx[-1]]), np.array([mesh.vectorCCy[yslice],mesh.vectorCCy[yslice]]),c='w',linestyle = '--')
 plt.title('Z: ' + str(mesh.vectorCCz[-8]) + ' m')
 plt.xlabel('x');plt.ylabel('z')
 plt.gca().set_aspect('equal', adjustable='box')
 
-
 ax = plt.subplot(212)
-mesh.plotSlice(m_out, ax = ax, normal = 'Y', ind=yslice, clim = (-1e-3, mrec.max()))
+mesh.plotSlice(m_out, ax = ax, normal = 'Y', ind=yslice, clim = (vmin,vmax))
 plt.title('Cross Section')
 plt.xlabel('x');plt.ylabel('z')
 plt.gca().set_aspect('equal', adjustable='box')
 
+plt.figure()
+ax = plt.subplot(121)
+plt.hist(mrec,100)
+plt.yscale('log', nonposy='clip')
+plt.title('Histogram of model values - Smooth')
+ax = plt.subplot(122)
+plt.hist(reg.Wsmooth*mrec,100)
+plt.yscale('log', nonposy='clip')
+plt.title('Histogram of model gradient values - Smooth')
 #%% Run one more round for sparsity
 phim = invProb.phi_m_last
 phid =  invProb.phi_d
@@ -181,11 +189,12 @@ reg = Regularization.Sparse(mesh, indActive = actv, mapping = idenMap)
 reg.recModel = mrec
 reg.mref = mref
 reg.wght = wr
-reg.eps = 1e-2
+reg.eps_p = eps_p
+reg.eps_q = eps_q
 reg.p   = lpnorms[0]
 reg.qx  = lpnorms[1]
-reg.qz  = lpnorms[2]
+reg.qy  = lpnorms[2]
-reg.qy  = lpnorms[3]
+reg.qz  = lpnorms[3]
 
 diagA = np.sum(prob.G**2.,axis=0) + beta_in*(reg.W.T*reg.W).diagonal()*(wr)
 PC     = Utils.sdiag(diagA**-1.)
@@ -193,8 +202,8 @@ PC     = Utils.sdiag(diagA**-1.)
 #reg.alpha_s = 1.
 
 dmis = DataMisfit.l2_DataMisfit(survey)
-dmis.Wd = 1/wd
+dmis.Wd = 1./wd
-opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=0.,upper=1., maxIterCG= 10, tolCG = 1e-4)
+opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=bounds[0],upper=bounds[1], maxIterCG= 25, tolCG = 1e-4)
 opt.approxHinv = PC
 #opt.phim_last = reg.eval(mrec)
 
@@ -203,7 +212,7 @@ invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta)
 beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1)
 #betaest = Directives.BetaEstimate_ByEig()
 target = Directives.TargetMisfit()
-IRLS =Directives.update_IRLS( phi_m_last = phim, phi_d_last = phid )
+IRLS =Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid )
 
 inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS])
 
@@ -219,29 +228,39 @@ Mesh.TensorMesh.writeModelUBC(mesh,'SimPEG_inv_l0l2.sus',m_out)
 pred = prob.fields(mrec)
 
 #%% Plot obs data
-PF.Magnetics.plot_obs_2D(rxLoc,pred,wd,'Predicted Data')
+PF.Magnetics.plot_obs_2D(rxLoc,pred,'Predicted Data', vmin = np.min(d), vmax = np.max(d))
-PF.Magnetics.plot_obs_2D(rxLoc,d,wd,'Observed Data')
+PF.Magnetics.plot_obs_2D(rxLoc,d,'Observed Data')
 print "Final misfit:" + str(np.sum( ((d-pred)/wd)**2. ) ) 
 #%% Plot out a section of the model
 
 yslice = midx
 plt.figure()
 ax = plt.subplot(221)
-mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-5, clim = (-1e-2, mrec.max()))
+mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-5, clim = (vmin,vmax))
 plt.plot(np.array([mesh.vectorCCx[0],mesh.vectorCCx[-1]]), np.array([mesh.vectorCCy[yslice],mesh.vectorCCy[yslice]]),c='w',linestyle = '--')
 plt.title('Z: ' + str(mesh.vectorCCz[-5]) + ' m')
 plt.xlabel('x');plt.ylabel('z')
 plt.gca().set_aspect('equal', adjustable='box')
 
 ax = plt.subplot(222)
-mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-8, clim = (-1e-2, mrec.max()))
+mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-8, clim = (vmin,vmax))
 plt.plot(np.array([mesh.vectorCCx[0],mesh.vectorCCx[-1]]), np.array([mesh.vectorCCy[yslice],mesh.vectorCCy[yslice]]),c='w',linestyle = '--')
 plt.title('Z: ' + str(mesh.vectorCCz[-8]) + ' m')
 plt.xlabel('x');plt.ylabel('z')
 plt.gca().set_aspect('equal', adjustable='box')
 
 ax = plt.subplot(212)
-mesh.plotSlice(m_out, ax = ax, normal = 'Y', ind=yslice, clim = (-1e-2, mrec.max()))
+mesh.plotSlice(m_out, ax = ax, normal = 'Y', ind=yslice, clim = (vmin,vmax))
 plt.title('Cross Section')
 plt.xlabel('x');plt.ylabel('z')
 plt.gca().set_aspect('equal', adjustable='box')
+
+plt.figure()
+ax = plt.subplot(121)
+plt.hist(mrec,100)
+plt.yscale('log', nonposy='clip')
+plt.title('Histogram of model values - Sparse lp:'+str(lpnorms[0]))
+ax = plt.subplot(122)
+plt.hist(reg.Wsmooth*mrec,100)
+plt.yscale('log', nonposy='clip')
+plt.title('Histogram of model gradient values - Sparse lqx: ' + str(lpnorms[1]) + ' lqy:'+ str(lpnorms[2]) + ' lqz:' + str(lpnorms[3]))

simpegPF/Dev/GRAV/PYGRAV3D_inv.inp

@@ -1,10 +1,10 @@
 Mesh_10m.msh			! Mesh file
-gzfor3d.dat	! Obsfile
+Grav.dat			! Obsfile
 Gaussian.topo			! Topofile | null
 VALUE 1e-4			! Starting model
-VALUE 0			! Reference model
+VALUE 0.0			! Reference model
 DEFAULT				! Cell based weight file
 1				! target chi factor | DEFAULT=1
 1 1 1 1			! alpha s, x ,y ,z
-VALUE 0 1			! Lower and Upper Bounds for p-component
+VALUE -1 1		! Lower and Upper Bounds for p-component
-VALUE 0 2 2 2 1		! lp-norm for amplitude inversion FILE pqxqyqzr.dat  ! Norms VALUE p, qx, qy, qz, r | FILE m-by-5 matrix
+VALUE 0 1 1 1 1		! lp-norm for amplitude inversion FILE pqxqyqzr.dat  ! Norms VALUE p, qx, qy, qz, r | FILE m-by-5 matrix

simpegPF/Gravity.py

@@ -34,11 +34,11 @@ class GravityIntegral(Problem.BaseProblem):
 
         # return self.G.dot(self.mapping*(m))
             
-    def Jvec(self, m, v, u=None):
+    def Jvec(self, m, v, f=None):
         dmudm = self.mapping.deriv(m)
         return self.G.dot(dmudm*v)
 
-    def Jtvec(self, m, v, u=None):
+    def Jtvec(self, m, v, f=None):
         dmudm = self.mapping.deriv(m)
         return dmudm.T * (self.G.T.dot(v))