Adapt MAG code for refactor IRLS iterations

simpegPF/Dev/MAG/Intgrl_MAG_Inv_Driver.py

@@ -11,12 +11,9 @@ survey = driver.survey
 rxLoc = survey.srcField.rxList[0].locs
 d = survey.dobs
 wd = survey.std
-
+    
 ndata = survey.srcField.rxList[0].locs.shape[0]
 
-eps_p = 1e-4
-eps_q = 1e-4
-
 actv = driver.activeCells
 nC = len(actv)
 
@@ -41,11 +38,14 @@ prob = PF.Magnetics.Problem3D_Integral(mesh, mapping=idenMap, actInd=actv)
 prob.solverOpts['accuracyTol'] = 1e-4
 survey.pair(prob)
 
+dmis = DataMisfit.l2_DataMisfit(survey)
+dmis.Wd = 1./wd
+
 # Write out the predicted
 pred = prob.fields(driver.m0)
 PF.Magnetics.writeUBCobs('Pred.dat', survey, pred)
 
-wr = np.sum(prob.G**2.,axis=0)**0.5 / mesh.vol[actv]
+wr = np.sum(prob.G**2.,axis=0)**0.5
 wr = ( wr/np.max(wr) )
 wr_out = actvMap * wr
 
@@ -55,117 +55,82 @@ mesh.plotSlice(wr_out, ax=ax, normal='Y', ind=midx ,clim=(-1e-3, wr.max()))
 plt.title('Distance weighting')
 plt.xlabel('x');plt.ylabel('z')
 plt.gca().set_aspect('equal', adjustable='box')
-
-reg = Regularization.Simple(mesh, indActive=actv, mapping=idenMap)
+reg = Regularization.Sparse(mesh, indActive=actv, mapping=idenMap)
 reg.mref = driver.mref
-reg.wght = wr
-
-dmis = DataMisfit.l2_DataMisfit(survey)
-dmis.Wd = 1/wd
-opt = Optimization.ProjectedGNCG(maxIter=10,lower=0.,upper=1., maxIterCG= 20, tolCG=1e-3)
+reg.cell_weights = wr
+    
+#reg.mref = np.zeros(mesh.nC)
+eps_p = 2e-4
+eps_q = 5e-5
+norms   = [0., 2., 2., 2.]
 
+opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3)
 invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
-
-beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
-beta_init = Directives.BetaEstimate_ByEig()
-target = Directives.TargetMisfit()
+#beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1)
+#update_beta = Directives.Scale_Beta(tol = 0.05, coolingRate=5)
+betaest = Directives.BetaEstimate_ByEig()
+IRLS = Directives.Update_IRLS( norms=norms,  eps_p=eps_p, eps_q=eps_q, f_min_change = 1e-2)
 update_Jacobi = Directives.Update_lin_PreCond(onlyOnStart=True)
-
-inv = Inversion.BaseInversion(invProb, directiveList=[beta,target,beta_init,update_Jacobi])
+inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi])
 
 # Run inversion
 mrec = inv.run(driver.m0)
 
-m_out = actvMap*mrec
-
-# Write result
-Mesh.TensorMesh.writeModelUBC(mesh,'SimPEG_inv_l2l2.sus',m_out)
-
-# Plot predicted
-pred = prob.fields(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. ) )
-
-#%% Plot out a section of the model
-
-yslice = midx
-m_out[m_out==-100] = np.nan
-
-plt.figure()
-ax = plt.subplot(221)
-mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-5, clim = (mrec.min(), mrec.max()))
-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 = (mrec.min(), mrec.max()))
-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 = (mrec.min(), mrec.max()))
-plt.title('Cross Section')
-plt.xlabel('x');plt.ylabel('z')
-plt.gca().set_aspect('equal', adjustable='box')
-
-#%% Run one more round for sparsity
-phim = invProb.phi_m_last
-phid =  invProb.phi_d
-
-# Set parameters for sparsity
-reg = Regularization.Sparse(mesh, indActive = actv, mapping = idenMap)
-reg.recModel = mrec
-reg.mref = driver.mref
-reg.wght = wr
-reg.eps_p = eps_p
-reg.eps_q = eps_q
-reg.norms   = driver.lpnorms
-
-
-dmis = DataMisfit.l2_DataMisfit(survey)
-dmis.Wd = wd
-opt = Optimization.ProjectedGNCG(maxIter=20 ,lower=0.,upper=1., maxIterCG= 10, tolCG = 1e-4)
-
-invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta)
-
-# Create inversion directives
-beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1)
-update_beta = Directives.Scale_Beta(tol = 0.05)
-target = Directives.TargetMisfit()
-IRLS =Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid )
-update_Jacobi = Directives.Update_lin_PreCond(onlyOnStart=False)
-save_log = Directives.SaveOutputEveryIteration()
-save_log.fileName = 'LogName_blabla'
-
-inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS,update_beta,update_Jacobi,save_log])
-
-m0 = mrec
-
-# Run inversion
-mrec = inv.run(m0)
-
-m_out = actvMap*mrec
-
-# Write final model out.
-Mesh.TensorMesh.writeModelUBC(mesh,'SimPEG_inv_l0l2.sus',m_out)
-
 pred = prob.fields(mrec)
 
 #%% Plot obs data
 PF.Magnetics.plot_obs_2D(rxLoc,pred,'Predicted 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
 
+
+m_out = actvMap*reg.l2model
+m_out[m_out==-100] = np.nan
+
+# Write result
+Mesh.TensorMesh.writeModelUBC(mesh,'SimPEG_inv_l2l2.sus',m_out)
+
+plt.figure()
+ax = plt.subplot(221)
+mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-5, clim = (mrec.min(), mrec.max()))
+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 = (mrec.min(), mrec.max()))
+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 = (mrec.min(), mrec.max()))
+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(reg.l2model,100)
+plt.yscale('log', nonposy='clip')
+plt.title('Histogram of model values - Smooth')
+ax = plt.subplot(122)
+plt.hist(reg.regmesh.cellDiffxStencil*reg.l2model,100)
+plt.yscale('log', nonposy='clip')
+plt.title('Histogram of model gradient values - Smooth')
+
+#%% Plot out a section of the model
+
+yslice = midx
+
+m_out = actvMap*mrec
 m_out[m_out==-100] = np.nan
 
 plt.figure()
@@ -189,5 +154,14 @@ 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 - Compact')
+ax = plt.subplot(122)
+plt.hist(reg.regmesh.cellDiffxStencil*mrec,100)
+plt.yscale('log', nonposy='clip')
+plt.title('Histogram of model gradient values - Smooth')
 
 plt.show()

simpegPF/Dev/MAG/LogName_blabla.txt

@@ -0,0 +1,21 @@
+  #     beta     phi_d     phi_m       f
+   1 1.4081e+07 1.3177e+01 9.9671e-06 2.8314e+02
+   2 1.4081e+07 1.2333e+01 9.2396e-06 1.5353e+02
+   3 1.4081e+07 1.0202e+01 9.1982e-06 1.4244e+02
+   4 1.4081e+07 1.0662e+01 9.0940e-06 1.3972e+02
+   5 4.4406e+08 9.7668e+00 9.1243e-06 1.3872e+02
+   6 4.4406e+08 3.2344e+02 3.1445e-06 4.0480e+03
+   7 4.4406e+08 5.1839e+02 2.2545e-06 1.7198e+03
+   8 4.4406e+08 5.1121e+02 2.1668e-06 1.5195e+03
+   9 4.4406e+08 5.4030e+02 2.0596e-06 1.4734e+03
+  10 2.4887e+08 5.4956e+02 2.0182e-06 1.4549e+03
+  11 2.4887e+08 3.1821e+02 2.2401e-06 1.0621e+03
+  12 2.4887e+08 2.8746e+02 2.0863e-06 8.7570e+02
+  13 2.4887e+08 2.3223e+02 2.1657e-06 8.0668e+02
+  14 2.4887e+08 2.2477e+02 2.1052e-06 7.7121e+02
+  15 3.9177e+08 1.9565e+02 2.1619e-06 7.4869e+02
+  16 3.9177e+08 2.9301e+02 1.4206e-06 1.0204e+03
+  17 3.9177e+08 2.7202e+02 1.3931e-06 8.4957e+02
+  18 3.9177e+08 2.7366e+02 1.3448e-06 8.1779e+02
+  19 3.9177e+08 2.4196e+02 1.3947e-06 8.0051e+02
+  20 4.9390e+08 2.4431e+02 1.3653e-06 7.8835e+02