From 27900498c65e3d0374e30d4c8cad06130a8a1a25 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Mon, 11 Jan 2016 16:40:28 -0800 Subject: [PATCH] Add surface to mesh test Fix forward function call --- .../Dev/Intgrl_MAG_Fwr_DualSphereTopo_Test.py | 59 ++++--- .../Dev/Intgrl_MAG_Fwr_InterpData_Test.py | 158 +++++++----------- simpegPF/Dev/Intgrl_MAG_Fwr_SphereTest.py | 4 +- simpegPF/Dev/interpFFT.py | 85 ++++++++++ 4 files changed, 184 insertions(+), 122 deletions(-) create mode 100644 simpegPF/Dev/interpFFT.py diff --git a/simpegPF/Dev/Intgrl_MAG_Fwr_DualSphereTopo_Test.py b/simpegPF/Dev/Intgrl_MAG_Fwr_DualSphereTopo_Test.py index 62d52c2e..115b5e4a 100644 --- a/simpegPF/Dev/Intgrl_MAG_Fwr_DualSphereTopo_Test.py +++ b/simpegPF/Dev/Intgrl_MAG_Fwr_DualSphereTopo_Test.py @@ -9,6 +9,7 @@ from SimPEG import * import matplotlib.pyplot as plt import simpegPF as PF import scipy.interpolate as interpolation +import time #from fwr_MAG_data import fwr_MAG_data @@ -39,24 +40,33 @@ xr = np.linspace(-99., 99., 40) yr = np.linspace(-49., 49., 20) X, Y = np.meshgrid(xr, yr) -d_iter = 1 -lrl = np.zeros(d_iter) + + sclx = 100. -dx = 5 +dx = np.asarray([15., 10., 5., 2.5]) + +d_iter = len(dx) +l1_r = np.zeros(d_iter) +l2_r = np.zeros(d_iter) +linf_r = np.zeros(d_iter) +timer = np.zeros(d_iter) +mcell = np.zeros(d_iter) #%% Loop through decreasing meshes and measure the residual # Create mesh using simpeg and write out in GIF format for ii in range(d_iter): - nc = int(sclx/dx) + nc = int(sclx/dx[ii]) - hxind = [(dx, 2*nc)] - hyind = [(dx, nc)] - hzind = [(dx, nc)] + hxind = [(dx[ii], 2*nc)] + hyind = [(dx[ii], nc)] + hzind = [(dx[ii], nc)] mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCN') + mcell[ii] = mesh.nC + actv = PF.Magnetics.getActiveTopo(mesh,topo,'N') # Drape observations on topo + offset @@ -74,15 +84,13 @@ for ii in range(d_iter): xn = mesh.vectorNx yn = mesh.vectorNy zn = mesh.vectorNz - - mcell = mesh.nC - - print 'Mesh size: ' + str(mcell) + + print 'Mesh size: ' + str(mcell[ii]) #%% Create model chibkg = 0. chiblk = 0.01 - model = np.ones(mcell)*chibkg + model = np.ones(mcell[ii])*chibkg # Do a three sphere problem for more frequencies sph_ind = PF.MagAnalytics.spheremodel(mesh, 0., 0., -sclx/3, R) @@ -98,7 +106,13 @@ for ii in range(d_iter): Utils.writeUBCTensorModel('Model.sus',mesh,model) actv = np.ones(mesh.nC) #%% Forward mode ldata + + start_time = time.time() + d = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,'tmi') + + timer[ii] = (time.time() - start_time) + #fwr_tmi = d[0:ndata] #fwr_y = d[ndata:2*ndata] #fwr_z = d[2*ndata:] @@ -130,8 +144,9 @@ for ii in range(d_iter): #r_By = fwr_y - bya #r_Bz = fwr_z - bza - lrl[ii] = sum( r_tmi**2 ) **0.5 - + l2_r[ii] = np.sum( r_tmi**2 ) **0.5 + l1_r[ii] = np.sum( np.abs( r_tmi ) ) + linf_r[ii] = np.max( np.abs( r_tmi ) ) #%% Write predicted to file @@ -139,34 +154,34 @@ PF.Magnetics.writeUBCobs('Obsloc.loc',B,M,rxLoc,d,np.ones(len(d))) #%% Plot results print 'Residual between analytical sphere and integral forward' +print "dx \t nc \t l1 \t l2 \t linf \t Runtime" for ii in range(d_iter): - nc = 3**(ii+1) - print "||r||= " + str(lrl[ii]) + "\t dx= " + str(1./nc) + print str(dx[ii]) + "\t" + str(mcell[ii]) + "\t" + str(l1_r[ii]) + "\t" + str(l2_r[ii]) + "\t" + str(linf_r[ii]) + "\t" + str(timer[ii]) #%% Plot fields plt.figure(1) -#ax = plt.subplot(221) +ax = plt.subplot() plt.imshow(np.reshape(b_tmi,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower') -plt.colorbar(fraction=0.04) +plt.colorbar(fraction=0.02) plt.contour(X,Y, np.reshape(b_tmi,X.shape),10) plt.scatter(X,Y, c=np.reshape(b_tmi,X.shape), s=20) ax.set_title('Analytical') #%% Plot the forward solution from integral plt.figure(2) -#ax = plt.subplot(222) +ax = plt.subplot() plt.imshow(np.reshape(d,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max() ], origin = 'lower') -plt.colorbar(fraction=0.04) +plt.colorbar(fraction=0.02) plt.contour(X,Y, np.reshape(d,X.shape),10) plt.scatter(X,Y, c=np.reshape(d,X.shape), s=20) ax.set_title('Numerical') #%% Plot residual data plt.figure(3) -#ax = plt.subplot(212) +ax = plt.subplot() plt.imshow(np.reshape(r_tmi,X.shape), interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max()], origin = 'lower') -plt.colorbar(fraction=0.04) +plt.colorbar(fraction=0.02) plt.contour(X,Y, np.reshape(r_tmi,X.shape),10) plt.scatter(X,Y, c=np.reshape(r_tmi,X.shape), s=20) ax.set_title('Sphere Ana Bx') \ No newline at end of file diff --git a/simpegPF/Dev/Intgrl_MAG_Fwr_InterpData_Test.py b/simpegPF/Dev/Intgrl_MAG_Fwr_InterpData_Test.py index 3a9525ca..6c747f4a 100644 --- a/simpegPF/Dev/Intgrl_MAG_Fwr_InterpData_Test.py +++ b/simpegPF/Dev/Intgrl_MAG_Fwr_InterpData_Test.py @@ -1,19 +1,18 @@ -import os - -home_dir = 'C:\\LC\\Private\\dominiquef\\Projects\\4414_Minsim\\Modeling\\MAG' - -os.chdir(home_dir) - #%% from SimPEG import * import matplotlib.pyplot as plt import simpegPF as PF import scipy.interpolate as interpolation -import scipy.signal as sn import time - +from interpFFT import interpFFT #from fwr_MAG_data import fwr_MAG_data +import os + +home_dir = 'C:\\LC\\Private\\dominiquef\\Projects\\4414_Minsim\\Modeling\\MAG' + +os.chdir(home_dir) + plt.close('all') topofile = 'Gaussian.topo' @@ -78,7 +77,7 @@ mcell = mesh.nC print 'Mesh size: ' + str(mcell) #%% Create model -chibkg = 0.000 +chibkg = 0.0001 chiblk = 0.01 model = np.ones(mcell)*chibkg @@ -92,9 +91,12 @@ model[sph_ind] = 4.*chiblk sph_ind = PF.MagAnalytics.spheremodel(mesh, sclx/2., 0., -sclx/2.5, R/2.5) model[sph_ind] = 2.5*chiblk +# Zero out +model[actv==0] = -100 Utils.writeUBCTensorMesh('Mesh.msh',mesh) Utils.writeUBCTensorModel('Model.sus',mesh,model) -actv = np.ones(mesh.nC) +Utils.writeUBCTensorModel('nullcell.dat',mesh,actv) +#actv = np.ones(mesh.nC) #%% Forward mode ldata d = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,'tmi') #fwr_tmi = d[0:ndata] @@ -103,7 +105,7 @@ d = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,'tmi') #%% Compute data on a line -xx = np.linspace(-101., 99., 200) +xx = np.linspace(xr.min(), xr.max(), 200) yy = np.zeros(len(xx)) zz = F(xx,0.) + zoffset @@ -118,20 +120,21 @@ l2_r = np.zeros((d_iter,5)) linf_r = np.zeros((d_iter,5)) timer = np.zeros((d_iter,5)) +d2d = np.reshape(d, (len(yr),len(xr))) #%% Try different interpolation schemes for ii in range(d_iter): indx = ii+1 - d2d = np.reshape(d, (len(yr),len(xr))) - d2d = mkvc(d2d[::indx,::indx].T) + + dsub = d2d[::indx,::indx] xsub = xr[::indx] ysub = yr[::indx] # Nearest Neighbourg start_time = time.time() - F = interpolation.NearestNDInterpolator(np.c_[mkvc(Y[::indx,::indx].T),mkvc(X[::indx,::indx].T)],d2d) + F = interpolation.NearestNDInterpolator(np.c_[mkvc(Y[::indx,::indx].T),mkvc(X[::indx,::indx].T)],mkvc(dsub.T)) d_i2d_nnb = mkvc( F(0.,xx) ) l1_r[ii,0] = np.sum( np.abs(d_line - d_i2d_nnb) )**0.5 l2_r[ii,0] = np.sum( (d_line - d_i2d_nnb)**2. ) @@ -140,7 +143,7 @@ for ii in range(d_iter): # Linear interpolation start_time = time.time() - F = interpolation.interp2d(ysub,xsub,d2d) + F = interpolation.interp2d(ysub,xsub,mkvc(dsub.T)) d_i2d_lin = mkvc( F(0.,xx) ) l1_r[ii,1] = np.sum( np.abs(d_line - d_i2d_lin) )**0.5 l2_r[ii,1] = np.sum( (d_line - d_i2d_lin)**2. ) @@ -150,7 +153,7 @@ for ii in range(d_iter): # Cubic interpolation start_time = time.time() - F = interpolation.interp2d(ysub,xsub,d2d,kind='cubic') + F = interpolation.interp2d(ysub,xsub,mkvc(dsub.T),kind='cubic') d_i2d_cub = mkvc( F(0.,xx) ) l1_r[ii,2] = np.sum( np.abs(d_line - d_i2d_cub) )**0.5 l2_r[ii,2] = np.sum( (d_line - d_i2d_cub)**2. ) @@ -160,95 +163,48 @@ for ii in range(d_iter): # Quintic interpolation start_time = time.time() - F = interpolation.interp2d(ysub,xsub,d2d,kind='quintic') + F = interpolation.interp2d(ysub,xsub,mkvc(dsub.T),kind='quintic') d_i2d_qui = mkvc( F(0.,xx) ) l1_r[ii,3] = np.sum( np.abs(d_line - d_i2d_qui) )**0.5 l2_r[ii,3] = np.sum( (d_line - d_i2d_qui)**2. ) linf_r[ii,3] = np.max( np.abs(d_line - d_i2d_qui) ) timer[ii,3] = (time.time() - start_time) + # CloughTocher interpolation + start_time = time.time() + F = interpolation.CloughTocher2DInterpolator(np.c_[mkvc(Y[::indx,::indx].T),mkvc(X[::indx,::indx].T)],mkvc(dsub.T)) + d_i2d_CTI = mkvc( F(0.,xx) ) + l1_r[ii,4] = np.sum( np.abs(d_line - d_i2d_CTI) )**0.5 + l2_r[ii,4] = np.sum( (d_line - d_i2d_CTI)**2. ) + linf_r[ii,4] = np.max( np.abs(d_line - d_i2d_CTI) ) + timer[ii,4] = (time.time() - start_time) + - -#%% FFT interpolation - -d2d = np.reshape(d, (len(yr),len(xr))) - -# Add padding values by reflection (2**n) -lenx = np.floor( np.log2( 2*len(xr) ) ) -npadx = int(np.floor( ( 2**lenx - len(xr) ) /2. )) - -#Create hemming taper -if np.mod(npadx*2+len(xr),2) != 0: - oddx = 1 - -else: - oddx = 0 - -tap0 = sn.hamming(npadx*2) -tapl = sp.spdiags(tap0[0:npadx],0,npadx,npadx) -tapr = sp.spdiags(tap0[npadx:],0,npadx,npadx+oddx) - -# Mirror the 2d data over the half lenght and apply 0-taper -d2dpad = np.hstack([np.fliplr(d2d[:,0:npadx]) * tapl, d2d, np.fliplr(d2d[:,-npadx:]) * tapr]) - -# Repeat along the second dimension -leny = np.floor( np.log2( 2*len(yr) ) ) -npady = int(np.floor( ( 2**leny - len(yr) ) /2. )) - -#Create hemming taper -if np.mod(npady*2+len(yr),2) != 0: - oddy = 1 - -else: - oddy = 0 - -tap0 = sn.hamming(npady*2) -tapu = sp.spdiags(tap0[0:npady],0,npady,npady) -tapd = sp.spdiags(tap0[npady:],0,npady+oddy,npady) - -d2dpad = np.vstack([tapu*np.flipud(d2dpad[0:npady,:]), d2dpad, tapd*np.flipud(d2dpad[-npady:,:])]) - -# Compute FFT -FFTd2d = np.fft.fft2(d2dpad) - -# Compute IFFT at a given location -# Do an FFT shift -#FFTshift = np.fft.fftshift(FFTd2d) -FFTd2d = np.hstack([FFTd2d[:,0:31],np.zeros((32,64)),FFTd2d[:,31:]]) -FFTd2d = np.vstack([FFTd2d[0:15,:],np.zeros((32,128)),FFTd2d[15:,:]]) -# Pad with zeros -#temp = np.zeros((FFTd2d.shape[0]*2,FFTd2d.shape[1]*2)) - - -# Compute inverse FFT -IFFTd2d = np.fft.ifft2(FFTd2d)*FFTd2d.size/d2dpad.size - - - -# Extract core -#d2d_out = np.real(IFFTd2d[npady*2:-(npady*2+oddy+1),npadx*2:-(npadx*2+oddx+1)]) -d2d_out = np.real(IFFTd2d[:-1,:-1]) -d2d_out = d2d_out[npady*2:-(npady*2+oddy),npadx*2:-(npadx*2+oddx)] - -# Create new distance vector -XX = np.linspace(np.min(xr),np.max(xr),d2d_out.shape[1]) -YY = np.linspace(np.min(yr),np.max(yr),d2d_out.shape[0]) - -start_time = time.time() -F = interpolation.interp2d(XX,YY,d2d_out) -d_i2d_fft = mkvc( F(xx,0.) ) -l1_r[ii,4] = np.sum( np.abs(d_line - d_i2d_nnb) )**0.5 -l2_r[ii,4] = np.sum( (d_line - d_i2d_nnb)**2. ) -linf_r[ii,4] = np.max( np.abs(d_line - d_i2d_nnb) ) -timer[ii,4] = (time.time() - start_time) - -print("--- FFT completed in %s seconds ---" % (time.time() - start_time)) - -plt.figure() -ax = plt.subplot() -plt.imshow(d2d_out, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max() ], origin = 'lower') -plt.colorbar(fraction=0.04) +#============================================================================== +# #%% FFT interpolation +# d2d_out = interpFFT(xsub,ysub,dsub) +# +# # Create new distance vector +# XX = np.linspace(np.min(xsub),np.max(xsub),d2d_out.shape[1]) +# YY = np.linspace(np.min(ysub),np.max(ysub),d2d_out.shape[0]) +# +# start_time = time.time() +# F = interpolation.interp2d(XX,YY,d2d_out) +# d_i2d_fft = mkvc( F(xx,0.) ) +# l1_r[ii,4] = np.sum( np.abs(d_line - d_i2d_nnb) )**0.5 +# l2_r[ii,4] = np.sum( (d_line - d_i2d_nnb)**2. ) +# linf_r[ii,4] = np.max( np.abs(d_line - d_i2d_nnb) ) +# timer[ii,4] = (time.time() - start_time) +# +# print("--- FFT completed in %s seconds ---" % (time.time() - start_time)) +# +# +# plt.figure() +# ax = plt.subplot() +# plt.imshow(d2d_out, interpolation="bicubic", extent=[xr.min(), xr.max(), yr.min(), yr.max() ], origin = 'lower') +# plt.colorbar(fraction=0.04) +#============================================================================== #%% Write predicted to file @@ -287,8 +243,12 @@ plt.plot(xx,d_i2d_lin,c='b') plt.plot(xx,d_i2d_cub,c='g') plt.plot(xx,d_i2d_qui,c='m') plt.plot(xx,d_i2d_nnb,c='k') -plt.plot(xx,d_i2d_fft,c='c',marker='o') -plt.plot(xr[::indx],np.zeros(len(xr[::indx])),c='r',linewidth=0.,marker='o') +plt.plot(xx,d_i2d_CTI,c='c') + +# Plot interpolation from true value on line +F = interpolation.interp1d(xx,d_line) +dtrue = F(xr[::indx]) +plt.plot(xr[::indx],dtrue,c='r',linewidth=0.,marker='o') ax.set_title('Analytical') #%% Write result to file diff --git a/simpegPF/Dev/Intgrl_MAG_Fwr_SphereTest.py b/simpegPF/Dev/Intgrl_MAG_Fwr_SphereTest.py index ddea3760..fbf48ae7 100644 --- a/simpegPF/Dev/Intgrl_MAG_Fwr_SphereTest.py +++ b/simpegPF/Dev/Intgrl_MAG_Fwr_SphereTest.py @@ -61,8 +61,10 @@ for ii in range(d_iter): model = np.ones(mcell)*chibkg model[sph_ind] = chiblk + actv = np.ones(mcell) + #%% Forward mode ldata - d = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,'xyz') + d = PF.Magnetics.Intgrl_Fwr_Data(mesh,B,M,rxLoc,model,actv,'xyz') fwr_x = d[0:ndata] fwr_y = d[ndata:2*ndata] fwr_z = d[2*ndata:] diff --git a/simpegPF/Dev/interpFFT.py b/simpegPF/Dev/interpFFT.py new file mode 100644 index 00000000..76269f96 --- /dev/null +++ b/simpegPF/Dev/interpFFT.py @@ -0,0 +1,85 @@ +def interpFFT(x,y,m): + """ + Load in a 2D grid and resample + + OUTPUT: + m_out + + + """ + + from SimPEG import np, sp + import scipy.signal as sn + + + # Add padding values by reflection (2**n) + lenx = np.round( np.log2( 2*len(x) ) ) + npadx = int(np.floor( ( 2**lenx - len(x) ) /2. )) + + #Create hemming taper + if np.mod(npadx*2+len(x),2) != 0: + oddx = 1 + + else: + oddx = 0 + + tap0 = sn.hamming(npadx*2) + tapl = sp.spdiags(tap0[0:npadx],0,npadx,npadx) + tapr = sp.spdiags(tap0[npadx:],0,npadx,npadx+oddx) + + # Mirror the 2d data over the half lenght and apply 0-taper + mpad = np.hstack([np.fliplr(m[:,0:npadx]) * tapl, m, np.fliplr(m[:,-npadx:]) * tapr]) + + # Repeat along the second dimension + leny = np.round( np.log2( 2*len(y) ) ) + npady = int(np.floor( ( 2**leny - len(y) ) /2. )) + + #Create hemming taper + if np.mod(npady*2+len(y),2) != 0: + oddy = 1 + + else: + oddy = 0 + + tap0 = sn.hamming(npady*2) + tapu = sp.spdiags(tap0[0:npady],0,npady,npady) + tapd = sp.spdiags(tap0[npady:],0,npady+oddy,npady) + + mpad = np.vstack([tapu*np.flipud(mpad[0:npady,:]), mpad, tapd*np.flipud(mpad[-npady:,:])]) + + # Compute FFT + FFTm = np.fft.fft2(mpad) + + # Do an FFT shift + FFTshift = np.fft.fftshift(FFTm) + + # Pad high frequencies with zeros to increase the sampling rate + py = int(FFTm.shape[0]/2) + px = int(FFTm.shape[1]/2) + + FFTshift = np.hstack([np.zeros((FFTshift.shape[0],px)),FFTshift,np.zeros((FFTshift.shape[0],px))]) + FFTshift = np.vstack([np.zeros((py,FFTshift.shape[1])),FFTshift,np.zeros((py,FFTshift.shape[1]))]) + + # Inverse shift + FFTm = np.fft.ifftshift(FFTshift) + + + # Compute inverse FFT + IFFTm = np.fft.ifft2(FFTm)*FFTm.size/mpad.size + + + m_out = np.real(IFFTm) + # Extract core + #m_out = np.real(IFFTm[npady*2:-(npady*2+oddy+1),npadx*2:-(npadx*2+oddx+1)]) + + + + m_out = m_out[npady*2:-(npady+oddy)*2,npadx*2:-(npadx+oddx)*2] + + if np.mod(m.shape[0],2) != 0: + m_out = m_out[:-1,:] + + if np.mod(m.shape[1],2) != 0: + m_out = m_out[:,:-1] + + return m_out