Implement sparse norm and test on nutcracker

simpegPF/Dev/Intgrl_MAG_Inv_Driver.py

@@ -94,7 +94,7 @@ PF.Magnetics.writeUBCobs(home_dir + dsep + 'Pred.dat',B,M,rxLoc,pred,wd)
 #%%
 plt.figure()
 ax = plt.subplot()
-mesh.plotSlice(wr_out, 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')
@@ -115,12 +115,6 @@ survey.dobs=d
 diagA = np.sum(F**2.,axis=0) + beta_in*np.ones(nC)
 PC     = sp.spdiags(diagA**-1., 0, nC, nC);
 
-# Create mesh with unit cells to remove dimensions from regularization
-hx = np.ones(mesh.nCx)
-hy = np.ones(mesh.nCy)
-hz = np.ones(mesh.nCz)
-meshreg = Mesh.TensorMesh([hx,hy,hz], mesh.x0)
-
 reg = Regularization.Simple(mesh, mapping=wrMap)
 reg.mref = mref
 #reg.alpha_s = 1.
@@ -132,7 +126,7 @@ opt.approxHinv = PC
 
 # opt = Optimization.InexactGaussNewton(maxIter=6)
 invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = beta_in)
-beta = Directives.BetaSchedule(coolingFactor=8, coolingRate=2)
+beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
 #betaest = Directives.BetaEstimate_ByEig()
 target = Directives.TargetMisfit()
 
@@ -143,6 +137,7 @@ m0 = mstart
 # Run inversion
 mrec = inv.run(m0)
 
+
 m_out = np.ones(mesh.nC)
 m_out[actv==1] = mrec
 
@@ -151,22 +146,72 @@ 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')
+#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-7
 plt.figure()
 ax = plt.subplot(211)
-mesh.plotSlice(m_out, ax = ax, normal = 'Z', ind=-5)
+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')
 plt.title('Inverted model')
 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=midx-12)
+mesh.plotSlice(m_out, ax = ax, normal = 'Y', ind=yslice, clim = (-mrec.min(), mrec.max()))
 plt.title('Inverted model')
 plt.xlabel('x');plt.ylabel('z')
 plt.gca().set_aspect('equal', adjustable='box')
+
+#%% Run one more round for sparsity
+
+
+reg = Regularization.SparseRegularization(mesh, mapping=wrMap)
+#reg.m = mrec
+reg.mref = mref
+reg.eps = 1e-3
+#reg.alpha_s = 1.
+
+dmis = DataMisfit.l2_DataMisfit(survey)
+dmis.Wd = wd
+opt = Optimization.ProjectedGNCG(maxIter=20,tolX = 1e-2, tolF = 1e-2,lower=0.,upper=1.)
+opt.approxHinv = PC
+
+# opt = Optimization.InexactGaussNewton(maxIter=6)
+invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = 1e-2)
+beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=10)
+#betaest = Directives.BetaEstimate_ByEig()
+target = Directives.TargetMisfit()
+
+inv = Inversion.BaseInversion(invProb, directiveList=[beta])
+
+m0 = mrec
+
+# Run inversion
+mrec = inv.run(m0)
+
+m_out[actv==1] = mrec
+
+#%% Plot out a section of the model
+
+yslice = midx-7
+plt.figure()
+ax = plt.subplot(211)
+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')
+plt.title('Inverted model')
+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('Inverted model')
+plt.xlabel('x');plt.ylabel('z')
+plt.gca().set_aspect('equal', adjustable='box')