From 606488d152591605876c1f767b0c0afcac355034 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Thu, 21 Apr 2016 21:58:40 -0700 Subject: [PATCH 01/10] Major fix to IRLS. --- SimPEG/DCIP/DCIPUtils.py | 152 +++++++++++++++++++----------- SimPEG/Directives.py | 4 +- SimPEG/Examples/Inversion_IRLS.py | 6 +- SimPEG/Regularization.py | 12 +-- 4 files changed, 109 insertions(+), 65 deletions(-) diff --git a/SimPEG/DCIP/DCIPUtils.py b/SimPEG/DCIP/DCIPUtils.py index e94b930f..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='dpdp', dtype="appc", clim=None): +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. @@ -192,9 +192,6 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None): from scipy.interpolate import griddata import pylab as plt - # Set depth to 0 for now - z0 = 0. - # Pre-allocate midx = [] midz = [] @@ -259,38 +256,53 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None): midx = np.hstack([midx, ( Cmid + Pmid )/2 ]) midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + (Tx[0][2] + Tx[1][2])/2 ]) - ax = axs - # 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, clim=(vmin, vmax)) - cbar = plt.colorbar(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) + 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 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) + 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") - # Plot apparent resistivity - ax.scatter(midx,midz,s=10,c=rho.T, vmin =vmin, vmax = vmax, clim=(vmin, vmax)) + 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) - #ax.set_xticklabels([]) - #ax.set_yticklabels([]) + + + if not axlabel: + axs.set_xticklabels([]) + axs.set_yticklabels([]) plt.gca().set_aspect('equal', adjustable='box') @@ -448,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 @@ -471,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): @@ -498,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') @@ -506,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': @@ -529,11 +554,12 @@ 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() @@ -640,51 +666,59 @@ def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'): 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 survey + Read UBC GIF IP 3D observation file and generate survey Input: :param fileName, path to the UBC GIF 3D obs file Output: - :param DCIPsurvey + :param IPsurvey :return - Created on Mon April 6th, 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] @@ -692,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]) @@ -703,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]) @@ -711,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:]) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index b694a2f1..193911fe 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -300,6 +300,8 @@ class Update_IRLS(InversionDirective): self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d + self.reg._W = None + class Update_lin_PreCond(InversionDirective): @@ -308,7 +310,7 @@ class Update_lin_PreCond(InversionDirective): 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. + 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(diagA**-1.) self.opt.approxHinv = PC print 'Updated pre-cond' diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index c236860d..06ef4be4 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -86,12 +86,12 @@ def run(N=200, plotIt=True): #reg.recModel = mrec reg.wght = np.ones(mesh.nC) reg.mref = np.zeros(mesh.nC) - reg.eps_p = 2e-3 - reg.eps_q = 2e-3 + reg.eps_p = 5e-2 + reg.eps_q = 1e-2 reg.norms = [0., 0., 2., 2.] reg.wght = wr - opt = Optimization.ProjectedGNCG(maxIter=5 ,lower=-2.,upper=2., maxIterCG= 100, tolCG = 1e-3) + 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() diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 195e6203..ed1039ec 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -728,14 +728,14 @@ class Sparse(Simple): @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 sp.vstack(wlist) + if getattr(self, '_W', None) is None: + wlist = (self.Wsmall, self.Wsmooth) + self._W = sp.vstack(wlist) + return self._W def R(self, f_m , eps, exponent): - eta = (eps**(1-exponent/2.))**0.5 - r = eta / (f_m**2.+ eps**2.)**((1-exponent/2.)/2.) + eta = (eps**(1.-exponent/2.))**0.5 + r = eta / (f_m**2.+ eps**2.)**((1.-exponent/2.)/2.) return r From 79183ae9fbfc06d0e0575d1b90a12f35dfc5b91c Mon Sep 17 00:00:00 2001 From: D Fournier Date: Fri, 22 Apr 2016 16:05:43 -0700 Subject: [PATCH 02/10] fIX MESH io --- SimPEG/Mesh/MeshIO.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) 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='!') From d8bfb27415ad6f34a38e8097f32dae8c473caddf Mon Sep 17 00:00:00 2001 From: D Fournier Date: Sat, 23 Apr 2016 15:25:44 -0700 Subject: [PATCH 03/10] Quick fix to MeshIO --- SimPEG/Directives.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 193911fe..3cdb5000 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -301,7 +301,7 @@ class Update_IRLS(InversionDirective): self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d self.reg._W = None - + class Update_lin_PreCond(InversionDirective): @@ -313,7 +313,7 @@ class Update_lin_PreCond(InversionDirective): 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(diagA**-1.) self.opt.approxHinv = PC - print 'Updated pre-cond' + class Update_Wj(InversionDirective): """ From 225394f74e8c03ba5d938cd9eed70f726d4f27b7 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Fri, 29 Apr 2016 11:10:04 -0700 Subject: [PATCH 04/10] Latest commit --- SimPEG/Directives.py | 28 +++++++++++++++++++++++----- 1 file changed, 23 insertions(+), 5 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 3cdb5000..ee29e770 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -258,6 +258,7 @@ class Update_IRLS(InversionDirective): phi_m_last = None phi_d_last = None + def initialize(self): # Scale the regularization for changes in norm @@ -275,7 +276,7 @@ class Update_IRLS(InversionDirective): self.phi_d_last = self.invProb.phi_d def endIter(self): - # Cool the threshold parameter + # Cool the threshold parameter if required if getattr(self, 'factor', None) is not None: eps = self.reg.eps / self.factor @@ -290,16 +291,17 @@ class Update_IRLS(InversionDirective): # Update the model used for the IRLS weights self.reg.curModel = self.invProb.curModel - # Temporarely set gamma to 1. + # Temporarely set gamma to 1. to get raw phi_m self.reg.gamma = 1. - # Compute change in model objective function and update scaling + # 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 - self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d - + # 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): @@ -340,3 +342,19 @@ class Update_Wj(InversionDirective): 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 From c83b460672df464598600f7b793927d26e4a4605 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 29 Apr 2016 11:43:31 -0700 Subject: [PATCH 05/10] Surface to Indices (GoCAD and VTK) --- SimPEG/Utils/io_utils.py | 132 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 132 insertions(+) create mode 100644 SimPEG/Utils/io_utils.py diff --git a/SimPEG/Utils/io_utils.py b/SimPEG/Utils/io_utils.py new file mode 100644 index 00000000..f68ee68c --- /dev/null +++ b/SimPEG/Utils/io_utils.py @@ -0,0 +1,132 @@ +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 + + """ + + + 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=True 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 From 028a16a45a3f0306e4911f675e29f5acee120de1 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 29 Apr 2016 11:44:42 -0700 Subject: [PATCH 06/10] Syntax bug. --- SimPEG/Utils/io_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Utils/io_utils.py b/SimPEG/Utils/io_utils.py index f68ee68c..e46c13db 100644 --- a/SimPEG/Utils/io_utils.py +++ b/SimPEG/Utils/io_utils.py @@ -113,7 +113,7 @@ def surface2inds(vrtx, trgl, mesh, boundaries=True, internal=True): else: extractImpDistRectGridFilt.ExtractBoundaryCellsOff() - if internal=True is True: + if internal is True: extractImpDistRectGridFilt.ExtractInsideOn() else: From 00db6746d4272dd953773300dc815c0a0ee0008e Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Fri, 29 Apr 2016 11:50:56 -0700 Subject: [PATCH 07/10] Add a warnign about mesh attributes --- SimPEG/Utils/io_utils.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/SimPEG/Utils/io_utils.py b/SimPEG/Utils/io_utils.py index e46c13db..323feb0c 100644 --- a/SimPEG/Utils/io_utils.py +++ b/SimPEG/Utils/io_utils.py @@ -18,6 +18,11 @@ def read_GOCAD_ts(tsfile): Author: @fourndo + + .. note:: + + Remove all attributes from the GoCAD surface before exporting it! + """ From 056dc09fa63a3cc0e6582a915b5605f1f4206338 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Fri, 29 Apr 2016 15:10:30 -0700 Subject: [PATCH 08/10] Fix Update_Precondition directive --- SimPEG/Directives.py | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index ee29e770..b32cd798 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -305,11 +305,24 @@ class Update_IRLS(InversionDirective): self.reg._W = None class Update_lin_PreCond(InversionDirective): - - + """ + Create a Jacobi preconditioner for the linear problem + """ + onlyOnStart=True + + 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(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. From 3d1dfc13d7a2d5aee0f0859c1df28fb3be617fcd Mon Sep 17 00:00:00 2001 From: D Fournier Date: Fri, 29 Apr 2016 15:49:44 -0700 Subject: [PATCH 09/10] Change Update_PreConditioner to default False --- SimPEG/Directives.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index b32cd798..e51c904e 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -308,7 +308,7 @@ class Update_lin_PreCond(InversionDirective): """ Create a Jacobi preconditioner for the linear problem """ - onlyOnStart=True + onlyOnStart=False def initialize(self): From 4e296c4cd5194495b9b397fc87a724925ca608a6 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Wed, 4 May 2016 16:01:29 -0700 Subject: [PATCH 10/10] Update PreCond Directive to allow inactive cells mapping --- SimPEG/Directives.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index e51c904e..4fc1ffc5 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -315,7 +315,7 @@ class Update_lin_PreCond(InversionDirective): 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(diagA**-1.) + PC = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.) self.opt.approxHinv = PC def endIter(self): @@ -326,7 +326,7 @@ class Update_lin_PreCond(InversionDirective): 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(diagA**-1.) + PC = Utils.sdiag((self.prob.mapping.deriv(None).T *diagA)**-1.) self.opt.approxHinv = PC