wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-08-06 13:30:16 +08:00
Code Issues Packages Projects Releases Wiki Activity

Implement sparse norm on mag problem.

Browse Source 
Add Nutcracker example to folder.
This commit is contained in:
D Fournier
2016-01-31 15:33:57 -08:00
parent 90838aa74b
commit bb9624f4ff
8 changed files with 29104 additions and 8104 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegPF/Dev/Intgrl_MAG_Inv_Driver.py
+58 -34
View File
@@ -1,19 +1,18 @@
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\\Nut_Cracker\\Induced_MAG3C'
inpfile = 'PYMAG3D_inv.inp'
dsep = '\\'
os.chdir(home_dir)
#%%
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\\Block_Gaussian_topo'
home_dir = '.\\'
inpfile = 'PYMAG3D_inv.inp'
dsep = '\\'
os.chdir(home_dir)
## New scripts to be added to basecode
#from fwr_MAG_data import fwr_MAG_data
#from read_MAGfwr_inp import read_MAGfwr_inp
@@ -23,8 +22,9 @@ import pylab as plt
[mshfile, obsfile, topofile, mstart, mref, magfile, wgtfile, chi, alphas, bounds, lpnorms] = PF.BaseMag.read_MAGinv_inp(home_dir + dsep + inpfile)
# Load mesh file
mesh = Utils.meshutils.readUBCTensorMesh(mshfile)
mesh = Mesh.TensorMesh.readUBC(mshfile)
#mesh = Utils.meshutils.readUBCTensorMesh(mshfile)
# Load in observation file
[B,M,dobs] = PF.BaseMag.readUBCmagObs(obsfile)
@@ -85,8 +85,8 @@ wrMap = PF.BaseMag.WeightMap(mesh, wr)
wr_out = np.zeros(mesh.nC)
wr_out[actv==1] = wr
Utils.writeUBCTensorModel(home_dir+dsep+'wr.dat',mesh,wr_out)
Mesh.TensorMesh.writeModelUBC(mesh,home_dir+dsep+'wr.dat',wr_out)
#Utils.meshutils.writeUBCTensorModel(home_dir+dsep+'wr.dat',mesh,wr_out)
# Write out the predicted
pred = F.dot(mstart)
PF.Magnetics.writeUBCobs(home_dir + dsep + 'Pred.dat',B,M,rxLoc,pred,wd)
@@ -111,14 +111,16 @@ survey.pair(prob)
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.Simple(mesh, mapping=wrMap)
reg.mref = mref
#reg.alpha_s = 1.
# Create pre-conditioner
diagA = np.sum(F**2.,axis=0) + beta_in*(reg.W.T*reg.W).diagonal()*(wr**2.0)
PC = Utils.sdiag(diagA**-1.)
dmis = DataMisfit.l2_DataMisfit(survey)
dmis.Wd = wd
opt = Optimization.ProjectedGNCG(maxIter=10,lower=0.,upper=1.)
@@ -142,7 +144,8 @@ m_out = np.ones(mesh.nC)
m_out[actv==1] = mrec
# Write result
Utils.meshutils.writeUBCTensorModel('SimPEG_inv.sus',mesh,m_out)
Mesh.TensorMesh.writeModelUBC(mesh,'SimPEG_inv.sus',m_out)
#Utils.meshutils.writeUBCTensorModel(home_dir+dsep+'wr.dat',mesh,wr_out)
# Plot predicted
pred = F.dot(mrec)
@@ -155,41 +158,55 @@ print "Final misfit:" + str(np.sum( ((d-pred)/wd)**2. ) )
yslice = midx-7
plt.figure()
ax = plt.subplot(211)
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')
plt.title('Inverted model')
plt.plot(np.array([mesh.vectorCCx[0],mesh.vectorCCx[-1]]), np.array([mesh.vectorCCy[yslice],mesh.vectorCCy[yslice]]),c='w',linestyle = '--')
plt.title('Z Section')
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=-1, 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('Top')
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.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
reg = Regularization.SparseRegularization(mesh, mapping=wrMap)
#reg.m = mrec
reg.m = mrec
reg.mref = mref
reg.eps = 1e-3
reg.eps = 1e-1
diagA = np.sum(F**2.,axis=0) + beta_in*(reg.W.T*reg.W).diagonal()*(wr**2.0)
PC = Utils.sdiag(diagA**-1.)
#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 = Optimization.ProjectedGNCG(maxIter=20 ,lower=0.,upper=1.)
opt.approxHinv = PC
#opt.phim_last = reg.eval(mrec)
# opt = Optimization.InexactGaussNewton(maxIter=6)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = 1e-2)
beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=10)
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(factor = 2, eps_min=1e-3, phi_m_last = phim )
inv = Inversion.BaseInversion(invProb, directiveList=[beta])
inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS])
m0 = mrec
@@ -202,16 +219,23 @@ m_out[actv==1] = mrec
yslice = midx-7
plt.figure()
ax = plt.subplot(211)
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')
plt.title('Inverted model')
plt.plot(np.array([mesh.vectorCCx[0],mesh.vectorCCx[-1]]), np.array([mesh.vectorCCy[yslice],mesh.vectorCCy[yslice]]),c='w',linestyle = '--')
plt.title('Z Section')
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=-1, 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('Top')
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.title('Cross Section')
plt.xlabel('x');plt.ylabel('z')
plt.gca().set_aspect('equal', adjustable='box')