diff --git a/simpegPF/Dev/Intgrl_MAG_Fwr_InterpData_Test.py b/simpegPF/Dev/Intgrl_MAG_Fwr_InterpData_Test.py index aa5e9680..7b2ff2dd 100644 --- a/simpegPF/Dev/Intgrl_MAG_Fwr_InterpData_Test.py +++ b/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 diff --git a/simpegPF/Dev/Intgrl_MAG_Inv_Driver.py b/simpegPF/Dev/Intgrl_MAG_Inv_Driver.py index b49f38dc..c7d12546 100644 --- a/simpegPF/Dev/Intgrl_MAG_Inv_Driver.py +++ b/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. ) ) \ No newline at end of file diff --git a/simpegPF/Magnetics.py b/simpegPF/Magnetics.py index f2408d95..5227be73 100644 --- a/simpegPF/Magnetics.py +++ b/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') + \ No newline at end of file diff --git a/simpegPF/notebooks/SimPEG Tutorial - MAG Linear Problem.ipynb b/simpegPF/notebooks/SimPEG Tutorial - MAG Linear Problem.ipynb index 63e15e34..619729be 100644 --- a/simpegPF/notebooks/SimPEG Tutorial - MAG Linear Problem.ipynb +++ b/simpegPF/notebooks/SimPEG Tutorial - MAG Linear Problem.ipynb @@ -100,7 +100,7 @@ }, { "cell_type": "code", - "execution_count": 59, + "execution_count": 13, "metadata": { "collapsed": false, "scrolled": true @@ -119,9 +119,9 @@ "# Create a mesh\n", "dx = 5.\n", "\n", - "hxind = [(dx,5,-1.3), (dx, 20), (dx,5,1.3)]\n", - "hyind = [(dx,5,-1.3), (dx, 20), (dx,5,1.3)]\n", - "hzind = [(dx,5,-1.3),(5, 10)]\n", + "hxind = [(dx, 20)]\n", + "hyind = [(dx, 20)]\n", + "hzind = [(5, 10)]\n", "\n", "mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCC')\n", "\n", @@ -145,7 +145,7 @@ }, { "cell_type": "code", - "execution_count": 60, + "execution_count": 14, "metadata": { "collapsed": false }, @@ -153,10 +153,10 @@ { "data": { "text/plain": [ - "59.390075000000003" + "30.0" ] }, - "execution_count": 60, + "execution_count": 14, "metadata": {}, "output_type": "execute_result" } @@ -200,7 +200,7 @@ }, { "cell_type": "code", - "execution_count": 61, + "execution_count": 15, "metadata": { "collapsed": false }, @@ -237,7 +237,7 @@ }, { "cell_type": "code", - "execution_count": 91, + "execution_count": 16, "metadata": { "collapsed": false }, @@ -1000,7 +1000,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -1012,7 +1012,7 @@ ], "source": [ "# Generate a distance weighting\n", - "wr = PF.Magnetics.get_dist_wgt(mesh,rxLoc,3.,np.min(mesh.hx)/4)\n", + "wr = PF.Magnetics.get_dist_wgt(mesh,rxLoc,actv,3.,np.min(mesh.hx)/4)\n", "wrMap = PF.BaseMag.WeightMap(mesh, wr)\n", "\n", "plt.figure()\n", @@ -1025,7 +1025,7 @@ }, { "cell_type": "code", - "execution_count": 92, + "execution_count": 17, "metadata": { "collapsed": false }, @@ -1769,7 +1769,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -1809,7 +1809,7 @@ }, { "cell_type": "code", - "execution_count": 93, + "execution_count": 18, "metadata": { "collapsed": false }, @@ -1850,7 +1850,7 @@ }, { "cell_type": "code", - "execution_count": 97, + "execution_count": 19, "metadata": { "collapsed": false }, @@ -1869,12 +1869,13 @@ "dmis = DataMisfit.l2_DataMisfit(survey)\n", "dmis.Wd = np.ones(len(data))*1.\n", "opt = Optimization.ProjectedGNCG(maxIter=6,lower=0.,upper=1.)\n", + "opt.tolX = 1e-10\n", "# opt = Optimization.InexactGaussNewton(maxIter=6)\n", - "invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = 1e+4)\n", + "invProb = InvProblem.BaseInvProblem(dmis, reg, opt)\n", "beta = Directives.BetaSchedule()\n", - "#betaest = Directives.BetaEstimate_ByEig()\n", + "betaest = Directives.BetaEstimate_ByEig(beta0_ratio = 10**-10)\n", "target = Directives.TargetMisfit()\n", - "inv = Inversion.BaseInversion(invProb, directiveList=[beta, target])\n", + "inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])\n", "reg.alpha_s =0.0025\n", "m0 = np.ones_like(survey.mtrue)*1e-4" ] @@ -1888,7 +1889,7 @@ }, { "cell_type": "code", - "execution_count": 98, + "execution_count": 20, "metadata": { "collapsed": false }, @@ -1903,12 +1904,12 @@ "=============================== Projected GNCG ===============================\n", " # beta phi_d phi_m f |proj(x-g)-x| LS Comment \n", "-----------------------------------------------------------------------------\n", - " 0 1.00e+04 5.58e+04 2.22e-07 5.58e+04 8.41e+01 0 \n", + " 0 1.13e-02 5.62e+04 2.15e-07 5.62e+04 5.06e+01 0 \n", "------------------------- STOP! -------------------------\n", - "1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 5.5847e+03\n", - "1 : |xc-x_last| = 4.1645e-02 <= tolX*(1+|x0|) = 1.0116e-01\n", - "0 : |proj(x-g)-x| = 8.4094e+01 <= tolG = 1.0000e-01\n", - "0 : |proj(x-g)-x| = 8.4094e+01 <= 1e3*eps = 1.0000e-02\n", + "1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 5.6238e+03\n", + "0 : |xc-x_last| = 3.9811e-02 <= tolX*(1+|x0|) = 1.0063e-10\n", + "0 : |proj(x-g)-x| = 5.0612e+01 <= tolG = 1.0000e-01\n", + "0 : |proj(x-g)-x| = 5.0612e+01 <= 1e3*eps = 1.0000e-02\n", "0 : maxIter = 6 <= iter = 1\n", "------------------------- DONE! -------------------------\n" ] @@ -1920,7 +1921,29 @@ }, { "cell_type": "code", - "execution_count": 99, + "execution_count": 21, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "text/plain": [ + "False" + ] + }, + "execution_count": 21, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.nan <= 0.1" + ] + }, + { + "cell_type": "code", + "execution_count": 22, "metadata": { "collapsed": false }, @@ -2664,7 +2687,7 @@ { "data": { "text/html": [ - "" + "" ], "text/plain": [ "" @@ -2676,10 +2699,10 @@ { "data": { "text/plain": [ - "" + "" ] }, - "execution_count": 99, + "execution_count": 22, "metadata": {}, "output_type": "execute_result" } @@ -2722,17 +2745,6 @@ "Hopefully you now have an idea of how to create a Problem class in SimPEG, and how this can be used with the other tools available." ] }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "collapsed": false - }, - "outputs": [], - "source": [ - "reg.W" - ] - }, { "cell_type": "code", "execution_count": null,