diff --git a/README.rst b/README.rst index 1c645334..c9b30a78 100644 --- a/README.rst +++ b/README.rst @@ -25,6 +25,10 @@ SimPEG :target: https://coveralls.io/r/simpeg/simpeg?branch=master :alt: Coverage status +.. image:: http://img.shields.io/badge/GITTER-JOIN_CHAT-brightgreen.svg?style=flat-square + :alt: gitter chat room at https://gitter.im/simpeg/simpeg + :target: https://gitter.im/simpeg/simpeg + Simulation and Parameter Estimation in Geophysics - A python package for simulation and gradient based parameter estimation in the context of geophysical applications. The vision is to create a package for finite volume simulation with applications to geophysical imaging and subsurface flow. To enable the understanding of the many different components, this package has the following features: diff --git a/SimPEG/DCIP/DCIPUtils.py b/SimPEG/DCIP/DCIPUtils.py index 90bbbca0..97c30e89 100644 --- a/SimPEG/DCIP/DCIPUtils.py +++ b/SimPEG/DCIP/DCIPUtils.py @@ -169,7 +169,7 @@ def readUBC_DC2DModel(fileName): return model -def plot_pseudoSection(DCsurvey, axs, stype): +def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None, cblabel=True, axlabel = True, colorbar = True, contour = None): """ Read list of 2D tx-rx location and plot a speudo-section of apparent resistivity. @@ -179,7 +179,7 @@ def plot_pseudoSection(DCsurvey, axs, stype): Input: :param d2D, z0 :switch stype -> Either 'pdp' (pole-dipole) | 'dpdp' (dipole-dipole) - + :switch dtype=-> Either 'appr' (app. res) | 'appc' (app. con) | 'volt' (potential) Output: :figure scatter plot overlayed on image @@ -192,9 +192,6 @@ def plot_pseudoSection(DCsurvey, axs, stype): from scipy.interpolate import griddata import pylab as plt - # Set depth to 0 for now - z0 = 0. - # Pre-allocate midx = [] midz = [] @@ -221,47 +218,97 @@ def plot_pseudoSection(DCsurvey, axs, stype): Cmid = (Tx[0][0] + Tx[1][0])/2 Pmid = (Rx[0][:,0] + Rx[1][:,0])/2 - # Compute pant leg of apparent rho - if stype == 'pdp': - leg = data * 2*np.pi * MA * ( MA + MN ) / MN + # Change output for dtype + if dtype == 'volt': - leg = np.log10(abs(1/leg)) + rho = np.hstack([rho,data]) - elif stype == 'dpdp': - leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA ) + else: + # Compute pant leg of apparent rho + if stype == 'pdp': + + leg = data * 2*np.pi * MA * ( MA + MN ) / MN + + elif stype == 'dpdp': + + leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA ) + + else: + print """dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """ + break + + + if dtype == 'appc': + + leg = np.log10(abs(1./leg)) + rho = np.hstack([rho,leg]) + + elif dtype == 'appr': + + leg = np.log10(abs(leg)) + rho = np.hstack([rho,leg]) + + else: + print """dtype must be 'appr' | 'appc' | 'volt' """ + break midx = np.hstack([midx, ( Cmid + Pmid )/2 ]) - midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + z0 ]) - rho = np.hstack([rho,leg]) - - - ax = axs + midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2 ]) # Grid points grid_x, grid_z = np.mgrid[np.min(midx):np.max(midx), np.min(midz):np.max(midz)] grid_rho = griddata(np.c_[midx,midz], rho.T, (grid_x, grid_z), method='linear') + + # Scale the color scheme + if clim == None: + vmin, vmax = rho.min(), rho.max() + else: + vmin, vmax = clim[0], clim[1] + + # Plot data + grid_rho = np.ma.masked_where(np.isnan(grid_rho), grid_rho) + ph = plt.pcolormesh(grid_x[:,0],grid_z[0,:],grid_rho.T, vmin = vmin, vmax = vmax) + plt.gca().tick_params(axis='both', which='major', labelsize=8) + + if contour is not None: + plt.contour(grid_x,grid_z,grid_rho,levels = contour,colors = 'r', vmin = vmin, vmax = vmax) + + # Add scatter points + axs.scatter(midx,midz,s=10,c=rho.T, vmin = vmin, vmax = vmax) + + if colorbar: + + if dtype == 'volt': + cbar = plt.colorbar(ph, ax = axs, format="%4.1f",fraction=0.04,orientation="horizontal") + + else: + cbar = plt.colorbar(ph, ax = axs, format="$10^{%.1f}$",fraction=0.04,orientation="horizontal") + + cmin,cmax = cbar.get_clim() + ticks = np.linspace(cmin,cmax,3) + cbar.set_ticks(ticks) + cbar.ax.tick_params(labelsize=10) + + if cblabel: + if dtype == 'appc': + cbar.set_label("App.Cond",size=12) + elif dtype == 'appr': + cbar.set_label("App.Res.",size=12) + elif dtype == 'volt': + cbar.set_label("Potential (V)",size=12) - plt.imshow(grid_rho.T, extent = (np.min(midx),np.max(midx),np.min(midz),np.max(midz)), origin='lower', alpha=0.8, vmin = np.min(rho), vmax = np.max(rho)) - cbar = plt.colorbar(format = '%.2f',fraction=0.04,orientation="horizontal") - cmin,cmax = cbar.get_clim() - ticks = np.linspace(cmin,cmax,3) - cbar.set_ticks(ticks) + if not axlabel: + axs.set_xticklabels([]) + axs.set_yticklabels([]) - # Plot apparent resistivity - plt.scatter(midx,midz,s=50,c=rho.T) - - ax.set_xticklabels([]) - - ax.set_ylabel('Z') - ax.yaxis.tick_right() - ax.yaxis.set_label_position('right') plt.gca().set_aspect('equal', adjustable='box') - return ax + + return ph def gen_DCIPsurvey(endl, mesh, stype, a, b, n): """ @@ -361,16 +408,6 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n): srcClass = DC.SrcDipole([rxClass], M[ii,:],M[ii,:]) SrcList.append(srcClass) -#============================================================================== -# elif re.match(stype,'dpdp'): -# -# for ii in range(0, int(nstn)-2): -# -# indx = np.min([ii+n+1,nstn]) -# Tx.append(np.c_[M[ii,:],N[ii,:]]) -# Rx.append(np.c_[M[ii+2:indx,:],N[ii+2:indx,:]]) -#============================================================================== - elif stype == 'gradient': # Gradient survey only requires Tx at end of line and creates a square @@ -423,15 +460,15 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n): survey = DC.SurveyDC(SrcList) return survey, Tx, Rx -def writeUBC_DCobs(fileName, DCsurvey, dtype, stype): +def writeUBC_DCobs(fileName, DCsurvey, dtype='3D', stype='SURFACE', iptype = 0): """ Write UBC GIF DCIP 2D or 3D observation file Input: - :string fileName -> including path where the file is written out - :DCsurvey -> DC survey class object - :string dtype -> either '2D' | '3D' - :string stype -> either 'SURFACE' | 'GENERAL' + :string fileName -> including path where the file is written out + :DCsurvey DC survey class object + :string dtype -> either '2D' | '3D' + :string stype -> either 'SURFACE' | 'GENERAL' Output: :param UBC2D-Data file @@ -446,10 +483,16 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype): assert (dtype=='2D') | (dtype=='3D'), "Data must be either '2D' | '3D'" assert (stype=='SURFACE') | (stype=='GENERAL') | (stype=='SIMPLE'), "Data must be either 'SURFACE' | 'GENERAL' | 'SIMPLE'" - + fid = open(fileName,'w') - fid.write('! ' + stype + ' FORMAT\n') - + + + if iptype!=0: + fid.write('IPTYPE=%i\n'%iptype) + + else: + fid.write('! ' + stype + ' FORMAT\n') + count = 0 for ii in range(DCsurvey.nSrc): @@ -473,7 +516,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype): B = np.repeat(tx[0,1],M.shape[0],axis=0) M = M[:,0] N = N[:,0] - + np.savetxt(fid, np.c_[A, B, M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n') @@ -481,18 +524,25 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype): if stype == 'SURFACE': - fid.writelines("%e " % ii for ii in mkvc(tx[0,:])) + fid.writelines("%f " % ii for ii in mkvc(tx[0,:])) M = M[:,0] N = N[:,0] if stype == 'GENERAL': + # Flip sign for z-elevation to depth + tx[2::2,:] = -tx[2::2,:] + fid.writelines("%e " % ii for ii in mkvc(tx[::2,:])) M = M[:,0::2] N = N[:,0::2] + # Flip sign for z-elevation to depth + M[:,1::2] = -M[:,1::2] + N[:,1::2] = -N[:,1::2] + fid.write('%i\n'% nD) - np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n') + np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%f',delimiter=' ',newline='\n') if dtype=='3D': @@ -504,31 +554,32 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype): if stype == 'GENERAL': - fid.writelines("%e " % ii for ii in mkvc(tx)) + fid.writelines("%e " % ii for ii in mkvc(tx[0:3,:])) fid.write('%i\n'% nD) np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n') - + fid.write('\n') + count += nD fid.close() -def convertObs_DC3D_to_2D(DCsurvey,lineID): +def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'): """ - Read DC survey and data and change - coordinate system to distance along line assuming - all data is acquired along line. - First transmitter pole is assumed to be at the origin + Read DC survey and projects the coordinate system + according to the flag = 'Xloc' | 'Yloc' | 'local' (default) + In the 'local' system, station coordinates are referenced + to distance from the first srcLoc[0].loc[0] - Assumes flat topo for now... + The Z value is preserved, but Y coordinates zeroed. Input: - :param Tx, Rx + :param survey3D Output: - :figure Tx2d, Rx2d + :figure survey2D - Edited Feb 17th, 2016 + Edited April 6th, 2016 @author: dominiquef @@ -570,25 +621,39 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID): Rx = DCsurvey.srcList[indx[ii]].rxList[0].locs nrx = Rx[0].shape[0] - # Find A electrode along line - vec, r = r_unit(x0,Tx[ii][0,0:2]) - A = stn_id(vecTx,vec,r) + if flag == 'local': + # Find A electrode along line + vec, r = r_unit(x0,Tx[ii][0,0:2]) + A = stn_id(vecTx,vec,r) - # Find B electrode along line - vec, r = r_unit(x0,Tx[ii][1,0:2]) - B = stn_id(vecTx,vec,r) + # Find B electrode along line + vec, r = r_unit(x0,Tx[ii][1,0:2]) + B = stn_id(vecTx,vec,r) - M = np.zeros(nrx) - N = np.zeros(nrx) - for kk in range(nrx): + M = np.zeros(nrx) + N = np.zeros(nrx) + for kk in range(nrx): - # Find all M electrodes along line - vec, r = r_unit(x0,Rx[0][kk,0:2]) - M[kk] = stn_id(vecTx,vec,r) + # Find all M electrodes along line + vec, r = r_unit(x0,Rx[0][kk,0:2]) + M[kk] = stn_id(vecTx,vec,r) - # Find all N electrodes along line - vec, r = r_unit(x0,Rx[1][kk,0:2]) - N[kk] = stn_id(vecTx,vec,r) + # Find all N electrodes along line + vec, r = r_unit(x0,Rx[1][kk,0:2]) + N[kk] = stn_id(vecTx,vec,r) + elif flag == 'Yloc': + """ Flip the XY axis locs""" + A = Tx[ii][0,1] + B = Tx[ii][1,1] + M = Rx[0][:,1] + N = Rx[1][:,1] + + elif flag == 'Xloc': + """ Copy the rx-tx locs""" + A = Tx[ii][0,0] + B = Tx[ii][1,0] + M = Rx[0][:,0] + N = Rx[1][:,0] Rx = DC.RxDipole(np.c_[M,np.zeros(nrx),Rx[0][:,2]],np.c_[N,np.zeros(nrx),Rx[1][:,2]]) @@ -601,51 +666,59 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID): DCsurvey2D.std = np.asarray(DCsurvey.std) return DCsurvey2D - -def readUBC_DC3Dobs(fileName): + +def readUBC_DC3Dobs(fileName, dtype = 'DC'): """ - Read UBC GIF DCIP 3D observation file and generate arrays for tx-rx location + Read UBC GIF IP 3D observation file and generate survey Input: :param fileName, path to the UBC GIF 3D obs file Output: - :param rx, tx, d, wd + :param IPsurvey :return - Created on Mon December 7th, 2015 - @author: dominiquef """ - + zflag = True # Flag for z value provided + # Load file - obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!') - + if dtype == 'IP': + obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='IPTYPE') + + elif dtype == 'DC': + obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!') + + else: + print "dtype must be 'DC'(default) | 'IP'" + # Pre-allocate srcLists = [] Rx = [] d = [] wd = [] - zflag = True # Flag for z value provided + # Countdown for number of obs/tx count = 0 for ii in range(obsfile.shape[0]): + # Skip if blank line if not obsfile[ii]: continue - # First line is transmitter with number of receivers + # First line or end of a transmitter block, read transmitter info if count==0: - - temp = (np.fromstring(obsfile[ii], dtype=float,sep=' ').T) + # Read the line + temp = (np.fromstring(obsfile[ii], dtype=float, sep=' ').T) count = int(temp[-1]) # Check if z value is provided, if False -> nan if len(temp)==5: - tx = np.r_[temp[0:2],np.nan,temp[0:2],np.nan] - zflag = False + tx = np.r_[temp[0:2],np.nan,temp[2:4],np.nan] + + zflag = False # Pass on the flag to the receiver loc else: tx = temp[:-1] @@ -653,8 +726,16 @@ def readUBC_DC3Dobs(fileName): rx = [] continue - temp = np.fromstring(obsfile[ii], dtype=float,sep=' ') + temp = np.fromstring(obsfile[ii], dtype=float,sep=' ') # Get the string + # Filter out negative IP +# if temp[-2] < 0: +# count = count -1 +# print "Negative!" +# +# else: + + # If the Z-location is provided, otherwise put nan if zflag: rx.append(temp[:-2]) @@ -664,7 +745,7 @@ def readUBC_DC3Dobs(fileName): wd.append(temp[-1]) else: - rx.append(np.r_[temp[0:2],np.nan,temp[0:2],np.nan] ) + rx.append(np.r_[temp[0:2],np.nan,temp[2:4],np.nan] ) # Check if there is data with the location if len(temp)==6: d.append(temp[-2]) @@ -672,7 +753,7 @@ def readUBC_DC3Dobs(fileName): count = count -1 - # Reach the end of transmitter block + # Reach the end of transmitter block, append the src, rx and continue if count == 0: rx = np.asarray(rx) Rx = DC.RxDipole(rx[:,:3],rx[:,3:]) @@ -688,6 +769,7 @@ def readUBC_DC3Dobs(fileName): def readUBC_DC2Dobs(fileName): """ + ------- NEEDS TO BE UPDATED ------ Read UBC GIF 2D observation file and generate arrays for tx-rx location Input: @@ -735,6 +817,73 @@ def readUBC_DC2Dobs(fileName): return tx, rx, d, wd +def readUBC_DC2Dpre(fileName): + """ + Read UBC GIF DCIP 2D observation file and generate arrays for tx-rx location + + Input: + :param fileName, path to the UBC GIF 3D obs file + + Output: + DCsurvey + :return + + Created on Mon March 9th, 2016 << Doug's 70th Birthday !! >> + + @author: dominiquef + + """ + + # Load file + obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!') + + # Pre-allocate + srcLists = [] + Rx = [] + d = [] + zflag = True # Flag for z value provided + + for ii in range(obsfile.shape[0]): + + if not obsfile[ii]: + continue + + # First line is transmitter with number of receivers + + + temp = (np.fromstring(obsfile[ii], dtype=float,sep=' ').T) + + + # Check if z value is provided, if False -> nan + if len(temp)==5: + tx = np.r_[temp[0],np.nan,np.nan,temp[1],np.nan,np.nan] + zflag = False + + else: + tx = np.r_[temp[0],np.nan,temp[1],temp[2],np.nan,temp[3]] + + + if zflag: + rx = np.c_[temp[4],np.nan,temp[5],temp[6],np.nan,temp[7]] + + + else: + rx = np.c_[temp[2],np.nan,np.nan,temp[3],np.nan,np.nan] + # Check if there is data with the location + + d.append(temp[-1]) + + + Rx = DC.RxDipole(rx[:,:3],rx[:,3:]) + srcLists.append( DC.SrcDipole( [Rx], tx[:3],tx[3:]) ) + + # Create survey class + survey = DC.SurveyDC(srcLists) + + survey.dobs = np.asarray(d) + + return {'DCsurvey':survey} + def readUBC_DC2DMesh(fileName): """ Read UBC GIF 2DTensor mesh and generate 2D Tensor mesh in simpeg @@ -928,7 +1077,6 @@ def getSrc_locs(DCsurvey): srcMat = np.zeros((DCsurvey.nSrc,2,3)) for ii in range(DCsurvey.nSrc): - print np.asarray(DCsurvey.srcList[ii].loc).shape srcMat[ii,:,:] = np.asarray(DCsurvey.srcList[ii].loc) return srcMat diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 3d726712..4fc1ffc5 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -216,13 +216,13 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration): # Save the data. ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) ) phi_ms = 0.5*ms.dot(ms) - if self.reg.smoothModel == True: + if self.reg.mrefInSmooth == True: mref = self.reg.mref else: mref = 0 mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) ) phi_mx = 0.5 * mx.dot(mx) - if self.prob.mesh.dim==2: + if self.prob.mesh.dim >= 2: my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) ) phi_my = 0.5 * my.dot(my) else: @@ -237,40 +237,6 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration): # Save the file as a npz np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred) -class SaveOutputDictEveryIteration(_SaveEveryIteration): - """SaveOutputDictEveryIteration - A directive that saves some relevant information from the inversion run to a numpy .npz dictionary file (see numpy.savez function for further info). - """ - - def initialize(self): - print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '%s-###.npz'"%self.fileName - - def endIter(self): - # Save the data. - ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) ) - phi_ms = 0.5*ms.dot(ms) - if self.reg.smoothModel == True: - mref = self.reg.mref - else: - mref = 0 - mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) ) - phi_mx = 0.5 * mx.dot(mx) - if self.prob.mesh.dim==2: - my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) ) - phi_my = 0.5 * my.dot(my) - else: - phi_my = 'NaN' - if self.prob.mesh.dim==3 and 'CYL' not in self.prob.mesh._meshType: - mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) ) - phi_mz = 0.5 * mz.dot(mz) - else: - phi_mz = 'NaN' - - - # Save the file as a npz - np.savez('{:s}-{:03d}'.format(self.fileName,self.opt.iter), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred) - - # class UpdateReferenceModel(Parameter): @@ -283,3 +249,125 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration): # mref = self.mref0 # self.m_prev = self.invProb.m_current # return mref + +class Update_IRLS(InversionDirective): + + eps_min = None + factor = None + gamma = None + phi_m_last = None + phi_d_last = None + + + def initialize(self): + + # Scale the regularization for changes in norm + if getattr(self, 'phi_m_last', None) is not None: + + self.reg.curModel = self.invProb.curModel + self.reg.gamma = 1. + phim_new = self.reg.eval(self.invProb.curModel) + self.gamma = self.phi_m_last / phim_new + + self.reg.curModel = self.invProb.curModel + self.reg.gamma = self.gamma + + if getattr(self, 'phi_d_last', None) is None: + self.phi_d_last = self.invProb.phi_d + + def endIter(self): + # Cool the threshold parameter if required + if getattr(self, 'factor', None) is not None: + eps = self.reg.eps / self.factor + + if getattr(self, 'eps_min', None) is not None: + self.reg.eps = np.max([self.eps_min,eps]) + else: + self.reg.eps = eps + + # Get phi_m at the end of current iteration + self.phi_m_last = self.invProb.phi_m_last + + # Update the model used for the IRLS weights + self.reg.curModel = self.invProb.curModel + + # Temporarely set gamma to 1. to get raw phi_m + self.reg.gamma = 1. + + # Compute new model objective function value + phim_new = self.reg.eval(self.invProb.curModel) + + # Update gamma to scale the regularization between IRLS iterations + self.reg.gamma = self.phi_m_last / phim_new + + # Set the weighting matrix to None so that it is recomputed next time + # it is called in the inversion + self.reg._W = None + +class Update_lin_PreCond(InversionDirective): + """ + Create a Jacobi preconditioner for the linear problem + """ + onlyOnStart=False + + def initialize(self): + + if getattr(self.opt, 'approxHinv', None) is None: + # Update the pre-conditioner + diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() #* (self.reg.mapping * np.ones(self.reg.curModel.size))**2. + PC = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.) + self.opt.approxHinv = PC + + def endIter(self): + # Cool the threshold parameter + if self.onlyOnStart==True: + return + + if getattr(self.opt, 'approxHinv', None) is not None: + # Update the pre-conditioner + diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() #* (self.reg.mapping * np.ones(self.reg.curModel.size))**2. + PC = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.) + self.opt.approxHinv = PC + + +class Update_Wj(InversionDirective): + """ + Create approx-sensitivity base weighting using the probing method + """ + k = None # Number of probing cycles + itr = None # Iteration number to update Wj, or always update if None + + def endIter(self): + + if self.itr is None or self.itr == self.opt.iter: + + m = self.invProb.curModel + if self.k is None: + self.k = int(self.survey.nD/10) + + def JtJv(v): + + Jv = self.prob.Jvec(m, v) + + return self.prob.Jtvec(m,Jv) + + JtJdiag = Utils.diagEst(JtJv,len(m),k=self.k) + JtJdiag = JtJdiag / max(JtJdiag) + + self.reg.wght = JtJdiag + +class Scale_Beta(InversionDirective): + """ + Instead of a linear cooling schedule, beta is allowed to change based + on the ratio between the target misfit and the current data misfit. The + update is done only if the misfit is outside some threshold bounds. + """ + tol = 0.05 + + def endIter(self): + + # Check if misfit is within the tolerance, otherwise adjust beta + val = self.invProb.phi_d / (self.survey.nD*0.5) + + if np.abs(1.-val) > self.tol: + self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d diff --git a/SimPEG/Examples/DC_Forward_PseudoSection.py b/SimPEG/Examples/DC_Forward_PseudoSection.py index 2f198238..53467826 100644 --- a/SimPEG/Examples/DC_Forward_PseudoSection.py +++ b/SimPEG/Examples/DC_Forward_PseudoSection.py @@ -2,19 +2,27 @@ from SimPEG import Mesh, Utils, np, sp import SimPEG.DCIP as DC import time -def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True): +def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', plotIt=True): """ DC Forward Simulation ===================== - Forward model conductive spheres in a half-space and plot a pseudo-section + Forward model two conductive spheres in a half-space and plot a + pseudo-section. Assumes an infinite line source and measures along the + center of the spheres. - Created by @fourndo on Mon Feb 01 19:28:06 2016 + INPUT: + loc = Location of spheres [[x1,y1,z1],[x2,y2,z2]] + radi = Radius of spheres [r1,r2] + param = Conductivity of background and two spheres [m0,m1,m2] + stype = survey type "pdp" (pole dipole) or "dpdp" (dipole dipole) + dtype = Data type "appr" (app res) | "appc" (app cond) | "volt" (potential) + Created by @fourndo """ assert stype in ['pdp', 'dpdp'], "Source type (stype) must be pdp or dpdp (pole dipole or dipole dipole)" - + assert dtype in ['appr', 'appc', 'volt'], "Data type (dtype) must be appr (app res) or appc (app cond) or volt (potential)" if loc is None: loc = np.c_[[-50.,0.,-50.],[50.,0.,-50.]] @@ -27,7 +35,6 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True): # First we need to create a mesh and a model. - # This is our mesh dx = 5. @@ -52,14 +59,10 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True): # Get index of the center indy = int(mesh.nCy/2) - # Plot the model for reference # Define core mesh extent xlim = 200 - zlim = 125 - - # Specify the survey type: "pdp" | "dpdp" - + zlim = 100 # Then specify the end points of the survey. Let's keep it simple for now and survey above the anomalies, top of the mesh ends = [(-175,0),(175,0)] @@ -77,12 +80,13 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True): dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) ) dl_x = ( Tx[-1][0,1] - Tx[0][0,0] ) / dl_len dl_y = ( Tx[-1][1,1] - Tx[0][1,0] ) / dl_len - azm = np.arctan(dl_y/dl_x) + #azm = np.arctan(dl_y/dl_x) #Set boundary conditions mesh.setCellGradBC('neumann') - # Define the differential operators needed for the DC problem + # Define the linear system needed for the DC problem. We assume an infitite + # line source for simplicity. Div = mesh.faceDiv Grad = mesh.cellGrad Msig = Utils.sdiag(1./(mesh.aveF2CC.T*(1./model))) @@ -145,16 +149,23 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True): print 'Forward completed' # Let's just convert the 3D format into 2D (distance along line) and plot - # [Tx2d, Rx2d] = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc)) - survey2D = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc)) + survey2D = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc) , 'Xloc') survey2D.dobs =np.hstack(data) - # Here is an example for the first tx-rx array + if plotIt: import matplotlib.pyplot as plt - fig = plt.figure() + fig = plt.figure(figsize=(7,7)) ax = plt.subplot(2,1,1, aspect='equal') - mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y', ind = indy,grid=True) - ax.set_title('E-W section at '+str(mesh.vectorCCy[indy])+' m') + # Plot the location of the spheres for reference + circle1=plt.Circle((loc[0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3) + circle2=plt.Circle((loc[0,1],loc[2,1]),radi[1],color='k',fill=False, lw=3) + ax.add_artist(circle1) + ax.add_artist(circle2) + + dat = mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y', + ind = indy,grid=True, clim = np.log10([sig.min(),sig.max()])) + + ax.set_title('3-D model') plt.gca().set_aspect('equal', adjustable='box') plt.scatter(Tx[0][0,:],Tx[0][2,:],s=40,c='g', marker='v') @@ -163,22 +174,34 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True): plt.ylim([-zlim,mesh.vectorNz[-1]+dx]) - ax = plt.subplot(2,1,2, aspect='equal') + pos = ax.get_position() + ax.set_position([pos.x0 , pos.y0 + 0.025 , pos.width, pos.height]) + pos = ax.get_position() + cbarax = fig.add_axes([pos.x0 , pos.y0 + 0.025 , pos.width, pos.height * 0.04]) ## the parameters are the specified position you set + cb = fig.colorbar(dat[0],cax=cbarax, orientation="horizontal", + ax = ax, ticks=np.linspace(np.log10(sig.min()), + np.log10(sig.max()), 3), format="$10^{%.1f}$") + cb.set_label("Conductivity (S/m)",size=12) + cb.ax.tick_params(labelsize=12) + + # Second plot for the predicted apparent resistivity data + ax2 = plt.subplot(2,1,2, aspect='equal') # Plot the location of the spheres for reference - circle1=plt.Circle((loc[0,0]-Tx[0][0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3) - circle2=plt.Circle((loc[0,1]-Tx[0][0,0],loc[2,1]),radi[1],color='k',fill=False, lw=3) - ax.add_artist(circle1) - ax.add_artist(circle2) + circle1=plt.Circle((loc[0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3) + circle2=plt.Circle((loc[0,1],loc[2,1]),radi[1],color='k',fill=False, lw=3) + ax2.add_artist(circle1) + ax2.add_artist(circle2) # Add the speudo section - DC.plot_pseudoSection(survey2D,ax,stype) + dat = DC.plot_pseudoSection(survey2D,ax2,stype=stype, dtype = dtype) # plt.scatter(Tx2d[0][:],Tx[0][2,:],s=40,c='g', marker='v') # plt.scatter(Rx2d[0][:],Rx[0][:,2::3],s=40,c='y') # plt.plot(np.r_[Tx2d[0][0],Rx2d[-1][-1,-1]],np.ones(2)*mesh.vectorNz[-1], color='k') - plt.ylim([-zlim,mesh.vectorNz[-1]+dx]) + ax2.set_title('Apparent Conductivity data') + plt.ylim([-zlim,mesh.vectorNz[-1]+dx]) plt.show() return fig, ax diff --git a/SimPEG/Examples/EM_FDEM_1D_Inversion.py b/SimPEG/Examples/EM_FDEM_1D_Inversion.py index 1764f28c..0a5c59a0 100644 --- a/SimPEG/Examples/EM_FDEM_1D_Inversion.py +++ b/SimPEG/Examples/EM_FDEM_1D_Inversion.py @@ -48,8 +48,7 @@ def run(plotIt=True): freqs = np.logspace(1,3,10) srcLoc = np.array([0., 0., 10.]) - srcList = [] - [srcList.append(EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z')) for freq in freqs] + srcList = [EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z') for freq in freqs] survey = EM.FDEM.Survey(srcList) prb = EM.FDEM.Problem3D_b(mesh, mapping=mapping) diff --git a/SimPEG/Examples/EM_Schenkel_Morrison_Casing.py b/SimPEG/Examples/EM_Schenkel_Morrison_Casing.py new file mode 100644 index 00000000..76af4a3d --- /dev/null +++ b/SimPEG/Examples/EM_Schenkel_Morrison_Casing.py @@ -0,0 +1,275 @@ +from SimPEG import * +from SimPEG.EM import FDEM, Analytics, mu_0 +import time + +try: + from pymatsolver import MumpsSolver + solver = MumpsSolver +except Exception: + solver = SolverLU + pass + +def run(plotIt=True): + """ + EM: Schenkel and Morrison Casing Model + ====================================== + + Here we create and run a FDEM forward simulation to calculate the vertical + current inside a steel-cased. The model is based on the Schenkel and + Morrison Casing Model, and the results are used in a 2016 SEG abstract by + Yang et al. + + - Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on potential field measurements using downhole current sources: Geophysical prospecting, 38, 663-686. + + + The model consists of: + - Air: Conductivity 1e-8 S/m, above z = 0 + - Background: conductivity 1e-2 S/m, below z = 0 + - Casing: conductivity 1e6 S/m + - 300m long + - radius of 0.1m + - thickness of 6e-3m + + Inside the casing, we take the same conductivity as the background. + + We are using an EM code to simulate DC, so we use frequency low enough + that the skin depth inside the casing is longer than the casing length (f + = 1e-6 Hz). The plot produced is of the current inside the casing. + + These results are shown in the SEG abstract by Yang et al., 2016: 3D DC + resistivity modeling of steel casing for reservoir monitoring using + equivalent resistor network. The solver used to produce these results and + achieve the CPU time of ~30s is Mumps, which was installed using pymatsolver_ + + .. _pymatsolver: https://github.com/rowanc1/pymatsolver + + This example is on figshare: https://dx.doi.org/10.6084/m9.figshare.3126961.v1 + + If you would use this example for a code comparison, or build upon it, a + citation would be much appreciated! + + """ + + if plotIt: + import matplotlib.pylab as plt + + # ------------------ MODEL ------------------ + sigmaair = 1e-8 # air + sigmaback = 1e-2 # background + sigmacasing = 1e6 # casing + sigmainside = sigmaback # inside the casing + + + casing_t = 0.006 # 1cm thickness + casing_l = 300 # length of the casing + + casing_r = 0.1 + casing_a = casing_r - casing_t/2. # inner radius + casing_b = casing_r + casing_t/2. # outer radius + casing_z = np.r_[-casing_l,0.] + + + # ------------------ SURVEY PARAMETERS ------------------ + freqs = np.r_[1e-6] #[1e-1, 1, 5] # frequencies + dsz = -300 # down-hole z source location + src_loc = np.r_[0.,0.,dsz] + inf_loc = np.r_[0.,0.,1e4] + + print 'Skin Depth: ', [(500./np.sqrt(sigmaback*_)) for _ in freqs] + + + # ------------------ MESH ------------------ + # fine cells near well bore + csx1, csx2 = 2e-3, 60. + pfx1, pfx2 = 1.3, 1.3 + ncx1 = np.ceil(casing_b/csx1+2) + + # pad nicely to second cell size + npadx1 = np.floor(np.log(csx2/csx1) / np.log(pfx1)) + hx1a,hx1b = Utils.meshTensor([(csx1,ncx1)]),Utils.meshTensor([(csx1,npadx1,pfx1)]) + dx1 = sum(hx1a)+sum(hx1b) + dx1 = np.floor(dx1/csx2) + hx1b *= (dx1*csx2 - sum(hx1a))/sum(hx1b) + + # second chunk of mesh + dx2 = 300. # uniform mesh out to here + ncx2 = np.ceil((dx2 - dx1)/csx2) + npadx2 = 45 + hx2a, hx2b = Utils.meshTensor([(csx2,ncx2)]), Utils.meshTensor([(csx2,npadx2,pfx2)]) + hx = np.hstack([hx1a,hx1b,hx2a,hx2b]) + + # z-direction + csz = 0.05 + nza = 10 + ncz, npadzu, npadzd = np.int(np.ceil(np.diff(casing_z)[0]/csz))+10, 68, 68 # cell size, number of core cells, number of padding cells in the x- direction + hz = Utils.meshTensor([(csz,npadzd,-1.3), (csz,ncz), (csz,npadzu,1.3)]) # vector of cell widths in the z-direction + + # Mesh + mesh = Mesh.CylMesh([hx,1.,hz], [0.,0.,-np.sum(hz[:npadzu+ncz-nza])]) + + print 'Mesh Extent xmax: %f,: zmin: %f, zmax: %f'%(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max()) + print 'Number of cells', mesh.nC + + if plotIt is True: + fig, ax = plt.subplots(1, 1, figsize=(6, 4)) + ax.set_title('Simulation Mesh') + mesh.plotGrid(ax=ax) + plt.show() + + # Put the model on the mesh + sigWholespace = sigmaback*np.ones((mesh.nC)) + + sigBack = sigWholespace.copy() + sigBack[mesh.gridCC[:,2] > 0.] = sigmaair + + sigCasing = sigBack.copy() + iCasingZ = (mesh.gridCC[:,2] <= casing_z[1]) & (mesh.gridCC[:,2] >= casing_z[0]) + iCasingX = (mesh.gridCC[:,0] >= casing_a) & (mesh.gridCC[:,0] <= casing_b) + iCasing = iCasingX & iCasingZ + sigCasing[iCasing] = sigmacasing + + + if plotIt is True: + + # plotting parameters + xlim = np.r_[0., 0.2] + zlim = np.r_[-350., 10.] + clim_sig = np.r_[-8,6] + + # plot models + fig, ax = plt.subplots(1,1,figsize=(4,4)) + + f = plt.colorbar(mesh.plotImage(np.log10(sigCasing),ax=ax)[0], ax=ax) + ax.grid(which='both') + ax.set_title('Log_10 (Sigma)') + ax.set_xlim(xlim) + ax.set_ylim(zlim) + f.set_clim(clim_sig) + + plt.show() + + + # -------------- Sources -------------------- + # Define Custom Current Sources + + # surface source + sg_x = np.zeros(mesh.vnF[0],dtype=complex) + sg_y = np.zeros(mesh.vnF[1],dtype=complex) + sg_z = np.zeros(mesh.vnF[2],dtype=complex) + + nza = 2 # put the wire two cells above the surface + ncin = 2 + + # vertically directed wire + sgv_indx = (mesh.gridFz[:,0] > casing_a) & (mesh.gridFz[:,0] < casing_a + csx1) # hook it up to casing at the surface + sgv_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] >= -csz*2) + sgv_ind = sgv_indx & sgv_indz + sg_z[sgv_ind] = -1. + + # horizontally directed wire + sgh_indx = (mesh.gridFx[:,0] > casing_a) & (mesh.gridFx[:,0] <= inf_loc[2]) + sgh_indz = (mesh.gridFx[:,2] > csz*(nza-0.5)) & (mesh.gridFx[:,2] < csz*(nza+0.5)) + sgh_ind = sgh_indx & sgh_indz + sg_x[sgh_ind] = -1. + + sgv2_indx = (mesh.gridFz[:,0] >= mesh.gridFx[sgh_ind,0].max()) & (mesh.gridFz[:,0] <= inf_loc[2]*1.2) # hook it up to casing at the surface + sgv2_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] >= -csz*2) + sgv2_ind = sgv2_indx & sgv2_indz + sg_z[sgv2_ind] = 1. + + # assemble the source + sg = np.hstack([sg_x,sg_y,sg_z]) + sg_p = [FDEM.Src.RawVec_e([],_,sg/mesh.area) for _ in freqs] + + # downhole source + dg_x = np.zeros(mesh.vnF[0],dtype=complex) + dg_y = np.zeros(mesh.vnF[1],dtype=complex) + dg_z = np.zeros(mesh.vnF[2],dtype=complex) + + # vertically directed wire + dgv_indx = (mesh.gridFz[:,0] < csx1) # go through the center of the well + dgv_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] > dsz + csz/2.) + dgv_ind = dgv_indx & dgv_indz + dg_z[dgv_ind] = -1. + + # couple to the casing downhole + dgh_indx = mesh.gridFx[:,0] < casing_a + csx1 + dgh_indz = (mesh.gridFx[:,2] < dsz + csz) & (mesh.gridFx[:,2] >= dsz) + dgh_ind = dgh_indx & dgh_indz + dg_x[dgh_ind] = 1. + + # horizontal part at surface + dgh2_indx = mesh.gridFx[:,0] <= inf_loc[2]*1.2 + dgh2_indz = sgh_indz.copy() + dgh2_ind = dgh2_indx & dgh2_indz + dg_x[dgh2_ind] = -1. + + # vertical part at surface + dgv2_ind = sgv2_ind.copy() + dg_z[dgv2_ind] = 1. + + # assemble the source + dg = np.hstack([dg_x,dg_y,dg_z]) + dg_p = [FDEM.Src.RawVec_e([],_,dg/mesh.area) for _ in freqs] + + # ------------ Problem and Survey --------------- + survey = FDEM.Survey(sg_p + dg_p) + mapping = [('sigma', Maps.IdentityMap(mesh))] + problem = FDEM.Problem_h(mesh, mapping=mapping) + problem.pair(survey) + + # ------------- Solve --------------------------- + t0 = time.time() + fieldsCasing = problem.fields(sigCasing) + print 'Time to solve 2 sources', time.time() - t0 + + # Plot current + + # current density + jn0 = fieldsCasing[dg_p,'j'] + jn1 = fieldsCasing[sg_p,'j'] + + # current + in0 = [mesh.area*fieldsCasing[dg_p,'j'][:,i] for i in range(len(freqs))] + in1 = [mesh.area*fieldsCasing[sg_p,'j'][:,i] for i in range(len(freqs))] + + in0 = np.vstack(in0).T + in1 = np.vstack(in1).T + + # integrate to get z-current inside casing + inds_inx = (mesh.gridFz[:,0] >= casing_a) & (mesh.gridFz[:,0] <= casing_b) + inds_inz = (mesh.gridFz[:,2] >= dsz ) & (mesh.gridFz[:,2] <= 0) + inds_fz = inds_inx & inds_inz + + indsx = [False]*mesh.nFx + inds = list(indsx) + list(inds_fz) + + in0_in = in0[np.r_[inds]] + in1_in = in1[np.r_[inds]] + z_in = mesh.gridFz[inds_fz,2] + + in0_in = in0_in.reshape([in0_in.shape[0]/3,3]) + in1_in = in1_in.reshape([in1_in.shape[0]/3,3]) + z_in = z_in.reshape([z_in.shape[0]/3,3]) + + I0 = in0_in.sum(1).real + I1 = in1_in.sum(1).real + z_in = z_in[:,0] + + if plotIt is True: + fig, ax = plt.subplots(1,2,figsize=(12,4)) + + ax[0].plot(z_in,np.absolute(I0), z_in,np.absolute(I1)) + ax[0].legend(['top casing', 'bottom casing'],loc='best') + ax[0].set_title('Magnitude of Vertical Current in Casing') + + ax[1].semilogy(z_in,np.absolute(I0), z_in,np.absolute(I1)) + ax[1].legend(['top casing', 'bottom casing'],loc='best') + ax[1].set_title('Magnitude of Vertical Current in Casing') + ax[1].set_ylim([1e-2, 1.]) + + plt.show() + +if __name__ == '__main__': + run() + diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py new file mode 100644 index 00000000..06ef4be4 --- /dev/null +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -0,0 +1,132 @@ +from SimPEG import * + + +def run(N=200, plotIt=True): + """ + Inversion: Linear Problem + ========================= + + Here we go over the basics of creating a linear problem and inversion. + + """ + + + np.random.seed(1) + + std_noise = 1e-2 + + mesh = Mesh.TensorMesh([N]) + + m0 = np.ones(mesh.nC) * 1e-4 + nk = 10 + jk = np.linspace(1.,nk,nk) + p = -2. + q = 1. + + g = lambda k: np.exp(p*jk[k]*mesh.vectorCCx)*np.cos(np.pi*q*jk[k]*mesh.vectorCCx) + + G = np.empty((nk, mesh.nC)) + + for i in range(nk): + G[i,:] = g(i) + + mtrue = np.zeros(mesh.nC) + mtrue[mesh.vectorCCx > 0.3] = 1. + mtrue[mesh.vectorCCx > 0.45] = -0.5 + mtrue[mesh.vectorCCx > 0.6] = 0 + + + prob = Problem.LinearProblem(mesh, G) + survey = Survey.LinearSurvey() + survey.pair(prob) + survey.dobs = prob.fields(mtrue) + std_noise * np.random.randn(nk) + #survey.makeSyntheticData(mtrue, std=std_noise) + + wd = np.ones(nk) * std_noise + + #print survey.std[0] + #M = prob.mesh + # Distance weighting + wr = np.sum(prob.G**2.,axis=0)**0.5 + wr = ( wr/np.max(wr) ) + + reg = Regularization.Simple(mesh) + reg.wght = wr + + dmis = DataMisfit.l2_DataMisfit(survey) + dmis.Wd = 1./wd + + opt = Optimization.ProjectedGNCG(maxIter=30,lower=-2.,upper=2., maxIterCG= 20, tolCG = 1e-4) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt) + invProb.curModel = m0 + + beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1) + target = Directives.TargetMisfit() + + betaest = Directives.BetaEstimate_ByEig() + inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target]) + + + mrec = inv.run(m0) + ml2 = mrec + print "Final misfit:" + str(invProb.dmisfit.eval(mrec)) + + # Switch regularization to sparse + phim = invProb.phi_m_last + phid = invProb.phi_d + + reg = Regularization.Sparse(mesh) + +#============================================================================== +# fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2)) +# dmdx = reg.mesh.cellDiffxStencil * mrec +# plt.plot(np.sort(dmdx)) +#============================================================================== + + #reg.recModel = mrec + reg.wght = np.ones(mesh.nC) + reg.mref = np.zeros(mesh.nC) + reg.eps_p = 5e-2 + reg.eps_q = 1e-2 + reg.norms = [0., 0., 2., 2.] + reg.wght = wr + + opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 20, tolCG = 1e-3) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.) + beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1) + #betaest = Directives.BetaEstimate_ByEig() + target = Directives.TargetMisfit() + IRLS =Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid ) + + inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS]) + + m0 = mrec + + # Run inversion + mrec = inv.run(m0) + + print "Final misfit:" + str(invProb.dmisfit.eval(mrec)) + + + if plotIt: + import matplotlib.pyplot as plt + + fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2)) + for i in range(prob.G.shape[0]): + axes[0].plot(prob.G[i,:]) + axes[0].set_title('Columns of matrix G') + + axes[1].plot(mesh.vectorCCx, mtrue, 'b-') + axes[1].plot(mesh.vectorCCx, ml2, 'r-') + #axes[1].legend(('True Model', 'Recovered Model')) + axes[1].set_ylim(-1.0,1.25) + + axes[1].plot(mesh.vectorCCx, mrec, 'k-',lw = 2) + axes[1].legend(('True Model', 'Smooth l2-l2', + 'Sparse lp:' + str(reg.norms[0]) + ', lqx:' + str(reg.norms[1]) ), fontsize = 12) + plt.show() + + return prob, survey, mesh, mrec + +if __name__ == '__main__': + run() diff --git a/SimPEG/Examples/MT_1D_ForwardAndInversion.py b/SimPEG/Examples/MT_1D_ForwardAndInversion.py index cf17e77b..69cedd98 100644 --- a/SimPEG/Examples/MT_1D_ForwardAndInversion.py +++ b/SimPEG/Examples/MT_1D_ForwardAndInversion.py @@ -100,7 +100,7 @@ def run(plotIt=True): # Regularization - with a regularization mesh regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]],m1d.x0) reg = simpeg.Regularization.Tikhonov(regMesh) - reg.smoothModel = True + reg.mrefInSmooth = True reg.alpha_s = 1e-7 reg.alpha_x = 1. # Inversion problem diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index cce22296..c67652a2 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -5,9 +5,11 @@ import DC_Analytic_Dipole import DC_Forward_PseudoSection import EM_FDEM_1D_Inversion import EM_FDEM_Analytic_MagDipoleWholespace +import EM_Schenkel_Morrison_Casing import EM_TDEM_1D_Inversion import FLOW_Richards_1D_Celia1990 import Forward_BasicDirectCurrent +import Inversion_IRLS import Inversion_Linear import Mesh_Basic_PlotImage import Mesh_Basic_Types @@ -19,7 +21,7 @@ import Mesh_Tensor_Creation import MT_1D_ForwardAndInversion import MT_3D_Foward -__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward"] +__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward"] ##### AUTOIMPORTS ##### diff --git a/SimPEG/Inversion.py b/SimPEG/Inversion.py index 3fe9cfd0..3e92374f 100644 --- a/SimPEG/Inversion.py +++ b/SimPEG/Inversion.py @@ -33,7 +33,9 @@ class BaseInversion(object): self._directiveList = value self._directiveList.inversion = self - def __init__(self, invProb, directiveList=[], **kwargs): + def __init__(self, invProb, directiveList=None, **kwargs): + if directiveList is None: + directiveList = [] self.directiveList = directiveList Utils.setKwargs(self, **kwargs) diff --git a/SimPEG/MT/Utils/ediFilesUtils.py b/SimPEG/MT/Utils/ediFilesUtils.py index 55818a87..d8f5e0ad 100644 --- a/SimPEG/MT/Utils/ediFilesUtils.py +++ b/SimPEG/MT/Utils/ediFilesUtils.py @@ -7,17 +7,16 @@ from SimPEG.MT.Utils.dataUtils import rec2ndarr # Import modules import numpy as np import os, sys, re -try: - import osr -except ImportError as e: - print 'Could not import osr, missing the gdal package' - pass + class EDIimporter: """ A class to import EDIfiles. """ + + + # Define data converters _impUnitEDI2SI = 4*np.pi*1e-4 # Convert Z[mV/km/nT] (as in EDI)to Z[V/A] SI unit _impUnitSI2EDI = 1./_impUnitEDI2SI # ConvertZ[V/A] SI unit to Z[mV/km/nT] (as in EDI) @@ -26,8 +25,8 @@ class EDIimporter: comps = None # Hidden properties - _outEPSG = None - _2out = None + _outEPSG = None # Project info + _2out = None # The projection operator def __init__(self, EDIfilesList, compList=None, outEPSG=None): @@ -113,6 +112,12 @@ class EDIimporter: # nOutData=length(obj.data); # obj.data(nOutData+1:nOutData+length(TEMP.data),:) = TEMP.data; def _transfromPoints(self,longD,latD): + # Import the coordinate projections + try: + import osr + except ImportError as e: + print 'Could not import osr, missing the gdal package\nCan not project coordinates' + raise e # Coordinates convertor if self._2out is None: src = osr.SpatialReference() diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 5f71d87c..f17d2314 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -759,15 +759,29 @@ class PolyMap(IdentityMap): m = [\sigma_1, \sigma_2, c] + Can take in an actInd vector to account for topography. + """ - def __init__(self, mesh, order, logSigma=True, normal='X'): + def __init__(self, mesh, order, logSigma=True, normal='X', actInd = None): IdentityMap.__init__(self, mesh) self.logSigma = logSigma self.order = order self.normal = normal + self.actInd = actInd + + if getattr(self, 'actInd', None) is None: + self.actInd = range(self.mesh.nC) + self.nC = self.mesh.nC + + else: + self.nC = len(self.actInd) slope = 1e4 + @property + def shape(self): + return (self.nC, self.nP) + @property def nP(self): if np.isscalar(self.order): @@ -785,8 +799,8 @@ class PolyMap(IdentityMap): sig1, sig2 = np.exp(sig1), np.exp(sig2) #2D if self.mesh.dim == 2: - X = self.mesh.gridCC[:,0] - Y = self.mesh.gridCC[:,1] + X = self.mesh.gridCC[self.actInd,0] + Y = self.mesh.gridCC[self.actInd,1] if self.normal =='X': f = polynomial.polyval(Y, c) - X elif self.normal =='Y': @@ -795,9 +809,9 @@ class PolyMap(IdentityMap): raise(Exception("Input for normal = X or Y or Z")) #3D elif self.mesh.dim == 3: - X = self.mesh.gridCC[:,0] - Y = self.mesh.gridCC[:,1] - Z = self.mesh.gridCC[:,2] + X = self.mesh.gridCC[self.actInd,0] + Y = self.mesh.gridCC[self.actInd,1] + Z = self.mesh.gridCC[self.actInd,2] if self.normal =='X': f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X elif self.normal =='Y': @@ -806,6 +820,7 @@ class PolyMap(IdentityMap): f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z else: raise(Exception("Input for normal = X or Y or Z")) + else: raise(Exception("Only supports 2D")) @@ -819,8 +834,8 @@ class PolyMap(IdentityMap): sig1, sig2 = np.exp(sig1), np.exp(sig2) #2D if self.mesh.dim == 2: - X = self.mesh.gridCC[:,0] - Y = self.mesh.gridCC[:,1] + X = self.mesh.gridCC[self.actInd,0] + Y = self.mesh.gridCC[self.actInd,1] if self.normal =='X': f = polynomial.polyval(Y, c) - X @@ -832,9 +847,9 @@ class PolyMap(IdentityMap): raise(Exception("Input for normal = X or Y or Z")) #3D elif self.mesh.dim == 3: - X = self.mesh.gridCC[:,0] - Y = self.mesh.gridCC[:,1] - Z = self.mesh.gridCC[:,2] + X = self.mesh.gridCC[self.actInd,0] + Y = self.mesh.gridCC[self.actInd,1] + Z = self.mesh.gridCC[self.actInd,2] if self.normal =='X': f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X diff --git a/SimPEG/Mesh/CylMesh.py b/SimPEG/Mesh/CylMesh.py index ecdf36ac..f0b6751e 100644 --- a/SimPEG/Mesh/CylMesh.py +++ b/SimPEG/Mesh/CylMesh.py @@ -330,7 +330,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView): raise NotImplementedError('wrapping in the averaging is not yet implemented') return self._aveF2CCV - def getInterpolationMatCartMesh(self, Mrect, locType='CC'): + def getInterpolationMatCartMesh(self, Mrect, locType='CC', locTypeTo=None): """ Takes a cartesian mesh and returns a projection to translate onto the cartesian grid. """ @@ -338,19 +338,22 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView): assert self.isSymmetric, "Currently we have not taken into account other projections for more complicated CylMeshes" + if locTypeTo is None: + locTypeTo = locType + if locType == 'F': # do this three times for each component - X = self.getInterpolationMatCartMesh(Mrect, locType='Fx') - Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy') - Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz') + X = self.getInterpolationMatCartMesh(Mrect, locType='Fx', locTypeTo=locTypeTo+'x') + Y = self.getInterpolationMatCartMesh(Mrect, locType='Fy', locTypeTo=locTypeTo+'y') + Z = self.getInterpolationMatCartMesh(Mrect, locType='Fz', locTypeTo=locTypeTo+'z') return sp.vstack((X,Y,Z)) if locType == 'E': - X = self.getInterpolationMatCartMesh(Mrect, locType='Ex') - Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey') - Z = spzeros(Mrect.nEz, self.nE) + X = self.getInterpolationMatCartMesh(Mrect, locType='Ex', locTypeTo=locTypeTo+'x') + Y = self.getInterpolationMatCartMesh(Mrect, locType='Ey', locTypeTo=locTypeTo+'y') + Z = spzeros(getattr(Mrect, 'n' + locTypeTo + 'z'), self.nE) return sp.vstack((X,Y,Z)) - grid = getattr(Mrect, 'grid' + locType) + grid = getattr(Mrect, 'grid' + locTypeTo) # This is unit circle stuff, 0 to 2*pi, starting at x-axis, rotating counter clockwise in an x-y slice theta = - np.arctan2(grid[:,0] - self.cartesianOrigin[0], grid[:,1] - self.cartesianOrigin[1]) + np.pi/2 theta[theta < 0] += np.pi*2.0 @@ -366,7 +369,7 @@ class CylMesh(BaseTensorMesh, BaseRectangularMesh, InnerProducts, CylView): 'Ex': Mrect.tangents[:Mrect.nEx,:], 'Ey': Mrect.tangents[Mrect.nEx:(Mrect.nEx+Mrect.nEy),:], 'Ez': Mrect.tangents[-Mrect.nEz:,:], - }[locType] + }[locTypeTo] if 'F' in locType: normals = np.c_[np.cos(theta), np.sin(theta), np.zeros(theta.size)] proj = ( normals * dotMe ).sum(axis=1) diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index 793b2581..d0363001 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -307,24 +307,28 @@ class DiffOperators(object): return BC _cellGradBC_list = 'neumann' + def _cellGradStencil(self): + BC = self.setCellGradBC(self._cellGradBC_list) + n = self.vnC + if(self.dim == 1): + G = ddxCellGrad(n[0], BC[0]) + elif(self.dim == 2): + G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0])) + G2 = sp.kron(ddxCellGrad(n[1], BC[1]), speye(n[0])) + G = sp.vstack((G1, G2), format="csr") + elif(self.dim == 3): + G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC[0])) + G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0])) + G3 = kron3(ddxCellGrad(n[2], BC[2]), speye(n[1]), speye(n[0])) + G = sp.vstack((G1, G2, G3), format="csr") + return G + def cellGrad(): doc = "The cell centered Gradient, takes you to cell faces." def fget(self): if(self._cellGrad is None): - BC = self.setCellGradBC(self._cellGradBC_list) - n = self.vnC - if(self.dim == 1): - G = ddxCellGrad(n[0], BC[0]) - elif(self.dim == 2): - G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0])) - G2 = sp.kron(ddxCellGrad(n[1], BC[1]), speye(n[0])) - G = sp.vstack((G1, G2), format="csr") - elif(self.dim == 3): - G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC[0])) - G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0])) - G3 = kron3(ddxCellGrad(n[2], BC[2]), speye(n[1]), speye(n[0])) - G = sp.vstack((G1, G2, G3), format="csr") + G = self._cellGradStencil() # Compute areas of cell faces & volumes S = self.area V = self.aveCC2F*self.vol # Average volume between adjacent cells @@ -361,19 +365,24 @@ class DiffOperators(object): _cellGradBC = None cellGradBC = property(**cellGradBC()) + def _cellGradxStencil(self): + BC = ['neumann', 'neumann'] + n = self.vnC + if(self.dim == 1): + G1 = ddxCellGrad(n[0], BC) + elif(self.dim == 2): + G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC)) + elif(self.dim == 3): + G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC)) + return G1 + + def cellGradx(): doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." def fget(self): if getattr(self, '_cellGradx', None) is None: - BC = ['neumann', 'neumann'] - n = self.vnC - if(self.dim == 1): - G1 = ddxCellGrad(n[0], BC) - elif(self.dim == 2): - G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC)) - elif(self.dim == 3): - G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC)) + G1 = self._cellGradxStencil() # Compute areas of cell faces & volumes V = self.aveCC2F*self.vol L = self.r(self.area/V, 'F','Fx', 'V') @@ -382,17 +391,22 @@ class DiffOperators(object): return locals() cellGradx = property(**cellGradx()) + def _cellGradyStencil(self): + if self.dim < 2: return None + BC = ['neumann', 'neumann'] + n = self.vnC + if(self.dim == 2): + G2 = sp.kron(ddxCellGrad(n[1], BC), speye(n[0])) + elif(self.dim == 3): + G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0])) + return G2 + def cellGrady(): doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." def fget(self): if self.dim < 2: return None if getattr(self, '_cellGrady', None) is None: - BC = ['neumann', 'neumann'] - n = self.vnC - if(self.dim == 2): - G2 = sp.kron(ddxCellGrad(n[1], BC), speye(n[0])) - elif(self.dim == 3): - G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0])) + G2 = self._cellGradyStencil() # Compute areas of cell faces & volumes V = self.aveCC2F*self.vol L = self.r(self.area/V, 'F','Fy', 'V') @@ -401,14 +415,19 @@ class DiffOperators(object): return locals() cellGrady = property(**cellGrady()) + def _cellGradzStencil(self): + if self.dim < 3: return None + BC = ['neumann', 'neumann'] + n = self.vnC + G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0])) + return G3 + def cellGradz(): doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." def fget(self): if self.dim < 3: return None if getattr(self, '_cellGradz', None) is None: - BC = ['neumann', 'neumann'] - n = self.vnC - G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0])) + G3 = self._cellGradzStencil() # Compute areas of cell faces & volumes V = self.aveCC2F*self.vol L = self.r(self.area/V, 'F','Fz', 'V') diff --git a/SimPEG/Mesh/MeshIO.py b/SimPEG/Mesh/MeshIO.py index 7501a66f..1c042237 100644 --- a/SimPEG/Mesh/MeshIO.py +++ b/SimPEG/Mesh/MeshIO.py @@ -21,10 +21,9 @@ class TensorMeshIO(object): if '*' in seg: st = seg sp = seg.split('*') - re = np.array(sp[0],dtype=int)*(' ' + sp[1]) + re = int(sp[0])*(' ' + sp[1]) line = line.replace(st,re.strip()) return np.array(line.split(),dtype=float) - # Read the file as line strings, remove lines with comment = ! msh = np.genfromtxt(fileName,delimiter='\n',dtype=np.str,comments='!') diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index 02c23dee..2c75fda1 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -2131,10 +2131,16 @@ class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO): def plotSlice(self, v, vType='CC', normal='Z', ind=None, grid=True, view='real', ax=None, clim=None, showIt=False, - pcolorOpts={}, - streamOpts={'color':'k'}, - gridOpts={'color':'k', 'alpha':0.5}): + pcolorOpts=None, + streamOpts=None, + gridOpts=None): + if pcolorOpts is None: + pcolorOpts = {} + if streamOpts is None: + streamOpts = {'color':'k'} + if gridOpts is None: + gridOpts = {'color':'k', 'alpha':0.5} assert vType in ['CC','F','E'] assert self.dim == 3 diff --git a/SimPEG/Mesh/View.py b/SimPEG/Mesh/View.py index 089d7d9a..8eb22098 100644 --- a/SimPEG/Mesh/View.py +++ b/SimPEG/Mesh/View.py @@ -42,9 +42,9 @@ class TensorView(object): def plotImage(self, v, vType='CC', grid=False, view='real', ax=None, clim=None, showIt=False, - pcolorOpts={}, - streamOpts={'color':'k'}, - gridOpts={'color':'k'}, + pcolorOpts=None, + streamOpts=None, + gridOpts=None, numbering=True, annotationColor='w' ): """ @@ -84,6 +84,12 @@ class TensorView(object): M.plotImage(v, annotationColor='k', showIt=True) """ + if pcolorOpts is None: + pcolorOpts = {} + if streamOpts is None: + streamOpts = {'color':'k'} + if gridOpts is None: + gridOpts = {'color':'k'} if ax is None: fig = plt.figure() @@ -174,9 +180,9 @@ class TensorView(object): def plotSlice(self, v, vType='CC', normal='Z', ind=None, grid=False, view='real', ax=None, clim=None, showIt=False, - pcolorOpts={}, - streamOpts={'color':'k'}, - gridOpts={'color':'k', 'alpha':0.5} + pcolorOpts=None, + streamOpts=None, + gridOpts=None ): """ @@ -197,6 +203,12 @@ class TensorView(object): M.plotSlice(M.cellGrad*b, 'F', view='vec', grid=True, showIt=True, pcolorOpts={'alpha':0.8}) """ + if pcolorOpts is None: + pcolorOpts = {} + if streamOpts is None: + streamOpts = {'color':'k'} + if gridOpts is None: + gridOpts = {'color':'k', 'alpha':0.5} if type(vType) in [list, tuple]: assert ax is None, "cannot specify an axis to plot on with this function." fig, axs = plt.subplots(1,len(vType)) @@ -289,11 +301,17 @@ class TensorView(object): def _plotImage2D(self, v, vType='CC', grid=False, view='real', ax=None, clim=None, showIt=False, - pcolorOpts={}, - streamOpts={'color':'k'}, - gridOpts={'color':'k'} + pcolorOpts=None, + streamOpts=None, + gridOpts=None ): + if pcolorOpts is None: + pcolorOpts = {} + if streamOpts is None: + streamOpts = {'color':'k'} + if gridOpts is None: + gridOpts = {'color':'k'} vTypeOptsCC = ['N','CC','Fx','Fy','Ex','Ey'] vTypeOptsV = ['CCv','F','E'] vTypeOpts = vTypeOptsCC + vTypeOptsV diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py index 4f2cb062..77704733 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -888,6 +888,8 @@ class ProjectedGNCG(BFGS, Minimize, Remember): maxIterCG = 5 tolCG = 1e-1 + stepOffBoundsFact = 0.1 # perturbation of the inactive set off the bounds + lower = -np.inf upper = np.inf @@ -990,4 +992,20 @@ class ProjectedGNCG(BFGS, Minimize, Remember): cgFlag = 1 # End CG Iterations + # Take a gradient step on the active cells if exist + if temp != self.xc.size: + + rhs_a = (Active) * -self.g + + dm_i = max( abs( delx ) ) + dm_a = max( abs(rhs_a) ) + + # perturb inactive set off of bounds so that they are included in the step + delx = delx + self.stepOffBoundsFact * (rhs_a * dm_i / dm_a) + + + # Only keep gradients going in the right direction on the active set + indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0)) + delx[indx] = 0. + return delx diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 5b41b91b..fc101a61 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -1,5 +1,289 @@ import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp +class RegularizationMesh(object): + """ + **Regularization Mesh** + + This contains the operators used in the regularization. Note that these + are not necessarily true differential operators, but are constructed from + a SimPEG Mesh. + + :param Mesh mesh: problem mesh + :param numpy.array indActive: bool array, size nC, that is True where we have active cells. Used to reduce the operators so we regularize only on active cells + """ + + def __init__(self, mesh, indActive=None): + self.mesh = mesh + assert indActive is None or indActive.dtype == 'bool', 'indActive needs to be None or a bool' + self.indActive = indActive + + @property + def vol(self): + """ + reduced volume vector + :rtype: numpy.array + :return: reduced cell volume + """ + if getattr(self, '_vol', None) is None: + self._vol = self._Pac.T * self.mesh.vol + return self._vol + + @property + def nC(self): + """ + reduced number of cells + :rtype: int + :return: number of cells being regularized + """ + if getattr(self, '_nC', None) is None: + if self.indActive is None: + self._nC = self.mesh.nC + else: + self._nC = sum(self.indActive) + return self._nC + + @property + def dim(self): + """ + dimension of regularization mesh (1D, 2D, 3D) + :rtype: int + :return: dimension + """ + if getattr(self, '_dim', None) is None: + self._dim = self.mesh.dim + return self._dim + + + @property + def _Pac(self): + """ + projection matrix that takes from the reduced space of active cells to full modelling space (ie. nC x nindActive) + :rtype: scipy.sparse.csr_matrix + :return: active cell projection matrix + """ + if getattr(self, '__Pac', None) is None: + if self.indActive is None: + self.__Pac = Utils.speye(self.mesh.nC) + else: + self.__Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + return self.__Pac + + @property + def _Pafx(self): + """ + projection matrix that takes from the reduced space of active x-faces to full modelling space (ie. nFx x nindActive_Fx ) + :rtype: scipy.sparse.csr_matrix + :return: active face-x projection matrix + """ + if getattr(self, '__Pafx', None) is None: + if self.indActive is None: + self.__Pafx = Utils.speye(self.mesh.nFx) + else: + indActive_Fx = (self.mesh.aveFx2CC.T * self.indActive) == 1 + self.__Pafx = Utils.speye(self.mesh.nFx)[:,indActive_Fx] + return self.__Pafx + + @property + def _Pafy(self): + """ + projection matrix that takes from the reduced space of active y-faces to full modelling space (ie. nFy x nindActive_Fy ) + :rtype: scipy.sparse.csr_matrix + :return: active face-y projection matrix + """ + if getattr(self, '__Pafy', None) is None: + if self.indActive is None: + self.__Pafy = Utils.speye(self.mesh.nFy) + else: + indActive_Fy = (self.mesh.aveFy2CC.T * self.indActive) == 1 + self.__Pafy = Utils.speye(self.mesh.nFy)[:,indActive_Fy] + return self.__Pafy + + @property + def _Pafz(self): + """ + projection matrix that takes from the reduced space of active z-faces to full modelling space (ie. nFz x nindActive_Fz ) + :rtype: scipy.sparse.csr_matrix + :return: active face-z projection matrix + """ + if getattr(self, '__Pafz', None) is None: + if self.indActive is None: + self.__Pafz = Utils.speye(self.mesh.nFz) + else: + indActive_Fz = (self.mesh.aveFz2CC.T * self.indActive) == 1 + self.__Pafz = Utils.speye(self.mesh.nFz)[:,indActive_Fz] + return self.__Pafz + + @property + def aveFx2CC(self): + """ + averaging from active cell centers to active x-faces + :rtype: scipy.sparse.csr_matrix + :return: averaging from active cell centers to active x-faces + """ + if getattr(self, '_aveFx2CC', None) is None: + self._aveFx2CC = self._Pac.T * self.mesh.aveFx2CC * self._Pafx + return self._aveFx2CC + + @property + def aveCC2Fx(self): + """ + averaging from active x-faces to active cell centers + :rtype: scipy.sparse.csr_matrix + :return: averaging matrix from active x-faces to active cell centers + """ + if getattr(self, '_aveCC2Fx', None) is None: + self._aveCC2Fx = Utils.sdiag(1./(self.aveFx2CC.T).sum(1)) * self.aveFx2CC.T + return self._aveCC2Fx + + @property + def aveFy2CC(self): + """ + averaging from active cell centers to active y-faces + :rtype: scipy.sparse.csr_matrix + :return: averaging from active cell centers to active y-faces + """ + if getattr(self, '_aveFy2CC', None) is None: + self._aveFy2CC = self._Pac.T * self.mesh.aveFy2CC * self._Pafy + return self._aveFy2CC + + @property + def aveCC2Fy(self): + """ + averaging from active y-faces to active cell centers + :rtype: scipy.sparse.csr_matrix + :return: averaging matrix from active y-faces to active cell centers + """ + if getattr(self, '_aveCC2Fy', None) is None: + self._aveCC2Fy = Utils.sdiag(1./(self.aveFy2CC.T).sum(1)) * self.aveFy2CC.T + return self._aveCC2Fy + + @property + def aveFz2CC(self): + """ + averaging from active cell centers to active z-faces + :rtype: scipy.sparse.csr_matrix + :return: averaging from active cell centers to active z-faces + """ + if getattr(self, '_aveFz2CC', None) is None: + self._aveFz2CC = self._Pac.T * self.mesh.aveFz2CC * self._Pafz + return self._aveFz2CC + + @property + def aveCC2Fz(self): + """ + averaging from active z-faces to active cell centers + :rtype: scipy.sparse.csr_matrix + :return: averaging matrix from active z-faces to active cell centers + """ + if getattr(self, '_aveCC2Fz', None) is None: + self._aveCC2Fz = Utils.sdiag(1./(self.aveFz2CC.T).sum(1)) * self.aveFz2CC.T + return self._aveCC2Fz + + @property + def cellDiffx(self): + """ + cell centered difference in the x-direction + :rtype: scipy.sparse.csr_matrix + :return: differencing matrix for active cells in the x-direction + """ + if getattr(self, '_cellDiffx', None) is None: + self._cellDiffx = self._Pafx.T * self.mesh.cellGradx * self._Pac + return self._cellDiffx + + @property + def cellDiffy(self): + """ + cell centered difference in the y-direction + :rtype: scipy.sparse.csr_matrix + :return: differencing matrix for active cells in the y-direction + """ + if getattr(self, '_cellDiffy', None) is None: + self._cellDiffy = self._Pafy.T * self.mesh.cellGrady * self._Pac + return self._cellDiffy + + @property + def cellDiffz(self): + """ + cell centered difference in the z-direction + :rtype: scipy.sparse.csr_matrix + :return: differencing matrix for active cells in the z-direction + """ + if getattr(self, '_cellDiffz', None) is None: + self._cellDiffz = self._Pafz.T * self.mesh.cellGradz * self._Pac + return self._cellDiffz + + @property + def faceDiffx(self): + """ + x-face differences + :rtype: scipy.sparse.csr_matrix + :return: differencing matrix for active faces in the x-direction + """ + if getattr(self, '_faceDiffx', None) is None: + self._faceDiffx = self._Pac.T * self.mesh.faceDivx * self._Pafx + return self._faceDiffx + + @property + def faceDiffy(self): + """ + y-face differences + :rtype: scipy.sparse.csr_matrix + :return: differencing matrix for active faces in the y-direction + """ + if getattr(self, '_faceDiffy', None) is None: + self._faceDiffy = self._Pac.T * self.mesh.faceDivy * self._Pafy + return self._faceDiffy + + @property + def faceDiffz(self): + """ + z-face differences + :rtype: scipy.sparse.csr_matrix + :return: differencing matrix for active faces in the z-direction + """ + if getattr(self, '_faceDiffz', None) is None: + self._faceDiffz = self._Pac.T * self.mesh.faceDivz * self._Pafz + return self._faceDiffz + + @property + def cellDiffxStencil(self): + """ + cell centered difference stencil (no cell lengths include) in the x-direction + :rtype: scipy.sparse.csr_matrix + :return: differencing matrix for active cells in the x-direction + """ + if getattr(self, '_cellDiffxStencil', None) is None: + + self._cellDiffxStencil = self._Pafx.T * self.mesh._cellGradxStencil() * self._Pac + return self._cellDiffxStencil + + @property + def cellDiffyStencil(self): + """ + cell centered difference stencil (no cell lengths include) in the y-direction + :rtype: scipy.sparse.csr_matrix + :return: differencing matrix for active cells in the y-direction + """ + if self.dim < 2: return None + if getattr(self, '_cellDiffyStencil', None) is None: + + self._cellDiffyStencil = self._Pafy.T * self.mesh._cellGradyStencil() * self._Pac + return self._cellDiffyStencil + + @property + def cellDiffzStencil(self): + """ + cell centered difference stencil (no cell lengths include) in the y-direction + :rtype: scipy.sparse.csr_matrix + :return: differencing matrix for active cells in the y-direction + """ + if self.dim < 3: return None + if getattr(self, '_cellDiffzStencil', None) is None: + + self._cellDiffzStencil = self._Pafz.T * self.mesh._cellGradzStencil() * self._Pac + return self._cellDiffzStencil + + class BaseRegularization(object): """ **Base Regularization Class** @@ -18,12 +302,19 @@ class BaseRegularization(object): mapping = None #: A SimPEG.Map instance. mesh = None #: A SimPEG.Mesh instance. - mref = None #: Reference model. + mref = None #: Reference model. def __init__(self, mesh, mapping=None, indActive=None, **kwargs): Utils.setKwargs(self, **kwargs) - self.mesh = mesh assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object." + if indActive is not None and indActive.dtype != 'bool': + tmp = indActive + indActive = np.zeros(mesh.nC, dtype=bool) + indActive[tmp] = True + if indActive is not None and mapping is None: + mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size) + + self.regmesh = RegularizationMesh(mesh,indActive) self.mapping = mapping or self.mapPair(mesh) self.mapping._assertMatchesPair(self.mapPair) self.indActive = indActive @@ -55,8 +346,7 @@ class BaseRegularization(object): @property def W(self): """Full regularization weighting matrix W.""" - return sp.identity(self.mapping.nP) - + return sp.identity(self.regmesh.nC) @Utils.timeIt def eval(self, m): @@ -87,11 +377,12 @@ class BaseRegularization(object): @Utils.timeIt def eval2Deriv(self, m, v=None): """ + Second derivative - :param numpy.array m: geophysical model - :param numpy.array v: vector to multiply - :rtype: scipy.sparse.csr_matrix or numpy.ndarray - :return: WtW or WtW*v + :param numpy.array m: geophysical model + :param numpy.array v: vector to multiply + :rtype: scipy.sparse.csr_matrix or numpy.ndarray + :return: WtW or WtW*v The regularization is: @@ -112,112 +403,94 @@ class BaseRegularization(object): return mD.T * ( self.W.T * ( self.W * ( mD * v) ) ) - class Tikhonov(BaseRegularization): """ + L2 Tikhonov regularization with both smallness and smoothness (first order + derivative) contributions. + + .. math:: + \phi_m(\mathbf{m}) = \\alpha_s \| W_s (\mathbf{m} - \mathbf{m_{ref}} ) \|^2 + + \\alpha_x \| W_x \\frac{\partial}{\partial x} (\mathbf{m} - \mathbf{m_{ref}} ) \|^2 + + \\alpha_y \| W_y \\frac{\partial}{\partial y} (\mathbf{m} - \mathbf{m_{ref}} ) \|^2 + + \\alpha_z \| W_z \\frac{\partial}{\partial z} (\mathbf{m} - \mathbf{m_{ref}} ) \|^2 + + Note if the key word argument `mrefInSmooth` is False, then mref is not + included in the smoothness contribution. + + :param Mesh mesh: SimPEG mesh + :param Maps mapping: regularization mapping, takes the model from model space to the thing you want to regularize + :param numpy.ndarray indActive: active cell indices for reducing the size of differential operators in the definition of a regularization mesh + :param bool mrefInSmooth: (default = False) put mref in the smoothness component? + :param float alpha_s: (default 1e-6) smallness weight + :param float alpha_x: (default 1) smoothness weight for first derivative in the x-direction + :param float alpha_y: (default 1) smoothness weight for first derivative in the y-direction + :param float alpha_z: (default 1) smoothness weight for first derivative in the z-direction + :param float alpha_xx: (default 1) smoothness weight for second derivative in the x-direction + :param float alpha_yy: (default 1) smoothness weight for second derivative in the y-direction + :param float alpha_zz: (default 1) smoothness weight for second derivative in the z-direction """ - smoothModel = True #: SMOOTH and SMOOTH_MOD_DIF options - alpha_s = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Ws'], "Smallness weight") - alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction") - alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction") - alpha_z = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction") - alpha_xx = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction") - alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction") - alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction") + mrefInSmooth = False # put mref in the smoothness contribution + alpha_s = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Wsmall'], "Smallness weight") + alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction") + alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction") + alpha_z = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction") + alpha_xx = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction") + alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction") + alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction") def __init__(self, mesh, mapping=None, indActive = None, **kwargs): - BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs) - self.indActive = indActive + BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) @property - def Ws(self): - """Regularization matrix Ws""" - if getattr(self,'_Ws', None) is None: - self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5) - if self.indActive is not None: - Pac = Utils.speye(self.mesh.nC)[:,self.indActive] - self._Ws = Pac.T * self._Ws * Pac - return self._Ws + def Wsmall(self): + """Regularization matrix Wsmall""" + if getattr(self,'_Wsmall', None) is None: + self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s)**0.5) + return self._Wsmall @property def Wx(self): """Regularization matrix Wx""" if getattr(self, '_Wx', None) is None: - Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFx].T*self.mesh.vol - self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx - - if self.indActive is not None: - indActive_Fx = (self.mesh.aveFx2CC.T * self.indActive) == 1 - Pac = Utils.speye(self.mesh.nC)[:,self.indActive] - Pafx = Utils.speye(self.mesh.nFx)[:,indActive_Fx] - self._Wx = Pafx.T*self._Wx*Pac - + Ave_x_vol = self.regmesh.aveCC2Fx * self.regmesh.vol + self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.regmesh.cellDiffx return self._Wx @property def Wy(self): """Regularization matrix Wy""" if getattr(self, '_Wy', None) is None: - Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFx:np.sum(self.mesh.vnF[:2])].T*self.mesh.vol - self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady - - if self.indActive is not None: - indActive_Fy = (self.mesh.aveFy2CC.T * self.indActive) == 1 - Pac = Utils.speye(self.mesh.nC)[:,self.indActive] - Pafy = Utils.speye(self.mesh.nFy)[:,indActive_Fy] - self._Wy = Pafy.T*self._Wy*Pac - + Ave_y_vol = self.regmesh.aveCC2Fy * self.regmesh.vol + self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.regmesh.cellDiffy return self._Wy @property def Wz(self): """Regularization matrix Wz""" if getattr(self, '_Wz', None) is None: - Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.vnF[:2]):].T*self.mesh.vol - self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz - - if self.indActive is not None: - indActive_Fz = (self.mesh.aveFz2CC.T * self.indActive) == 1 - Pac = Utils.speye(self.mesh.nC)[:,self.indActive] - Pafz = Utils.speye(self.mesh.nFz)[:,indActive_Fz] - self._Wz = Pafz.T*self._Wz*Pac - + Ave_z_vol = self.regmesh.aveCC2Fz * self.regmesh.vol + self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.regmesh.cellDiffz return self._Wz @property def Wxx(self): """Regularization matrix Wxx""" if getattr(self, '_Wxx', None) is None: - self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx - - if self.indActive is not None: - Pac = Utils.speye(self.mesh.nC)[:,self.indActive] - self._Wxx = Pac.T*self._Wxx*Pac - + self._Wxx = Utils.sdiag((self.regmesh.vol*self.alpha_xx)**0.5)*self.regmesh.faceDiffx*self.regmesh.cellDiffx return self._Wxx @property def Wyy(self): """Regularization matrix Wyy""" if getattr(self, '_Wyy', None) is None: - self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady - - if self.indActive is not None: - Pac = Utils.speye(self.mesh.nC)[:,self.indActive] - self._Wyy = Pac.T*self._Wyy*Pac - + self._Wyy = Utils.sdiag((self.regmesh.vol*self.alpha_yy)**0.5)*self.regmesh.faceDiffy*self.regmesh.cellDiffy return self._Wyy @property def Wzz(self): """Regularization matrix Wzz""" if getattr(self, '_Wzz', None) is None: - self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz - - if self.indActive is not None: - Pac = Utils.speye(self.mesh.nC)[:,self.indActive] - self._Wzz = Pac.T*self._Wzz*Pac - + self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDiffz*self.regmesh.cellDiffz return self._Wzz @property @@ -225,9 +498,9 @@ class Tikhonov(BaseRegularization): """Full smoothness regularization matrix W""" if getattr(self, '_Wsmooth', None) is None: wlist = (self.Wx, self.Wxx) - if self.mesh.dim > 1: + if self.regmesh.dim > 1: wlist += (self.Wy, self.Wyy) - if self.mesh.dim > 2: + if self.regmesh.dim > 2: wlist += (self.Wz, self.Wzz) self._Wsmooth = sp.vstack(wlist) return self._Wsmooth @@ -236,25 +509,44 @@ class Tikhonov(BaseRegularization): def W(self): """Full regularization matrix W""" if getattr(self, '_W', None) is None: - wlist = (self.Ws, self.Wsmooth) + wlist = (self.Wsmall, self.Wsmooth) self._W = sp.vstack(wlist) return self._W @Utils.timeIt - def eval(self, m): - if self.smoothModel == True: - r1 = self.Wsmooth * ( self.mapping * (m) ) - r2 = self.Ws * ( self.mapping * (m - self.mref) ) - return 0.5*(r1.dot(r1)+r2.dot(r2)) - elif self.smoothModel == False: - r = self.W * ( self.mapping * (m - self.mref) ) - return 0.5*r.dot(r) + def _evalSmall(self, m): + r = self.Wsmall * ( self.mapping * (m - self.mref) ) + return 0.5 * r.dot(r) + @Utils.timeIt + def _evalSmooth(self, m): + if self.mrefInSmooth == True: + r = self.Wsmooth * ( self.mapping * (m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wsmooth * ( self.mapping * (m) ) + return 0.5 * r.dot(r) + + @Utils.timeIt + def eval(self, m): + return self._evalSmall(m) + self._evalSmooth(m) + + @Utils.timeIt + def _evalSmallDeriv(self,m): + r = self.Wsmall * ( self.mapping * (m - self.mref) ) + return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) ) + + @Utils.timeIt + def _evalSmoothDeriv(self,m): + if self.mrefInSmooth == True: + r = self.Wsmooth * ( self.mapping * ( m - self.mref ) ) + return r.T * ( self.Wsmooth * self.mapping.deriv(m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wsmooth * ( self.mapping * m ) + return r.T * ( self.Wsmooth * self.mapping.deriv(m) ) @Utils.timeIt def evalDeriv(self, m): """ - The regularization is: .. math:: @@ -268,257 +560,185 @@ class Tikhonov(BaseRegularization): R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} """ - if self.smoothModel == True: - mD1 = self.mapping.deriv(m) - mD2 = self.mapping.deriv(m - self.mref) - r1 = self.Wsmooth * ( self.mapping * (m)) - r2 = self.Ws * ( self.mapping * (m - self.mref) ) - out1 = mD1.T * ( self.Wsmooth.T * r1 ) - out2 = mD2.T * ( self.Ws.T * r2 ) - out = out1+out2 - elif self.smoothModel == False: - mD = self.mapping.deriv(m - self.mref) - r = self.W * ( self.mapping * (m - self.mref) ) - out = mD.T * ( self.W.T * r ) - return out - -# <<<<<<< HEAD - -# class Simple(BaseRegularization): -# """ -# Only for tensor mesh -# """ - -# smoothModel = True #: SMOOTH and SMOOTH_MOD_DIF options -# alpha_s = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Ws'], "Smallness weight") -# alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction") -# alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction") -# alpha_z = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction") -# alpha_xx = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction") -# alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction") -# alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction") - -# def __init__(self, mesh, mapping=None, **kwargs): -# BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs) + return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) +class Simple(Tikhonov): + """ + Simple regularization that does not include length scales in the derivatives. + """ -# @property -# def Ws(self): -# """Regularization matrix Ws""" -# if getattr(self,'_Ws', None) is None: -# self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5) -# return self._Ws + mrefInSmooth = False #: SMOOTH and SMOOTH_MOD_DIF options + alpha_s = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight") + alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction") + alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction") + alpha_z = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction") + wght = 1. -# @property -# def Wx(self): -# """Regularization matrix Wx""" -# if getattr(self, '_Wx', None) is None: -# self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x)**0.5)*self.mesh.unitCellGradx -# return self._Wx + def __init__(self, mesh, mapping=None, indActive=None, **kwargs): + BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) -# @property -# def Wy(self): -# """Regularization matrix Wy""" -# if getattr(self, '_Wy', None) is None: -# self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y)**0.5)*self.mesh.unitCellGrady -# return self._Wy + if isinstance(self.wght,float): + self.wght = np.ones(self.regmesh.nC) * self.wght -# @property -# def Wz(self): -# """Regularization matrix Wz""" -# if getattr(self, '_Wz', None) is None: -# self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z)**0.5)*self.mesh.unitCellGradz -# return self._Wz + @property + def Wsmall(self): + """Regularization matrix Wsmall""" + if getattr(self,'_Wsmall', None) is None: + self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5) + return self._Wsmall -# @property -# def Wxx(self): -# """Regularization matrix Wxx""" -# if getattr(self, '_Wxx', None) is None: -# self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx -# return self._Wxx + @property + def Wx(self): + """Regularization matrix Wx""" + if getattr(self, '_Wx', None) is None: + self._Wx = Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.regmesh.cellDiffxStencil + return self._Wx -# @property -# def Wyy(self): -# """Regularization matrix Wyy""" -# if getattr(self, '_Wyy', None) is None: -# self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady -# return self._Wyy + @property + def Wy(self): + """Regularization matrix Wy""" + if getattr(self, '_Wy', None) is None: + self._Wy = Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol * self.alpha_y*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.regmesh.cellDiffyStencil + return self._Wy -# @property -# def Wzz(self): -# """Regularization matrix Wzz""" -# if getattr(self, '_Wzz', None) is None: -# self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz -# return self._Wzz + @property + def Wz(self): + """Regularization matrix Wz""" + if getattr(self, '_Wz', None) is None: + self._Wz = Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.regmesh.cellDiffzStencil + return self._Wz -# @property -# def Wsmooth(self): -# """Full smoothness regularization matrix W""" -# if getattr(self, '_Wsmooth', None) is None: -# wlist = (self.Wx, self.Wxx) -# if self.mesh.dim > 1: -# wlist += (self.Wy, self.Wyy) -# if self.mesh.dim > 2: -# wlist += (self.Wz, self.Wzz) -# self._Wsmooth = sp.vstack(wlist) -# return self._Wsmooth + @property + def Wsmooth(self): + """Full smoothness regularization matrix W""" + if getattr(self, '_Wsmooth', None) is None: + wlist = (self.Wx,) + if self.regmesh.dim > 1: + wlist += (self.Wy,) + if self.regmesh.dim > 2: + wlist += (self.Wz,) + self._Wsmooth = sp.vstack(wlist) + return self._Wsmooth -# @property -# def W(self): -# """Full regularization matrix W""" -# if getattr(self, '_W', None) is None: -# wlist = (self.Ws, self.Wsmooth) -# self._W = sp.vstack(wlist) -# return self._W + @property + def W(self): + """Full regularization matrix W""" + if getattr(self, '_W', None) is None: + wlist = (self.Wsmall, self.Wsmooth) + self._W = sp.vstack(wlist) + return self._W -# @Utils.timeIt -# def eval(self, m): -# if self.smoothModel == True: -# r1 = self.Wsmooth * ( self.mapping * (m) ) -# r2 = self.Ws * ( self.mapping * (m - self.mref) ) -# return 0.5*(r1.dot(r1)+r2.dot(r2)) -# elif self.smoothModel == False: -# r = self.W * ( self.mapping * (m - self.mref) ) -# return 0.5*r.dot(r) + @Utils.timeIt + def _evalSmall(self, m): + r = self.Wsmall * ( self.mapping * (m - self.mref) ) + return 0.5 * r.dot(r) + + @Utils.timeIt + def _evalSmooth(self, m): + if self.mrefInSmooth == True: + r = self.Wsmooth * ( self.mapping * (m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wsmooth * ( self.mapping * m) + return 0.5 * r.dot(r) -# @Utils.timeIt -# def evalDeriv(self, m): -# """ +class Sparse(Simple): -# The regularization is: + # set default values + eps_p = 1e-1 + eps_q = 1e-1 + curModel = None # use a model to compute the weights + gamma = 1. + norms = [0., 2., 2., 2.] + wght = 1. -# .. math:: + def __init__(self, mesh, mapping=None, indActive=None, **kwargs): + Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) -# R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})} + if isinstance(self.wght,float): + self.wght = np.ones(self.regmesh.nC) * self.wght -# So the derivative is straight forward: + @property + def Wsmall(self): + """Regularization matrix Wsmall""" + if getattr(self, 'curModel', None) is None: + self.Rs = Utils.speye(self.regmesh.nC) -# .. math:: + else: + f_m = self.curModel - self.reg.mref + self.rs = self.R(f_m , self.eps_p, self.norms[0]) + #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs)) + self.Rs = Utils.sdiag( self.rs ) -# R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} - -# """ -# if self.smoothModel == True: -# mD1 = self.mapping.deriv(m) -# mD2 = self.mapping.deriv(m - self.mref) -# r1 = self.Wsmooth * ( self.mapping * (m)) -# r2 = self.Ws * ( self.mapping * (m - self.mref) ) -# out1 = mD1.T * ( self.Wsmooth.T * r1 ) -# out2 = mD2.T * ( self.Ws.T * r2 ) -# out = out1+out2 -# elif self.smoothModel == False: -# mD = self.mapping.deriv(m - self.mref) -# r = self.W * ( self.mapping * (m - self.mref) ) -# out = mD.T * ( self.W.T * r ) -# return out - -# class SparseRegularization(Simple): - -# eps = 1e-1 - -# m = None -# gamma = 1. -# p = 0. -# qx = 2. -# qy = 2. -# qz = 2. - -# def __init__(self, mesh, mapping=None, **kwargs): -# Simple.__init__(self, mesh, mapping=mapping, **kwargs) + return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.wght)**0.5)*self.Rs -# @property -# def Wsmooth(self): -# """Full smoothness regularization matrix W""" -# if getattr(self, '_Wsmooth', None) is None: -# wlist = (self.Wx, self.Wxx) -# if self.mesh.dim > 1: -# wlist += (self.Wy, self.Wyy) -# if self.mesh.dim > 2: -# wlist += (self.Wz, self.Wzz) -# self._Wsmooth = sp.vstack(wlist) -# return self._Wsmooth + @property + def Wx(self): + """Regularization matrix Wx""" -# @property -# def W(self): -# """Full regularization matrix W""" -# if getattr(self, '_W', None) is None: -# wlist = (self.Ws, self.Wsmooth) -# self._W = sp.vstack(wlist) -# return self._W + if getattr(self, 'curModel', None) is None: + self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0]) -# @property -# def Ws(self): -# """Regularization matrix Ws""" -# if getattr(self, 'm', None) is None: -# self.Rs = Utils.speye(self.mesh.nC) + else: + f_m = self.regmesh.cellDiffxStencil * self.curModel + self.rx = self.R( f_m , self.eps_q, self.norms[1]) + self.Rx = Utils.sdiag( self.rx ) -# else: -# f_m = self.m -# self.rs = self.R(f_m , self.p, self.eps) -# #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs)) -# self.Rs = Utils.sdiag( self.rs ) + return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.Rx*self.regmesh.cellDiffxStencil -# self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s*self.gamma)**0.5)*self.Rs + @property + def Wy(self): + """Regularization matrix Wy""" -# return self._Ws + if getattr(self, 'curModel', None) is None: + self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0]) -# @property -# def Wx(self): -# """Regularization matrix Wx""" + else: + f_m = self.regmesh.cellDiffyStencil * self.curModel + self.ry = self.R( f_m , self.eps_q, self.norms[2]) + self.Ry = Utils.sdiag( self.ry ) -# if getattr(self, 'm', None) is None: -# self.Rx = Utils.speye(self.mesh.unitCellGradx.shape[0]) + return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.Ry*self.regmesh.cellDiffyStencil -# else: -# f_m = self.mesh.unitCellGradx * self.m -# self.rx = self.R( f_m , self.qx, self.eps) -# self.Rx = Utils.sdiag( self.rx ) + @property + def Wz(self): + """Regularization matrix Wz""" -# if getattr(self, '_Wx', None) is None: -# self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x*self.gamma)**0.5)*self.Rx*self.mesh.unitCellGradx -# return self._Wx + if getattr(self, 'curModel', None) is None: + self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0]) -# @property -# def Wy(self): -# """Regularization matrix Wy""" + else: + f_m = self.regmesh.cellDiffzStencil * self.curModel + self.rz = self.R( f_m , self.eps_q, self.norms[3]) + self.Rz = Utils.sdiag( self.rz ) -# if getattr(self, 'm', None) is None: -# self.Ry = Utils.speye(self.mesh.unitCellGrady.shape[0]) + return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.Rz*self.regmesh.cellDiffzStencil -# else: -# f_m = self.mesh.unitCellGrady * self.m -# self.ry = self.R( f_m , self.qy, self.eps) -# self.Ry = Utils.sdiag( self.ry ) + @property + def Wsmooth(self): + """Full smoothness regularization matrix W""" + #if getattr(self, '_Wsmooth', None) is None: + wlist = (self.Wx,) + if self.regmesh.dim > 1: + wlist += (self.Wy,) + if self.regmesh.dim > 2: + wlist += (self.Wz,) + #self._Wsmooth = sp.vstack(wlist) + return sp.vstack(wlist) -# if getattr(self, '_Wy', None) is None: -# self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y*self.gamma)**0.5)*self.Ry*self.mesh.unitCellGrady -# return self._Wy + @property + def W(self): + """Full regularization matrix W""" + if getattr(self, '_W', None) is None: + wlist = (self.Wsmall, self.Wsmooth) + self._W = sp.vstack(wlist) + return self._W -# @property -# def Wz(self): -# """Regularization matrix Wz""" + def R(self, f_m , eps, exponent): -# if getattr(self, 'm', None) is None: -# self.Rz = Utils.speye(self.mesh.unitCellGradz.shape[0]) + eta = (eps**(1.-exponent/2.))**0.5 + r = eta / (f_m**2.+ eps**2.)**((1.-exponent/2.)/2.) -# else: -# f_m = self.mesh.unitCellGradz * self.m -# self.rz = self.R( f_m , self.qz, self.eps) -# self.Rz = Utils.sdiag( self.rz ) - -# if getattr(self, '_Wz', None) is None: -# self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z*self.gamma)**0.5)*self.Rz*self.mesh.unitCellGradz -# return self._Wz - - -# def R(self, f_m , p, dec): - -# eta = (self.eps**(1-p/2.))**0.5 -# r = eta / (f_m**2.+self.eps**2.)**((1-p/2.)/2.) - -# return r -# ======= -# >>>>>>> 834de582844e8e1eac95819fbe03eed55dbeb001 + return r diff --git a/SimPEG/Utils/ModelBuilder.py b/SimPEG/Utils/ModelBuilder.py index 435856f7..1b868fe5 100644 --- a/SimPEG/Utils/ModelBuilder.py +++ b/SimPEG/Utils/ModelBuilder.py @@ -88,12 +88,14 @@ def getIndicesBlock(p0,p1,ccMesh): # Return a tuple return ind -def defineBlock(ccMesh,p0,p1,vals=[0,1]): +def defineBlock(ccMesh,p0,p1,vals=None): """ Build a block with the conductivity specified by condVal. Returns an array. vals[0] conductivity of the block vals[1] conductivity of the ground """ + if vals is None: + vals = [0,1] sigma = np.zeros(ccMesh.shape[0]) + vals[1] ind = getIndicesBlock(p0,p1,ccMesh) @@ -101,7 +103,11 @@ def defineBlock(ccMesh,p0,p1,vals=[0,1]): return mkvc(sigma) -def defineElipse(ccMesh, center=[0,0,0], anisotropy=[1,1,1], slope=10., theta=0.): +def defineElipse(ccMesh, center=None, anisotropy=None, slope=10., theta=0.): + if center is None: + center = [0,0,0] + if anisotropy is None: + anisotropy = [1,1,1] G = ccMesh.copy() dim = ccMesh.shape[1] for i in range(dim): @@ -156,7 +162,7 @@ def getIndicesSphere(center,radius,ccMesh): # Return a tuple return ind -def defineTwoLayers(ccMesh,depth,vals=[0,1]): +def defineTwoLayers(ccMesh,depth,vals=None): """ Define a two layered model. Depth of the first layer must be specified. CondVals vector with the conductivity values of the layers. Eg: @@ -167,6 +173,8 @@ def defineTwoLayers(ccMesh,depth,vals=[0,1]): 0 depth zf 1st layer 2nd layer """ + if vals is None: + vals = [0,1] sigma = np.zeros(ccMesh.shape[0]) + vals[1] dim = np.size(ccMesh[0,:]) @@ -252,7 +260,7 @@ def layeredModel(ccMesh, layerTops, layerValues): -def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]): +def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=None): """ Create a random model by convolving a kernel with a uniformly distributed model. @@ -276,6 +284,8 @@ def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]): """ + if bounds is None: + bounds = [0,1] if seed is None: seed = np.random.randint(1e3) diff --git a/SimPEG/Utils/codeutils.py b/SimPEG/Utils/codeutils.py index bfd00889..6151c540 100644 --- a/SimPEG/Utils/codeutils.py +++ b/SimPEG/Utils/codeutils.py @@ -55,8 +55,10 @@ def hook(obj, method, name=None, overwrite=False, silent=False): print 'Method '+name+' was not overwritten.' -def setKwargs(obj, ignore=[], **kwargs): +def setKwargs(obj, ignore=None, **kwargs): """Sets key word arguments (kwargs) that are present in the object, throw an error if they don't exist.""" + if ignore is None: + ignore = [] for attr in kwargs: if attr in ignore: continue diff --git a/SimPEG/Utils/io_utils.py b/SimPEG/Utils/io_utils.py new file mode 100644 index 00000000..323feb0c --- /dev/null +++ b/SimPEG/Utils/io_utils.py @@ -0,0 +1,137 @@ +from SimPEG import np, Mesh +import time as tm +import vtk, vtk.util.numpy_support as npsup +import re + +def read_GOCAD_ts(tsfile): + """ + + Read GOCAD triangulated surface (*.ts) file + INPUT: + tsfile: Triangulated surface + + OUTPUT: + vrts : Array of vertices in XYZ coordinates [n x 3] + trgl : Array of index for triangles [m x 3]. The order of the vertices + is important and describes the normal + n = cross( (P2 - P1 ) , (P3 - P1) ) + + Author: @fourndo + + + .. note:: + + Remove all attributes from the GoCAD surface before exporting it! + + """ + + + fid = open(tsfile,'r') + line = fid.readline() + + # Skip all the lines until the vertices + while re.match('TFACE',line)==None: + line = fid.readline() + + line = fid.readline() + vrtx = [] + + # Run down all the vertices and save in array + while re.match('VRTX',line): + l_input = re.split('[\s*]',line) + temp = np.array(l_input[2:5]) + vrtx.append(temp.astype(np.float)) + + # Read next line + line = fid.readline() + + vrtx = np.asarray(vrtx) + + # Skip lines to the triangles + while re.match('TRGL',line)==None: + line = fid.readline() + + # Run down the list of triangles + trgl = [] + + # Run down all the vertices and save in array + while re.match('TRGL',line): + l_input = re.split('[\s*]',line) + temp = np.array(l_input[1:4]) + trgl.append(temp.astype(np.int)) + + # Read next line + line = fid.readline() + + trgl = np.asarray(trgl) + + return vrtx, trgl + +def surface2inds(vrtx, trgl, mesh, boundaries=True, internal=True): + """" + Function to read gocad polystructure file and output indexes of mesh with in the structure. + + """ + # Adjust the index + trgl = trgl - 1 + + # Make vtk pts + ptsvtk = vtk.vtkPoints() + ptsvtk.SetData(npsup.numpy_to_vtk(vrtx,deep=1)) + + # Make the polygon connection + polys = vtk.vtkCellArray() + for face in trgl: + poly = vtk.vtkPolygon() + poly.GetPointIds().SetNumberOfIds(len(face)) + for nrv, vert in enumerate(face): + poly.GetPointIds().SetId(nrv,vert) + polys.InsertNextCell(poly) + + # Make the polydata, structure of connections and vrtx + polyData = vtk.vtkPolyData() + polyData.SetPoints(ptsvtk) + polyData.SetPolys(polys) + + # Make implicit func + ImpDistFunc = vtk.vtkImplicitPolyDataDistance() + ImpDistFunc.SetInput(polyData) + + # Convert the mesh + vtkMesh = vtk.vtkRectilinearGrid() + vtkMesh.SetDimensions(mesh.nNx,mesh.nNy,mesh.nNz) + vtkMesh.SetXCoordinates(npsup.numpy_to_vtk(mesh.vectorNx, deep=1)) + vtkMesh.SetYCoordinates(npsup.numpy_to_vtk(mesh.vectorNy, deep=1)) + vtkMesh.SetZCoordinates(npsup.numpy_to_vtk(mesh.vectorNz, deep=1)) + # Add indexes + vtkInd = npsup.numpy_to_vtk(np.arange(mesh.nC), deep=1) + vtkInd.SetName('Index') + vtkMesh.GetCellData().AddArray(vtkInd) + + extractImpDistRectGridFilt = vtk.vtkExtractGeometry() # Object constructor + extractImpDistRectGridFilt.SetImplicitFunction(ImpDistFunc) # + extractImpDistRectGridFilt.SetInputData(vtkMesh) + + if boundaries is True: + extractImpDistRectGridFilt.ExtractBoundaryCellsOn() + + else: + extractImpDistRectGridFilt.ExtractBoundaryCellsOff() + + if internal is True: + extractImpDistRectGridFilt.ExtractInsideOn() + + else: + extractImpDistRectGridFilt.ExtractInsideOff() + + print "Extracting indices from grid..." + # Executing the pipe + extractImpDistRectGridFilt.Update() + + # Get index inside + insideGrid = extractImpDistRectGridFilt.GetOutput() + insideGrid = npsup.vtk_to_numpy(insideGrid.GetCellData().GetArray('Index')) + + + # Return the indexes inside + return insideGrid diff --git a/docs/examples/DC_Forward_PseudoSection.rst b/docs/examples/DC_Forward_PseudoSection.rst index 1a500cae..80cf0307 100644 --- a/docs/examples/DC_Forward_PseudoSection.rst +++ b/docs/examples/DC_Forward_PseudoSection.rst @@ -12,8 +12,16 @@ DC Forward Simulation ===================== -Forward model conductive spheres in a half-space and plot a pseudo-section +Forward model two conductive spheres in a half-space and plot a +pseudo-section. Assumes an infinite line source and measures along the +center of the spheres. +INPUT: +loc = Location of spheres [[x1,y1,z1],[x2,y2,z2]] +radi = Radius of spheres [r1,r2] +param = Conductivity of background and two spheres [m0,m1,m2] +stype = survey type "pdp" (pole dipole) or "dpdp" (dipole dipole) +dtype = Data type "appr" (app res) | "appc" (app cond) | "volt" (potential) Created by @fourndo on Mon Feb 01 19:28:06 2016 diff --git a/docs/examples/EM_Schenkel_Morrison_Casing.rst b/docs/examples/EM_Schenkel_Morrison_Casing.rst new file mode 100644 index 00000000..55f00168 --- /dev/null +++ b/docs/examples/EM_Schenkel_Morrison_Casing.rst @@ -0,0 +1,58 @@ +.. _examples_EM_Schenkel_Morrison_Casing: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +EM: Schenkel and Morrison Casing Model +====================================== + +Here we create and run a FDEM forward simulation to calculate the vertical +current inside a steel-cased. The model is based on the Schenkel and +Morrison Casing Model, and the results are used in a 2016 SEG abstract by +Yang et al. + +- Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on potential field measurements using downhole current sources: Geophysical prospecting, 38, 663-686. + + +The model consists of: +- Air: Conductivity 1e-8 S/m, above z = 0 +- Background: conductivity 1e-2 S/m, below z = 0 +- Casing: conductivity 1e6 S/m + - 300m long + - radius of 0.1m + - thickness of 6e-3m + +Inside the casing, we take the same conductivity as the background. + +We are using an EM code to simulate DC, so we use frequency low enough +that the skin depth inside the casing is longer than the casing length (f += 1e-6 Hz). The plot produced is of the current inside the casing. + +These results are shown in the SEG abstract by Yang et al., 2016: 3D DC +resistivity modeling of steel casing for reservoir monitoring using +equivalent resistor network. The solver used to produce these results and +achieve the CPU time of ~30s is Mumps, which was installed using pymatsolver_ + +.. _pymatsolver: https://github.com/rowanc1/pymatsolver + +This example is on figshare: https://dx.doi.org/10.6084/m9.figshare.3126961.v1 + +If you would use this example for a code comparison, or build upon it, a +citation would be much appreciated! + + + +.. plot:: + + from SimPEG import Examples + Examples.EM_Schenkel_Morrison_Casing.run() + +.. literalinclude:: ../../SimPEG/Examples/EM_Schenkel_Morrison_Casing.py + :language: python + :linenos: diff --git a/docs/examples/Inversion_IRLS.rst b/docs/examples/Inversion_IRLS.rst new file mode 100644 index 00000000..64a072f1 --- /dev/null +++ b/docs/examples/Inversion_IRLS.rst @@ -0,0 +1,26 @@ +.. _examples_Inversion_IRLS: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +Inversion: Linear Problem +========================= + +Here we go over the basics of creating a linear problem and inversion. + + + +.. plot:: + + from SimPEG import Examples + Examples.Inversion_IRLS.run() + +.. literalinclude:: ../../SimPEG/Examples/Inversion_IRLS.py + :language: python + :linenos: diff --git a/tests/base/test_regularization.py b/tests/base/test_regularization.py index 050c46ac..6dddfb74 100644 --- a/tests/base/test_regularization.py +++ b/tests/base/test_regularization.py @@ -5,6 +5,8 @@ from scipy.sparse.linalg import dsolve import inspect TOL = 1e-20 +testReg = True +testRegMesh = True class RegularizationTests(unittest.TestCase): @@ -16,44 +18,80 @@ class RegularizationTests(unittest.TestCase): mesh3 = Mesh.TensorMesh([hx, hy, hz]) self.meshlist = [mesh1,mesh2, mesh3] - def test_regularization(self): - for R in dir(Regularization): - r = getattr(Regularization, R) - if not inspect.isclass(r): continue - if not issubclass(r, Regularization.BaseRegularization): - continue + if testReg: + def test_regularization(self): + for R in dir(Regularization): + r = getattr(Regularization, R) + if not inspect.isclass(r): continue + if not issubclass(r, Regularization.BaseRegularization): + continue + + for i, mesh in enumerate(self.meshlist): + + print 'Testing %iD'%mesh.dim + + mapping = r.mapPair(mesh) + reg = r(mesh, mapping=mapping) + m = np.random.rand(mapping.nP) + reg.mref = np.ones_like(m)*np.mean(m) + + print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref) + passed = reg.eval(reg.mref) < TOL + self.assertTrue(passed) + + print 'Check:', R + passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False) + self.assertTrue(passed) + + print 'Check 2 Deriv:', R + passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False) + self.assertTrue(passed) + + def test_regularization_ActiveCells(self): + for R in dir(Regularization): + r = getattr(Regularization, R) + if not inspect.isclass(r): continue + if not issubclass(r, Regularization.BaseRegularization): + continue + + for i, mesh in enumerate(self.meshlist): + + print 'Testing Active Cells %iD'%(mesh.dim) + + if mesh.dim == 1: + indActive = Utils.mkvc(mesh.gridCC <= 0.8) + elif mesh.dim == 2: + indActive = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5) + elif mesh.dim == 3: + indActive = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5) + + for indAct in [indActive, indActive.nonzero()[0]]: # test both bool and integers + reg = r(mesh, indActive=indAct) + m = np.random.rand(mesh.nC)[indAct] + reg.mref = np.ones_like(m)*np.mean(m) + + print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref) + passed = reg.eval(reg.mref) < TOL + self.assertTrue(passed) + + print 'Check:', R + passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False) + self.assertTrue(passed) + + print 'Check 2 Deriv:', R + passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False) + self.assertTrue(passed) + + if testRegMesh: + def test_regularizationMesh(self): for i, mesh in enumerate(self.meshlist): print 'Testing %iD'%mesh.dim - mapping = r.mapPair(mesh) - reg = r(mesh, mapping=mapping) - m = np.random.rand(mapping.nP) - reg.mref = np.ones_like(m)*np.mean(m) - - print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref) - passed = reg.eval(reg.mref) < TOL - self.assertTrue(passed) - - print 'Check:', R - passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False) - self.assertTrue(passed) - - print 'Check 2 Deriv:', R - passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False) - self.assertTrue(passed) - - def test_regularization_ActiveCells(self): - for R in dir(Regularization): - r = getattr(Regularization, R) - if not inspect.isclass(r): continue - if not issubclass(r, Regularization.BaseRegularization): - continue - - for i, mesh in enumerate(self.meshlist): - - print 'Testing Active Cells %iD'%(mesh.dim) + # mapping = r.mapPair(mesh) + # reg = r(mesh, mapping=mapping) + # m = np.random.rand(mapping.nP) if mesh.dim == 1: indAct = Utils.mkvc(mesh.gridCC <= 0.8) @@ -62,23 +100,9 @@ class RegularizationTests(unittest.TestCase): elif mesh.dim == 3: indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5) - mapping = Maps.IdentityMap(nP=indAct.nonzero()[0].size) + regmesh = Regularization.RegularizationMesh(mesh, indActive=indAct) - reg = r(mesh, mapping=mapping, indActive=indAct) - m = np.random.rand(mesh.nC)[indAct] - reg.mref = np.ones_like(m)*np.mean(m) - - print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref) - passed = reg.eval(reg.mref) < TOL - self.assertTrue(passed) - - print 'Check:', R - passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False) - self.assertTrue(passed) - - print 'Check 2 Deriv:', R - passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False) - self.assertTrue(passed) + assert (regmesh.vol == mesh.vol[indAct]).all() if __name__ == '__main__': diff --git a/tests/em/tdem/test_TDEM_forward_Analytic.py b/tests/em/tdem/test_TDEM_forward_Analytic.py index dc3696ef..7748b5fc 100644 --- a/tests/em/tdem/test_TDEM_forward_Analytic.py +++ b/tests/em/tdem/test_TDEM_forward_Analytic.py @@ -10,7 +10,9 @@ except ImportError, e: MumpsSolver = SolverLU -def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,1e-3], showIt=False): +def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=None, showIt=False): + if bounds is None: + bounds = [1e-5,1e-3] if meshType == 'CYL': cs, ncx, ncz, npad = 5., 30, 10, 15 hx = [(cs,ncx), (cs,npad,1.3)] diff --git a/tests/mesh/test_cylMesh.py b/tests/mesh/test_cylMesh.py index 5a963398..c7cd44ae 100644 --- a/tests/mesh/test_cylMesh.py +++ b/tests/mesh/test_cylMesh.py @@ -146,6 +146,20 @@ class TestCyl2DMesh(unittest.TestCase): assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3 + def test_getInterpMatCartMesh_Cells2Nodes(self): + + Mr = Mesh.TensorMesh([100,100,2], x0='CC0') + Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0]) + + mc = np.arange(Mc.nC) + xr = np.linspace(0,0.4,50) + xc = np.linspace(0,0.4,50) + 0.2 + Pr = Mr.getInterpolationMat(np.c_[xr,np.ones(50)*-0.2,np.ones(50)*0.5],'N') + Pc = Mc.getInterpolationMat(np.c_[xc,np.zeros(50),np.ones(50)*0.5],'CC') + Pc2r = Mc.getInterpolationMatCartMesh(Mr, 'CC', locTypeTo='N') + + assert np.abs(Pr*(Pc2r*mc) - Pc*mc).max() < 1e-3 + def test_getInterpMatCartMesh_Faces(self): Mr = Mesh.TensorMesh([100,100,2], x0='CC0') @@ -177,6 +191,37 @@ class TestCyl2DMesh(unittest.TestCase): assert np.abs(mag[dist > 0.1].min() - 1) < TOL + def test_getInterpMatCartMesh_Faces2Edges(self): + + Mr = Mesh.TensorMesh([100,100,2], x0='CC0') + Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0]) + + Pf2e = Mc.getInterpolationMatCartMesh(Mr, 'F', locTypeTo='E') + mf = np.ones(Mc.nF) + + ecart = Pf2e * mf + + excc = Mr.aveEx2CC*Mr.r(ecart, 'E', 'Ex') + eycc = Mr.aveEy2CC*Mr.r(ecart, 'E', 'Ey') + ezcc = Mr.r(ecart, 'E', 'Ez') + + indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5]) + indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5]) + + TOL = 1e-2 + assert np.abs(float(excc[indX]) - 1) < TOL + assert np.abs(float(excc[indY]) - 0) < TOL + assert np.abs(float(eycc[indX]) - 0) < TOL + assert np.abs(float(eycc[indY]) - 1) < TOL + assert np.abs((ezcc - 1).sum()) < TOL + + mag = (excc**2 + eycc**2)**0.5 + dist = ((Mr.gridCC[:,0] + 0.2)**2 + (Mr.gridCC[:,1] + 0.2)**2)**0.5 + + assert np.abs(mag[dist > 0.1].max() - 1) < TOL + assert np.abs(mag[dist > 0.1].min() - 1) < TOL + + def test_getInterpMatCartMesh_Edges(self): Mr = Mesh.TensorMesh([100,100,2], x0='CC0') @@ -185,11 +230,42 @@ class TestCyl2DMesh(unittest.TestCase): Pe = Mc.getInterpolationMatCartMesh(Mr, 'E') me = np.ones(Mc.nE) - erect = Pe * me + ecart = Pe * me - excc = Mr.aveEx2CC*Mr.r(erect, 'E', 'Ex') - eycc = Mr.aveEy2CC*Mr.r(erect, 'E', 'Ey') - ezcc = Mr.r(erect, 'E', 'Ez') + excc = Mr.aveEx2CC*Mr.r(ecart, 'E', 'Ex') + eycc = Mr.aveEy2CC*Mr.r(ecart, 'E', 'Ey') + ezcc = Mr.aveEz2CC*Mr.r(ecart, 'E', 'Ez') + + indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5]) + indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5]) + + TOL = 1e-2 + assert np.abs(float(excc[indX]) - 0) < TOL + assert np.abs(float(excc[indY]) + 1) < TOL + assert np.abs(float(eycc[indX]) - 1) < TOL + assert np.abs(float(eycc[indY]) - 0) < TOL + assert np.abs(ezcc.sum()) < TOL + + mag = (excc**2 + eycc**2)**0.5 + dist = ((Mr.gridCC[:,0] + 0.2)**2 + (Mr.gridCC[:,1] + 0.2)**2)**0.5 + + assert np.abs(mag[dist > 0.1].max() - 1) < TOL + assert np.abs(mag[dist > 0.1].min() - 1) < TOL + + + def test_getInterpMatCartMesh_Edges2Faces(self): + + Mr = Mesh.TensorMesh([100,100,2], x0='CC0') + Mc = Mesh.CylMesh([np.ones(10)/5,1,10],x0='0C0',cartesianOrigin=[-0.2,-0.2,0]) + + Pe2f = Mc.getInterpolationMatCartMesh(Mr, 'E', locTypeTo='F') + me = np.ones(Mc.nE) + + frect = Pe2f * me + + excc = Mr.aveFx2CC*Mr.r(frect, 'F', 'Fx') + eycc = Mr.aveFy2CC*Mr.r(frect, 'F', 'Fy') + ezcc = Mr.r(frect, 'F', 'Fz') indX = Utils.closestPoints(Mr, [0.45, -0.2, 0.5]) indY = Utils.closestPoints(Mr, [-0.2, 0.45, 0.5]) diff --git a/tests/mt/test_Problem3D_againstAnalytic.py b/tests/mt/test_Problem3D_againstAnalytic.py index f68e515f..602dbaeb 100644 --- a/tests/mt/test_Problem3D_againstAnalytic.py +++ b/tests/mt/test_Problem3D_againstAnalytic.py @@ -242,9 +242,6 @@ class TestAnalytics(unittest.TestCase): def test_appRes1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3)) def test_appPhs1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3,False)) - # Do a derivative test - def test_derivProj1(self):self.assertTrue(DerivProjfieldsTest(halfSpace(1e-2))) - # Do a derivative test of Jvec # def test_derivJvec_zxxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxr',.1)) # def test_derivJvec_zxxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxi',.1))