From a0174e4f3043216cbc302be4fd775cda2674cf95 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 29 Apr 2016 12:52:45 -0700 Subject: [PATCH 01/27] kwarg name updates --- SimPEG/DCIP/DCIPUtils.py | 84 ++++++++++----------- SimPEG/Examples/DC_Forward_PseudoSection.py | 32 ++++---- 2 files changed, 55 insertions(+), 61 deletions(-) diff --git a/SimPEG/DCIP/DCIPUtils.py b/SimPEG/DCIP/DCIPUtils.py index e94b930f..91eab0b3 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, surveyType='dipole-dipole', unitType='volt', clim=None): """ Read list of 2D tx-rx location and plot a speudo-section of apparent resistivity. @@ -177,16 +177,12 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None): Assumes flat topo for now... Input: - :param d2D, z0 - :switch stype -> Either 'pdp' (pole-dipole) | 'dpdp' (dipole-dipole) - :switch dtype=-> Either 'appr' (app. res) | 'appc' (app. con) | 'volt' (potential) + :param DCsurvey, z0 + :switch surveyType -> Either 'pole-dipole' | 'dipole-dipole' + :switch unitType=-> Either 'appResistivity' | 'appConductivity' | 'volt' Output: :figure scatter plot overlayed on image - Edited Feb 17th, 2016 - - @author: dominiquef - """ from SimPEG import np from scipy.interpolate import griddata @@ -221,39 +217,39 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None): Cmid = (Tx[0][0] + Tx[1][0])/2 Pmid = (Rx[0][:,0] + Rx[1][:,0])/2 - # Change output for dtype - if dtype == 'volt': + # Change output for unitType + if unitType == 'volt': rho = np.hstack([rho,data]) else: # Compute pant leg of apparent rho - if stype == 'pdp': + if surveyType == 'pole-dipole': leg = data * 2*np.pi * MA * ( MA + MN ) / MN - elif stype == 'dpdp': + elif surveyType == 'dipole-dipole': leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA ) else: - print """dtype must be 'pdp'(pole-dipole) | 'dpdp' (dipole-dipole) """ + print """unitType must be 'pole-dipole' | 'dipole-dipole' """ break - if dtype == 'appc': + if unitType == 'appConductivity': leg = np.log10(abs(1./leg)) rho = np.hstack([rho,leg]) - elif dtype == 'appr': + elif unitType == 'appResistivity': leg = np.log10(abs(leg)) rho = np.hstack([rho,leg]) else: - print """dtype must be 'appr' | 'appc' | 'volt' """ + print """unitType must be 'appResistivity' | 'appConductivity' | 'volt' """ break midx = np.hstack([midx, ( Cmid + Pmid )/2 ]) @@ -278,12 +274,12 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None): ticks = np.linspace(cmin,cmax,3) cbar.set_ticks(ticks) cbar.ax.tick_params(labelsize=10) - - if dtype == 'appc': + + if unitType == 'appConductivity': cbar.set_label("App.Cond",size=12) - elif dtype == 'appr': + elif unitType == 'appResistivity': cbar.set_label("App.Res.",size=12) - elif dtype == 'volt': + elif unitType == 'volt': cbar.set_label("Potential (V)",size=12) # Plot apparent resistivity @@ -298,7 +294,7 @@ def plot_pseudoSection(DCsurvey, axs, stype='dpdp', dtype="appc", clim=None): return ph -def gen_DCIPsurvey(endl, mesh, stype, a, b, n): +def gen_DCIPsurvey(endl, mesh, surveyType, a, b, n): """ Load in endpoints and survey specifications to generate Tx, Rx location stations. @@ -308,7 +304,7 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n): Input: :param endl -> input endpoints [x1, y1, z1, x2, y2, z2] :object mesh -> SimPEG mesh object - :switch stype -> "dpdp" (dipole-dipole) | "pdp" (pole-dipole) | 'gradient' + :switch surveyType -> 'dipole-dipole' | 'pole-dipole' | 'gradient' : param a, n -> pole seperation, number of rx dipoles per tx Output: @@ -354,14 +350,14 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n): SrcList = [] - if stype != 'gradient': + if surveyType != 'gradient': for ii in range(0, int(nstn)-1): - if stype == 'dpdp': + if surveyType == 'dipole-dipole': tx = np.c_[M[ii,:],N[ii,:]] - elif stype == 'pdp': + elif surveyType == 'pole-dipole': tx = np.c_[M[ii,:],M[ii,:]] # Rx.append(np.c_[M[ii+1:indx,:],N[ii+1:indx,:]]) @@ -390,13 +386,13 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n): Rx.append(np.c_[P1,P2]) rxClass = DC.RxDipole(P1, P2) Tx.append(tx) - if stype == 'dpdp': + if surveyType == 'dipole-dipole': srcClass = DC.SrcDipole([rxClass], M[ii,:],N[ii,:]) - elif stype == 'pdp': + elif surveyType == 'pole-dipole': srcClass = DC.SrcDipole([rxClass], M[ii,:],M[ii,:]) SrcList.append(srcClass) - elif stype == 'gradient': + elif surveyType == 'gradient': # Gradient survey only requires Tx at end of line and creates a square # grid of receivers at in the middle at a pre-set minimum distance @@ -443,20 +439,20 @@ def gen_DCIPsurvey(endl, mesh, stype, a, b, n): srcClass = DC.SrcDipole([rxClass], M[0,:], N[-1,:]) SrcList.append(srcClass) else: - print """stype must be either 'pdp', 'dpdp' or 'gradient'. """ + print """surveyType must be either 'pole-dipole', 'dipole-dipole' or 'gradient'. """ survey = DC.SurveyDC(SrcList) return survey, Tx, Rx -def writeUBC_DCobs(fileName, DCsurvey, dtype, stype): +def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType): """ 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 dim -> either '2D' | '3D' + :string surveyType -> either 'SURFACE' | 'GENERAL' Output: :param UBC2D-Data file @@ -469,11 +465,11 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype): """ from SimPEG import mkvc - 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'" + assert (dim=='2D') | (dim=='3D'), "Data must be either '2D' | '3D'" + assert (surveyType=='SURFACE') | (surveyType=='GENERAL') | (surveyType=='SIMPLE'), "Data must be either 'SURFACE' | 'GENERAL' | 'SIMPLE'" fid = open(fileName,'w') - fid.write('! ' + stype + ' FORMAT\n') + fid.write('! ' + surveyType + ' FORMAT\n') count = 0 @@ -488,10 +484,10 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype): M = rx[0] N = rx[1] - # Adapt source-receiver location for dtype and stype - if dtype=='2D': + # Adapt source-receiver location for dim and surveyType + if dim=='2D': - if stype == 'SIMPLE': + if surveyType == 'SIMPLE': #fid.writelines("%e " % ii for ii in mkvc(tx[0,:])) A = np.repeat(tx[0,0],M.shape[0],axis=0) @@ -504,13 +500,13 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype): else: - if stype == 'SURFACE': + if surveyType == 'SURFACE': fid.writelines("%e " % ii for ii in mkvc(tx[0,:])) M = M[:,0] N = N[:,0] - if stype == 'GENERAL': + if surveyType == 'GENERAL': fid.writelines("%e " % ii for ii in mkvc(tx[::2,:])) M = M[:,0::2] @@ -519,15 +515,15 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype): 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') - if dtype=='3D': + if dim=='3D': - if stype == 'SURFACE': + if surveyType == 'SURFACE': fid.writelines("%e " % ii for ii in mkvc(tx[0:2,:])) M = M[:,0:2] N = N[:,0:2] - if stype == 'GENERAL': + if surveyType == 'GENERAL': fid.writelines("%e " % ii for ii in mkvc(tx)) @@ -538,7 +534,7 @@ def writeUBC_DCobs(fileName, DCsurvey, dtype, stype): fid.close() -def convertObs_DC3D_to_2D(DCsurvey,lineID, flag = 'local'): +def convertObs_DC3D_to_2D(DCsurvey, lineID, flag='local'): """ Read DC survey and projects the coordinate system according to the flag = 'Xloc' | 'Yloc' | 'local' (default) diff --git a/SimPEG/Examples/DC_Forward_PseudoSection.py b/SimPEG/Examples/DC_Forward_PseudoSection.py index 53467826..240cbb33 100644 --- a/SimPEG/Examples/DC_Forward_PseudoSection.py +++ b/SimPEG/Examples/DC_Forward_PseudoSection.py @@ -2,7 +2,7 @@ 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', dtype='appc', plotIt=True): +def run(loc=None, sig=None, radi=None, param=None, surveyType='dipole-dipole', unitType='appConductivity', plotIt=True): """ DC Forward Simulation ===================== @@ -15,14 +15,14 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p 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) + surveyType = survey type 'pole-dipole' or 'dipole-dipole' + unitType = Data type "appResistivity" | "appConductivity" | "volt" 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)" + assert surveyType in ['pole-dipole', 'dipole-dipole'], "Source type (surveyType) must be pdp or dpdp (pole dipole or dipole dipole)" + assert unitType in ['appResistivity', 'appConductivity', 'volt'], "Unit type (unitType) must be appResistivity or appConductivity or volt (potential)" if loc is None: loc = np.c_[[-50.,0.,-50.],[50.,0.,-50.]] @@ -73,8 +73,8 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p locs = np.c_[mesh.gridCC[indx,0],mesh.gridCC[indx,1],np.ones(2).T*mesh.vectorNz[-1]] # We will handle the geometry of the survey for you and create all the combination of tx-rx along line - # [Tx, Rx] = DC.gen_DCIPsurvey(locs, mesh, stype, param[0], param[1], param[2]) - survey, Tx, Rx = DC.gen_DCIPsurvey(locs, mesh, stype, param[0], param[1], param[2]) + # [Tx, Rx] = DC.gen_DCIPsurvey(locs, mesh, surveyType, param[0], param[1], param[2]) + survey, Tx, Rx = DC.gen_DCIPsurvey(locs, mesh, surveyType, param[0], param[1], param[2]) # Define some global geometry dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) ) @@ -118,8 +118,8 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p rxloc_N = np.asarray(Rx[ii][:,3:]) - # For usual cases "dpdp" or "gradient" - if stype == 'pdp': + # For usual cases 'dipole-dipole' or "gradient" + if surveyType == 'pole-dipole': # Create an "inifinity" pole tx = np.squeeze(Tx[ii][:,0:1]) tinf = tx + np.array([dl_x,dl_y,0])*dl_len*2 @@ -157,12 +157,12 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p fig = plt.figure(figsize=(7,7)) ax = plt.subplot(2,1,1, aspect='equal') # 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) + 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', + 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') @@ -188,15 +188,13 @@ def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', dtype='appc', p ax2 = plt.subplot(2,1,2, aspect='equal') # 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) + 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 - dat = DC.plot_pseudoSection(survey2D,ax2,stype=stype, dtype = dtype) - - # plt.scatter(Tx2d[0][:],Tx[0][2,:],s=40,c='g', marker='v') + dat = DC.plot_pseudoSection(survey2D, ax2, surveyType=surveyType, unitType=unitType) # 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') ax2.set_title('Apparent Conductivity data') From 4257ea77b357bdab1811ac8f9adfb82845a6cf55 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 29 Apr 2016 15:09:04 -0700 Subject: [PATCH 02/27] remove InjectActiveCellsTopo. you should use InjectActiveCells --- SimPEG/Maps.py | 77 -------------------------------------------------- 1 file changed, 77 deletions(-) diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index f17d2314..3e97499f 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -533,83 +533,6 @@ class ActiveCells(InjectActiveCells): FutureWarning) InjectActiveCells.__init__(self, mesh, indActive, valInactive, nC) -class InjectActiveCellsTopo(IdentityMap): - """ - Active model parameters. Extend for cells on topography to air cell (only works for tensor mesh) - - """ - - indActive = None #: Active Cells - valInactive = None #: Values of inactive Cells - nC = None #: Number of cells in the full model - - def __init__(self, mesh, indActive, nC=None): - self.mesh = mesh - - self.nC = nC or mesh.nC - - if indActive.dtype is not bool: - z = np.zeros(self.nC,dtype=bool) - z[indActive] = True - indActive = z - self.indActive = indActive - - self.indInactive = np.logical_not(indActive) - inds = np.nonzero(self.indActive)[0] - self.P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(self.nC, self.nP)) - - @property - def shape(self): - return (self.nC, self.nP) - - @property - def nP(self): - """Number of parameters in the model.""" - return self.indActive.sum() - - def _transform(self, m): - val_temp = np.zeros(self.mesh.nC) - val_temp[self.indActive] = m - valInactive = np.zeros(self.mesh.nC) - #1D - if self.mesh.dim == 1: - z_temp = self.mesh.gridCC - val_temp[~self.indActive] = val_temp[np.argmax(z_temp[self.indActive])] - #2D - elif self.mesh.dim == 2: - act_temp = self.indActive.reshape((self.mesh.nCx, self.mesh.nCy), order = 'F') - val_temp = val_temp.reshape((self.mesh.nCx, self.mesh.nCy), order = 'F') - y_temp = self.mesh.gridCC[:,1].reshape((self.mesh.nCx, self.mesh.nCy), order = 'F') - for i in range(self.mesh.nCx): - act_tempx = act_temp[i,:] == 1 - val_temp[i,~act_tempx] = val_temp[i,np.argmax(y_temp[i,act_tempx])] - valInactive[~self.indActive] = Utils.mkvc(val_temp)[~self.indActive] - #3D - elif self.mesh.dim == 3: - act_temp = self.indActive.reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F') - val_temp = val_temp.reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F') - z_temp = self.mesh.gridCC[:,2].reshape((self.mesh.nCx*self.mesh.nCy, self.mesh.nCz), order = 'F') - for i in range(self.mesh.nCx*self.mesh.nCy): - act_tempxy = act_temp[i,:] == 1 - val_temp[i,~act_tempxy] = val_temp[i,np.argmax(z_temp[i,act_tempxy])] - valInactive[~self.indActive] = Utils.mkvc(val_temp)[~self.indActive] - - self.valInactive = valInactive - - return self.P*m + self.valInactive - - def inverse(self, D): - return self.P.T*D - - def deriv(self, m): - return self.P - -class ActiveCellsTopo(InjectActiveCellsTopo): - def __init__(self, mesh, indActive, valInactive, nC=None): - warnings.warn( - "`ActiveCellsTopo` is deprecated and will be removed in future versions. Use `InjectActiveCellsTopo` instead", - FutureWarning) - InjectActiveCellsTopo.__init__(self, mesh, indActive, valInactive, nC) class Weighting(IdentityMap): """ From ba8f270b3adb7d9aa6553b3ffac739ed7fef0cba Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 1 May 2016 13:17:16 -0700 Subject: [PATCH 03/27] start of surface2ind_topo --- SimPEG/Examples/Utils_surface2ind_topo.py | 42 +++++++++++++++ SimPEG/Utils/__init__.py | 1 + SimPEG/Utils/modelutils.py | 63 +++++++++++++++++++++++ 3 files changed, 106 insertions(+) create mode 100644 SimPEG/Examples/Utils_surface2ind_topo.py create mode 100644 SimPEG/Utils/modelutils.py diff --git a/SimPEG/Examples/Utils_surface2ind_topo.py b/SimPEG/Examples/Utils_surface2ind_topo.py new file mode 100644 index 00000000..61a0f937 --- /dev/null +++ b/SimPEG/Examples/Utils_surface2ind_topo.py @@ -0,0 +1,42 @@ +from SimPEG import * +from SimPEG.Utils import surface2ind_topo + + +def run(plotIt=False, nx = 5, ny = 5): + """ + Here we show how to use :code:`Utils.surface2ind_topo` to identify cells below + a topographic surface. + + """ + + mesh = Mesh.TensorMesh([nx,ny], x0='CC') # 2D mesh + xtopo = np.linspace(mesh.gridN[:,0].min(), mesh.gridN[:,0].max()) + topo = 0.4*np.sin(xtopo*5) # define a topographic surface + + Topo = np.hstack([Utils.mkvc(xtopo,2),Utils.mkvc(topo,2)]) #make it an array + + indcc = surface2ind_topo(mesh, Topo,'CC') + + if plotIt: + from matplotlib.pylab import plt + from scipy.interpolate import interp1d + fig, ax = plt.subplots(1,1,figsize=(6,6)) + mesh.plotGrid(ax=ax, nodes=True, centers=True) + ax.plot(xtopo,topo,'k',linewidth=1) + # ax.plot(mesh.vectorNx, interp1d(xtopo,topo)(mesh.vectorNx),'--k',linewidth=3) + ax.plot(mesh.vectorCCx, interp1d(xtopo,topo)(mesh.vectorCCx),'--k',linewidth=3) + + + aveN2CC = Utils.sdiag(mesh.aveN2CC.T.sum(1))*mesh.aveN2CC.T + a = aveN2CC * indcc + a[a > 0] = 1. + a[a < 0.25] = np.nan + a = a.reshape(mesh.vnN, order='F') + masked_array = np.ma.array(a, mask=np.isnan(a)) + ax.pcolor(mesh.vectorNx,mesh.vectorNy,masked_array.T, cmap = plt.cm.gray,alpha=0.2) + plt.show() + + + +if __name__ == '__main__': + run(plotIt=True) diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 18c1994f..c7597fee 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -7,3 +7,4 @@ from CounterUtils import * import ModelBuilder import SolverUtils from coordutils import * +from modelutils import * diff --git a/SimPEG/Utils/modelutils.py b/SimPEG/Utils/modelutils.py new file mode 100644 index 00000000..dad92fae --- /dev/null +++ b/SimPEG/Utils/modelutils.py @@ -0,0 +1,63 @@ +from matutils import mkvc, ndgrid +import numpy as np + +def surface2ind_topo(mesh, topo, gridLoc='CC'): +# def genActiveindfromTopo(mesh, topo): + """ + Get active indices from topography + """ + + + if mesh.dim == 3: + from scipy.interpolate import NearestNDInterpolator + Ftopo = NearestNDInterpolator(topo[:,:2], topo[:,2]) + + if gridLoc == 'CC': + XY = ndgrid(mesh.vectorCCx, mesh.vectorCCy) + Zcc = mesh.gridCC[:,2].reshape((np.prod(mesh.vnC[:2]), mesh.nCz), order='F') + + gridTopo = Ftopo(XY) + actind = [gridTopo[ixy] <= Zcc[ixy,:] for ixy in range(np.prod(mesh.vnC[0]))] + actind = np.hstack(actind) + + elif gridLoc == 'N': + + XY = ndgrid(mesh.vectorNx, mesh.vectorNy) + gridTopo = Ftopo(XY).reshape(mesh.vnN[:2], order='F') + + if mesh._meshType not in ['TENSOR', 'CYL', 'BASETENSOR']: + raise NotImplementedError('Nodal surface2ind_topo not implemented for %s mesh'%mesh._meshType) + + Nz = mesh.vectorNz[1:] # TODO: this will only work for tensor meshes + actind = np.array([False]*mesh.nC).reshape(mesh.vnC, order='F') + + for ii in range(mesh.nCx): + for jj in range(mesh.nCy): + actind[ii,jj,:] = [np.all(gridTopo[ii:ii+2, jj:jj+2] >= Nz[kk]) for kk in range(len(Nz)) ] + + elif mesh.dim == 2: + from scipy.interpolate import interp1d + Ftopo = interp1d(topo[:,0], topo[:,1]) + + if gridLoc == 'CC': + gridTopo = Ftopo(mesh.gridCC[:,0]) + actind = mesh.gridCC[:,1] <= gridTopo + + elif gridLoc == 'N': + + gridTopo = Ftopo(mesh.vectorNx) + if mesh._meshType not in ['TENSOR', 'CYL', 'BASETENSOR']: + raise NotImplementedError('Nodal surface2ind_topo not implemented for %s mesh'%mesh._meshType) + + Ny = mesh.vectorNy[1:] # TODO: this will only work for tensor meshes + actind = np.array([False]*mesh.nC).reshape(mesh.vnC, order='F') + + for ii in range(mesh.nCx): + actind[ii,:] = [np.all(gridTopo[ii:ii+2] > Ny[kk]) for kk in range(len(Ny)) ] + + else: + raise NotImplementedError('surface2ind_topo not implemented for 1D mesh') + + return mkvc(actind) + + From dd45a6a0855a90bc1ef382293661a5240f700463 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Mon, 2 May 2016 11:40:02 -0700 Subject: [PATCH 04/27] name updates in DC_Forward_PseudoSection, DC_Utils, example for Utils_surface2ind_topo --- SimPEG/DCIP/DCIPUtils.py | 113 ++++++++++----------- SimPEG/EM/Base.py | 3 +- SimPEG/Examples/Utils_surface2ind_topo.py | 1 - SimPEG/Examples/__init__.py | 3 +- docs/examples/DC_Forward_PseudoSection.rst | 6 +- docs/examples/Utils_surface2ind_topo.rst | 24 +++++ 6 files changed, 82 insertions(+), 68 deletions(-) create mode 100644 docs/examples/Utils_surface2ind_topo.rst diff --git a/SimPEG/DCIP/DCIPUtils.py b/SimPEG/DCIP/DCIPUtils.py index 91eab0b3..800b478c 100644 --- a/SimPEG/DCIP/DCIPUtils.py +++ b/SimPEG/DCIP/DCIPUtils.py @@ -1,4 +1,4 @@ -from SimPEG import np +from SimPEG import np, Utils import BaseDC as DC import BaseDC as IP @@ -118,34 +118,27 @@ def readUBC_DC3Dobstopo(filename,mesh,topo,probType="CC"): def readUBC_DC2DModel(fileName): """ - Read UBC GIF 2DTensor model and generate 2D Tensor model in simpeg + Read UBC GIF 2DTensor model and generate 2D Tensor model in simpeg - Input: - :param fileName, path to the UBC GIF 2D model file - - Output: - :param SimPEG TensorMesh 2D object - :return - - Created on Thu Nov 12 13:14:10 2015 - - @author: dominiquef + :param string fileName: path to the UBC GIF 2D model file + :rtype: TensorMesh + :return: SimPEG TensorMesh 2D object """ from SimPEG import np, mkvc # Open fileand skip header... assume that we know the mesh already - obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!') + obsfile = np.genfromtxt(fileName, delimiter=' \n', dtype=np.str, comments='!') - dim = np.array(obsfile[0].split(),dtype=float) + dim = np.array(obsfile[0].split(), dtype=float) - temp = np.array(obsfile[1].split(),dtype=float) + temp = np.array(obsfile[1].split(), dtype=float) if len(temp) > 1: model = np.zeros(dim) for ii in range(len(obsfile)-1): - mm = np.array(obsfile[ii+1].split(),dtype=float) + mm = np.array(obsfile[ii+1].split(), dtype=float) model[:,ii] = mm model = model[:,::-1] @@ -153,10 +146,10 @@ def readUBC_DC2DModel(fileName): else: if len(obsfile[1:])==1: - mm = np.array(obsfile[1:].split(),dtype=float) + mm = np.array(obsfile[1:].split(), dtype=float) else: - mm = np.array(obsfile[1:],dtype=float) + mm = np.array(obsfile[1:], dtype=float) # Permute the second dimension to flip the order model = mm.reshape(dim[1],dim[0]) @@ -171,17 +164,16 @@ def readUBC_DC2DModel(fileName): def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt', clim=None): """ - Read list of 2D tx-rx location and plot a speudo-section of apparent - resistivity. + Read list of 2D tx-rx location and plot a speudo-section of apparent + resistivity. - Assumes flat topo for now... + Assumes flat topo for now... - Input: - :param DCsurvey, z0 - :switch surveyType -> Either 'pole-dipole' | 'dipole-dipole' - :switch unitType=-> Either 'appResistivity' | 'appConductivity' | 'volt' - Output: - :figure scatter plot overlayed on image + :param SurveyDC DCsurvey: + :param string surveyType: Either 'pole-dipole' | 'dipole-dipole' + :param string unitType: Either 'appResistivity' | 'appConductivity' | 'volt' + :rtype: matplotlib.plt + :return: figure scatter plot overlayed on image """ from SimPEG import np @@ -294,27 +286,24 @@ def plot_pseudoSection(DCsurvey, axs, surveyType='dipole-dipole', unitType='volt return ph -def gen_DCIPsurvey(endl, mesh, surveyType, a, b, n): +def gen_DCIPsurvey(endl, mesh, surveyType, AM_sep, MN_sep, nrx): """ - Load in endpoints and survey specifications to generate Tx, Rx location - stations. + Load in endpoints and survey specifications to generate Tx, Rx location + stations. - Assumes flat topo for now... + Assumes flat topo for now... - Input: - :param endl -> input endpoints [x1, y1, z1, x2, y2, z2] - :object mesh -> SimPEG mesh object - :switch surveyType -> 'dipole-dipole' | 'pole-dipole' | 'gradient' - : param a, n -> pole seperation, number of rx dipoles per tx + :param numpy.array endl: input endpoints [[x1, y1] , [x2, y2]] + :param Mesh mesh: SimPEG mesh object + :param string surveyType: 'dipole-dipole' | 'pole-dipole' | 'gradient' + :param float AM_sep: transmitter (A) - receiver (M) seperation + :param float b: receiver dipole seperation + :param float nrx: pole seperation, number of rx dipoles per tx - Output: - :param Tx, Rx -> List objects for each tx location - Lines: P1x, P1y, P1z, P2x, P2y, P2z + :rtype: DC.Survey, Src, Rx + :returns: DC survey, Source - Created on Wed December 9th, 2015 - - @author: dominiquef - !! Require clean up to deal with DCsurvey + !! Require clean up to deal with DCsurvey """ from SimPEG import np @@ -330,17 +319,17 @@ def gen_DCIPsurvey(endl, mesh, surveyType, a, b, n): dl_x = ( endl[1,0] - endl[0,0] ) / dl_len dl_y = ( endl[1,1] - endl[0,1] ) / dl_len - nstn = np.floor( dl_len / a ) + nstn = np.floor( dl_len / AM_sep ) # Compute discrete pole location along line - stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*a - stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*a + stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*AM_sep + stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*AM_sep # Create line of P1 locations M = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]] # Create line of P2 locations - N = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] + N = np.c_[stn_x+AM_sep*dl_x, stn_y+AM_sep*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] ## Build list of Tx-Rx locations depending on survey type # Dipole-dipole: Moving tx with [a] spacing -> [AB a MN1 a MN2 ... a MNn] @@ -366,22 +355,22 @@ def gen_DCIPsurvey(endl, mesh, surveyType, a, b, n): AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1]) # Number of receivers to fit - nstn = np.min([np.floor( (AB - b) / a ) , n]) + nstn = np.min([np.floor( (AB - MN_sep) / AM_sep ) , nrx]) # Check if there is enough space, else break the loop if nstn <= 0: continue # Compute discrete pole location along line - stn_x = N[ii,0] + dl_x*b + np.array(range(int(nstn)))*dl_x*a - stn_y = N[ii,1] + dl_y*b + np.array(range(int(nstn)))*dl_y*a + stn_x = N[ii,0] + dl_x*MN_sep + np.array(range(int(nstn)))*dl_x*AM_sep + stn_y = N[ii,1] + dl_y*MN_sep + np.array(range(int(nstn)))*dl_y*AM_sep # Create receiver poles # Create line of P1 locations P1 = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]] # Create line of P2 locations - P2 = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] + P2 = np.c_[stn_x+AM_sep*dl_x, stn_y+AM_sep*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] Rx.append(np.c_[P1,P2]) rxClass = DC.RxDipole(P1, P2) @@ -400,23 +389,23 @@ def gen_DCIPsurvey(endl, mesh, surveyType, a, b, n): Tx.append(np.c_[M[0,:],N[-1,:]]) # Get the edge limit of survey area - min_x = endl[0,0] + dl_x * b - min_y = endl[0,1] + dl_y * b + min_x = endl[0,0] + dl_x * MN_sep + min_y = endl[0,1] + dl_y * MN_sep - max_x = endl[1,0] - dl_x * b - max_y = endl[1,1] - dl_y * b + max_x = endl[1,0] - dl_x * MN_sep + max_y = endl[1,1] - dl_y * MN_sep box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 ) box_w = box_l/2. - nstn = np.floor( box_l / a ) + nstn = np.floor( box_l / AM_sep ) # Compute discrete pole location along line - stn_x = min_x + np.array(range(int(nstn)))*dl_x*a - stn_y = min_y + np.array(range(int(nstn)))*dl_y*a + stn_x = min_x + np.array(range(int(nstn)))*dl_x*AM_sep + stn_y = min_y + np.array(range(int(nstn)))*dl_y*AM_sep # Define number of cross lines - nlin = int(np.floor( box_w / a )) + nlin = int(np.floor( box_w / AM_sep )) lind = range(-nlin,nlin+1) ngrad = nstn * len(lind) @@ -425,12 +414,12 @@ def gen_DCIPsurvey(endl, mesh, surveyType, a, b, n): for ii in range( len(lind) ): # Move line in perpendicular direction by dipole spacing - lxx = stn_x - lind[ii]*a*dl_y - lyy = stn_y + lind[ii]*a*dl_x + lxx = stn_x - lind[ii]*AM_sep*dl_y + lyy = stn_y + lind[ii]*AM_sep*dl_x M = np.c_[ lxx, lyy , np.ones(nstn).T*mesh.vectorNz[-1]] - N = np.c_[ lxx+a*dl_x, lyy+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] + N = np.c_[ lxx+AM_sep*dl_x, lyy+AM_sep*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]] rx[(ii*nstn):((ii+1)*nstn),:] = np.c_[M,N] diff --git a/SimPEG/EM/Base.py b/SimPEG/EM/Base.py index a16cdb91..032b4429 100644 --- a/SimPEG/EM/Base.py +++ b/SimPEG/EM/Base.py @@ -165,7 +165,8 @@ class BaseEMProblem(Problem.BaseProblem): """ Derivative of :code:`MfRho` with respect to the model. """ - return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv) + return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * self.curModel.rhoDeriv + # (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv) # self.curModel.rhoDeriv @property diff --git a/SimPEG/Examples/Utils_surface2ind_topo.py b/SimPEG/Examples/Utils_surface2ind_topo.py index 61a0f937..4ed1cdaf 100644 --- a/SimPEG/Examples/Utils_surface2ind_topo.py +++ b/SimPEG/Examples/Utils_surface2ind_topo.py @@ -37,6 +37,5 @@ def run(plotIt=False, nx = 5, ny = 5): plt.show() - if __name__ == '__main__': run(plotIt=True) diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index c67652a2..24828466 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -20,8 +20,9 @@ import Mesh_QuadTree_HangingNodes import Mesh_Tensor_Creation import MT_1D_ForwardAndInversion import MT_3D_Foward +import Utils_surface2ind_topo -__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"] +__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", "Utils_surface2ind_topo"] ##### AUTOIMPORTS ##### diff --git a/docs/examples/DC_Forward_PseudoSection.rst b/docs/examples/DC_Forward_PseudoSection.rst index 80cf0307..4231e944 100644 --- a/docs/examples/DC_Forward_PseudoSection.rst +++ b/docs/examples/DC_Forward_PseudoSection.rst @@ -20,9 +20,9 @@ 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 +surveyType = survey type 'pole-dipole' or 'dipole-dipole' +unitType = Data type "appResistivity" | "appConductivity" | "volt" +Created by @fourndo diff --git a/docs/examples/Utils_surface2ind_topo.rst b/docs/examples/Utils_surface2ind_topo.rst new file mode 100644 index 00000000..04b27d5d --- /dev/null +++ b/docs/examples/Utils_surface2ind_topo.rst @@ -0,0 +1,24 @@ +.. _examples_Utils_surface2ind_topo: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +Here we show how to use :code:`Utils.surface2ind_topo` to identify cells below +a topographic surface. + + + +.. plot:: + + from SimPEG import Examples + Examples.Utils_surface2ind_topo.run() + +.. literalinclude:: ../../SimPEG/Examples/Utils_surface2ind_topo.py + :language: python + :linenos: From dbdcc3cefb28cfaa7802b28468d5f0cc60f70e9f Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 4 May 2016 22:06:32 -0700 Subject: [PATCH 05/27] use sigma in MfRhoDeriv - due to propmap bug --- SimPEG/EM/Base.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/SimPEG/EM/Base.py b/SimPEG/EM/Base.py index 032b4429..a16cdb91 100644 --- a/SimPEG/EM/Base.py +++ b/SimPEG/EM/Base.py @@ -165,8 +165,7 @@ class BaseEMProblem(Problem.BaseProblem): """ Derivative of :code:`MfRho` with respect to the model. """ - return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * self.curModel.rhoDeriv - # (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv) + return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv) # self.curModel.rhoDeriv @property From 66440b0478b434adaa9bbb07c13ce056fd3f45f3 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 4 May 2016 22:14:41 -0700 Subject: [PATCH 06/27] add depreciation warnings to DCIP utils for activeind from topo --- SimPEG/DCIP/DCIPUtils.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/SimPEG/DCIP/DCIPUtils.py b/SimPEG/DCIP/DCIPUtils.py index 800b478c..f4e643eb 100644 --- a/SimPEG/DCIP/DCIPUtils.py +++ b/SimPEG/DCIP/DCIPUtils.py @@ -1,12 +1,16 @@ from SimPEG import np, Utils import BaseDC as DC import BaseDC as IP +import warnings def getActiveindfromTopo(mesh, topo): # def genActiveindfromTopo(mesh, topo): """ Get active indices from topography """ + warnings.warn( + "`getActiveindfromTopo` is deprecated and will be removed in future versions. Use `SimPEG.Utils.surface2ind_topo` instead", + FutureWarning) from scipy.interpolate import NearestNDInterpolator if mesh.dim==3: nCxy = mesh.nCx*mesh.nCy @@ -28,6 +32,9 @@ def gettopoCC(mesh, airind): """ Get topography from active indices of mesh. """ + warnings.warn( + "`gettopoCC` is deprecated and will be removed in future versions. Use `SimPEG.Utils.surface2ind_topo` instead", + FutureWarning) mesh2D = Mesh.TensorMesh([mesh.hx, mesh.hy], mesh.x0[:2]) zc = mesh.gridCC[:,2] AIRIND = airind.reshape((mesh.vnC[0]*mesh.vnC[1],mesh.vnC[2]), order='F') @@ -630,16 +637,9 @@ def readUBC_DC3Dobs(fileName): """ Read UBC GIF DCIP 3D observation file and generate survey - Input: - :param fileName, path to the UBC GIF 3D obs file - - Output: - :param DCIPsurvey - :return - - Created on Mon April 6th, 2015 - - @author: dominiquef + :param string fileName:, path to the UBC GIF 3D obs file + :rtype: Survey + :return: DCIPsurvey """ From 0379df2bf22bed20151ed29ece5e41e95e6e98cf Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 4 May 2016 22:27:02 -0700 Subject: [PATCH 07/27] attempt to clean up docs in DCIP utils --- SimPEG/DCIP/DCIPUtils.py | 73 ++++++++++++---------------------------- 1 file changed, 22 insertions(+), 51 deletions(-) diff --git a/SimPEG/DCIP/DCIPUtils.py b/SimPEG/DCIP/DCIPUtils.py index f4e643eb..49f9ba2b 100644 --- a/SimPEG/DCIP/DCIPUtils.py +++ b/SimPEG/DCIP/DCIPUtils.py @@ -444,21 +444,14 @@ def writeUBC_DCobs(fileName, DCsurvey, dim, surveyType): """ 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 dim -> either '2D' | '3D' - :string surveyType -> either 'SURFACE' | 'GENERAL' - - Output: - :param UBC2D-Data file - :return - - Last edit: February 16th, 2016 - - @author: dominiquef - + :param string fileName: including path where the file is written out + :param Survey DCsurvey: DC survey class object + :param string dim: either '2D' | '3D' + :param string surveyType: either 'SURFACE' | 'GENERAL' + :rtype: file + :return: UBC2D-Data file """ + from SimPEG import mkvc assert (dim=='2D') | (dim=='3D'), "Data must be either '2D' | '3D'" @@ -539,15 +532,9 @@ def convertObs_DC3D_to_2D(DCsurvey, lineID, flag='local'): The Z value is preserved, but Y coordinates zeroed. - Input: - :param survey3D - - Output: - :figure survey2D - - Edited April 6th, 2016 - - @author: dominiquef + :param DC.Survey survey3D: 3D simpeg DC survey + :rtype: DC.Survey + :return: survey2D """ from SimPEG import np @@ -715,17 +702,9 @@ def readUBC_DC2Dobs(fileName): ------- NEEDS TO BE UPDATED ------ Read UBC GIF 2D observation file and generate arrays for tx-rx location - Input: - :param fileName, path to the UBC GIF 2D model file - - Output: - :param rx, tx - :return - - Created on Thu Nov 12 13:14:10 2015 - - @author: dominiquef - + :param string fileName: path to the UBC GIF 2D model file + :rtype: (DC.Src, DC.Rx, ??, ??) + :return: source_locs, rx_locs, ??, ?? """ from SimPEG import np @@ -765,11 +744,9 @@ 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 + :param string fileName: path to the UBC GIF 3D obs file + :rtype: DC.Survey + :return: DCsurvey Created on Mon March 9th, 2016 << Doug's 70th Birthday !! >> @@ -831,12 +808,9 @@ def readUBC_DC2DMesh(fileName): """ Read UBC GIF 2DTensor mesh and generate 2D Tensor mesh in simpeg - Input: - :param fileName, path to the UBC GIF mesh file - - Output: - :param SimPEG TensorMesh 2D object - :return + :param string fileName: path to the UBC GIF mesh file + :rtype: Mesh.TensorMesh + :return: SimPEG TensorMesh 2D object Created on Thu Nov 12 13:14:10 2015 @@ -902,12 +876,9 @@ def xy_2_lineID(DCsurvey): they were collected. May need to generalize for random point locations, but will be more expensive - Input: - :param DCdict Vectors of station location - - Output: - :param LineID Vector of integers - :return + :param numpy.array DCdict: Vectors of station location + :rtype: numpy.array + :return: LineID Vector of integers Created on Thu Feb 11, 2015 From fbb8cf2731e10000f71e2056215a363ffc0d4023 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 4 May 2016 23:17:01 -0700 Subject: [PATCH 08/27] modularizing regularization --- SimPEG/Regularization.py | 404 ++++++++++++++++++++++++++------------- 1 file changed, 269 insertions(+), 135 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index fc101a61..fd974130 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -403,7 +403,203 @@ class BaseRegularization(object): return mD.T * ( self.W.T * ( self.W * ( mD * v) ) ) -class Tikhonov(BaseRegularization): +class Simple(BaseRegularization): + """ + Simple regularization that does not include length scales in the derivatives. + """ + + mrefInSmooth = False #: SMOOTH and SMOOTH_MOD_DIF options + 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") + wght = 1. + + def __init__(self, mesh, mapping=None, indActive=None, **kwargs): + BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) + + if isinstance(self.wght,float): + self.wght = np.ones(self.regmesh.nC) * self.wght + + @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 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 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 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,) + 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.Wsmall, self.Wsmooth) + self._W = sp.vstack(wlist) + return self._W + + @Utils.timeIt + def _evalSmall(self, m): + r = self.Wsmall * ( self.mapping * (m - self.mref) ) + return 0.5 * r.dot(r) + + @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 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 + + + @Utils.timeIt + def _evalSmall(self, m): + r = self.Wsmall * ( self.mapping * (m - self.mref) ) + return 0.5 * r.dot(r) + + @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 _evalSmoothx(self, m): + if self.mrefInSmooth == True: + r = self.Wx * ( self.mapping * (m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wx * ( self.mapping * (m) ) + return 0.5 * r.dot(r) + + @Utils.timeIt + def _evalSmoothy(self, m): + if self.mrefInSmooth == True: + r = self.Wy * ( self.mapping * (m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wy * ( self.mapping * (m) ) + return 0.5 * r.dot(r) + + @Utils.timeIt + def _evalSmoothz(self, m): + if self.mrefInSmooth == True: + r = self.Wz * ( self.mapping * (m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wz * ( self.mapping * (m) ) + return 0.5 * r.dot(r) + + @Utils.timeIt + def _evalSmooth(self, m): + phiSmooth = self._evalSmoothx(m) + if self.regmesh.dim > 1: + phiSmooth += self._evalSmoothy(m) + if self.regmesh.dim > 2: + phiSmooth += self._evalSmoothz(m) + return phiSmooth + + @Utils.timeIt + def _evalSmoothxDeriv(self,m): + if self.mrefInSmooth == True: + r = self.Wx * ( self.mapping * ( m - self.mref ) ) + return r.T * ( self.Wx * self.mapping.deriv(m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wx * ( self.mapping * m ) + return r.T * ( self.Wx * self.mapping.deriv(m) ) + + @Utils.timeIt + def _evalSmoothyDeriv(self,m): + if self.mrefInSmooth == True: + r = self.Wy * ( self.mapping * ( m - self.mref ) ) + return r.T * ( self.Wy * self.mapping.deriv(m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wy * ( self.mapping * m ) + return r.T * ( self.Wy * self.mapping.deriv(m) ) + + @Utils.timeIt + def _evalSmoothzDeriv(self,m): + if self.mrefInSmooth == True: + r = self.Wz * ( self.mapping * ( m - self.mref ) ) + return r.T * ( self.Wz * self.mapping.deriv(m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wz * ( self.mapping * m ) + return r.T * ( self.Wz * self.mapping.deriv(m) ) + + @Utils.timeIt + def _evalSmoothDeriv(self,m): + deriv = self._evalSmoothxDeriv(m) + if self.regmesh.dim > 1: + deriv += self._evalSmoothyDeriv(m) + if self.regmesh.dim > 2: + deriv += self._evalSmoothzDeriv(m) + return deriv + + + @Utils.timeIt + def eval(self, m): + return self._evalSmall(m) + self._evalSmooth(m) + + @Utils.timeIt + def evalDeriv(self, m): + """ + The regularization is: + + .. math:: + + R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})} + + So the derivative is straight forward: + + .. math:: + + R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} + + """ + return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + + + +class Tikhonov(Simple): """ L2 Tikhonov regularization with both smallness and smoothness (first order derivative) contributions. @@ -493,56 +689,92 @@ class Tikhonov(BaseRegularization): self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDiffz*self.regmesh.cellDiffz return self._Wzz + @property - def Wsmooth(self): + def Wsmooth2(self): """Full smoothness regularization matrix W""" if getattr(self, '_Wsmooth', None) is None: - wlist = (self.Wx, self.Wxx) + wlist = (self.Wxx) if self.regmesh.dim > 1: - wlist += (self.Wy, self.Wyy) + wlist += (self.Wyy) if self.regmesh.dim > 2: - wlist += (self.Wz, self.Wzz) + wlist += (self.Wzz) 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.Wsmall, self.Wsmooth) - self._W = sp.vstack(wlist) - return self._W - @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): + def _evalSmoothxx(self, m): if self.mrefInSmooth == True: - r = self.Wsmooth * ( self.mapping * (m - self.mref) ) + r = self.Wxx * ( self.mapping * (m - self.mref) ) elif self.mrefInSmooth == False: - r = self.Wsmooth * ( self.mapping * (m) ) + r = self.Wxx * ( self.mapping * (m) ) return 0.5 * r.dot(r) + @Utils.timeIt + def _evalSmoothyy(self, m): + if self.mrefInSmooth == True: + r = self.Wyy * ( self.mapping * (m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wyy * ( self.mapping * (m) ) + return 0.5 * r.dot(r) + + @Utils.timeIt + def _evalSmoothzz(self, m): + if self.mrefInSmooth == True: + r = self.Wzz * ( self.mapping * (m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wzz * ( self.mapping * (m) ) + return 0.5 * r.dot(r) + + @Utils.timeIt + def _evalSmooth2(self, m): + phiSmooth2 = self._evalSmoothxx(m) + if self.regmesh.dim > 1: + phiSmooth2 += self._evalSmoothyy(m) + if self.regmesh.dim > 2: + phiSmooth2 += self._evalSmoothzz(m) + return phiSmooth2 + + @Utils.timeIt + def _evalSmoothxxDeriv(self,m): + if self.mrefInSmooth == True: + r = self.Wxx * ( self.mapping * ( m - self.mref ) ) + return r.T * ( self.Wxx * self.mapping.deriv(m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wxx * ( self.mapping * m ) + return r.T * ( self.Wxx * self.mapping.deriv(m) ) + + @Utils.timeIt + def _evalSmoothyyDeriv(self,m): + if self.mrefInSmooth == True: + r = self.Wyy * ( self.mapping * ( m - self.mref ) ) + return r.T * ( self.Wyy * self.mapping.deriv(m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wyy * ( self.mapping * m ) + return r.T * ( self.Wyy * self.mapping.deriv(m) ) + + @Utils.timeIt + def _evalSmoothzzDeriv(self,m): + if self.mrefInSmooth == True: + r = self.Wzz * ( self.mapping * ( m - self.mref ) ) + return r.T * ( self.Wzz * self.mapping.deriv(m - self.mref) ) + elif self.mrefInSmooth == False: + r = self.Wzz * ( self.mapping * m ) + return r.T * ( self.Wzz * self.mapping.deriv(m) ) + + @Utils.timeIt + def _evalSmooth2Deriv(self,m): + deriv = self._evalSmoothxxDeriv(m) + if self.regmesh.dim > 1: + deriv += self._evalSmoothyyDeriv(m) + if self.regmesh.dim > 2: + deriv += self._evalSmoothzzDeriv(m) + return deriv + + @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) ) + return self._evalSmall(m) + self._evalSmooth(m) + self._evalSmooth2(m) @Utils.timeIt def evalDeriv(self, m): @@ -560,88 +792,9 @@ class Tikhonov(BaseRegularization): R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} """ - return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + self._evalSmooth2Deriv(m) -class Simple(Tikhonov): - """ - Simple regularization that does not include length scales in the derivatives. - """ - - 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. - - def __init__(self, mesh, mapping=None, indActive=None, **kwargs): - BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) - - if isinstance(self.wght,float): - self.wght = np.ones(self.regmesh.nC) * self.wght - - @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 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 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 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,) - 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.Wsmall, self.Wsmooth) - self._W = sp.vstack(wlist) - return self._W - - @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) - class Sparse(Simple): @@ -716,25 +869,6 @@ class Sparse(Simple): 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 - @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) - - @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 def R(self, f_m , eps, exponent): From b4ab60c2601b6cb2dd64a65d58863f66f265f458 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Thu, 5 May 2016 11:55:56 -0700 Subject: [PATCH 09/27] Add model mapping to sparse regularization --- SimPEG/Regularization.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index ed1039ec..1b0d180f 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -663,7 +663,7 @@ class Sparse(Simple): self.Rs = Utils.speye(self.regmesh.nC) else: - f_m = self.curModel - self.reg.mref + f_m = self.mapping * (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 ) @@ -679,7 +679,7 @@ class Sparse(Simple): self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0]) else: - f_m = self.regmesh.cellDiffxStencil * self.curModel + f_m = self.regmesh.cellDiffxStencil * (self.mapping * self.curModel) self.rx = self.R( f_m , self.eps_q, self.norms[1]) self.Rx = Utils.sdiag( self.rx ) @@ -693,7 +693,7 @@ class Sparse(Simple): self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0]) else: - f_m = self.regmesh.cellDiffyStencil * self.curModel + f_m = self.regmesh.cellDiffyStencil * (self.mapping * self.curModel) self.ry = self.R( f_m , self.eps_q, self.norms[2]) self.Ry = Utils.sdiag( self.ry ) @@ -707,7 +707,7 @@ class Sparse(Simple): self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0]) else: - f_m = self.regmesh.cellDiffzStencil * self.curModel + f_m = self.regmesh.cellDiffzStencil * (self.mapping * self.curModel) self.rz = self.R( f_m , self.eps_q, self.norms[3]) self.Rz = Utils.sdiag( self.rz ) From fb5434695f13388790df41b6f4b293a5425d053d Mon Sep 17 00:00:00 2001 From: D Fournier Date: Tue, 10 May 2016 14:31:45 -0700 Subject: [PATCH 10/27] Alpha_s default to 1.0 --- SimPEG/Regularization.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 7f50ea59..8fdb9fd5 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -409,7 +409,7 @@ class Simple(BaseRegularization): """ mrefInSmooth = False #: SMOOTH and SMOOTH_MOD_DIF options - alpha_s = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Wsmall'], "Smallness weight") + 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") From eaa37f42e47b3e1d01ed10fd50d9b7cd8d29796f Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 10 May 2016 14:51:56 -0700 Subject: [PATCH 11/27] remove duplicate evalSmall --- SimPEG/Regularization.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 8fdb9fd5..1276cbb6 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -469,10 +469,6 @@ class Simple(BaseRegularization): self._W = sp.vstack(wlist) return self._W - @Utils.timeIt - def _evalSmall(self, m): - r = self.Wsmall * ( self.mapping * (m - self.mref) ) - return 0.5 * r.dot(r) @property def W(self): From 3f0c89f10b3783e63b14ccb530b565613d0cbff0 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 10 May 2016 14:53:43 -0700 Subject: [PATCH 12/27] remove extra Ws --- SimPEG/Regularization.py | 21 --------------------- 1 file changed, 21 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 1276cbb6..361a3206 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -470,27 +470,6 @@ class Simple(BaseRegularization): 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 - - @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 - - @Utils.timeIt def _evalSmall(self, m): r = self.Wsmall * ( self.mapping * (m - self.mref) ) From 2a802c1aa3ae90b120cb828194c36dfd2d57ccbf Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 10 May 2016 16:39:47 -0700 Subject: [PATCH 13/27] weights --> cell_weights, removed vol term from simple regularization --- SimPEG/Regularization.py | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 361a3206..5d40adf8 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -413,40 +413,40 @@ class Simple(BaseRegularization): 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. + cell_weights = 1. def __init__(self, mesh, mapping=None, indActive=None, **kwargs): BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) - if isinstance(self.wght,float): - self.wght = np.ones(self.regmesh.nC) * self.wght + if isinstance(self.cell_weights,float): + self.cell_weights = np.ones(self.regmesh.nC) * self.cell_weights @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) + self._Wsmall = Utils.sdiag((self.alpha_s*self.cell_weights)**0.5) return self._Wsmall @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 + self._Wx = Utils.sdiag((self.alpha_x * (self.regmesh.aveCC2Fx*self.cell_weights))**0.5)*self.regmesh.cellDiffxStencil return self._Wx @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 + self._Wy = Utils.sdiag((self.alpha_y * (self.regmesh.aveCC2Fy*self.cell_weights))**0.5)*self.regmesh.cellDiffyStencil return self._Wy @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 + self._Wz = Utils.sdiag((self.alpha_z * (self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.regmesh.cellDiffzStencil return self._Wz @property @@ -779,13 +779,13 @@ class Sparse(Simple): curModel = None # use a model to compute the weights gamma = 1. norms = [0., 2., 2., 2.] - wght = 1. + cell_weights = 1. def __init__(self, mesh, mapping=None, indActive=None, **kwargs): Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) - if isinstance(self.wght,float): - self.wght = np.ones(self.regmesh.nC) * self.wght + if isinstance(self.cell_weights,float): + self.cell_weights = np.ones(self.regmesh.nC) * self.cell_weights @property def Wsmall(self): @@ -799,7 +799,7 @@ class Sparse(Simple): #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.vol*self.alpha_s*self.gamma*self.wght)**0.5)*self.Rs + return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.cell_weights)**0.5)*self.Rs @property @@ -814,7 +814,7 @@ class Sparse(Simple): self.rx = self.R( f_m , self.eps_q, self.norms[1]) self.Rx = Utils.sdiag( self.rx ) - 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 + return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.cell_weights))**0.5)*self.Rx*self.regmesh.cellDiffxStencil @property def Wy(self): @@ -828,7 +828,7 @@ class Sparse(Simple): self.ry = self.R( f_m , self.eps_q, self.norms[2]) self.Ry = Utils.sdiag( self.ry ) - 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 + return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.cell_weights))**0.5)*self.Ry*self.regmesh.cellDiffyStencil @property def Wz(self): @@ -842,7 +842,7 @@ class Sparse(Simple): self.rz = self.R( f_m , self.eps_q, self.norms[3]) self.Rz = Utils.sdiag( self.rz ) - 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 + return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.Rz*self.regmesh.cellDiffzStencil def R(self, f_m , eps, exponent): From 955bd54019b9e6fccf8f766ce83d4fdde9ac61f0 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 10 May 2016 16:46:11 -0700 Subject: [PATCH 14/27] notation cleanup in Regularization --- SimPEG/Regularization.py | 32 +++++++++++++++++--------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 5d40adf8..16bb4b9f 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -1,4 +1,6 @@ -import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp +import Utils, Maps, Mesh +import numpy as np +import scipy.sparse as sp class RegularizationMesh(object): """ @@ -476,7 +478,7 @@ class Simple(BaseRegularization): return 0.5 * r.dot(r) @Utils.timeIt - def _evalSmallDeriv(self,m): + def _evalSmallDeriv(self, m): r = self.Wsmall * ( self.mapping * (m - self.mref) ) return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) ) @@ -514,7 +516,7 @@ class Simple(BaseRegularization): return phiSmooth @Utils.timeIt - def _evalSmoothxDeriv(self,m): + def _evalSmoothxDeriv(self, m): if self.mrefInSmooth == True: r = self.Wx * ( self.mapping * ( m - self.mref ) ) return r.T * ( self.Wx * self.mapping.deriv(m - self.mref) ) @@ -523,7 +525,7 @@ class Simple(BaseRegularization): return r.T * ( self.Wx * self.mapping.deriv(m) ) @Utils.timeIt - def _evalSmoothyDeriv(self,m): + def _evalSmoothyDeriv(self, m): if self.mrefInSmooth == True: r = self.Wy * ( self.mapping * ( m - self.mref ) ) return r.T * ( self.Wy * self.mapping.deriv(m - self.mref) ) @@ -532,7 +534,7 @@ class Simple(BaseRegularization): return r.T * ( self.Wy * self.mapping.deriv(m) ) @Utils.timeIt - def _evalSmoothzDeriv(self,m): + def _evalSmoothzDeriv(self, m): if self.mrefInSmooth == True: r = self.Wz * ( self.mapping * ( m - self.mref ) ) return r.T * ( self.Wz * self.mapping.deriv(m - self.mref) ) @@ -541,7 +543,7 @@ class Simple(BaseRegularization): return r.T * ( self.Wz * self.mapping.deriv(m) ) @Utils.timeIt - def _evalSmoothDeriv(self,m): + def _evalSmoothDeriv(self, m): deriv = self._evalSmoothxDeriv(m) if self.regmesh.dim > 1: deriv += self._evalSmoothyDeriv(m) @@ -711,7 +713,7 @@ class Tikhonov(Simple): return phiSmooth2 @Utils.timeIt - def _evalSmoothxxDeriv(self,m): + def _evalSmoothxxDeriv(self, m): if self.mrefInSmooth == True: r = self.Wxx * ( self.mapping * ( m - self.mref ) ) return r.T * ( self.Wxx * self.mapping.deriv(m - self.mref) ) @@ -720,7 +722,7 @@ class Tikhonov(Simple): return r.T * ( self.Wxx * self.mapping.deriv(m) ) @Utils.timeIt - def _evalSmoothyyDeriv(self,m): + def _evalSmoothyyDeriv(self, m): if self.mrefInSmooth == True: r = self.Wyy * ( self.mapping * ( m - self.mref ) ) return r.T * ( self.Wyy * self.mapping.deriv(m - self.mref) ) @@ -729,7 +731,7 @@ class Tikhonov(Simple): return r.T * ( self.Wyy * self.mapping.deriv(m) ) @Utils.timeIt - def _evalSmoothzzDeriv(self,m): + def _evalSmoothzzDeriv(self, m): if self.mrefInSmooth == True: r = self.Wzz * ( self.mapping * ( m - self.mref ) ) return r.T * ( self.Wzz * self.mapping.deriv(m - self.mref) ) @@ -738,7 +740,7 @@ class Tikhonov(Simple): return r.T * ( self.Wzz * self.mapping.deriv(m) ) @Utils.timeIt - def _evalSmooth2Deriv(self,m): + def _evalSmooth2Deriv(self, m): deriv = self._evalSmoothxxDeriv(m) if self.regmesh.dim > 1: deriv += self._evalSmoothyyDeriv(m) @@ -774,12 +776,12 @@ class Tikhonov(Simple): class Sparse(Simple): # set default values - eps_p = 1e-1 - eps_q = 1e-1 + eps_p = 1e-1 + eps_q = 1e-1 curModel = None # use a model to compute the weights - gamma = 1. - norms = [0., 2., 2., 2.] - cell_weights = 1. + gamma = 1. + norms = [0., 2., 2., 2.] + cell_weights = 1. def __init__(self, mesh, mapping=None, indActive=None, **kwargs): Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) From 7964ebce50f86b87f11873cfb0e40632f4eb930f Mon Sep 17 00:00:00 2001 From: D Fournier Date: Tue, 10 May 2016 17:20:46 -0700 Subject: [PATCH 15/27] Update directive to None the Wsmooth after iteration. --- SimPEG/Directives.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 4fc1ffc5..008cc57b 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -303,26 +303,27 @@ class Update_IRLS(InversionDirective): # Set the weighting matrix to None so that it is recomputed next time # it is called in the inversion self.reg._W = None + self.reg._Wsmooth = 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. From 90a3030796234798c3755745810415f52ef5a77e Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 10 May 2016 21:39:15 -0700 Subject: [PATCH 16/27] fixed 2 deriv --- SimPEG/Examples/Inversion_IRLS.py | 4 +- SimPEG/Regularization.py | 123 ++++++++++++++++++++++++++++-- 2 files changed, 119 insertions(+), 8 deletions(-) diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index 06ef4be4..c925cce2 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -84,12 +84,12 @@ def run(N=200, plotIt=True): #============================================================================== #reg.recModel = mrec - reg.wght = np.ones(mesh.nC) + # reg.cell_weight = 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 + reg.cell_weight = 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.) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 16bb4b9f..5947ed03 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -410,11 +410,11 @@ class Simple(BaseRegularization): Simple regularization that does not include length scales in the derivatives. """ - mrefInSmooth = False #: SMOOTH and SMOOTH_MOD_DIF options + mrefInSmooth = False #: include mref in the smoothness? 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") + 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") cell_weights = 1. def __init__(self, mesh, mapping=None, indActive=None, **kwargs): @@ -454,6 +454,8 @@ class Simple(BaseRegularization): @property def Wsmooth(self): """Full smoothness regularization matrix W""" + print 'wtf why are we using Wsmooth' + raise NotImplementedError if getattr(self, '_Wsmooth', None) is None: wlist = (self.Wx,) if self.regmesh.dim > 1: @@ -482,6 +484,13 @@ class Simple(BaseRegularization): r = self.Wsmall * ( self.mapping * (m - self.mref) ) return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) ) + @Utils.timeIt + def _evalSmall2Deriv(self, m, v = None): + rDeriv = self.Wsmall * ( self.mapping.deriv(m - self.mref) ) + if v is not None: + return rDeriv.T * (rDeriv * v) + return rDeriv.T * rDeriv + @Utils.timeIt def _evalSmoothx(self, m): if self.mrefInSmooth == True: @@ -524,6 +533,17 @@ class Simple(BaseRegularization): r = self.Wx * ( self.mapping * m ) return r.T * ( self.Wx * self.mapping.deriv(m) ) + @Utils.timeIt + def _evalSmoothx2Deriv(self, m, v=None): + if self.mrefInSmooth == True: + rDeriv = self.Wx * ( self.mapping.deriv( m - self.mref ) ) + elif self.mrefInSmooth == False: + rDeriv = self.Wx * ( self.mapping.deriv(m) ) + + if v is not None: + return rDeriv.T * ( rDeriv * v ) + return rDeriv.T * rDeriv + @Utils.timeIt def _evalSmoothyDeriv(self, m): if self.mrefInSmooth == True: @@ -533,6 +553,17 @@ class Simple(BaseRegularization): r = self.Wy * ( self.mapping * m ) return r.T * ( self.Wy * self.mapping.deriv(m) ) + @Utils.timeIt + def _evalSmoothy2Deriv(self, m, v=None): + if self.mrefInSmooth == True: + rDeriv = self.Wy * ( self.mapping.deriv( m - self.mref ) ) + elif self.mrefInSmooth == False: + rDeriv = self.Wy * ( self.mapping.deriv(m) ) + + if v is not None: + return rDeriv.T * ( rDeriv * v ) + return rDeriv.T * rDeriv + @Utils.timeIt def _evalSmoothzDeriv(self, m): if self.mrefInSmooth == True: @@ -542,6 +573,17 @@ class Simple(BaseRegularization): r = self.Wz * ( self.mapping * m ) return r.T * ( self.Wz * self.mapping.deriv(m) ) + @Utils.timeIt + def _evalSmoothz2Deriv(self, m, v=None): + if self.mrefInSmooth == True: + rDeriv = self.Wz * ( self.mapping.deriv( m - self.mref ) ) + elif self.mrefInSmooth == False: + rDeriv = self.Wz * ( self.mapping.deriv(m) ) + + if v is not None: + return rDeriv.T * ( rDeriv * v ) + return rDeriv.T * rDeriv + @Utils.timeIt def _evalSmoothDeriv(self, m): deriv = self._evalSmoothxDeriv(m) @@ -551,6 +593,15 @@ class Simple(BaseRegularization): deriv += self._evalSmoothzDeriv(m) return deriv + @Utils.timeIt + def _evalSmooth2Deriv(self, m, v=None): + deriv = self._evalSmoothx2Deriv(m, v) + if self.regmesh.dim > 1: + deriv += self._evalSmoothy2Deriv(m, v) + if self.regmesh.dim > 2: + deriv += self._evalSmoothz2Deriv(m, v) + return deriv + @Utils.timeIt def eval(self, m): @@ -574,6 +625,10 @@ class Simple(BaseRegularization): """ return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + @Utils.timeIt + def eval2Deriv(self, m, v=None): + return self._evalSmall2Deriv(m, v) + self._evalSmooth2Deriv(m, v) + class Tikhonov(Simple): @@ -740,7 +795,37 @@ class Tikhonov(Simple): return r.T * ( self.Wzz * self.mapping.deriv(m) ) @Utils.timeIt - def _evalSmooth2Deriv(self, m): + def _evalSmoothxx2Deriv(self, m, v=None): + if self.mrefInSmooth == True: + rDeriv = self.Wxx * ( self.mapping.deriv( m - self.mref ) ) + elif self.mrefInSmooth == False: + rDeriv = self.Wxx * self.mapping.deriv(m) + if v is not None: + return rDeriv.T * (rDeriv * v) + return rDeriv.T * rDeriv + + @Utils.timeIt + def _evalSmoothyy2Deriv(self, m, v=None): + if self.mrefInSmooth == True: + rDeriv = self.Wyy * ( self.mapping.deriv( m - self.mref ) ) + elif self.mrefInSmooth == False: + rDeriv = self.Wyy * self.mapping.deriv(m) + if v is not None: + return rDeriv.T * (rDeriv * v) + return rDeriv.T * rDeriv + + @Utils.timeIt + def _evalSmoothzz2Deriv(self, m, v=None): + if self.mrefInSmooth == True: + rDeriv = self.Wzz * ( self.mapping.deriv( m - self.mref ) ) + elif self.mrefInSmooth == False: + rDeriv = self.Wzz * self.mapping.deriv(m) + if v is not None: + return rDeriv.T * (rDeriv * v) + return rDeriv.T * rDeriv + + @Utils.timeIt + def _evalSmoothDeriv2(self, m): deriv = self._evalSmoothxxDeriv(m) if self.regmesh.dim > 1: deriv += self._evalSmoothyyDeriv(m) @@ -748,6 +833,15 @@ class Tikhonov(Simple): deriv += self._evalSmoothzzDeriv(m) return deriv + @Utils.timeIt + def _evalSmooth2Deriv2(self, m, v=None): + deriv = self._evalSmoothxx2Deriv(m, v) + if self.regmesh.dim > 1: + deriv += self._evalSmoothyy2Deriv(m, v) + if self.regmesh.dim > 2: + deriv += self._evalSmoothzz2Deriv(m, v) + return deriv + @Utils.timeIt def eval(self, m): @@ -769,7 +863,24 @@ class Tikhonov(Simple): R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} """ - return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + self._evalSmooth2Deriv(m) + return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + self._evalSmoothDeriv2(m) + + def eval2Deriv(self, m): + """ + The regularization is: + + .. math:: + + R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})} + + So the derivative is straight forward: + + .. math:: + + R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} + + """ + return self._evalSmall2Deriv(m) + self._evalSmooth2Deriv(m) + self._evalSmooth2Deriv2(m) From 3dd9ecc9cd7f95b596d9027595013f9ca4647035 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 10 May 2016 22:10:19 -0700 Subject: [PATCH 17/27] fix tikhonov 2Deriv --- SimPEG/Regularization.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 5947ed03..8f1b2d2c 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -865,7 +865,7 @@ class Tikhonov(Simple): """ return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + self._evalSmoothDeriv2(m) - def eval2Deriv(self, m): + def eval2Deriv(self, m, v=None): """ The regularization is: @@ -880,7 +880,7 @@ class Tikhonov(Simple): R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} """ - return self._evalSmall2Deriv(m) + self._evalSmooth2Deriv(m) + self._evalSmooth2Deriv2(m) + return self._evalSmall2Deriv(m, v) + self._evalSmooth2Deriv(m, v) + self._evalSmooth2Deriv2(m, v) From e10d6878fb86b129953e68c3749307321a7d8eea Mon Sep 17 00:00:00 2001 From: D Fournier Date: Wed, 11 May 2016 07:58:06 -0700 Subject: [PATCH 18/27] Remove Wsmooth from def W and replace by parts --- SimPEG/Regularization.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 8f1b2d2c..3263ba9a 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -469,7 +469,11 @@ class Simple(BaseRegularization): def W(self): """Full regularization matrix W""" if getattr(self, '_W', None) is None: - wlist = (self.Wsmall, self.Wsmooth) + wlist = (self.Wsmall, self.Wx) + if self.regmesh.dim > 1: + wlist += (self.Wy,) + if self.regmesh.dim > 2: + wlist += (self.Wz,) self._W = sp.vstack(wlist) return self._W From cd2360b8153a9b4f343cfdb81d084e5f36c49132 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Wed, 11 May 2016 23:04:14 -0700 Subject: [PATCH 19/27] Stash the regularization between each beta --- SimPEG/Directives.py | 24 ++++++++++++-- SimPEG/Regularization.py | 67 ++++++++++++++++++++++------------------ 2 files changed, 58 insertions(+), 33 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 008cc57b..a9b74276 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -271,6 +271,14 @@ class Update_IRLS(InversionDirective): self.reg.curModel = self.invProb.curModel self.reg.gamma = self.gamma + print "Initial gamma ", np.linalg.norm(self.reg.gamma) + # Reset the regularization matrices so that it is + # recalculated with new gamma + self.reg._Wsmall = None + self.reg._Wx = None + self.reg._Wy = None + self.reg._Wz = None + self.reg._W = None if getattr(self, 'phi_d_last', None) is None: self.phi_d_last = self.invProb.phi_d @@ -294,16 +302,26 @@ class Update_IRLS(InversionDirective): # Temporarely set gamma to 1. to get raw phi_m self.reg.gamma = 1. + # Reset the regularization matrices so that it is + # recalculated for current model + self.reg._Wsmall = None + self.reg._Wx = None + self.reg._Wy = None + self.reg._Wz = None + self.reg._W = None + # 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 + # Reset the regularization matrices again for new gamma + self.reg._Wsmall = None + self.reg._Wx = None + self.reg._Wy = None + self.reg._Wz = None self.reg._W = None - self.reg._Wsmooth = None class Update_lin_PreCond(InversionDirective): """ diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 3263ba9a..5dc5aba0 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -468,6 +468,7 @@ class Simple(BaseRegularization): @property def W(self): """Full regularization matrix W""" + print 'wtf why are we using W' if getattr(self, '_W', None) is None: wlist = (self.Wsmall, self.Wx) if self.regmesh.dim > 1: @@ -907,60 +908,66 @@ class Sparse(Simple): @property def Wsmall(self): """Regularization matrix Wsmall""" - if getattr(self, 'curModel', None) is None: - self.Rs = Utils.speye(self.regmesh.nC) + if getattr(self,'_Wsmall', None) is None: + if getattr(self, 'curModel', None) is None: + self.Rs = Utils.speye(self.regmesh.nC) - else: - f_m = self.mapping * (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 ) + else: + f_m = self.mapping * (self.curModel - self.reg.mref) + self.rs = self.R(f_m , self.eps_p, self.norms[0]) + self.Rs = Utils.sdiag( self.rs ) - return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.cell_weights)**0.5)*self.Rs + self._Wsmall = Utils.sdiag((self.alpha_s*self.gamma*self.cell_weights)**0.5)*self.Rs + return self._Wsmall @property def Wx(self): """Regularization matrix Wx""" + if getattr(self,'_Wx', None) is None: + if getattr(self, 'curModel', None) is None: + self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0]) - if getattr(self, 'curModel', None) is None: - self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0]) + else: + f_m = self.regmesh.cellDiffxStencil * (self.mapping * self.curModel) + self.rx = self.R( f_m , self.eps_q, self.norms[1]) + self.Rx = Utils.sdiag( self.rx ) - else: - f_m = self.regmesh.cellDiffxStencil * (self.mapping * self.curModel) - self.rx = self.R( f_m , self.eps_q, self.norms[1]) - self.Rx = Utils.sdiag( self.rx ) + self._Wx = Utils.sdiag(( self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.cell_weights))**0.5)*self.Rx*self.regmesh.cellDiffxStencil - return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.cell_weights))**0.5)*self.Rx*self.regmesh.cellDiffxStencil + return self._Wx @property def Wy(self): """Regularization matrix Wy""" + if getattr(self,'_Wy', None) is None: + if getattr(self, 'curModel', None) is None: + self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0]) - if getattr(self, 'curModel', None) is None: - self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0]) + else: + f_m = self.regmesh.cellDiffyStencil * (self.mapping * self.curModel) + self.ry = self.R( f_m , self.eps_q, self.norms[2]) + self.Ry = Utils.sdiag( self.ry ) - else: - f_m = self.regmesh.cellDiffyStencil * (self.mapping * self.curModel) - self.ry = self.R( f_m , self.eps_q, self.norms[2]) - self.Ry = Utils.sdiag( self.ry ) + self._Wy = Utils.sdiag((self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.cell_weights))**0.5)*self.Ry*self.regmesh.cellDiffyStencil - return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.cell_weights))**0.5)*self.Ry*self.regmesh.cellDiffyStencil + return self._Wy @property def Wz(self): """Regularization matrix Wz""" + if getattr(self,'_Wz', None) is None: + if getattr(self, 'curModel', None) is None: + self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0]) - if getattr(self, 'curModel', None) is None: - self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0]) + else: + f_m = self.regmesh.cellDiffzStencil * (self.mapping * self.curModel) + self.rz = self.R( f_m , self.eps_q, self.norms[3]) + self.Rz = Utils.sdiag( self.rz ) - else: - f_m = self.regmesh.cellDiffzStencil * (self.mapping * self.curModel) - self.rz = self.R( f_m , self.eps_q, self.norms[3]) - self.Rz = Utils.sdiag( self.rz ) - - return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.Rz*self.regmesh.cellDiffzStencil + self._Wz = Utils.sdiag((self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.Rz*self.regmesh.cellDiffzStencil + return self._Wz def R(self, f_m , eps, exponent): From 3cc46131a3138096441fb4c423c5f4f21543d2b4 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Thu, 12 May 2016 08:31:01 -0700 Subject: [PATCH 20/27] Temporary change ... comment out W and Wsmooth --- SimPEG/Directives.py | 6 ++--- SimPEG/Regularization.py | 52 ++++++++++++++++++++-------------------- 2 files changed, 28 insertions(+), 30 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index a9b74276..6b15be9f 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -278,7 +278,6 @@ class Update_IRLS(InversionDirective): self.reg._Wx = None self.reg._Wy = None self.reg._Wz = None - self.reg._W = None if getattr(self, 'phi_d_last', None) is None: self.phi_d_last = self.invProb.phi_d @@ -308,20 +307,19 @@ class Update_IRLS(InversionDirective): self.reg._Wx = None self.reg._Wy = None self.reg._Wz = None - self.reg._W = None # 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 - + print "New gamma ", np.linalg.norm(self.reg.gamma) + # Reset the regularization matrices again for new gamma self.reg._Wsmall = None self.reg._Wx = None self.reg._Wy = None self.reg._Wz = None - self.reg._W = None class Update_lin_PreCond(InversionDirective): """ diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 5dc5aba0..422c8f9c 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -451,32 +451,32 @@ class Simple(BaseRegularization): self._Wz = Utils.sdiag((self.alpha_z * (self.regmesh.aveCC2Fz*self.cell_weights))**0.5)*self.regmesh.cellDiffzStencil return self._Wz - @property - def Wsmooth(self): - """Full smoothness regularization matrix W""" - print 'wtf why are we using Wsmooth' - raise NotImplementedError - 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""" - print 'wtf why are we using W' - if getattr(self, '_W', None) is None: - wlist = (self.Wsmall, self.Wx) - if self.regmesh.dim > 1: - wlist += (self.Wy,) - if self.regmesh.dim > 2: - wlist += (self.Wz,) - self._W = sp.vstack(wlist) - return self._W +# @property +# def Wsmooth(self): +# """Full smoothness regularization matrix W""" +# print 'wtf why are we using Wsmooth' +# raise NotImplementedError +# 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""" +# print 'wtf why are we using W' +# if getattr(self, '_W', None) is None: +# wlist = (self.Wsmall, self.Wx) +# if self.regmesh.dim > 1: +# wlist += (self.Wy,) +# if self.regmesh.dim > 2: +# wlist += (self.Wz,) +# self._W = sp.vstack(wlist) +# return self._W @Utils.timeIt From fd3bde787f3a9bc570ad8386e36eb16a7d3da1d3 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Thu, 12 May 2016 14:58:16 -0700 Subject: [PATCH 21/27] Propose change to the Projected_GNCG solver. Add inner GN iterations. Nice improvement to the convergence of IRLS --- SimPEG/Directives.py | 49 ++++++++++++++++++---- SimPEG/Examples/Inversion_IRLS.py | 20 ++++----- SimPEG/Optimization.py | 70 +++++++++++++++++++++++++++++++ 3 files changed, 118 insertions(+), 21 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 6b15be9f..f5c25249 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -144,6 +144,35 @@ class BetaSchedule(InversionDirective): if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter self.invProb.beta /= self.coolingFactor +#class BetaSchedule_PGN_CG(InversionDirective): +# """BetaSchedule""" +# +# coolingFactor = 5. +# coolingRate = 1 +# GN_step_last = None +# GN_step_c = None +# +# def endIter(self): +# +# """ Compute the change in GN step, and proceed with cooling if below tol""" +# if self.opt.iter == 1: +# self.GN_step_last = np.linalg.norm(self.opt.xc - self.opt.x_last) +# d_GN_step = 1. +# +# else: +# self.GN_step_c = np.linalg.norm(self.opt.xc - self.opt.x_last) +# d_GN_step = self.GN_step_c / self.GN_step_last +# +# # Re-initiate last GN step +# self.GN_step_last = self.GN_step_c +# +# print "GN_step_last: ", self.GN_step_last +# print "d_GN_step: ", d_GN_step +# +# if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0: +# if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter +# self.invProb.beta /= self.coolingFactor + class TargetMisfit(InversionDirective): @property @@ -265,13 +294,16 @@ class Update_IRLS(InversionDirective): if getattr(self, 'phi_m_last', None) is not None: self.reg.curModel = self.invProb.curModel + self.reg._Wsmall = None + self.reg._Wx = None + self.reg._Wy = None + self.reg._Wz = None 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 - print "Initial gamma ", np.linalg.norm(self.reg.gamma) # Reset the regularization matrices so that it is # recalculated with new gamma self.reg._Wsmall = None @@ -295,12 +327,6 @@ class Update_IRLS(InversionDirective): # 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. - # Reset the regularization matrices so that it is # recalculated for current model self.reg._Wsmall = None @@ -308,13 +334,18 @@ class Update_IRLS(InversionDirective): self.reg._Wy = None self.reg._Wz = None + # 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 - print "New gamma ", np.linalg.norm(self.reg.gamma) - + # Reset the regularization matrices again for new gamma self.reg._Wsmall = None self.reg._Wx = None diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index c925cce2..6551bf21 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -18,6 +18,8 @@ def run(N=200, plotIt=True): mesh = Mesh.TensorMesh([N]) m0 = np.ones(mesh.nC) * 1e-4 + mref = np.zeros(mesh.nC) + nk = 10 jk = np.linspace(1.,nk,nk) p = -2. @@ -51,12 +53,13 @@ def run(N=200, plotIt=True): wr = ( wr/np.max(wr) ) reg = Regularization.Simple(mesh) - reg.wght = wr + reg.mref = mref + reg.cell_weights = wr dmis = DataMisfit.l2_DataMisfit(survey) dmis.Wd = 1./wd - opt = Optimization.ProjectedGNCG(maxIter=30,lower=-2.,upper=2., maxIterCG= 20, tolCG = 1e-4) + opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, maxIterGN=1, tolCG = 1e-4) invProb = InvProblem.BaseInvProblem(dmis, reg, opt) invProb.curModel = m0 @@ -76,22 +79,15 @@ def run(N=200, plotIt=True): phid = invProb.phi_d reg = Regularization.Sparse(mesh) + reg.mref = mref + reg.cell_weights = wr -#============================================================================== -# 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.cell_weight = 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.cell_weight = wr - opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 20, tolCG = 1e-3) + opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, 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/Optimization.py b/SimPEG/Optimization.py index 77704733..95f53320 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -893,6 +893,10 @@ class ProjectedGNCG(BFGS, Minimize, Remember): lower = -np.inf upper = np.inf + # Variables to control inner GN iterations + rdm_tol = 1e-2 # Tolerance for largest change in step (Default 1%) + maxIterGN = 3 # Maximum number of GN inner iterations + def _startup(self, x0): # ensure bound vectors are the same size as the model if type(self.lower) is not np.ndarray: @@ -938,6 +942,72 @@ class ProjectedGNCG(BFGS, Minimize, Remember): def approxHinv(self, value): self._approxHinv = value + @Utils.timeIt + def minimize(self, evalFunction, x0): + """minimize(evalFunction, x0) + + Minimizes the function (evalFunction) starting at the location x0. + + :param def evalFunction: function handle that evaluates: f, g, H = F(x) + :param numpy.ndarray x0: starting location + :rtype: numpy.ndarray + :return: x, the last iterate of the optimization algorithm + + The GN newton steps are repeated until it the maxIterGN is reached or the + relative step change falls below some tolrerance. + + """ + self.evalFunction = evalFunction + self.startup(x0) + self.printInit() + + + while True: + self.doStartIteration() + self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True) + self.printIter() + if self.stoppingCriteria(): break + + # Inner GN iterations, stop on maximum number of iterations + # or on tolerance for step length change + GN_count = 0 + dm0 = None # Initial GN step length + dmc = None # Current GN step length + rdm = 1. # Relative change in step length + + while rdm > self.rdm_tol and GN_count < self.maxIterGN: + + GN_count += 1 + self.searchDirection = self.findSearchDirection() + p = self.scaleSearchDirection(self.searchDirection) + xt, passLS = self.modifySearchDirection(p) + if not passLS: + xt, caught = self.modifySearchDirectionBreak(p) + if not caught: return self.xc + + if GN_count == 1: + dm0 = np.linalg.norm(self.xc - xt) + dmc = dm0 + + else: + dmc = np.linalg.norm(self.xc - xt) + + rdm = dmc / dm0 + self.xc = xt # Update current model + + # Form system for next iteration + self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True) + + print "GN iter: %i,\t dm: %8.5e,\t rdm: %8.5e"% (GN_count, dmc, rdm) + + self.doEndIteration(xt) + if self.stopNextIteration: break + + self.printDone() + self.finish() + + return self.xc + @Utils.timeIt def findSearchDirection(self): From 022e1f7660dd0b20e44b00ebb919a0a2cf77b036 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Fri, 27 May 2016 13:11:31 -0700 Subject: [PATCH 22/27] Update IRLS directive to allow multiple GN iterations. Remove modifications to the ProjGN solver. Update IRLS example. --- SimPEG/Directives.py | 156 ++++++++++++++++-------------- SimPEG/Examples/Inversion_IRLS.py | 14 +-- SimPEG/Optimization.py | 72 +------------- SimPEG/Regularization.py | 38 ++++++-- 4 files changed, 121 insertions(+), 159 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index f5c25249..b84d78d1 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -144,34 +144,6 @@ class BetaSchedule(InversionDirective): if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter self.invProb.beta /= self.coolingFactor -#class BetaSchedule_PGN_CG(InversionDirective): -# """BetaSchedule""" -# -# coolingFactor = 5. -# coolingRate = 1 -# GN_step_last = None -# GN_step_c = None -# -# def endIter(self): -# -# """ Compute the change in GN step, and proceed with cooling if below tol""" -# if self.opt.iter == 1: -# self.GN_step_last = np.linalg.norm(self.opt.xc - self.opt.x_last) -# d_GN_step = 1. -# -# else: -# self.GN_step_c = np.linalg.norm(self.opt.xc - self.opt.x_last) -# d_GN_step = self.GN_step_c / self.GN_step_last -# -# # Re-initiate last GN step -# self.GN_step_last = self.GN_step_c -# -# print "GN_step_last: ", self.GN_step_last -# print "d_GN_step: ", d_GN_step -# -# if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0: -# if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter -# self.invProb.beta /= self.coolingFactor class TargetMisfit(InversionDirective): @@ -286,10 +258,16 @@ class Update_IRLS(InversionDirective): gamma = None phi_m_last = None phi_d_last = None - - + f_old = None + f_min_change = 1e-1 + coolingRate = 3 + maxIRLSiter = 10 + + def initialize(self): - + + self.IRLSiter = 0 + # Scale the regularization for changes in norm if getattr(self, 'phi_m_last', None) is not None: @@ -310,48 +288,75 @@ class Update_IRLS(InversionDirective): self.reg._Wx = None self.reg._Wy = None self.reg._Wz = None - + if getattr(self, 'phi_d_last', None) is None: self.phi_d_last = self.invProb.phi_d + if getattr(self, 'f_last', None) is None: + self.f_old = self.invProb.evalFunction(self.reg.curModel, return_g=False, return_H=False) + print self.f_old def endIter(self): - # Cool the threshold parameter if required - if getattr(self, 'factor', None) is not None: - eps = self.reg.eps / self.factor + + + # Only update after GN iterations + if self.opt.iter % self.coolingRate == 0: + + self.IRLSiter += 1 - 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 - - # Reset the regularization matrices so that it is - # recalculated for current model - self.reg._Wsmall = None - self.reg._Wx = None - self.reg._Wy = None - self.reg._Wz = None - - # 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 - - # Reset the regularization matrices again for new gamma - self.reg._Wsmall = None - self.reg._Wx = None - self.reg._Wy = None - self.reg._Wz = None + self.f_change = np.abs(self.f_old - self.opt.f_last) / self.f_old + + print "Function decrease" + str(self.f_change) + + # Check for maximum number of IRLS cycles + if self.IRLSiter == self.maxIRLSiter: + self.opt.stopNextIteration = True + return + + # Check if the function has changed enough + if self.f_change < self.f_min_change: + self.opt.stopNextIteration = True + return + else: + self.f_old = self.opt.f_last + + # 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 + + # Reset the regularization matrices so that it is + # recalculated for current model + self.reg._Wsmall = None + self.reg._Wx = None + self.reg._Wy = None + self.reg._Wz = None + + # 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 + + # Reset the regularization matrices again for new gamma + self.reg._Wsmall = None + self.reg._Wx = None + self.reg._Wy = None + self.reg._Wz = None + # Compute the change in class Update_lin_PreCond(InversionDirective): """ Create a Jacobi preconditioner for the linear problem @@ -411,11 +416,14 @@ class Scale_Beta(InversionDirective): update is done only if the misfit is outside some threshold bounds. """ tol = 0.05 - + coolingRate=5 + 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 + + # Only update after GN iterations + if self.opt.iter % self.coolingRate == 0: + # 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/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index 6551bf21..d1da048f 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -20,7 +20,7 @@ def run(N=200, plotIt=True): m0 = np.ones(mesh.nC) * 1e-4 mref = np.zeros(mesh.nC) - nk = 10 + nk = 15 jk = np.linspace(1.,nk,nk) p = -2. q = 1. @@ -59,7 +59,7 @@ def run(N=200, plotIt=True): dmis = DataMisfit.l2_DataMisfit(survey) dmis.Wd = 1./wd - opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, maxIterGN=1, tolCG = 1e-4) + opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, tolCG = 1e-4) invProb = InvProblem.BaseInvProblem(dmis, reg, opt) invProb.curModel = m0 @@ -83,18 +83,18 @@ def run(N=200, plotIt=True): reg.cell_weights = wr reg.mref = np.zeros(mesh.nC) - reg.eps_p = 5e-2 + reg.eps_p = 1e-3 reg.eps_q = 1e-2 reg.norms = [0., 0., 2., 2.] - opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3) + opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3) invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.) beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1) - #betaest = Directives.BetaEstimate_ByEig() + update_beta = Directives.Scale_Beta(tol = 0.05, coolingRate=5) target = Directives.TargetMisfit() - IRLS =Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid ) + IRLS = Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid, coolingRate=5 ) - inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS]) + inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS,update_beta]) m0 = mrec diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py index 95f53320..05214dd5 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -893,10 +893,6 @@ class ProjectedGNCG(BFGS, Minimize, Remember): lower = -np.inf upper = np.inf - # Variables to control inner GN iterations - rdm_tol = 1e-2 # Tolerance for largest change in step (Default 1%) - maxIterGN = 3 # Maximum number of GN inner iterations - def _startup(self, x0): # ensure bound vectors are the same size as the model if type(self.lower) is not np.ndarray: @@ -942,72 +938,6 @@ class ProjectedGNCG(BFGS, Minimize, Remember): def approxHinv(self, value): self._approxHinv = value - @Utils.timeIt - def minimize(self, evalFunction, x0): - """minimize(evalFunction, x0) - - Minimizes the function (evalFunction) starting at the location x0. - - :param def evalFunction: function handle that evaluates: f, g, H = F(x) - :param numpy.ndarray x0: starting location - :rtype: numpy.ndarray - :return: x, the last iterate of the optimization algorithm - - The GN newton steps are repeated until it the maxIterGN is reached or the - relative step change falls below some tolrerance. - - """ - self.evalFunction = evalFunction - self.startup(x0) - self.printInit() - - - while True: - self.doStartIteration() - self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True) - self.printIter() - if self.stoppingCriteria(): break - - # Inner GN iterations, stop on maximum number of iterations - # or on tolerance for step length change - GN_count = 0 - dm0 = None # Initial GN step length - dmc = None # Current GN step length - rdm = 1. # Relative change in step length - - while rdm > self.rdm_tol and GN_count < self.maxIterGN: - - GN_count += 1 - self.searchDirection = self.findSearchDirection() - p = self.scaleSearchDirection(self.searchDirection) - xt, passLS = self.modifySearchDirection(p) - if not passLS: - xt, caught = self.modifySearchDirectionBreak(p) - if not caught: return self.xc - - if GN_count == 1: - dm0 = np.linalg.norm(self.xc - xt) - dmc = dm0 - - else: - dmc = np.linalg.norm(self.xc - xt) - - rdm = dmc / dm0 - self.xc = xt # Update current model - - # Form system for next iteration - self.f, self.g, self.H = evalFunction(self.xc, return_g=True, return_H=True) - - print "GN iter: %i,\t dm: %8.5e,\t rdm: %8.5e"% (GN_count, dmc, rdm) - - self.doEndIteration(xt) - if self.stopNextIteration: break - - self.printDone() - self.finish() - - return self.xc - @Utils.timeIt def findSearchDirection(self): @@ -1078,4 +1008,4 @@ class ProjectedGNCG(BFGS, Minimize, Remember): indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0)) delx[indx] = 0. - return delx + return delx \ No newline at end of file diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 422c8f9c..c3b1c9e5 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -890,14 +890,37 @@ class Tikhonov(Simple): class Sparse(Simple): - + """ + The regularization is: + + .. math:: + + R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top R^\\top R W(m-m_\\text{ref})} + + where the IRLS weight + + .. math:: + + R = \eta TO FINISH LATER!!! + + So the derivative is straight forward: + + .. math:: + + R(m) = \mathbf{W^\\top R^\\top R W (m-m_\\text{ref})} + + The IRLS weights are recomputed after each beta solves. + It is strongly recommended to do a few Gauss-Newton iterations + before updating. + """ + # 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.] - cell_weights = 1. + eps_p = 1e-1 # Threshold value for the model norm + eps_q = 1e-1 # Threshold value for the model gradient norm + curModel = None # Requires model to compute the weights + gamma = 1. # Model norm scaling to smooth out convergence + norms = [0., 2., 2., 2.] # Values for norm on (m, dmdx, dmdy, dmdz) + cell_weights = 1. # Consider overwriting with sensitivity weights def __init__(self, mesh, mapping=None, indActive=None, **kwargs): Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) @@ -971,6 +994,7 @@ class Sparse(Simple): def R(self, f_m , eps, exponent): + # Eta scaling is important for mix-norms...do not mess with it eta = (eps**(1.-exponent/2.))**0.5 r = eta / (f_m**2.+ eps**2.)**((1.-exponent/2.)/2.) From 3b4bec9c0b40b0ce949cea8422adef1f29b81bb7 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Sat, 28 May 2016 11:27:09 -0700 Subject: [PATCH 23/27] Refactor IRLS iterations, full solves from l2->lp Adapt Example --- SimPEG/Directives.py | 206 ++++++++++++++++++------------ SimPEG/Examples/Inversion_IRLS.py | 66 +++++----- SimPEG/Regularization.py | 1 + 3 files changed, 158 insertions(+), 115 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index b84d78d1..8727d925 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -144,7 +144,7 @@ class BetaSchedule(InversionDirective): if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter self.invProb.beta /= self.coolingFactor - + class TargetMisfit(InversionDirective): @property @@ -238,12 +238,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 UpdateReferenceModel(Parameter): - -# mref0 = None - -# def nextIter(self): # mref = getattr(self, 'm_prev', None) # if mref is None: # if self.debug: print 'UpdateReferenceModel is using mref0' @@ -254,109 +248,165 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration): class Update_IRLS(InversionDirective): eps_min = None + eps_p = None + eps_q = None + norms = [2.,2.,2.,2.] factor = None gamma = None phi_m_last = None phi_d_last = None f_old = None - f_min_change = 1e-1 - coolingRate = 3 - maxIRLSiter = 10 - + f_min_change = 1e-2 + beta_tol = 5e-2 + # Solving parameter for IRLS (mode:2) + IRLSiter = 0 + minGNiter = 6 + maxIRLSiter = 10 + iterStart = 0 + + # Beta schedule + coolingFactor = 2. + coolingRate = 1 + + mode = 1 + + @property + def target(self): + if getattr(self, '_target', None) is None: + self._target = self.survey.nD + return self._target + @target.setter + def target(self, val): + self._target = val + def initialize(self): - - self.IRLSiter = 0 - - # 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._Wsmall = None - self.reg._Wx = None - self.reg._Wy = None - self.reg._Wz = None - self.reg.gamma = 1. - phim_new = self.reg.eval(self.invProb.curModel) - self.gamma = self.phi_m_last / phim_new + if self.mode == 1: + self.reg.norms = [2., 2., 2., 2.] +# # 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._Wsmall = None +# self.reg._Wx = None +# self.reg._Wy = None +# self.reg._Wz = None +# 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 +# # Reset the regularization matrices so that it is +# # recalculated with new gamma +# self.reg._Wsmall = None +# self.reg._Wx = None +# self.reg._Wy = None +# self.reg._Wz = None - self.reg.curModel = self.invProb.curModel - self.reg.gamma = self.gamma - # Reset the regularization matrices so that it is - # recalculated with new gamma - self.reg._Wsmall = None - self.reg._Wx = None - self.reg._Wy = None - self.reg._Wz = None - - if getattr(self, 'phi_d_last', None) is None: - self.phi_d_last = self.invProb.phi_d - - if getattr(self, 'f_last', None) is None: - self.f_old = self.invProb.evalFunction(self.reg.curModel, return_g=False, return_H=False) - print self.f_old - def endIter(self): - - - # Only update after GN iterations - if self.opt.iter % self.coolingRate == 0: +# if getattr(self, 'phi_d_last', None) is None: +# self.phi_d_last = self.invProb.phi_d +# +# if getattr(self, 'f_last', None) is None: +# self.f_old = self.invProb.evalFunction(self.reg.curModel, return_g=False, return_H=False) +# print self.f_old + def endIter(self): + + # After reaching target misfit with l2-norm, switch to IRLS (mode:2) + if self.invProb.phi_d < self.target and self.mode == 1: + print "Convergence with smooth l2-norm regularization: Start IRLS steps..." + + self.mode = 2 + print self.eps_p, self.eps_q, self.norms + self.reg.eps_p = self.eps_p + self.reg.eps_q = self.eps_q + self.reg.norms = self.norms + self.coolingFactor = 1. + self.coolingRate = 1 + self.iterStart = self.opt.iter + self.phi_d_last = self.invProb.phi_d + self.phi_m_last = self.invProb.phi_m_last + + self.reg.l2model = self.invProb.curModel + self.reg.curModel = self.invProb.curModel + + if getattr(self, 'f_old', None) is None: + self.f_old = self.reg.eval(self.invProb.curModel)#self.invProb.evalFunction(self.invProb.curModel, return_g=False, return_H=False) + + if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0: + if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter + self.invProb.beta /= self.coolingFactor + + + # Only update after GN iterations + if (self.opt.iter-self.iterStart) % self.minGNiter == 0 and self.mode==2: + self.IRLSiter += 1 - self.f_change = np.abs(self.f_old - self.opt.f_last) / self.f_old - - print "Function decrease" + str(self.f_change) - + phim_new = self.reg.eval(self.invProb.curModel) + self.f_change = np.abs(self.f_old - phim_new) / self.f_old + + print "Regularization decrease: %6.3e" % (self.f_change) + # Check for maximum number of IRLS cycles if self.IRLSiter == self.maxIRLSiter: + print "Reach maximum number of IRLS cycles: %i" % self.maxIRLSiter self.opt.stopNextIteration = True return - + # Check if the function has changed enough - if self.f_change < self.f_min_change: + if self.f_change < self.f_min_change and self.IRLSiter > 1: + print "Minimum decrease in regularization. End of IRLS" self.opt.stopNextIteration = True - return - else: - self.f_old = self.opt.f_last - + return + else: + self.f_old = phim_new + # 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 - + # Reset the regularization matrices so that it is # recalculated for current model self.reg._Wsmall = None self.reg._Wx = None self.reg._Wy = None self.reg._Wz = None - + # 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 - + # Reset the regularization matrices again for new gamma self.reg._Wsmall = None self.reg._Wx = None self.reg._Wy = None self.reg._Wz = None - # Compute the change in + # 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.beta_tol: + self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d + class Update_lin_PreCond(InversionDirective): """ Create a Jacobi preconditioner for the linear problem @@ -409,21 +459,15 @@ class Update_Wj(InversionDirective): 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 - coolingRate=5 - - def endIter(self): - - # Only update after GN iterations - if self.opt.iter % self.coolingRate == 0: - # 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 +#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 +# coolingRate=5 +# +# def endIter(self): +# +# diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index d1da048f..5c1a8770 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -52,51 +52,49 @@ def run(N=200, plotIt=True): wr = np.sum(prob.G**2.,axis=0)**0.5 wr = ( wr/np.max(wr) ) - reg = Regularization.Simple(mesh) - reg.mref = mref - reg.cell_weights = wr - +# reg = Regularization.Simple(mesh) +# reg.mref = mref +# reg.cell_weights = wr +# dmis = DataMisfit.l2_DataMisfit(survey) dmis.Wd = 1./wd - - opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, tolCG = 1e-4) - invProb = InvProblem.BaseInvProblem(dmis, reg, opt) - invProb.curModel = m0 - - beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1) - target = Directives.TargetMisfit() - +# +# opt = Optimization.ProjectedGNCG(maxIter=20,lower=-2.,upper=2., maxIterCG= 10, 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 +# 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) reg.mref = mref reg.cell_weights = wr reg.mref = np.zeros(mesh.nC) - reg.eps_p = 1e-3 - reg.eps_q = 1e-2 - reg.norms = [0., 0., 2., 2.] + eps_p = 1e-3 + eps_q = 1e-2 + norms = [0., 0., 2., 2.] opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3) - invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.) - beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1) - update_beta = Directives.Scale_Beta(tol = 0.05, coolingRate=5) - target = Directives.TargetMisfit() - IRLS = Directives.Update_IRLS( phi_m_last = phim, phi_d_last = phid, coolingRate=5 ) + invProb = InvProblem.BaseInvProblem(dmis, reg, opt) + #beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1) + #update_beta = Directives.Scale_Beta(tol = 0.05, coolingRate=5) +# target = Directives.TargetMisfit() + IRLS = Directives.Update_IRLS( norms=norms, eps_p=eps_p, eps_q=eps_q) - inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS,update_beta]) - - m0 = mrec + inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest]) # Run inversion mrec = inv.run(m0) @@ -113,7 +111,7 @@ def run(N=200, plotIt=True): axes[0].set_title('Columns of matrix G') axes[1].plot(mesh.vectorCCx, mtrue, 'b-') - axes[1].plot(mesh.vectorCCx, ml2, 'r-') + axes[1].plot(mesh.vectorCCx, reg.l2model, 'r-') #axes[1].legend(('True Model', 'Recovered Model')) axes[1].set_ylim(-1.0,1.25) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index c3b1c9e5..3304022a 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -918,6 +918,7 @@ class Sparse(Simple): eps_p = 1e-1 # Threshold value for the model norm eps_q = 1e-1 # Threshold value for the model gradient norm curModel = None # Requires model to compute the weights + l2model = None gamma = 1. # Model norm scaling to smooth out convergence norms = [0., 2., 2., 2.] # Values for norm on (m, dmdx, dmdy, dmdz) cell_weights = 1. # Consider overwriting with sensitivity weights From 825511e9d3756db1cd92985b31a43f1d18c681b1 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Sun, 29 May 2016 13:17:22 -0700 Subject: [PATCH 24/27] Add a curvilinear plotImage function, update example. --- SimPEG/Examples/Forward_BasicDirectCurrent.py | 33 +---- SimPEG/Mesh/CurvilinearMesh.py | 99 +-------------- SimPEG/Mesh/View.py | 118 ++++++++++++------ 3 files changed, 84 insertions(+), 166 deletions(-) diff --git a/SimPEG/Examples/Forward_BasicDirectCurrent.py b/SimPEG/Examples/Forward_BasicDirectCurrent.py index efbb287c..30cf0ed1 100644 --- a/SimPEG/Examples/Forward_BasicDirectCurrent.py +++ b/SimPEG/Examples/Forward_BasicDirectCurrent.py @@ -5,18 +5,15 @@ from SimPEG import Mesh, Utils, np, SolverLU def run(plotIt=True): # Step1: Generate Tensor and Curvilinear Mesh sz = [40,40] - # Tensor Mesh tM = Mesh.TensorMesh(sz) - # Curvilinear Mesh rM = Mesh.CurvilinearMesh(Utils.meshutils.exampleLrmGrid(sz,'rotate')) + # Step2: Direct Current (DC) operator def DCfun(mesh, pts): D = mesh.faceDiv - G = D.T sigma = 1e-2*np.ones(mesh.nC) - Msigi = mesh.getFaceInnerProduct(1./sigma) - MsigI = Utils.sdInv(Msigi) - A = D*MsigI*G + MsigI = mesh.getFaceInnerProduct(sigma, invProp=True, invMat=True) + A = -D*MsigI*D.T A[-1,-1] /= mesh.vol[-1] # Remove null space rhs = np.zeros(mesh.nC) txind = Utils.meshutils.closestPoints(mesh, pts) @@ -37,39 +34,17 @@ def run(plotIt=True): if not plotIt: return import matplotlib.pyplot as plt - import matplotlib - from matplotlib.mlab import griddata #Step4: Making Figure fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2)) - label = ["(a)", "(b)"] - opts = {} vmin, vmax = phitM.min(), phitM.max() dat = tM.plotImage(phitM, ax=axes[0], clim=(vmin, vmax), grid=True) - - #TODO: At the moment Curvilinear Mesh do not have plotimage - - Xi = tM.gridCC[:,0].reshape(sz[0], sz[1], order='F') - Yi = tM.gridCC[:,1].reshape(sz[0], sz[1], order='F') - PHIrM = griddata(rM.gridCC[:,0], rM.gridCC[:,1], phirM, Xi, Yi, interp='linear') - axes[1].contourf(Xi, Yi, PHIrM, 100, vmin=vmin, vmax=vmax) - + dat = rM.plotImage(phirM, ax=axes[1], clim=(vmin, vmax), grid=True) cb = plt.colorbar(dat[0], ax=axes[0]); cb.set_label("Voltage (V)") cb = plt.colorbar(dat[0], ax=axes[1]); cb.set_label("Voltage (V)") - tM.plotGrid(ax=axes[0], **opts) axes[0].set_title('TensorMesh') - rM.plotGrid(ax=axes[1], **opts) axes[1].set_title('CurvilinearMesh') - for i in range(2): - axes[i].set_xlim(0.025, 0.975) - axes[i].set_ylim(0.025, 0.975) - axes[i].text(0., 1.0, label[i], fontsize=20) - if i==0: - axes[i].set_ylabel("y") - else: - axes[i].set_ylabel(" ") - axes[i].set_xlabel("x") plt.show() diff --git a/SimPEG/Mesh/CurvilinearMesh.py b/SimPEG/Mesh/CurvilinearMesh.py index f8b0bcd2..abc9ff3a 100644 --- a/SimPEG/Mesh/CurvilinearMesh.py +++ b/SimPEG/Mesh/CurvilinearMesh.py @@ -2,6 +2,7 @@ from SimPEG import Utils, np from BaseMesh import BaseRectangularMesh from DiffOperators import DiffOperators from InnerProducts import InnerProducts +from View import CurvView # Some helper functions. length2D = lambda x: (x[:, 0]**2 + x[:, 1]**2)**0.5 @@ -10,7 +11,7 @@ normalize2D = lambda x: x/np.kron(np.ones((1, 2)), Utils.mkvc(length2D(x), 2)) normalize3D = lambda x: x/np.kron(np.ones((1, 3)), Utils.mkvc(length3D(x), 2)) -class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts): +class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts, CurvView): """ CurvilinearMesh is a mesh class that deals with curvilinear meshes. @@ -330,102 +331,6 @@ class CurvilinearMesh(BaseRectangularMesh, DiffOperators, InnerProducts): - ############################################# - # Plotting Functions # - ############################################# - - def plotGrid(self, ax=None, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False): - """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions. - - - .. plot:: - :include-source: - - from SimPEG import Mesh, Utils - X, Y = Utils.exampleLrmGrid([3,3],'rotate') - M = Mesh.CurvilinearMesh([X, Y]) - M.plotGrid(showIt=True) - - """ - import matplotlib.pyplot as plt - import matplotlib - from mpl_toolkits.mplot3d import Axes3D - mkvc = Utils.mkvc - - axOpts = {'projection':'3d'} if self.dim == 3 else {} - if ax is None: ax = plt.subplot(111, **axOpts) - - NN = self.r(self.gridN, 'N', 'N', 'M') - if self.dim == 2: - - if lines: - X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten() - Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten() - - X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten() - Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten() - - X = np.r_[X1, X2] - Y = np.r_[Y1, Y2] - - ax.plot(X, Y, 'b-') - if centers: - ax.plot(self.gridCC[:,0],self.gridCC[:,1],'ro') - - # Nx = self.r(self.normals, 'F', 'Fx', 'V') - # Ny = self.r(self.normals, 'F', 'Fy', 'V') - # Tx = self.r(self.tangents, 'E', 'Ex', 'V') - # Ty = self.r(self.tangents, 'E', 'Ey', 'V') - - # ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo') - - # nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten() - # nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten() - # ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs') - # ax.plot(nX, nY, 'r-') - - # nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten() - # nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten() - # #ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs') - # ax.plot(nX, nY, 'g-') - - # tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten() - # tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten() - # ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^') - # ax.plot(tX, tY, 'r-') - - # nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten() - # nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten() - # #ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^') - # ax.plot(nX, nY, 'g-') - - elif self.dim == 3: - X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten() - Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten() - Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten() - - X2 = np.c_[mkvc(NN[0][:, :-1, :]), mkvc(NN[0][:, 1:, :]), mkvc(NN[0][:, :-1, :])*np.nan].flatten() - Y2 = np.c_[mkvc(NN[1][:, :-1, :]), mkvc(NN[1][:, 1:, :]), mkvc(NN[1][:, :-1, :])*np.nan].flatten() - Z2 = np.c_[mkvc(NN[2][:, :-1, :]), mkvc(NN[2][:, 1:, :]), mkvc(NN[2][:, :-1, :])*np.nan].flatten() - - X3 = np.c_[mkvc(NN[0][:, :, :-1]), mkvc(NN[0][:, :, 1:]), mkvc(NN[0][:, :, :-1])*np.nan].flatten() - Y3 = np.c_[mkvc(NN[1][:, :, :-1]), mkvc(NN[1][:, :, 1:]), mkvc(NN[1][:, :, :-1])*np.nan].flatten() - Z3 = np.c_[mkvc(NN[2][:, :, :-1]), mkvc(NN[2][:, :, 1:]), mkvc(NN[2][:, :, :-1])*np.nan].flatten() - - X = np.r_[X1, X2, X3] - Y = np.r_[Y1, Y2, Y3] - Z = np.r_[Z1, Z2, Z3] - - ax.plot(X, Y, 'b', zs=Z) - ax.set_zlabel('x3') - - ax.grid(True) - ax.set_xlabel('x1') - ax.set_ylabel('x2') - - if showIt: plt.show() - - if __name__ == '__main__': nc = 5 h1 = np.cumsum(np.r_[0, np.ones(nc)/(nc)]) diff --git a/SimPEG/Mesh/View.py b/SimPEG/Mesh/View.py index 8eb22098..72225040 100644 --- a/SimPEG/Mesh/View.py +++ b/SimPEG/Mesh/View.py @@ -552,7 +552,8 @@ class CurvView(object): def __init__(self): pass - def plotGrid(self, length=0.05, showIt=False): + + def plotGrid(self, ax=None, nodes=False, faces=False, centers=False, edges=False, lines=True, showIt=False): """Plot the nodal, cell-centered and staggered grids for 1,2 and 3 dimensions. @@ -560,60 +561,63 @@ class CurvView(object): :include-source: from SimPEG import Mesh, Utils - X, Y = Utils.exampleCurvGird([3,3],'rotate') + X, Y = Utils.exampleLrmGrid([3,3],'rotate') M = Mesh.CurvilinearMesh([X, Y]) M.plotGrid(showIt=True) """ + import matplotlib.pyplot as plt + import matplotlib + from mpl_toolkits.mplot3d import Axes3D + + axOpts = {'projection':'3d'} if self.dim == 3 else {} + if ax is None: ax = plt.subplot(111, **axOpts) + NN = self.r(self.gridN, 'N', 'N', 'M') if self.dim == 2: - fig = plt.figure(2) - fig.clf() - ax = plt.subplot(111) - X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten() - Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten() - X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten() - Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten() + if lines: + X1 = np.c_[mkvc(NN[0][:-1, :]), mkvc(NN[0][1:, :]), mkvc(NN[0][:-1, :])*np.nan].flatten() + Y1 = np.c_[mkvc(NN[1][:-1, :]), mkvc(NN[1][1:, :]), mkvc(NN[1][:-1, :])*np.nan].flatten() - X = np.r_[X1, X2] - Y = np.r_[Y1, Y2] + X2 = np.c_[mkvc(NN[0][:, :-1]), mkvc(NN[0][:, 1:]), mkvc(NN[0][:, :-1])*np.nan].flatten() + Y2 = np.c_[mkvc(NN[1][:, :-1]), mkvc(NN[1][:, 1:]), mkvc(NN[1][:, :-1])*np.nan].flatten() - plt.plot(X, Y) + X = np.r_[X1, X2] + Y = np.r_[Y1, Y2] - plt.hold(True) - Nx = self.r(self.normals, 'F', 'Fx', 'V') - Ny = self.r(self.normals, 'F', 'Fy', 'V') - Tx = self.r(self.tangents, 'E', 'Ex', 'V') - Ty = self.r(self.tangents, 'E', 'Ey', 'V') + ax.plot(X, Y, 'b-') + if centers: + ax.plot(self.gridCC[:,0],self.gridCC[:,1],'ro') - plt.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo') + # Nx = self.r(self.normals, 'F', 'Fx', 'V') + # Ny = self.r(self.normals, 'F', 'Fy', 'V') + # Tx = self.r(self.tangents, 'E', 'Ex', 'V') + # Ty = self.r(self.tangents, 'E', 'Ey', 'V') - nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten() - nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten() - plt.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs') - plt.plot(nX, nY, 'r-') + # ax.plot(self.gridN[:, 0], self.gridN[:, 1], 'bo') - nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten() - nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten() - #plt.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs') - plt.plot(nX, nY, 'g-') + # nX = np.c_[self.gridFx[:, 0], self.gridFx[:, 0] + Nx[0]*length, self.gridFx[:, 0]*np.nan].flatten() + # nY = np.c_[self.gridFx[:, 1], self.gridFx[:, 1] + Nx[1]*length, self.gridFx[:, 1]*np.nan].flatten() + # ax.plot(self.gridFx[:, 0], self.gridFx[:, 1], 'rs') + # ax.plot(nX, nY, 'r-') - tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten() - tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten() - plt.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^') - plt.plot(tX, tY, 'r-') + # nX = np.c_[self.gridFy[:, 0], self.gridFy[:, 0] + Ny[0]*length, self.gridFy[:, 0]*np.nan].flatten() + # nY = np.c_[self.gridFy[:, 1], self.gridFy[:, 1] + Ny[1]*length, self.gridFy[:, 1]*np.nan].flatten() + # #ax.plot(self.gridFy[:, 0], self.gridFy[:, 1], 'gs') + # ax.plot(nX, nY, 'g-') - nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten() - nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten() - #plt.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^') - plt.plot(nX, nY, 'g-') - plt.axis('equal') + # tX = np.c_[self.gridEx[:, 0], self.gridEx[:, 0] + Tx[0]*length, self.gridEx[:, 0]*np.nan].flatten() + # tY = np.c_[self.gridEx[:, 1], self.gridEx[:, 1] + Tx[1]*length, self.gridEx[:, 1]*np.nan].flatten() + # ax.plot(self.gridEx[:, 0], self.gridEx[:, 1], 'r^') + # ax.plot(tX, tY, 'r-') + + # nX = np.c_[self.gridEy[:, 0], self.gridEy[:, 0] + Ty[0]*length, self.gridEy[:, 0]*np.nan].flatten() + # nY = np.c_[self.gridEy[:, 1], self.gridEy[:, 1] + Ty[1]*length, self.gridEy[:, 1]*np.nan].flatten() + # #ax.plot(self.gridEy[:, 0], self.gridEy[:, 1], 'g^') + # ax.plot(nX, nY, 'g-') elif self.dim == 3: - fig = plt.figure(3) - fig.clf() - ax = fig.add_subplot(111, projection='3d') X1 = np.c_[mkvc(NN[0][:-1, :, :]), mkvc(NN[0][1:, :, :]), mkvc(NN[0][:-1, :, :])*np.nan].flatten() Y1 = np.c_[mkvc(NN[1][:-1, :, :]), mkvc(NN[1][1:, :, :]), mkvc(NN[1][:-1, :, :])*np.nan].flatten() Z1 = np.c_[mkvc(NN[2][:-1, :, :]), mkvc(NN[2][1:, :, :]), mkvc(NN[2][:-1, :, :])*np.nan].flatten() @@ -630,16 +634,50 @@ class CurvView(object): Y = np.r_[Y1, Y2, Y3] Z = np.r_[Z1, Z2, Z3] - plt.plot(X, Y, 'b', zs=Z) + ax.plot(X, Y, 'b', zs=Z) ax.set_zlabel('x3') ax.grid(True) - ax.hold(False) ax.set_xlabel('x1') ax.set_ylabel('x2') if showIt: plt.show() + def plotImage(self, I, ax=None, showIt=False, grid=False, clim=None): + if self.dim == 3: raise NotImplementedError('This is not yet done!') + + import matplotlib.pyplot as plt + import matplotlib + from mpl_toolkits.mplot3d import Axes3D + import matplotlib.colors as colors + import matplotlib.cm as cmx + + if ax is None: ax = plt.subplot(111) + jet = cm = plt.get_cmap('jet') + cNorm = colors.Normalize( + vmin=I.min() if clim is None else clim[0], + vmax=I.max() if clim is None else clim[1]) + + scalarMap = cmx.ScalarMappable(norm=cNorm, cmap=jet) + # ax.set_xlim((self.x0[0], self.h[0].sum())) + # ax.set_ylim((self.x0[1], self.h[1].sum())) + + Nx = self.r(self.gridN[:,0],'N','N','M') + Ny = self.r(self.gridN[:,1],'N','N','M') + cell = self.r(I,'CC','CC','M') + + for ii in range(self.nCx): + for jj in range(self.nCy): + I = [ii,ii+1,ii+1,ii] + J = [jj,jj,jj+1,jj+1] + ax.add_patch(plt.Polygon(np.c_[Nx[I,J],Ny[I,J]], facecolor=scalarMap.to_rgba(cell[ii,jj]), edgecolor='k' if grid else 'none')) + + scalarMap._A = [] # http://stackoverflow.com/questions/8342549/matplotlib-add-colorbar-to-a-sequence-of-line-plots + ax.set_xlabel('x') + ax.set_ylabel('y') + if showIt: plt.show() + return [scalarMap] + if __name__ == '__main__': from SimPEG import * From e476bf0059b3c6f8e75665bffa1e699b586919a1 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Sun, 29 May 2016 14:09:29 -0700 Subject: [PATCH 25/27] Cleanup Sparse Reg and Directives --- SimPEG/Directives.py | 71 +++++++------------------------ SimPEG/Examples/Inversion_IRLS.py | 16 +++---- 2 files changed, 23 insertions(+), 64 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 8727d925..98770666 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -258,66 +258,39 @@ class Update_IRLS(InversionDirective): f_old = None f_min_change = 1e-2 beta_tol = 5e-2 - + # Solving parameter for IRLS (mode:2) IRLSiter = 0 - minGNiter = 6 + minGNiter = 5 maxIRLSiter = 10 iterStart = 0 - + # Beta schedule coolingFactor = 2. coolingRate = 1 - + mode = 1 - + @property def target(self): if getattr(self, '_target', None) is None: - self._target = self.survey.nD + self._target = self.survey.nD*0.5 return self._target @target.setter def target(self, val): self._target = val - + def initialize(self): if self.mode == 1: self.reg.norms = [2., 2., 2., 2.] -# # 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._Wsmall = None -# self.reg._Wx = None -# self.reg._Wy = None -# self.reg._Wz = None -# 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 -# # Reset the regularization matrices so that it is -# # recalculated with new gamma -# self.reg._Wsmall = None -# self.reg._Wx = None -# self.reg._Wy = None -# self.reg._Wz = None -# if getattr(self, 'phi_d_last', None) is None: -# self.phi_d_last = self.invProb.phi_d -# -# if getattr(self, 'f_last', None) is None: -# self.f_old = self.invProb.evalFunction(self.reg.curModel, return_g=False, return_H=False) -# print self.f_old - def endIter(self): # After reaching target misfit with l2-norm, switch to IRLS (mode:2) if self.invProb.phi_d < self.target and self.mode == 1: print "Convergence with smooth l2-norm regularization: Start IRLS steps..." - + self.mode = 2 print self.eps_p, self.eps_q, self.norms self.reg.eps_p = self.eps_p @@ -328,17 +301,18 @@ class Update_IRLS(InversionDirective): self.iterStart = self.opt.iter self.phi_d_last = self.invProb.phi_d self.phi_m_last = self.invProb.phi_m_last - + self.reg.l2model = self.invProb.curModel self.reg.curModel = self.invProb.curModel - + if getattr(self, 'f_old', None) is None: self.f_old = self.reg.eval(self.invProb.curModel)#self.invProb.evalFunction(self.invProb.curModel, return_g=False, return_H=False) - + + # Beta Schedule if self.opt.iter > 0 and self.opt.iter % self.coolingRate == 0: if self.debug: print 'BetaSchedule is cooling Beta. Iteration: %d' % self.opt.iter self.invProb.beta /= self.coolingFactor - + # Only update after GN iterations if (self.opt.iter-self.iterStart) % self.minGNiter == 0 and self.mode==2: @@ -401,12 +375,12 @@ class Update_IRLS(InversionDirective): self.reg._Wy = None self.reg._Wz = None - # Check if misfit is within the tolerance, otherwise adjust beta + # Check if misfit is within the tolerance, otherwise scale beta val = self.invProb.phi_d / (self.survey.nD*0.5) - + if np.abs(1.-val) > self.beta_tol: self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d - + class Update_lin_PreCond(InversionDirective): """ Create a Jacobi preconditioner for the linear problem @@ -458,16 +432,3 @@ 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 -# coolingRate=5 -# -# def endIter(self): -# -# diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index 5c1a8770..afd90525 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -1,7 +1,7 @@ from SimPEG import * -def run(N=200, plotIt=True): +def run(N=100, plotIt=True): """ Inversion: Linear Problem ========================= @@ -19,8 +19,8 @@ def run(N=200, plotIt=True): m0 = np.ones(mesh.nC) * 1e-4 mref = np.zeros(mesh.nC) - - nk = 15 + + nk = 10 jk = np.linspace(1.,nk,nk) p = -2. q = 1. @@ -83,18 +83,16 @@ def run(N=200, plotIt=True): reg.cell_weights = wr reg.mref = np.zeros(mesh.nC) - eps_p = 1e-3 - eps_q = 1e-2 + eps_p = 5e-2 + eps_q = 5e-2 norms = [0., 0., 2., 2.] opt = Optimization.ProjectedGNCG(maxIter=100 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 10, tolCG = 1e-3) invProb = InvProblem.BaseInvProblem(dmis, reg, opt) - #beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1) - #update_beta = Directives.Scale_Beta(tol = 0.05, coolingRate=5) -# target = Directives.TargetMisfit() + update_Jacobi = Directives.Update_lin_PreCond() IRLS = Directives.Update_IRLS( norms=norms, eps_p=eps_p, eps_q=eps_q) - inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest]) + inv = Inversion.BaseInversion(invProb, directiveList=[IRLS,betaest,update_Jacobi]) # Run inversion mrec = inv.run(m0) From 62eb4541cbdb086023870a1df0564a74498dafe2 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 29 May 2016 15:48:15 -0700 Subject: [PATCH 26/27] update example name --> based on mesh --- ...sicDirectCurrent.py => Mesh_Basic_ForwardDC.py} | 10 ++++++++-- SimPEG/Examples/__init__.py | 4 ++-- ...cDirectCurrent.rst => Mesh_Basic_ForwardDC.rst} | 14 +++++++++----- 3 files changed, 19 insertions(+), 9 deletions(-) rename SimPEG/Examples/{Forward_BasicDirectCurrent.py => Mesh_Basic_ForwardDC.py} (89%) rename docs/examples/{Forward_BasicDirectCurrent.rst => Mesh_Basic_ForwardDC.rst} (55%) diff --git a/SimPEG/Examples/Forward_BasicDirectCurrent.py b/SimPEG/Examples/Mesh_Basic_ForwardDC.py similarity index 89% rename from SimPEG/Examples/Forward_BasicDirectCurrent.py rename to SimPEG/Examples/Mesh_Basic_ForwardDC.py index 30cf0ed1..bd4ea18a 100644 --- a/SimPEG/Examples/Forward_BasicDirectCurrent.py +++ b/SimPEG/Examples/Mesh_Basic_ForwardDC.py @@ -1,8 +1,14 @@ from SimPEG import Mesh, Utils, np, SolverLU -## 2D DC forward modeling example with Tensor and Curvilinear Meshes - def run(plotIt=True): + + """ + Mesh: Basic Forward 2D DC Resistivity + ===================================== + + 2D DC forward modeling example with Tensor and Curvilinear Meshes + """ + # Step1: Generate Tensor and Curvilinear Mesh sz = [40,40] tM = Mesh.TensorMesh(sz) diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index c67652a2..c592e454 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -8,9 +8,9 @@ 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_ForwardDC import Mesh_Basic_PlotImage import Mesh_Basic_Types import Mesh_Operators_CahnHilliard @@ -21,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_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"] +__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", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_ForwardDC", "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/docs/examples/Forward_BasicDirectCurrent.rst b/docs/examples/Mesh_Basic_ForwardDC.rst similarity index 55% rename from docs/examples/Forward_BasicDirectCurrent.rst rename to docs/examples/Mesh_Basic_ForwardDC.rst index 20b39eb8..df18050e 100644 --- a/docs/examples/Forward_BasicDirectCurrent.rst +++ b/docs/examples/Mesh_Basic_ForwardDC.rst @@ -1,4 +1,4 @@ -.. _examples_Forward_BasicDirectCurrent: +.. _examples_Mesh_Basic_ForwardDC: .. --------------------------------- .. .. .. @@ -8,14 +8,18 @@ .. .. .. --------------------------------- .. -Forward BasicDirectCurrent -========================== + +Mesh: Basic Forward 2D DC Resistivity +===================================== + +2D DC forward modeling example with Tensor and Curvilinear Meshes + .. plot:: from SimPEG import Examples - Examples.Forward_BasicDirectCurrent.run() + Examples.Mesh_Basic_ForwardDC.run() -.. literalinclude:: ../../SimPEG/Examples/Forward_BasicDirectCurrent.py +.. literalinclude:: ../../SimPEG/Examples/Mesh_Basic_ForwardDC.py :language: python :linenos: From 279bd49b4c8cb401023aa8a62091575c8919993a Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 31 May 2016 15:38:05 -0700 Subject: [PATCH 27/27] =?UTF-8?q?Bump=20version:=200.1.10=20=E2=86=92=200.?= =?UTF-8?q?1.11?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .bumpversion.cfg | 2 +- SimPEG/__init__.py | 2 +- docs/conf.py | 4 ++-- setup.py | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/.bumpversion.cfg b/.bumpversion.cfg index 84469e3f..3ca3f780 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,4 +1,4 @@ [bumpversion] -current_version = 0.1.10 +current_version = 0.1.11 files = setup.py SimPEG/__init__.py docs/conf.py diff --git a/SimPEG/__init__.py b/SimPEG/__init__.py index a1e989b6..23d4de4a 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -15,7 +15,7 @@ import Directives import Inversion import Tests -__version__ = '0.1.10' +__version__ = '0.1.11' __author__ = 'Rowan Cockett' __license__ = 'MIT' __copyright__ = 'Copyright 2014 Rowan Cockett' diff --git a/docs/conf.py b/docs/conf.py index 45407435..67351319 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -51,9 +51,9 @@ copyright = u'2013, SimPEG Developers' # built documents. # # The short X.Y version. -version = '0.1.10' +version = '0.1.11' # The full version, including alpha/beta/rc tags. -release = '0.1.10' +release = '0.1.11' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. diff --git a/setup.py b/setup.py index 5308db47..06cd9b27 100644 --- a/setup.py +++ b/setup.py @@ -83,7 +83,7 @@ with open("README.rst") as f: setup( name = "SimPEG", - version = "0.1.10", + version = "0.1.11", packages = find_packages(), install_requires = ['numpy>=1.7', 'scipy>=0.13',