From 05a85018de792c16c8e728f4ab893c715bcd52de Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 16 Feb 2016 21:59:08 -0800 Subject: [PATCH 01/18] Create a Simple and a SparseRegularization class. The SparseRegularization class allows implementation of p-q norms. --- SimPEG/Regularization.py | 210 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 210 insertions(+) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 677fba5c..6a65860d 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -282,3 +282,213 @@ class Tikhonov(BaseRegularization): out = mD.T * ( self.W.T * r ) return out + +class Simple(BaseRegularization): + """ + Only for tensor mesh + """ + + smoothModel = True #: SMOOTH and SMOOTH_MOD_DIF options + alpha_s = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Ws'], "Smallness weight") + alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction") + alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction") + alpha_z = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction") + + def __init__(self, mesh, mapping=None, **kwargs): + BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs) + + @property + def Ws(self): + """Regularization matrix Ws""" + if getattr(self,'_Ws', None) is None: + self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5) + return self._Ws + + @property + def Wx(self): + """Regularization matrix Wx""" + if getattr(self, '_Wx', None) is None: + self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x)**0.5)*self.mesh.unitCellGradx + return self._Wx + + @property + def Wy(self): + """Regularization matrix Wy""" + if getattr(self, '_Wy', None) is None: + self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y)**0.5)*self.mesh.unitCellGrady + return self._Wy + + @property + def Wz(self): + """Regularization matrix Wz""" + if getattr(self, '_Wz', None) is None: + self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z)**0.5)*self.mesh.unitCellGradz + return self._Wz + + @property + def Wsmooth(self): + """Full smoothness regularization matrix W""" + if getattr(self, '_Wsmooth', None) is None: + wlist = (self.Wx,) + if self.mesh.dim > 1: + wlist += (self.Wy,) + if self.mesh.dim > 2: + wlist += (self.Wz,) + self._Wsmooth = sp.vstack(wlist) + return self._Wsmooth + + @property + def W(self): + """Full regularization matrix W""" + if getattr(self, '_W', None) is None: + wlist = (self.Ws, self.Wsmooth) + self._W = sp.vstack(wlist) + return self._W + + @Utils.timeIt + def eval(self, m): + if self.smoothModel == True: + r1 = self.Wsmooth * ( self.mapping * (m) ) + r2 = self.Ws * ( self.mapping * (m - self.mref) ) + return 0.5*(r1.dot(r1)+r2.dot(r2)) + elif self.smoothModel == False: + r = self.W * ( self.mapping * (m - self.mref) ) + return 0.5*r.dot(r) + + + @Utils.timeIt + def 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})} + + """ + if self.smoothModel == True: + mD1 = self.mapping.deriv(m) + mD2 = self.mapping.deriv(m - self.mref) + r1 = self.Wsmooth * ( self.mapping * (m)) + r2 = self.Ws * ( self.mapping * (m - self.mref) ) + out1 = mD1.T * ( self.Wsmooth.T * r1 ) + out2 = mD2.T * ( self.Ws.T * r2 ) + out = out1+out2 + elif self.smoothModel == False: + mD = self.mapping.deriv(m - self.mref) + r = self.W * ( self.mapping * (m - self.mref) ) + out = mD.T * ( self.W.T * r ) + return out + + +class SparseRegularization(Simple): + + eps = 1e-1 + m = None + gamma = 1. + p = 0. + qx = 2. + qy = 2. + qz = 2. + + def __init__(self, mesh, mapping=None, **kwargs): + Simple.__init__(self, mesh, mapping=mapping, **kwargs) + + @property + def Wsmooth(self): + """Full smoothness regularization matrix W""" + if getattr(self, '_Wsmooth', None) is None: + wlist = (self.Wx, self.Wxx) + if self.mesh.dim > 1: + wlist += (self.Wy, self.Wyy) + if self.mesh.dim > 2: + wlist += (self.Wz, self.Wzz) + self._Wsmooth = sp.vstack(wlist) + return self._Wsmooth + + @property + def W(self): + """Full regularization matrix W""" + if getattr(self, '_W', None) is None: + wlist = (self.Ws, self.Wsmooth) + self._W = sp.vstack(wlist) + return self._W + + @property + def Ws(self): + """Regularization matrix Ws""" + if getattr(self, 'm', None) is None: + self.Rs = Utils.speye(self.mesh.nC) + + else: + f_m = self.m + self.rs = self.R(f_m , self.p, self.eps) + #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs)) + self.Rs = Utils.sdiag( self.rs ) + + self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s*self.gamma)**0.5)*self.Rs + + return self._Ws + + @property + def Wx(self): + """Regularization matrix Wx""" + + if getattr(self, 'm', None) is None: + self.Rx = Utils.speye(self.mesh.unitCellGradx.shape[0]) + + else: + f_m = self.mesh.unitCellGradx * self.m + self.rx = self.R( f_m , self.qx, self.eps) + self.Rx = Utils.sdiag( self.rx ) + + if getattr(self, '_Wx', None) is None: + self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x*self.gamma)**0.5)*self.Rx*self.mesh.unitCellGradx + return self._Wx + + @property + def Wy(self): + """Regularization matrix Wy""" + + if getattr(self, 'm', None) is None: + self.Ry = Utils.speye(self.mesh.unitCellGrady.shape[0]) + + else: + f_m = self.mesh.unitCellGrady * self.m + self.ry = self.R( f_m , self.qy, self.eps) + self.Ry = Utils.sdiag( self.ry ) + + if getattr(self, '_Wy', None) is None: + self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y*self.gamma)**0.5)*self.Ry*self.mesh.unitCellGrady + return self._Wy + + @property + def Wz(self): + """Regularization matrix Wz""" + + if getattr(self, 'm', None) is None: + self.Rz = Utils.speye(self.mesh.unitCellGradz.shape[0]) + + else: + f_m = self.mesh.unitCellGradz * self.m + self.rz = self.R( f_m , self.qz, self.eps) + self.Rz = Utils.sdiag( self.rz ) + + if getattr(self, '_Wz', None) is None: + self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z*self.gamma)**0.5)*self.Rz*self.mesh.unitCellGradz + return self._Wz + + + def R(self, f_m , p, dec): + + eta = (self.eps**(1-p/2.))**0.5 + r = eta / (f_m**2.+self.eps**2.)**((1-p/2.)/2.) + + return r From c10777a245510e95329dd020f6567ce26aa60f41 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 16 Feb 2016 22:00:12 -0800 Subject: [PATCH 02/18] Addition of unitCellGrad. Possibly rename to cellGradStencil? --- SimPEG/Mesh/DiffOperators.py | 50 ++++++++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index 00bf0259..22b87edd 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -565,6 +565,56 @@ class DiffOperators(object): return Pbc, Pin, Pout + def unitCellGradx(): + doc = """Cell centered Gradient in the x dimension used for + regularization. The gradient operator is square (nC-by-nC)""" + def fget(self): + if self.dim < 3: return None + if getattr(self, '_unitCellGradx', None) is None: + + n = self.vnC + gx = ddx(n[0]-1) + gx_square = sp.vstack((gx,gx[-1,:]*-1), format="csr") + + self._unitCellGradx = kron3(speye(n[2]), speye(n[1]), gx_square) + + return self._unitCellGradx + return locals() + unitCellGradx = property(**unitCellGradx()) + + def unitCellGrady(): + doc = """Cell centered Gradient in they dimension used for + regularization. The gradient operator is square (nC-by-nC)""" + def fget(self): + if self.dim < 3: return None + if getattr(self, '_unitCellGrady', None) is None: + + n = self.vnC + gy = ddx(n[1]-1) + gy_square = sp.vstack((gy,gy[-1,:]*-1), format="csr") + + self._unitCellGrady = kron3(speye(n[2]), gy_square, speye(n[0])) + + return self._unitCellGrady + return locals() + unitCellGrady = property(**unitCellGrady()) + + def unitCellGradz(): + doc = """Cell centered Gradient in they dimension used for + regularization. The gradient operator is square (nC-by-nC)""" + def fget(self): + if self.dim < 3: return None + if getattr(self, '_unitCellGradz', None) is None: + + n = self.vnC + gz = ddx(n[2]-1) + gz_square = sp.vstack((gz,gz[-1,:]*-1), format="csr") + + self._unitCellGradz = kron3( gz_square , speye(n[1]), speye(n[0])) + + return self._unitCellGradz + return locals() + unitCellGradz = property(**unitCellGradz()) # --------------- Averaging --------------------- From 1c2fecf3a24e5a90918451fdd5a0ec3f013bff62 Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Tue, 16 Feb 2016 22:07:33 -0800 Subject: [PATCH 03/18] Add the IRLS Directive. --- SimPEG/Directives.py | 50 +++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 49 insertions(+), 1 deletion(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 2ed27c20..476c5916 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -271,7 +271,6 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration): np.savez('{:s}-{:03d}'.format(self.fileName,self.opt.iter), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred) - # class UpdateReferenceModel(Parameter): # mref0 = None @@ -283,3 +282,52 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration): # mref = self.mref0 # self.m_prev = self.invProb.m_current # return mref + +class update_IRLS(InversionDirective): + + m = None + eps_min = None + factor = None + gamma = None + phi_m_last = None + + def initialize(self): + + # Scale the regularization for changes in norm + if getattr(self, 'phi_m_last', None) is not None: + self.reg.gamma = 1. + phim_new = self.reg.eval(self.invProb.curModel) + self.gamma = self.phi_m_last / phim_new + + self.reg.gamma = self.gamma + + def endIter(self): + # Cool the threshold parameter + 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 + + + # Update the model used for the IRLS weights + if getattr(self, 'm', None) is None: + self.reg.m = self.invProb.curModel + + # 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.prob.mesh.nC))**2. + PC = Utils.sdiag(diagA**-1.) + + self.opt.approxHinv = PC + + phim_new = self.reg.eval(self.invProb.curModel) + self.reg.gamma = self.reg.gamma * self.invProb.phi_m_last / phim_new + +#============================================================================== +# import pylab as plt +# plt.figure() +# ax = plt.subplot(221) +# self.prob.mesh.plotSlice(self.invProb.curModel, ax = ax, normal = 'Z', ind=-5, clim = (0, 0.005)) +#============================================================================== From e4a3e0a16d3e07b9a2217fd30f0fc5274a8d9d8a Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 19 Feb 2016 16:23:26 -0800 Subject: [PATCH 04/18] break out the Pac, Pafx, ... and make part of base regularization --- SimPEG/Regularization.py | 121 ++++++++++++++++++++++----------------- 1 file changed, 70 insertions(+), 51 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 6a65860d..0a6f2622 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -55,8 +55,47 @@ class BaseRegularization(object): @property def W(self): """Full regularization weighting matrix W.""" - return sp.identity(self.mapping.nP) + return self._Pac.T * sp.identity(self.mesh.nC) * self. # or do we want sp.identity(self.mesh.nC) or even just Utils.Identity() ? + @property + def _Pac(self): + if getattr(self, '__Pac', None) is None: + if self.indActive is None: + self.__Pac = Utils.speye(self.mesh.nC) + else: + self.__Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + return self.__Pac + + @property + def _Pafx(self): + if getattr(self, '__Pafx', None) is None: + if self.indActive is None: + self.__Pafx = Utils.speye(self.mesh.nFx) + else: + indActive_Fx = (self.mesh.aveFx2CC.T * self.indActive) == 1 + self.__Pafx = Utils.speye(self.mesh.nFx)[:,indActive_Fx] + return self.__Pafx + + @property + def _Pafy(self): + if getattr(self, '__Pafy', None) is None: + if self.indActive is None: + self.__Pafy = Utils.speye(self.mesh.nFy) + else: + indActive_Fy = (self.mesh.aveFy2CC.T * self.indActive) == 1 + self.__Pafy = Utils.speye(self.mesh.nFy)[:,indActive_Fy] + return self.__Pafy + + @property + def _Pafz(self): + if getattr(self, '__Pafz', None) is None: + if self.indActive is None: + self.__Pafz = Utils.speye(self.mesh.nFz) + else: + indActive_Fz = (self.mesh.aveFz2CC.T * self.indActive) == 1 + self.__Pafz = Utils.speye(self.mesh.nFz)[:,indActive_Fz] + return self.__Pafz + @Utils.timeIt def eval(self, m): @@ -133,10 +172,8 @@ class Tikhonov(BaseRegularization): def Ws(self): """Regularization matrix Ws""" if getattr(self,'_Ws', None) is None: - self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5) - if self.indActive is not None: - Pac = Utils.speye(self.mesh.nC)[:,self.indActive] - self._Ws = Pac.T * self._Ws * Pac + Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5) + self._Ws = self._Pac.T * Ws * self._Pac return self._Ws @property @@ -144,14 +181,8 @@ class Tikhonov(BaseRegularization): """Regularization matrix Wx""" if getattr(self, '_Wx', None) is None: Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFx].T*self.mesh.vol - self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx - - if self.indActive is not None: - indActive_Fx = (self.mesh.aveFx2CC.T * self.indActive) == 1 - Pac = Utils.speye(self.mesh.nC)[:,self.indActive] - Pafx = Utils.speye(self.mesh.nFx)[:,indActive_Fx] - self._Wx = Pafx.T*self._Wx*Pac - + Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx + self._Wx = self._Pafx.T*Wx*self.self._Pac return self._Wx @property @@ -159,14 +190,8 @@ class Tikhonov(BaseRegularization): """Regularization matrix Wy""" if getattr(self, '_Wy', None) is None: Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFx:np.sum(self.mesh.vnF[:2])].T*self.mesh.vol - self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady - - if self.indActive is not None: - indActive_Fy = (self.mesh.aveFy2CC.T * self.indActive) == 1 - Pac = Utils.speye(self.mesh.nC)[:,self.indActive] - Pafy = Utils.speye(self.mesh.nFy)[:,indActive_Fy] - self._Wy = Pafy.T*self._Wy*Pac - + Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady + self._Wy = self._Pafy.T*Wy*self._Pac return self._Wy @property @@ -174,50 +199,32 @@ class Tikhonov(BaseRegularization): """Regularization matrix Wz""" if getattr(self, '_Wz', None) is None: Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.vnF[:2]):].T*self.mesh.vol - self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz - - if self.indActive is not None: - indActive_Fz = (self.mesh.aveFz2CC.T * self.indActive) == 1 - Pac = Utils.speye(self.mesh.nC)[:,self.indActive] - Pafz = Utils.speye(self.mesh.nFz)[:,indActive_Fz] - self._Wz = Pafz.T*self._Wz*Pac - + Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz + self._Wz = self._Pafz.T*Wz*self._Pac return self._Wz @property def Wxx(self): """Regularization matrix Wxx""" if getattr(self, '_Wxx', None) is None: - self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx - - if self.indActive is not None: - Pac = Utils.speye(self.mesh.nC)[:,self.indActive] - self._Wxx = Pac.T*self._Wxx*Pac - + Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx + self._Wxx = self._Pac.T*Wxx*self._Pac return self._Wxx @property def Wyy(self): """Regularization matrix Wyy""" if getattr(self, '_Wyy', None) is None: - self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady - - if self.indActive is not None: - Pac = Utils.speye(self.mesh.nC)[:,self.indActive] - self._Wyy = Pac.T*self._Wyy*Pac - + Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady + self._Wyy = self._Pac.T*self._Wyy*self._Pac return self._Wyy @property def Wzz(self): """Regularization matrix Wzz""" if getattr(self, '_Wzz', None) is None: - self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz - - if self.indActive is not None: - Pac = Utils.speye(self.mesh.nC)[:,self.indActive] - self._Wzz = Pac.T*self._Wzz*Pac - + Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz + self._Wzz = self._Pac.T*Wzz*self._Pac return self._Wzz @property @@ -346,14 +353,26 @@ class Simple(BaseRegularization): return self._W @Utils.timeIt - def eval(self, m): + def _evalSmall(self, m): + r = self.W * ( self.mapping * (m - self.mref) ) + return 0.5*r.dot(r) + + @Utils.timeIt + def _evalSmooth(self, m): if self.smoothModel == True: r1 = self.Wsmooth * ( self.mapping * (m) ) r2 = self.Ws * ( self.mapping * (m - self.mref) ) return 0.5*(r1.dot(r1)+r2.dot(r2)) - elif self.smoothModel == False: - r = self.W * ( self.mapping * (m - self.mref) ) - return 0.5*r.dot(r) + else: + return None + + @Utils.timeIt + def eval(self, m): + phim = self._evalSmall(m) + if self.smoothModel is True: + phim += self._evalSmooth(m) + return phim + @Utils.timeIt From b5f4d8e9997456b4bb6ef2ce84274d59bb8f49a6 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 19 Feb 2016 17:43:50 -0800 Subject: [PATCH 05/18] typo in Regularization.py --- SimPEG/Regularization.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 0a6f2622..d6780c45 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -55,7 +55,7 @@ class BaseRegularization(object): @property def W(self): """Full regularization weighting matrix W.""" - return self._Pac.T * sp.identity(self.mesh.nC) * self. # or do we want sp.identity(self.mesh.nC) or even just Utils.Identity() ? + return self._Pac.T * sp.identity(self.mesh.nC) * self._Pac # or do we want sp.identity(self.mesh.nC) or even just Utils.Identity() ? @property def _Pac(self): From 4e871a43a9f7e756e47ca4e7ba5e6d9b558aea01 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 24 Feb 2016 18:03:42 -0800 Subject: [PATCH 06/18] prototype of defining regularization mesh within Regularization.py for constructing operators for regularization that are not true differential operators --- SimPEG/Mesh/DiffOperators.py | 129 +++++------ SimPEG/Regularization.py | 344 +++++++++++++++++++----------- tests/base/test_regularization.py | 121 ++++++----- 3 files changed, 344 insertions(+), 250 deletions(-) diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index 22b87edd..10ed2268 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -307,24 +307,28 @@ class DiffOperators(object): return BC _cellGradBC_list = 'neumann' + def _cellGradStencil(self): + BC = self.setCellGradBC(self._cellGradBC_list) + n = self.vnC + if(self.dim == 1): + G = ddxCellGrad(n[0], BC[0]) + elif(self.dim == 2): + G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0])) + G2 = sp.kron(ddxCellGrad(n[1], BC[1]), speye(n[0])) + G = sp.vstack((G1, G2), format="csr") + elif(self.dim == 3): + G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC[0])) + G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0])) + G3 = kron3(ddxCellGrad(n[2], BC[2]), speye(n[1]), speye(n[0])) + G = sp.vstack((G1, G2, G3), format="csr") + return G + def cellGrad(): doc = "The cell centered Gradient, takes you to cell faces." def fget(self): if(self._cellGrad is None): - BC = self.setCellGradBC(self._cellGradBC_list) - n = self.vnC - if(self.dim == 1): - G = ddxCellGrad(n[0], BC[0]) - elif(self.dim == 2): - G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC[0])) - G2 = sp.kron(ddxCellGrad(n[1], BC[1]), speye(n[0])) - G = sp.vstack((G1, G2), format="csr") - elif(self.dim == 3): - G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC[0])) - G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC[1]), speye(n[0])) - G3 = kron3(ddxCellGrad(n[2], BC[2]), speye(n[1]), speye(n[0])) - G = sp.vstack((G1, G2, G3), format="csr") + G = self._cellGradStencil() # Compute areas of cell faces & volumes S = self.area V = self.aveCC2F*self.vol # Average volume between adjacent cells @@ -361,19 +365,24 @@ class DiffOperators(object): _cellGradBC = None cellGradBC = property(**cellGradBC()) + def _cellGradxStencil(self): + BC = ['neumann', 'neumann'] + n = self.vnC + if(self.dim == 1): + G1 = ddxCellGrad(n[0], BC) + elif(self.dim == 2): + G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC)) + elif(self.dim == 3): + G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC)) + return G1 + + def cellGradx(): doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." def fget(self): if getattr(self, '_cellGradx', None) is None: - BC = ['neumann', 'neumann'] - n = self.vnC - if(self.dim == 1): - G1 = ddxCellGrad(n[0], BC) - elif(self.dim == 2): - G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC)) - elif(self.dim == 3): - G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC)) + G1 = self._cellGradxStencil() # Compute areas of cell faces & volumes V = self.aveCC2F*self.vol L = self.r(self.area/V, 'F','Fx', 'V') @@ -382,17 +391,22 @@ class DiffOperators(object): return locals() cellGradx = property(**cellGradx()) + def _cellGradyStencil(self): + if self.dim < 2: return None + BC = ['neumann', 'neumann'] + n = self.vnC + if(self.dim == 2): + G2 = sp.kron(ddxCellGrad(n[1], BC), speye(n[0])) + elif(self.dim == 3): + G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0])) + return G2 + def cellGrady(): doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." def fget(self): if self.dim < 2: return None if getattr(self, '_cellGrady', None) is None: - BC = ['neumann', 'neumann'] - n = self.vnC - if(self.dim == 2): - G2 = sp.kron(ddxCellGrad(n[1], BC), speye(n[0])) - elif(self.dim == 3): - G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0])) + G2 = self._cellGradyStencil() # Compute areas of cell faces & volumes V = self.aveCC2F*self.vol L = self.r(self.area/V, 'F','Fy', 'V') @@ -401,14 +415,19 @@ class DiffOperators(object): return locals() cellGrady = property(**cellGrady()) + def _cellGradzStencil(self): + if self.dim < 3: return None + BC = ['neumann', 'neumann'] + n = self.vnC + G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0])) + return G3 + def cellGradz(): doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions." def fget(self): if self.dim < 3: return None if getattr(self, '_cellGradz', None) is None: - BC = ['neumann', 'neumann'] - n = self.vnC - G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0])) + G3 = self._cellGradzStencil() # Compute areas of cell faces & volumes V = self.aveCC2F*self.vol L = self.r(self.area/V, 'F','Fz', 'V') @@ -565,56 +584,6 @@ class DiffOperators(object): return Pbc, Pin, Pout - def unitCellGradx(): - doc = """Cell centered Gradient in the x dimension used for - regularization. The gradient operator is square (nC-by-nC)""" - def fget(self): - if self.dim < 3: return None - if getattr(self, '_unitCellGradx', None) is None: - - n = self.vnC - gx = ddx(n[0]-1) - gx_square = sp.vstack((gx,gx[-1,:]*-1), format="csr") - - self._unitCellGradx = kron3(speye(n[2]), speye(n[1]), gx_square) - - return self._unitCellGradx - return locals() - unitCellGradx = property(**unitCellGradx()) - - def unitCellGrady(): - doc = """Cell centered Gradient in they dimension used for - regularization. The gradient operator is square (nC-by-nC)""" - def fget(self): - if self.dim < 3: return None - if getattr(self, '_unitCellGrady', None) is None: - - n = self.vnC - gy = ddx(n[1]-1) - gy_square = sp.vstack((gy,gy[-1,:]*-1), format="csr") - - self._unitCellGrady = kron3(speye(n[2]), gy_square, speye(n[0])) - - return self._unitCellGrady - return locals() - unitCellGrady = property(**unitCellGrady()) - - def unitCellGradz(): - doc = """Cell centered Gradient in they dimension used for - regularization. The gradient operator is square (nC-by-nC)""" - def fget(self): - if self.dim < 3: return None - if getattr(self, '_unitCellGradz', None) is None: - - n = self.vnC - gz = ddx(n[2]-1) - gz_square = sp.vstack((gz,gz[-1,:]*-1), format="csr") - - self._unitCellGradz = kron3( gz_square , speye(n[1]), speye(n[0])) - - return self._unitCellGradz - return locals() - unitCellGradz = property(**unitCellGradz()) # --------------- Averaging --------------------- diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index d6780c45..31013dc5 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -1,61 +1,35 @@ import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp -class BaseRegularization(object): - """ - **Base Regularization Class** +class RegularizationMesh(object): - This is used to regularize the model space:: - - reg = Regularization(mesh) - - """ - - __metaclass__ = Utils.SimPEGMetaClass - - counter = None - - mapPair = Maps.IdentityMap #: A SimPEG.Map Class - - mapping = None #: A SimPEG.Map instance. - mesh = None #: A SimPEG.Mesh instance. - mref = None #: Reference model. - - def __init__(self, mesh, mapping=None, indActive=None, **kwargs): - Utils.setKwargs(self, **kwargs) + def __init__(self, mesh, indActive=None): self.mesh = mesh - assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object." - self.mapping = mapping or self.mapPair(mesh) - self.mapping._assertMatchesPair(self.mapPair) self.indActive = indActive @property - def parent(self): - """This is the parent of the regularization.""" - return getattr(self,'_parent',None) - @parent.setter - def parent(self, p): - if getattr(self,'_parent',None) is not None: - print 'Regularization has switched to a new parent!' - self._parent = p + def vol(self): + if getattr(self, '_vol', None) is None: + self._vol = self._Pac.T * self.mesh.vol + return self._vol @property - def inv(self): return self.parent.inv - @property - def invProb(self): return self.parent - @property - def reg(self): return self - @property - def opt(self): return self.parent.opt - @property - def prob(self): return self.parent.prob - @property - def survey(self): return self.parent.survey - + def nC(self): + if getattr(self, '_nC', None) is None: + if self.indActive is None: + self._nC = self.mesh.nC + else: + if self.indActive.dtype == 'bool': + self._nC = sum(self.indActive) + else: + self._nC = len(self.indActive) # you shouldn't pass a vector of int 0, 1 's + return self._nC @property - def W(self): - """Full regularization weighting matrix W.""" - return self._Pac.T * sp.identity(self.mesh.nC) * self._Pac # or do we want sp.identity(self.mesh.nC) or even just Utils.Identity() ? + def dim(self): + if getattr(self, '_dim', None) is None: + self._dim = self.mesh.dim + return self._dim + @property def _Pac(self): @@ -95,7 +69,170 @@ class BaseRegularization(object): indActive_Fz = (self.mesh.aveFz2CC.T * self.indActive) == 1 self.__Pafz = Utils.speye(self.mesh.nFz)[:,indActive_Fz] return self.__Pafz + + @property + def aveFx2CC(self): + if getattr(self, '_aveFx2CC', None) is None: + self._aveFx2CC = self._Pac.T * self.mesh.aveFx2CC * self._Pafx + return self._aveFx2CC + + @property + def aveCC2Fx(self): + if getattr(self, '_aveCC2Fx', None) is None: + self._aveCC2Fx = Utils.sdiag(1./(self.aveFx2CC.T).sum(1)) * self.aveFx2CC.T + return self._aveCC2Fx + + @property + def aveFy2CC(self): + if getattr(self, '_aveFy2CC', None) is None: + self._aveFy2CC = self._Pac.T * self.mesh.aveFy2CC * self._Pafy + return self._aveFy2CC + + @property + def aveCC2Fy(self): + if getattr(self, '_aveCC2Fy', None) is None: + self._aveCC2Fy = Utils.sdiag(1./(self.aveFy2CC.T).sum(1)) * self.aveFy2CC.T + return self._aveCC2Fy + + @property + def aveFz2CC(self): + if getattr(self, '_aveFz2CC', None) is None: + self._aveFz2CC = self._Pac.T * self.mesh.aveFz2CC * self._Pafz + return self._aveFz2CC + + @property + def aveCC2Fz(self): + if getattr(self, '_aveCC2Fz', None) is None: + self._aveCC2Fz = Utils.sdiag(1./(self.aveFz2CC.T).sum(1)) * self.aveFz2CC.T + return self._aveCC2Fz + + @property + def cellGradx(self): + if getattr(self, '_cellGradx', None) is None: + self._cellGradx = self._Pafx.T * self.mesh.cellGradx * self._Pac + return self._cellGradx + + @property + def cellGrady(self): + if getattr(self, '_cellGrady', None) is None: + self._cellGrady = self._Pafy.T * self.mesh.cellGrady * self._Pac + return self._cellGrady + + @property + def cellGradz(self): + if getattr(self, '_cellGradz', None) is None: + self._cellGradz = self._Pafz.T * self.mesh.cellGradz * self._Pac + return self._cellGradz + @property + def faceDivx(self): + if getattr(self, '_faceDivx', None) is None: + self._faceDivx = self._Pac.T * self.mesh.faceDivx * self._Pafx + return self._faceDivx + + @property + def faceDivy(self): + if getattr(self, '_faceDivy', None) is None: + self._faceDivy = self._Pac.T * self.mesh.faceDivy * self._Pafy + return self._faceDivy + + @property + def faceDivz(self): + if getattr(self, '_faceDivz', None) is None: + self._faceDivz = self._Pac.T * self.mesh.faceDivz * self._Pafz + return self._faceDivz + + @property + def cellGradxStencil(self): + """Cell centered Gradient in the x dimension used for + regularization. The gradient operator is square (nC-by-nC)""" + + # if self.dim < 3: return None + if getattr(self, '_cellGradxStencil', None) is None: + + self._cellGradxStencil = self._Pafx.T * self.mesh._cellGradxStencil() * self._Pac + return self._cellGradxStencil + + @property + def cellGradyStencil(self): + """Cell centered Gradient in the x dimension used for + regularization. The gradient operator is square (nC-by-nC)""" + + if self.dim < 2: return None + if getattr(self, '_cellGradyStencil', None) is None: + + self._cellGradyStencil = self._Pafy.T * self.mesh._cellGradyStencil() * self._Pac + return self._cellGradyStencil + + @property + def cellGradzStencil(self): + """Cell centered Gradient in the x dimension used for + regularization. The gradient operator is square (nC-by-nC)""" + + if self.dim < 3: return None + if getattr(self, '_cellGradzStencil', None) is None: + + self._cellGradzStencil = self._Pafz.T * self.mesh._cellGradzStencil() * self._Pac + return self._cellGradzStencil + + +class BaseRegularization(object): + """ + **Base Regularization Class** + + This is used to regularize the model space:: + + reg = Regularization(mesh) + + """ + + __metaclass__ = Utils.SimPEGMetaClass + + counter = None + + mapPair = Maps.IdentityMap #: A SimPEG.Map Class + + mapping = None #: A SimPEG.Map instance. + mesh = None #: A SimPEG.Mesh instance. + mref = None #: Reference model. + + def __init__(self, mesh, mapping=None, indActive=None, **kwargs): + Utils.setKwargs(self, **kwargs) + assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object." + self.regmesh = RegularizationMesh(mesh,indActive) + self.mapping = mapping or self.mapPair(mesh) + self.mapping._assertMatchesPair(self.mapPair) + self.indActive = indActive + + @property + def parent(self): + """This is the parent of the regularization.""" + return getattr(self,'_parent',None) + @parent.setter + def parent(self, p): + if getattr(self,'_parent',None) is not None: + print 'Regularization has switched to a new parent!' + self._parent = p + + @property + def inv(self): return self.parent.inv + @property + def invProb(self): return self.parent + @property + def reg(self): return self + @property + def opt(self): return self.parent.opt + @property + def prob(self): return self.parent.prob + @property + def survey(self): return self.parent.survey + + + @property + def W(self): + """Full regularization weighting matrix W.""" + return sp.identity(self.regmesh.nC) + # self.regmesh._Pac.T * sp.identity(self.regmesh.nC) * self.regmesh._Pac # or do we want sp.identity(self.mesh.nC) or even just Utils.Identity() ? @Utils.timeIt def eval(self, m): @@ -151,7 +288,6 @@ class BaseRegularization(object): return mD.T * ( self.W.T * ( self.W * ( mD * v) ) ) - class Tikhonov(BaseRegularization): """ """ @@ -165,66 +301,58 @@ class Tikhonov(BaseRegularization): alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction") def __init__(self, mesh, mapping=None, indActive = None, **kwargs): - BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs) - self.indActive = indActive + BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) @property def Ws(self): """Regularization matrix Ws""" if getattr(self,'_Ws', None) is None: - Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5) - self._Ws = self._Pac.T * Ws * self._Pac + self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s)**0.5) return self._Ws @property def Wx(self): """Regularization matrix Wx""" if getattr(self, '_Wx', None) is None: - Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFx].T*self.mesh.vol - Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx - self._Wx = self._Pafx.T*Wx*self.self._Pac + Ave_x_vol = self.regmesh.aveCC2Fx * self.regmesh.vol + self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.regmesh.cellGradx return self._Wx @property def Wy(self): """Regularization matrix Wy""" if getattr(self, '_Wy', None) is None: - Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFx:np.sum(self.mesh.vnF[:2])].T*self.mesh.vol - Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady - self._Wy = self._Pafy.T*Wy*self._Pac + Ave_y_vol = self.regmesh.aveCC2Fy * self.regmesh.vol + self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.regmesh.cellGrady return self._Wy @property def Wz(self): """Regularization matrix Wz""" if getattr(self, '_Wz', None) is None: - Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.vnF[:2]):].T*self.mesh.vol - Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz - self._Wz = self._Pafz.T*Wz*self._Pac + Ave_z_vol = self.regmesh.aveCC2Fz * self.regmesh.vol + self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.regmesh.cellGradz return self._Wz @property def Wxx(self): """Regularization matrix Wxx""" if getattr(self, '_Wxx', None) is None: - Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx - self._Wxx = self._Pac.T*Wxx*self._Pac + self._Wxx = Utils.sdiag((self.regmesh.vol*self.alpha_xx)**0.5)*self.regmesh.faceDivx*self.regmesh.cellGradx return self._Wxx @property def Wyy(self): """Regularization matrix Wyy""" if getattr(self, '_Wyy', None) is None: - Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady - self._Wyy = self._Pac.T*self._Wyy*self._Pac + self._Wyy = Utils.sdiag((self.regmesh.vol*self.alpha_yy)**0.5)*self.regmesh.faceDivy*self.regmesh.cellGrady return self._Wyy @property def Wzz(self): """Regularization matrix Wzz""" if getattr(self, '_Wzz', None) is None: - Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz - self._Wzz = self._Pac.T*Wzz*self._Pac + self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDivz*self.regmesh.cellGradz return self._Wzz @property @@ -232,9 +360,9 @@ class Tikhonov(BaseRegularization): """Full smoothness regularization matrix W""" if getattr(self, '_Wsmooth', None) is None: wlist = (self.Wx, self.Wxx) - if self.mesh.dim > 1: + if self.regmesh.dim > 1: wlist += (self.Wy, self.Wyy) - if self.mesh.dim > 2: + if self.regmesh.dim > 2: wlist += (self.Wz, self.Wzz) self._Wsmooth = sp.vstack(wlist) return self._Wsmooth @@ -301,35 +429,35 @@ class Simple(BaseRegularization): 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") - def __init__(self, mesh, mapping=None, **kwargs): - BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs) + def __init__(self, mesh, mapping=None, indActive=None, **kwargs): + BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) @property def Ws(self): """Regularization matrix Ws""" if getattr(self,'_Ws', None) is None: - self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5) + self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s)**0.5) return self._Ws @property def Wx(self): """Regularization matrix Wx""" if getattr(self, '_Wx', None) is None: - self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x)**0.5)*self.mesh.unitCellGradx + self._Wx = Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x)**0.5)*self.regmesh.cellGradxStencil return self._Wx @property def Wy(self): """Regularization matrix Wy""" if getattr(self, '_Wy', None) is None: - self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y)**0.5)*self.mesh.unitCellGrady + self._Wy = Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol * self.alpha_y)**0.5)*self.regmesh.cellGradyStencil return self._Wy @property def Wz(self): """Regularization matrix Wz""" if getattr(self, '_Wz', None) is None: - self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z)**0.5)*self.mesh.unitCellGradz + self._Wz = Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z)**0.5)*self.regmesh.cellGradzStencil return self._Wz @property @@ -337,9 +465,9 @@ class Simple(BaseRegularization): """Full smoothness regularization matrix W""" if getattr(self, '_Wsmooth', None) is None: wlist = (self.Wx,) - if self.mesh.dim > 1: + if self.regmesh.dim > 1: wlist += (self.Wy,) - if self.mesh.dim > 2: + if self.regmesh.dim > 2: wlist += (self.Wz,) self._Wsmooth = sp.vstack(wlist) return self._Wsmooth @@ -352,25 +480,16 @@ class Simple(BaseRegularization): self._W = sp.vstack(wlist) return self._W - @Utils.timeIt - def _evalSmall(self, m): - r = self.W * ( self.mapping * (m - self.mref) ) - return 0.5*r.dot(r) @Utils.timeIt - def _evalSmooth(self, m): + def eval(self, m): if self.smoothModel == True: r1 = self.Wsmooth * ( self.mapping * (m) ) r2 = self.Ws * ( self.mapping * (m - self.mref) ) return 0.5*(r1.dot(r1)+r2.dot(r2)) - else: - return None - - @Utils.timeIt - def eval(self, m): - phim = self._evalSmall(m) - if self.smoothModel is True: - phim += self._evalSmooth(m) + elif self.smoothModel == False: + r = self.W * ( self.mapping * (m - self.mref) ) + return 0.5*r.dot(r) return phim @@ -417,34 +536,15 @@ class SparseRegularization(Simple): qy = 2. qz = 2. - def __init__(self, mesh, mapping=None, **kwargs): - Simple.__init__(self, mesh, mapping=mapping, **kwargs) + def __init__(self, mesh, mapping=None, indActive=None, **kwargs): + Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) - @property - def Wsmooth(self): - """Full smoothness regularization matrix W""" - if getattr(self, '_Wsmooth', None) is None: - wlist = (self.Wx, self.Wxx) - if self.mesh.dim > 1: - wlist += (self.Wy, self.Wyy) - if self.mesh.dim > 2: - wlist += (self.Wz, self.Wzz) - self._Wsmooth = sp.vstack(wlist) - return self._Wsmooth - - @property - def W(self): - """Full regularization matrix W""" - if getattr(self, '_W', None) is None: - wlist = (self.Ws, self.Wsmooth) - self._W = sp.vstack(wlist) - return self._W @property def Ws(self): """Regularization matrix Ws""" if getattr(self, 'm', None) is None: - self.Rs = Utils.speye(self.mesh.nC) + self.Rs = Utils.speye(self.regmesh.nC) else: f_m = self.m @@ -452,7 +552,7 @@ class SparseRegularization(Simple): #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs)) self.Rs = Utils.sdiag( self.rs ) - self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s*self.gamma)**0.5)*self.Rs + self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma)**0.5)*self.Rs return self._Ws @@ -461,15 +561,15 @@ class SparseRegularization(Simple): """Regularization matrix Wx""" if getattr(self, 'm', None) is None: - self.Rx = Utils.speye(self.mesh.unitCellGradx.shape[0]) + self.Rx = Utils.speye(self.regmesh.cellGradxStencil.shape[0]) else: - f_m = self.mesh.unitCellGradx * self.m + f_m = self.regmesh.cellGradxStencil * self.m self.rx = self.R( f_m , self.qx, self.eps) self.Rx = Utils.sdiag( self.rx ) if getattr(self, '_Wx', None) is None: - self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x*self.gamma)**0.5)*self.Rx*self.mesh.unitCellGradx + self._Wx = Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma)**0.5)*self.Rx*self.regmesh.cellGradxStencil return self._Wx @property @@ -477,15 +577,15 @@ class SparseRegularization(Simple): """Regularization matrix Wy""" if getattr(self, 'm', None) is None: - self.Ry = Utils.speye(self.mesh.unitCellGrady.shape[0]) + self.Ry = Utils.speye(self.regmesh.cellGradyStencil.shape[0]) else: - f_m = self.mesh.unitCellGrady * self.m + f_m = self.regmesh.cellGradyStencil * self.m self.ry = self.R( f_m , self.qy, self.eps) self.Ry = Utils.sdiag( self.ry ) if getattr(self, '_Wy', None) is None: - self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y*self.gamma)**0.5)*self.Ry*self.mesh.unitCellGrady + self._Wy = Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma)**0.5)*self.Ry*self.regmesh.cellGradyStencil return self._Wy @property @@ -493,15 +593,15 @@ class SparseRegularization(Simple): """Regularization matrix Wz""" if getattr(self, 'm', None) is None: - self.Rz = Utils.speye(self.mesh.unitCellGradz.shape[0]) + self.Rz = Utils.speye(self.regmesh.cellGradzStencil.shape[0]) else: - f_m = self.mesh.unitCellGradz * self.m + f_m = self.regmesh.cellGradzStencil * self.m self.rz = self.R( f_m , self.qz, self.eps) self.Rz = Utils.sdiag( self.rz ) if getattr(self, '_Wz', None) is None: - self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z*self.gamma)**0.5)*self.Rz*self.mesh.unitCellGradz + self._Wz = Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma)**0.5)*self.Rz*self.regmesh.cellGradzStencil return self._Wz diff --git a/tests/base/test_regularization.py b/tests/base/test_regularization.py index 050c46ac..614d3158 100644 --- a/tests/base/test_regularization.py +++ b/tests/base/test_regularization.py @@ -5,6 +5,8 @@ from scipy.sparse.linalg import dsolve import inspect TOL = 1e-20 +testReg = True +testRegMesh = True class RegularizationTests(unittest.TestCase): @@ -16,69 +18,92 @@ class RegularizationTests(unittest.TestCase): mesh3 = Mesh.TensorMesh([hx, hy, hz]) self.meshlist = [mesh1,mesh2, mesh3] - def test_regularization(self): - for R in dir(Regularization): - r = getattr(Regularization, R) - if not inspect.isclass(r): continue - if not issubclass(r, Regularization.BaseRegularization): - continue + if testReg: + def test_regularization(self): + for R in dir(Regularization): + r = getattr(Regularization, R) + if not inspect.isclass(r): continue + if not issubclass(r, Regularization.BaseRegularization): + continue - for i, mesh in enumerate(self.meshlist): + for i, mesh in enumerate(self.meshlist): - print 'Testing %iD'%mesh.dim + print 'Testing %iD'%mesh.dim - mapping = r.mapPair(mesh) - reg = r(mesh, mapping=mapping) - m = np.random.rand(mapping.nP) - reg.mref = np.ones_like(m)*np.mean(m) + mapping = r.mapPair(mesh) + reg = r(mesh, mapping=mapping) + m = np.random.rand(mapping.nP) + reg.mref = np.ones_like(m)*np.mean(m) - print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref) - passed = reg.eval(reg.mref) < TOL - self.assertTrue(passed) + print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref) + passed = reg.eval(reg.mref) < TOL + self.assertTrue(passed) - print 'Check:', R - passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False) - self.assertTrue(passed) + print 'Check:', R + passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False) + self.assertTrue(passed) - print 'Check 2 Deriv:', R - passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False) - self.assertTrue(passed) + print 'Check 2 Deriv:', R + passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False) + self.assertTrue(passed) - def test_regularization_ActiveCells(self): - for R in dir(Regularization): - r = getattr(Regularization, R) - if not inspect.isclass(r): continue - if not issubclass(r, Regularization.BaseRegularization): - continue + def test_regularization_ActiveCells(self): + for R in dir(Regularization): + r = getattr(Regularization, R) + if not inspect.isclass(r): continue + if not issubclass(r, Regularization.BaseRegularization): + continue - for i, mesh in enumerate(self.meshlist): + for i, mesh in enumerate(self.meshlist): - print 'Testing Active Cells %iD'%(mesh.dim) + print 'Testing Active Cells %iD'%(mesh.dim) - if mesh.dim == 1: - indAct = Utils.mkvc(mesh.gridCC <= 0.8) - elif mesh.dim == 2: - indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5) - elif mesh.dim == 3: - indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5) + if mesh.dim == 1: + indAct = Utils.mkvc(mesh.gridCC <= 0.8) + elif mesh.dim == 2: + indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5) + elif mesh.dim == 3: + indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5) - mapping = Maps.IdentityMap(nP=indAct.nonzero()[0].size) + mapping = Maps.IdentityMap(nP=indAct.nonzero()[0].size) - reg = r(mesh, mapping=mapping, indActive=indAct) - m = np.random.rand(mesh.nC)[indAct] - reg.mref = np.ones_like(m)*np.mean(m) + reg = r(mesh, mapping=mapping, indActive=indAct) + m = np.random.rand(mesh.nC)[indAct] + reg.mref = np.ones_like(m)*np.mean(m) - print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref) - passed = reg.eval(reg.mref) < TOL - self.assertTrue(passed) + print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref) + passed = reg.eval(reg.mref) < TOL + self.assertTrue(passed) - print 'Check:', R - passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False) - self.assertTrue(passed) + print 'Check:', R + passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False) + self.assertTrue(passed) - print 'Check 2 Deriv:', R - passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False) - self.assertTrue(passed) + print 'Check 2 Deriv:', R + passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False) + self.assertTrue(passed) + + if testRegMesh: + def test_regularizationMesh(self): + + for i, mesh in enumerate(self.meshlist): + + print 'Testing %iD'%mesh.dim + + # mapping = r.mapPair(mesh) + # reg = r(mesh, mapping=mapping) + # m = np.random.rand(mapping.nP) + + if mesh.dim == 1: + indAct = Utils.mkvc(mesh.gridCC <= 0.8) + elif mesh.dim == 2: + indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5) + elif mesh.dim == 3: + indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5) + + regmesh = Regularization.RegularizationMesh(mesh, indActive=indAct) + + assert (regmesh.vol == mesh.vol[indAct]).all() if __name__ == '__main__': From e3af1fd94e9828caa57a931074fa01a91a3ebf1b Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 24 Feb 2016 20:28:09 -0800 Subject: [PATCH 07/18] convert indActive to a bool if an integer list is provided --- SimPEG/Regularization.py | 10 ++++++---- tests/base/test_regularization.py | 15 ++++++++------- 2 files changed, 14 insertions(+), 11 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 31013dc5..d0fc21b1 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -4,6 +4,7 @@ class RegularizationMesh(object): def __init__(self, mesh, indActive=None): self.mesh = mesh + assert indActive is None or indActive.dtype == 'bool', 'indActive needs to be None or a bool' self.indActive = indActive @property @@ -18,10 +19,7 @@ class RegularizationMesh(object): if self.indActive is None: self._nC = self.mesh.nC else: - if self.indActive.dtype == 'bool': - self._nC = sum(self.indActive) - else: - self._nC = len(self.indActive) # you shouldn't pass a vector of int 0, 1 's + self._nC = sum(self.indActive) return self._nC @property @@ -199,6 +197,10 @@ class BaseRegularization(object): def __init__(self, mesh, mapping=None, indActive=None, **kwargs): Utils.setKwargs(self, **kwargs) assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object." + if indActive is not None and indActive.dtype != 'bool': + tmp = indActive + indActive = np.zeros(mesh.nC, dtype=bool) + indActive[tmp] = True self.regmesh = RegularizationMesh(mesh,indActive) self.mapping = mapping or self.mapPair(mesh) self.mapping._assertMatchesPair(self.mapPair) diff --git a/tests/base/test_regularization.py b/tests/base/test_regularization.py index 614d3158..52cef349 100644 --- a/tests/base/test_regularization.py +++ b/tests/base/test_regularization.py @@ -59,17 +59,18 @@ class RegularizationTests(unittest.TestCase): print 'Testing Active Cells %iD'%(mesh.dim) if mesh.dim == 1: - indAct = Utils.mkvc(mesh.gridCC <= 0.8) + indActive = Utils.mkvc(mesh.gridCC <= 0.8) elif mesh.dim == 2: - indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5) + indActive = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5) elif mesh.dim == 3: - indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5) + indActive = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5) - mapping = Maps.IdentityMap(nP=indAct.nonzero()[0].size) + mapping = Maps.IdentityMap(nP=indActive.nonzero()[0].size) - reg = r(mesh, mapping=mapping, indActive=indAct) - m = np.random.rand(mesh.nC)[indAct] - reg.mref = np.ones_like(m)*np.mean(m) + for indAct in [indActive, indActive.nonzero()[0]]: # test both bool and integers + reg = r(mesh, mapping=mapping, indActive=indAct) + m = np.random.rand(mesh.nC)[indAct] + reg.mref = np.ones_like(m)*np.mean(m) print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref) passed = reg.eval(reg.mref) < TOL From 63bf8b9e4d6d16b7d3f3d9271246f86b337df626 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Wed, 24 Feb 2016 20:59:51 -0800 Subject: [PATCH 08/18] Add linear survey --- SimPEG/Survey.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 47a88ae2..cc877c69 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -375,3 +375,11 @@ class BaseSurvey(object): self.dobs = self.dtrue+noise self.std = self.dobs*0 + std return self.dobs + +class LinearSurvey(BaseSurvey): + def projectFields(self, u): + return u + + @property + def nD(self): + return self.prob.G.shape[0] \ No newline at end of file From 6c33455d15ee6d937029ecb7b7ccc3306f09b0c4 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Thu, 25 Feb 2016 08:56:18 -0800 Subject: [PATCH 09/18] update Directive for sparse norm --- SimPEG/Directives.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 476c5916..9ba38e4e 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -317,7 +317,7 @@ class update_IRLS(InversionDirective): self.reg.m = self.invProb.curModel # 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.prob.mesh.nC))**2. + diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() * (self.reg.mapping * np.ones(self.reg.m.size))**2. PC = Utils.sdiag(diagA**-1.) self.opt.approxHinv = PC From 5e5c7ba0fb408645eb5c2fa85340b3af75f01098 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 1 Mar 2016 17:31:37 -0800 Subject: [PATCH 10/18] docs for regmesh, cellGrad--> cellDiff, faceDiv--> faceDiff for regmesh --- SimPEG/Regularization.py | 238 ++++++++++++++++++++++++++++----------- 1 file changed, 174 insertions(+), 64 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index d0fc21b1..acd6c68d 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -1,6 +1,16 @@ import Utils, Maps, Mesh, numpy as np, scipy.sparse as sp class RegularizationMesh(object): + """ + **Regularization Mesh** + + This contains the operators used in the regularization. Note that these + are not necessarily true differential operators, but are constructed from + a SimPEG Mesh. + + :param Mesh mesh: problem mesh + :param numpy.array indActive: bool array, size nC, that is True where we have active cells. Used to reduce the operators so we regularize only on active cells + """ def __init__(self, mesh, indActive=None): self.mesh = mesh @@ -9,12 +19,22 @@ class RegularizationMesh(object): @property def vol(self): + """ + reduced volume vector + :rtype: numpy.array + :return: reduced cell volume + """ if getattr(self, '_vol', None) is None: self._vol = self._Pac.T * self.mesh.vol return self._vol @property def nC(self): + """ + reduced number of cells + :rtype: int + :return: number of cells being regularized + """ if getattr(self, '_nC', None) is None: if self.indActive is None: self._nC = self.mesh.nC @@ -24,6 +44,11 @@ class RegularizationMesh(object): @property def dim(self): + """ + dimension of regularization mesh (1D, 2D, 3D) + :rtype: int + :return: dimension + """ if getattr(self, '_dim', None) is None: self._dim = self.mesh.dim return self._dim @@ -31,6 +56,11 @@ class RegularizationMesh(object): @property def _Pac(self): + """ + projection matrix that takes from the reduced space of active cells to full modelling space (ie. nC x nindActive) + :rtype: scipy.sparse.csr_matrix + :return: active cell projection matrix + """ if getattr(self, '__Pac', None) is None: if self.indActive is None: self.__Pac = Utils.speye(self.mesh.nC) @@ -40,6 +70,11 @@ class RegularizationMesh(object): @property def _Pafx(self): + """ + projection matrix that takes from the reduced space of active x-faces to full modelling space (ie. nFx x nindActive_Fx ) + :rtype: scipy.sparse.csr_matrix + :return: active face-x projection matrix + """ if getattr(self, '__Pafx', None) is None: if self.indActive is None: self.__Pafx = Utils.speye(self.mesh.nFx) @@ -50,6 +85,11 @@ class RegularizationMesh(object): @property def _Pafy(self): + """ + projection matrix that takes from the reduced space of active y-faces to full modelling space (ie. nFy x nindActive_Fy ) + :rtype: scipy.sparse.csr_matrix + :return: active face-y projection matrix + """ if getattr(self, '__Pafy', None) is None: if self.indActive is None: self.__Pafy = Utils.speye(self.mesh.nFy) @@ -60,6 +100,11 @@ class RegularizationMesh(object): @property def _Pafz(self): + """ + projection matrix that takes from the reduced space of active z-faces to full modelling space (ie. nFz x nindActive_Fz ) + :rtype: scipy.sparse.csr_matrix + :return: active face-z projection matrix + """ if getattr(self, '__Pafz', None) is None: if self.indActive is None: self.__Pafz = Utils.speye(self.mesh.nFz) @@ -70,108 +115,173 @@ class RegularizationMesh(object): @property def aveFx2CC(self): + """ + averaging from active cell centers to active x-faces + :rtype: scipy.sparse.csr_matrix + :return: averaging from active cell centers to active x-faces + """ if getattr(self, '_aveFx2CC', None) is None: self._aveFx2CC = self._Pac.T * self.mesh.aveFx2CC * self._Pafx return self._aveFx2CC @property def aveCC2Fx(self): + """ + averaging from active x-faces to active cell centers + :rtype: scipy.sparse.csr_matrix + :return: averaging matrix from active x-faces to active cell centers + """ if getattr(self, '_aveCC2Fx', None) is None: self._aveCC2Fx = Utils.sdiag(1./(self.aveFx2CC.T).sum(1)) * self.aveFx2CC.T return self._aveCC2Fx @property def aveFy2CC(self): + """ + averaging from active cell centers to active y-faces + :rtype: scipy.sparse.csr_matrix + :return: averaging from active cell centers to active y-faces + """ if getattr(self, '_aveFy2CC', None) is None: self._aveFy2CC = self._Pac.T * self.mesh.aveFy2CC * self._Pafy return self._aveFy2CC @property def aveCC2Fy(self): + """ + averaging from active y-faces to active cell centers + :rtype: scipy.sparse.csr_matrix + :return: averaging matrix from active y-faces to active cell centers + """ if getattr(self, '_aveCC2Fy', None) is None: self._aveCC2Fy = Utils.sdiag(1./(self.aveFy2CC.T).sum(1)) * self.aveFy2CC.T return self._aveCC2Fy @property def aveFz2CC(self): + """ + averaging from active cell centers to active z-faces + :rtype: scipy.sparse.csr_matrix + :return: averaging from active cell centers to active z-faces + """ if getattr(self, '_aveFz2CC', None) is None: self._aveFz2CC = self._Pac.T * self.mesh.aveFz2CC * self._Pafz return self._aveFz2CC @property def aveCC2Fz(self): + """ + averaging from active z-faces to active cell centers + :rtype: scipy.sparse.csr_matrix + :return: averaging matrix from active z-faces to active cell centers + """ if getattr(self, '_aveCC2Fz', None) is None: self._aveCC2Fz = Utils.sdiag(1./(self.aveFz2CC.T).sum(1)) * self.aveFz2CC.T return self._aveCC2Fz @property - def cellGradx(self): - if getattr(self, '_cellGradx', None) is None: - self._cellGradx = self._Pafx.T * self.mesh.cellGradx * self._Pac - return self._cellGradx + def cellDiffx(self): + """ + cell centered difference in the x-direction + :rtype: scipy.sparse.csr_matrix + :return: differencing matrix for active cells in the x-direction + """ + if getattr(self, '_cellDiffx', None) is None: + self._cellDiffx = self._Pafx.T * self.mesh.cellGradx * self._Pac + return self._cellDiffx @property - def cellGrady(self): - if getattr(self, '_cellGrady', None) is None: - self._cellGrady = self._Pafy.T * self.mesh.cellGrady * self._Pac - return self._cellGrady + def cellDiffy(self): + """ + cell centered difference in the y-direction + :rtype: scipy.sparse.csr_matrix + :return: differencing matrix for active cells in the y-direction + """ + if getattr(self, '_cellDiffy', None) is None: + self._cellDiffy = self._Pafy.T * self.mesh.cellGrady * self._Pac + return self._cellDiffy @property - def cellGradz(self): - if getattr(self, '_cellGradz', None) is None: - self._cellGradz = self._Pafz.T * self.mesh.cellGradz * self._Pac - return self._cellGradz + def cellDiffz(self): + """ + cell centered difference in the z-direction + :rtype: scipy.sparse.csr_matrix + :return: differencing matrix for active cells in the z-direction + """ + if getattr(self, '_cellDiffz', None) is None: + self._cellDiffz = self._Pafz.T * self.mesh.cellGradz * self._Pac + return self._cellDiffz @property - def faceDivx(self): - if getattr(self, '_faceDivx', None) is None: - self._faceDivx = self._Pac.T * self.mesh.faceDivx * self._Pafx - return self._faceDivx + def faceDiffx(self): + """ + x-face differences + :rtype: scipy.sparse.csr_matrix + :return: differencing matrix for active faces in the x-direction + """ + if getattr(self, '_faceDiffx', None) is None: + self._faceDiffx = self._Pac.T * self.mesh.faceDivx * self._Pafx + return self._faceDiffx @property - def faceDivy(self): - if getattr(self, '_faceDivy', None) is None: - self._faceDivy = self._Pac.T * self.mesh.faceDivy * self._Pafy - return self._faceDivy + def faceDiffy(self): + """ + y-face differences + :rtype: scipy.sparse.csr_matrix + :return: differencing matrix for active faces in the y-direction + """ + if getattr(self, '_faceDiffy', None) is None: + self._faceDiffy = self._Pac.T * self.mesh.faceDivy * self._Pafy + return self._faceDiffy @property - def faceDivz(self): - if getattr(self, '_faceDivz', None) is None: - self._faceDivz = self._Pac.T * self.mesh.faceDivz * self._Pafz - return self._faceDivz + def faceDiffz(self): + """ + z-face differences + :rtype: scipy.sparse.csr_matrix + :return: differencing matrix for active faces in the z-direction + """ + if getattr(self, '_faceDiffz', None) is None: + self._faceDiffz = self._Pac.T * self.mesh.faceDivz * self._Pafz + return self._faceDiffz @property - def cellGradxStencil(self): - """Cell centered Gradient in the x dimension used for - regularization. The gradient operator is square (nC-by-nC)""" + def cellDiffxStencil(self): + """ + cell centered difference stencil (no cell lengths include) in the x-direction + :rtype: scipy.sparse.csr_matrix + :return: differencing matrix for active cells in the x-direction + """ + if getattr(self, '_cellDiffxStencil', None) is None: - # if self.dim < 3: return None - if getattr(self, '_cellGradxStencil', None) is None: - - self._cellGradxStencil = self._Pafx.T * self.mesh._cellGradxStencil() * self._Pac - return self._cellGradxStencil + self._cellDiffxStencil = self._Pafx.T * self.mesh._cellGradxStencil() * self._Pac + return self._cellDiffxStencil @property - def cellGradyStencil(self): - """Cell centered Gradient in the x dimension used for - regularization. The gradient operator is square (nC-by-nC)""" - + def cellDiffyStencil(self): + """ + cell centered difference stencil (no cell lengths include) in the y-direction + :rtype: scipy.sparse.csr_matrix + :return: differencing matrix for active cells in the y-direction + """ if self.dim < 2: return None - if getattr(self, '_cellGradyStencil', None) is None: + if getattr(self, '_cellDiffyStencil', None) is None: - self._cellGradyStencil = self._Pafy.T * self.mesh._cellGradyStencil() * self._Pac - return self._cellGradyStencil + self._cellDiffyStencil = self._Pafy.T * self.mesh._cellGradyStencil() * self._Pac + return self._cellDiffyStencil @property - def cellGradzStencil(self): - """Cell centered Gradient in the x dimension used for - regularization. The gradient operator is square (nC-by-nC)""" - + def cellDiffzStencil(self): + """ + cell centered difference stencil (no cell lengths include) in the y-direction + :rtype: scipy.sparse.csr_matrix + :return: differencing matrix for active cells in the y-direction + """ if self.dim < 3: return None - if getattr(self, '_cellGradzStencil', None) is None: + if getattr(self, '_cellDiffzStencil', None) is None: - self._cellGradzStencil = self._Pafz.T * self.mesh._cellGradzStencil() * self._Pac - return self._cellGradzStencil + self._cellDiffzStencil = self._Pafz.T * self.mesh._cellGradzStencil() * self._Pac + return self._cellDiffzStencil class BaseRegularization(object): @@ -317,7 +427,7 @@ class Tikhonov(BaseRegularization): """Regularization matrix Wx""" if getattr(self, '_Wx', None) is None: Ave_x_vol = self.regmesh.aveCC2Fx * self.regmesh.vol - self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.regmesh.cellGradx + self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.regmesh.cellDiffx return self._Wx @property @@ -325,7 +435,7 @@ class Tikhonov(BaseRegularization): """Regularization matrix Wy""" if getattr(self, '_Wy', None) is None: Ave_y_vol = self.regmesh.aveCC2Fy * self.regmesh.vol - self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.regmesh.cellGrady + self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.regmesh.cellDiffy return self._Wy @property @@ -333,28 +443,28 @@ class Tikhonov(BaseRegularization): """Regularization matrix Wz""" if getattr(self, '_Wz', None) is None: Ave_z_vol = self.regmesh.aveCC2Fz * self.regmesh.vol - self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.regmesh.cellGradz + self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.regmesh.cellDiffz return self._Wz @property def Wxx(self): """Regularization matrix Wxx""" if getattr(self, '_Wxx', None) is None: - self._Wxx = Utils.sdiag((self.regmesh.vol*self.alpha_xx)**0.5)*self.regmesh.faceDivx*self.regmesh.cellGradx + self._Wxx = Utils.sdiag((self.regmesh.vol*self.alpha_xx)**0.5)*self.regmesh.faceDiffx*self.regmesh.cellDiffx return self._Wxx @property def Wyy(self): """Regularization matrix Wyy""" if getattr(self, '_Wyy', None) is None: - self._Wyy = Utils.sdiag((self.regmesh.vol*self.alpha_yy)**0.5)*self.regmesh.faceDivy*self.regmesh.cellGrady + self._Wyy = Utils.sdiag((self.regmesh.vol*self.alpha_yy)**0.5)*self.regmesh.faceDiffy*self.regmesh.cellDiffy return self._Wyy @property def Wzz(self): """Regularization matrix Wzz""" if getattr(self, '_Wzz', None) is None: - self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDivz*self.regmesh.cellGradz + self._Wzz = Utils.sdiag((self.regmesh.vol*self.alpha_zz)**0.5)*self.regmesh.faceDiffz*self.regmesh.cellDiffz return self._Wzz @property @@ -445,21 +555,21 @@ class Simple(BaseRegularization): 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)**0.5)*self.regmesh.cellGradxStencil + self._Wx = Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x)**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)**0.5)*self.regmesh.cellGradyStencil + self._Wy = Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol * self.alpha_y)**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)**0.5)*self.regmesh.cellGradzStencil + self._Wz = Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z)**0.5)*self.regmesh.cellDiffzStencil return self._Wz @property @@ -563,15 +673,15 @@ class SparseRegularization(Simple): """Regularization matrix Wx""" if getattr(self, 'm', None) is None: - self.Rx = Utils.speye(self.regmesh.cellGradxStencil.shape[0]) + self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0]) else: - f_m = self.regmesh.cellGradxStencil * self.m + f_m = self.regmesh.cellDiffxStencil * self.m self.rx = self.R( f_m , self.qx, self.eps) self.Rx = Utils.sdiag( self.rx ) if getattr(self, '_Wx', None) is None: - self._Wx = Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma)**0.5)*self.Rx*self.regmesh.cellGradxStencil + self._Wx = Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma)**0.5)*self.Rx*self.regmesh.cellDiffxStencil return self._Wx @property @@ -579,15 +689,15 @@ class SparseRegularization(Simple): """Regularization matrix Wy""" if getattr(self, 'm', None) is None: - self.Ry = Utils.speye(self.regmesh.cellGradyStencil.shape[0]) + self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0]) else: - f_m = self.regmesh.cellGradyStencil * self.m + f_m = self.regmesh.cellDiffyStencil * self.m self.ry = self.R( f_m , self.qy, self.eps) self.Ry = Utils.sdiag( self.ry ) if getattr(self, '_Wy', None) is None: - self._Wy = Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma)**0.5)*self.Ry*self.regmesh.cellGradyStencil + self._Wy = Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma)**0.5)*self.Ry*self.regmesh.cellDiffyStencil return self._Wy @property @@ -595,15 +705,15 @@ class SparseRegularization(Simple): """Regularization matrix Wz""" if getattr(self, 'm', None) is None: - self.Rz = Utils.speye(self.regmesh.cellGradzStencil.shape[0]) + self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0]) else: - f_m = self.regmesh.cellGradzStencil * self.m + f_m = self.regmesh.cellDiffzStencil * self.m self.rz = self.R( f_m , self.qz, self.eps) self.Rz = Utils.sdiag( self.rz ) if getattr(self, '_Wz', None) is None: - self._Wz = Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma)**0.5)*self.Rz*self.regmesh.cellGradzStencil + self._Wz = Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma)**0.5)*self.Rz*self.regmesh.cellDiffzStencil return self._Wz From 2f8b8a36bf71e1ca2e501540b09ee539c4502e01 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 2 Mar 2016 09:46:50 -0800 Subject: [PATCH 11/18] smoothModel --> mrefInSmooth --- SimPEG/Directives.py | 4 ++-- SimPEG/Examples/MT_1D_ForwardAndInversion.py | 2 +- SimPEG/Regularization.py | 20 ++++++++++---------- 3 files changed, 13 insertions(+), 13 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 9ba38e4e..0a98bc09 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -216,7 +216,7 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration): # Save the data. ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) ) phi_ms = 0.5*ms.dot(ms) - if self.reg.smoothModel == True: + if self.reg.mrefInSmooth == True: mref = self.reg.mref else: mref = 0 @@ -249,7 +249,7 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration): # Save the data. ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) ) phi_ms = 0.5*ms.dot(ms) - if self.reg.smoothModel == True: + if self.reg.mrefInSmooth == True: mref = self.reg.mref else: mref = 0 diff --git a/SimPEG/Examples/MT_1D_ForwardAndInversion.py b/SimPEG/Examples/MT_1D_ForwardAndInversion.py index cf17e77b..69cedd98 100644 --- a/SimPEG/Examples/MT_1D_ForwardAndInversion.py +++ b/SimPEG/Examples/MT_1D_ForwardAndInversion.py @@ -100,7 +100,7 @@ def run(plotIt=True): # Regularization - with a regularization mesh regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]],m1d.x0) reg = simpeg.Regularization.Tikhonov(regMesh) - reg.smoothModel = True + reg.mrefInSmooth = True reg.alpha_s = 1e-7 reg.alpha_x = 1. # Inversion problem diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index acd6c68d..990b79c3 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -403,7 +403,7 @@ class BaseRegularization(object): class Tikhonov(BaseRegularization): """ """ - smoothModel = True #: SMOOTH and SMOOTH_MOD_DIF options + mrefInSmooth = True #: SMOOTH and SMOOTH_MOD_DIF options alpha_s = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Ws'], "Smallness weight") alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction") alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction") @@ -489,11 +489,11 @@ class Tikhonov(BaseRegularization): @Utils.timeIt def eval(self, m): - if self.smoothModel == True: + if self.mrefInSmooth == True: r1 = self.Wsmooth * ( self.mapping * (m) ) r2 = self.Ws * ( self.mapping * (m - self.mref) ) return 0.5*(r1.dot(r1)+r2.dot(r2)) - elif self.smoothModel == False: + elif self.mrefInSmooth == False: r = self.W * ( self.mapping * (m - self.mref) ) return 0.5*r.dot(r) @@ -515,7 +515,7 @@ class Tikhonov(BaseRegularization): R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} """ - if self.smoothModel == True: + if self.mrefInSmooth == True: mD1 = self.mapping.deriv(m) mD2 = self.mapping.deriv(m - self.mref) r1 = self.Wsmooth * ( self.mapping * (m)) @@ -523,7 +523,7 @@ class Tikhonov(BaseRegularization): out1 = mD1.T * ( self.Wsmooth.T * r1 ) out2 = mD2.T * ( self.Ws.T * r2 ) out = out1+out2 - elif self.smoothModel == False: + elif self.mrefInSmooth == False: mD = self.mapping.deriv(m - self.mref) r = self.W * ( self.mapping * (m - self.mref) ) out = mD.T * ( self.W.T * r ) @@ -535,7 +535,7 @@ class Simple(BaseRegularization): Only for tensor mesh """ - smoothModel = True #: SMOOTH and SMOOTH_MOD_DIF options + mrefInSmooth = True #: SMOOTH and SMOOTH_MOD_DIF options alpha_s = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Ws'], "Smallness weight") alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction") alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction") @@ -595,11 +595,11 @@ class Simple(BaseRegularization): @Utils.timeIt def eval(self, m): - if self.smoothModel == True: + if self.mrefInSmooth == True: r1 = self.Wsmooth * ( self.mapping * (m) ) r2 = self.Ws * ( self.mapping * (m - self.mref) ) return 0.5*(r1.dot(r1)+r2.dot(r2)) - elif self.smoothModel == False: + elif self.mrefInSmooth == False: r = self.W * ( self.mapping * (m - self.mref) ) return 0.5*r.dot(r) return phim @@ -623,7 +623,7 @@ class Simple(BaseRegularization): R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} """ - if self.smoothModel == True: + if self.mrefInSmooth == True: mD1 = self.mapping.deriv(m) mD2 = self.mapping.deriv(m - self.mref) r1 = self.Wsmooth * ( self.mapping * (m)) @@ -631,7 +631,7 @@ class Simple(BaseRegularization): out1 = mD1.T * ( self.Wsmooth.T * r1 ) out2 = mD2.T * ( self.Ws.T * r2 ) out = out1+out2 - elif self.smoothModel == False: + elif self.mrefInSmooth == False: mD = self.mapping.deriv(m - self.mref) r = self.W * ( self.mapping * (m - self.mref) ) out = mD.T * ( self.W.T * r ) From 33c9059e4e2b9fae8375923197d23d2dafaa62a7 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 8 Mar 2016 16:40:02 -0800 Subject: [PATCH 12/18] SparseRegularization --> Sparse --- SimPEG/Regularization.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 990b79c3..db1bc39b 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -638,7 +638,7 @@ class Simple(BaseRegularization): return out -class SparseRegularization(Simple): +class Sparse(Simple): eps = 1e-1 m = None From ef4513bcd42a0ddaba6b52f6c2fba6af4e70c6a0 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 10 Mar 2016 14:25:53 -0800 Subject: [PATCH 13/18] some cleanup inside of sparse regularization --- SimPEG/Directives.py | 54 ++++++++++++------------- SimPEG/Regularization.py | 66 ++++++++++++++----------------- tests/base/test_regularization.py | 2 +- 3 files changed, 57 insertions(+), 65 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 0a98bc09..46b0b087 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -285,45 +285,43 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration): class update_IRLS(InversionDirective): - m = None - eps_min = None - factor = None - gamma = None - phi_m_last = None + eps_min = None + factor = None + gamma = None + phi_m_last = None - def initialize(self): + def initialize(self): - # Scale the regularization for changes in norm - if getattr(self, 'phi_m_last', None) is not None: - self.reg.gamma = 1. - phim_new = self.reg.eval(self.invProb.curModel) - self.gamma = self.phi_m_last / phim_new + # Scale the regularization for changes in norm + if getattr(self, 'phi_m_last', None) is not None: + self.reg.gamma = 1. + phim_new = self.reg.eval(self.invProb.curModel) + self.gamma = self.phi_m_last / phim_new - self.reg.gamma = self.gamma + self.reg.gamma = self.gamma - def endIter(self): - # Cool the threshold parameter - if getattr(self, 'factor', None) is not None: - eps = self.reg.eps / self.factor + def endIter(self): + # Cool the threshold parameter + 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 + if getattr(self, 'eps_min', None) is not None: + self.reg.eps = np.max([self.eps_min,eps]) + else: + self.reg.eps = eps # Update the model used for the IRLS weights - if getattr(self, 'm', None) is None: - self.reg.m = self.invProb.curModel + self.reg.curModel = self.invProb.curModel - # 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.m.size))**2. - PC = Utils.sdiag(diagA**-1.) + # 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.m.size))**2. + PC = Utils.sdiag(diagA**-1.) - self.opt.approxHinv = PC + self.opt.approxHinv = PC - phim_new = self.reg.eval(self.invProb.curModel) - self.reg.gamma = self.reg.gamma * self.invProb.phi_m_last / phim_new + phim_new = self.reg.eval(self.invProb.curModel) + self.reg.gamma = self.reg.gamma * self.invProb.phi_m_last / phim_new #============================================================================== # import pylab as plt diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index db1bc39b..6ad94caf 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -6,7 +6,7 @@ class RegularizationMesh(object): This contains the operators used in the regularization. Note that these are not necessarily true differential operators, but are constructed from - a SimPEG Mesh. + a SimPEG Mesh. :param Mesh mesh: problem mesh :param numpy.array indActive: bool array, size nC, that is True where we have active cells. Used to reduce the operators so we regularize only on active cells @@ -52,7 +52,7 @@ class RegularizationMesh(object): if getattr(self, '_dim', None) is None: self._dim = self.mesh.dim return self._dim - + @property def _Pac(self): @@ -64,7 +64,7 @@ class RegularizationMesh(object): if getattr(self, '__Pac', None) is None: if self.indActive is None: self.__Pac = Utils.speye(self.mesh.nC) - else: + else: self.__Pac = Utils.speye(self.mesh.nC)[:,self.indActive] return self.__Pac @@ -211,7 +211,7 @@ class RegularizationMesh(object): if getattr(self, '_cellDiffz', None) is None: self._cellDiffz = self._Pafz.T * self.mesh.cellGradz * self._Pac return self._cellDiffz - + @property def faceDiffx(self): """ @@ -233,7 +233,7 @@ class RegularizationMesh(object): if getattr(self, '_faceDiffy', None) is None: self._faceDiffy = self._Pac.T * self.mesh.faceDivy * self._Pafy return self._faceDiffy - + @property def faceDiffz(self): """ @@ -310,7 +310,7 @@ class BaseRegularization(object): if indActive is not None and indActive.dtype != 'bool': tmp = indActive indActive = np.zeros(mesh.nC, dtype=bool) - indActive[tmp] = True + indActive[tmp] = True self.regmesh = RegularizationMesh(mesh,indActive) self.mapping = mapping or self.mapPair(mesh) self.mapping._assertMatchesPair(self.mapPair) @@ -427,7 +427,7 @@ class Tikhonov(BaseRegularization): """Regularization matrix Wx""" if getattr(self, '_Wx', None) is None: Ave_x_vol = self.regmesh.aveCC2Fx * self.regmesh.vol - self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.regmesh.cellDiffx + self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.regmesh.cellDiffx return self._Wx @property @@ -640,13 +640,14 @@ class Simple(BaseRegularization): class Sparse(Simple): - eps = 1e-1 - m = None - gamma = 1. - p = 0. - qx = 2. - qy = 2. - qz = 2. + # set default values + eps = 1e-1 + curModel = None # use a model to compute the weights + gamma = 1. + p = 0. + qx = 2. + qy = 2. + qz = 2. def __init__(self, mesh, mapping=None, indActive=None, **kwargs): Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs) @@ -655,71 +656,64 @@ class Sparse(Simple): @property def Ws(self): """Regularization matrix Ws""" - if getattr(self, 'm', None) is None: + if getattr(self, 'curModel', None) is None: self.Rs = Utils.speye(self.regmesh.nC) else: - f_m = self.m + f_m = self.curModel self.rs = self.R(f_m , self.p, self.eps) #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs)) self.Rs = Utils.sdiag( self.rs ) - self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma)**0.5)*self.Rs + return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma)**0.5)*self.Rs - return self._Ws @property def Wx(self): """Regularization matrix Wx""" - if getattr(self, 'm', None) is None: + if getattr(self, 'curModel', None) is None: self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0]) else: - f_m = self.regmesh.cellDiffxStencil * self.m + f_m = self.regmesh.cellDiffxStencil * self.curModel self.rx = self.R( f_m , self.qx, self.eps) self.Rx = Utils.sdiag( self.rx ) - if getattr(self, '_Wx', None) is None: - self._Wx = Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma)**0.5)*self.Rx*self.regmesh.cellDiffxStencil - return self._Wx + return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma)**0.5)*self.Rx*self.regmesh.cellDiffxStencil @property def Wy(self): """Regularization matrix Wy""" - if getattr(self, 'm', None) is None: + if getattr(self, 'curModel', None) is None: self.Ry = Utils.speye(self.regmesh.cellDiffyStencil.shape[0]) else: - f_m = self.regmesh.cellDiffyStencil * self.m + f_m = self.regmesh.cellDiffyStencil * self.curModel self.ry = self.R( f_m , self.qy, self.eps) self.Ry = Utils.sdiag( self.ry ) - if getattr(self, '_Wy', None) is None: - self._Wy = Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma)**0.5)*self.Ry*self.regmesh.cellDiffyStencil - return self._Wy + return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma)**0.5)*self.Ry*self.regmesh.cellDiffyStencil @property def Wz(self): """Regularization matrix Wz""" - if getattr(self, 'm', None) is None: + if getattr(self, 'curModel', None) is None: self.Rz = Utils.speye(self.regmesh.cellDiffzStencil.shape[0]) else: - f_m = self.regmesh.cellDiffzStencil * self.m + f_m = self.regmesh.cellDiffzStencil * self.curModel self.rz = self.R( f_m , self.qz, self.eps) self.Rz = Utils.sdiag( self.rz ) - if getattr(self, '_Wz', None) is None: - self._Wz = Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma)**0.5)*self.Rz*self.regmesh.cellDiffzStencil - return self._Wz + return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma)**0.5)*self.Rz*self.regmesh.cellDiffzStencil - def R(self, f_m , p, dec): + def R(self, f_m , exponent): - eta = (self.eps**(1-p/2.))**0.5 - r = eta / (f_m**2.+self.eps**2.)**((1-p/2.)/2.) + eta = (self.eps**(1-exponent/2.))**0.5 + r = eta / (f_m**2.+self.eps**2.)**((1-exponent/2.)/2.) return r diff --git a/tests/base/test_regularization.py b/tests/base/test_regularization.py index 52cef349..97223015 100644 --- a/tests/base/test_regularization.py +++ b/tests/base/test_regularization.py @@ -18,7 +18,7 @@ class RegularizationTests(unittest.TestCase): mesh3 = Mesh.TensorMesh([hx, hy, hz]) self.meshlist = [mesh1,mesh2, mesh3] - if testReg: + if testReg: def test_regularization(self): for R in dir(Regularization): r = getattr(Regularization, R) From 838035adaef7991dd48f8c6984b32951b032eea7 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 10 Mar 2016 14:30:43 -0800 Subject: [PATCH 14/18] fixed indentation level on test_regularization --- tests/base/test_regularization.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/tests/base/test_regularization.py b/tests/base/test_regularization.py index 97223015..5aa21ad5 100644 --- a/tests/base/test_regularization.py +++ b/tests/base/test_regularization.py @@ -87,24 +87,24 @@ class RegularizationTests(unittest.TestCase): if testRegMesh: def test_regularizationMesh(self): - for i, mesh in enumerate(self.meshlist): + for i, mesh in enumerate(self.meshlist): - print 'Testing %iD'%mesh.dim + print 'Testing %iD'%mesh.dim - # mapping = r.mapPair(mesh) - # reg = r(mesh, mapping=mapping) - # m = np.random.rand(mapping.nP) + # mapping = r.mapPair(mesh) + # reg = r(mesh, mapping=mapping) + # m = np.random.rand(mapping.nP) - if mesh.dim == 1: - indAct = Utils.mkvc(mesh.gridCC <= 0.8) - elif mesh.dim == 2: - indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5) - elif mesh.dim == 3: - indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5) + if mesh.dim == 1: + indAct = Utils.mkvc(mesh.gridCC <= 0.8) + elif mesh.dim == 2: + indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5) + elif mesh.dim == 3: + indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5) - regmesh = Regularization.RegularizationMesh(mesh, indActive=indAct) + regmesh = Regularization.RegularizationMesh(mesh, indActive=indAct) - assert (regmesh.vol == mesh.vol[indAct]).all() + assert (regmesh.vol == mesh.vol[indAct]).all() if __name__ == '__main__': From 38b4079f0bdb23eaa7b5810ed925e2a4f32102ce Mon Sep 17 00:00:00 2001 From: D Fournier Date: Fri, 11 Mar 2016 11:40:47 -0800 Subject: [PATCH 15/18] Move cell-based weights (i.e. distance weighting) inside regularization. Fix gamma parameter update TO DO: Check inversion print screen -> values don't match reality. --- SimPEG/Directives.py | 8 ++++-- SimPEG/Regularization.py | 60 ++++++++++++++++++++++++++++------------ 2 files changed, 48 insertions(+), 20 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 46b0b087..dbb525f7 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -298,6 +298,7 @@ class update_IRLS(InversionDirective): 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 def endIter(self): @@ -315,13 +316,14 @@ class update_IRLS(InversionDirective): self.reg.curModel = self.invProb.curModel # 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.m.size))**2. + diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() * (self.reg.mapping * np.ones(self.reg.curModel.size))**2. PC = Utils.sdiag(diagA**-1.) self.opt.approxHinv = PC - + + self.reg.gamma = 1. phim_new = self.reg.eval(self.invProb.curModel) - self.reg.gamma = self.reg.gamma * self.invProb.phi_m_last / phim_new + self.reg.gamma = self.invProb.phi_m_last / phim_new #============================================================================== # import pylab as plt diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 6ad94caf..6c6b4b4b 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -540,36 +540,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. + 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 Ws(self): """Regularization matrix Ws""" if getattr(self,'_Ws', None) is None: - self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s)**0.5) + self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5) return self._Ws @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)**0.5)*self.regmesh.cellDiffxStencil + 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)**0.5)*self.regmesh.cellDiffyStencil + 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)**0.5)*self.regmesh.cellDiffzStencil + 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 @@ -648,10 +652,13 @@ class Sparse(Simple): qx = 2. qy = 2. qz = 2. + wght = 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 @property def Ws(self): @@ -661,26 +668,26 @@ class Sparse(Simple): else: f_m = self.curModel - self.rs = self.R(f_m , self.p, self.eps) + self.rs = self.R(f_m , self.p) #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)**0.5)*self.Rs + + return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.wght)**0.5)*self.Rs @property def Wx(self): """Regularization matrix Wx""" - + if getattr(self, 'curModel', None) is None: self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0]) else: f_m = self.regmesh.cellDiffxStencil * self.curModel - self.rx = self.R( f_m , self.qx, self.eps) + self.rx = self.R( f_m , self.qx) self.Rx = Utils.sdiag( self.rx ) - return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma)**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.wght))**0.5)*self.Rx*self.regmesh.cellDiffxStencil @property def Wy(self): @@ -691,10 +698,10 @@ class Sparse(Simple): else: f_m = self.regmesh.cellDiffyStencil * self.curModel - self.ry = self.R( f_m , self.qy, self.eps) + self.ry = self.R( f_m , self.qy) self.Ry = Utils.sdiag( self.ry ) - - return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma)**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.wght))**0.5)*self.Ry*self.regmesh.cellDiffyStencil @property def Wz(self): @@ -705,12 +712,31 @@ class Sparse(Simple): else: f_m = self.regmesh.cellDiffzStencil * self.curModel - self.rz = self.R( f_m , self.qz, self.eps) + self.rz = self.R( f_m , self.qz) self.Rz = Utils.sdiag( self.rz ) - return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma)**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.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.Ws, self.Wsmooth) + #self._W = sp.vstack(wlist) + return sp.vstack(wlist) + def R(self, f_m , exponent): eta = (self.eps**(1-exponent/2.))**0.5 From d226186c8e68c33657a318db0de41ec2f9a0a16e Mon Sep 17 00:00:00 2001 From: D Fournier Date: Fri, 11 Mar 2016 15:09:31 -0800 Subject: [PATCH 16/18] Add auto-beta adjustment. --- SimPEG/Directives.py | 26 +++++++++++++++++++------- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index dbb525f7..de1c60c1 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -289,7 +289,8 @@ class update_IRLS(InversionDirective): factor = None gamma = None phi_m_last = None - + phi_d_last = None + def initialize(self): # Scale the regularization for changes in norm @@ -300,7 +301,10 @@ class update_IRLS(InversionDirective): self.reg.curModel = self.invProb.curModel self.reg.gamma = self.gamma - + + if getattr(self, 'phi_d_last', None) is None: + self.phi_d_last = self.invProb.phi_d + def endIter(self): # Cool the threshold parameter if getattr(self, 'factor', None) is not None: @@ -311,20 +315,28 @@ class update_IRLS(InversionDirective): else: self.reg.eps = eps - + # Get phi_m at the end of current iteration + self.phi_m_last = self.invProb.phi_m_last + # Update the model used for the IRLS weights self.reg.curModel = self.invProb.curModel # Update the pre-conditioner diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() * (self.reg.mapping * np.ones(self.reg.curModel.size))**2. PC = Utils.sdiag(diagA**-1.) - self.opt.approxHinv = PC - + + # Temporarely set gamma to 1. self.reg.gamma = 1. + + # Compute change in model objective function and update scaling phim_new = self.reg.eval(self.invProb.curModel) - self.reg.gamma = self.invProb.phi_m_last / phim_new - + self.reg.gamma = self.phi_m_last / phim_new + print self.reg.gamma + + # TO DO: Re-scale beta if too much change in misfit + self.invProb.beta = self.invProb.beta * self.phi_d_last / self.invProb.phi_d + #============================================================================== # import pylab as plt # plt.figure() From ef467efce03fb02025bc81712c689f3b3f2c8f11 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Tue, 15 Mar 2016 20:56:38 -0700 Subject: [PATCH 17/18] Small change to directive --- SimPEG/Directives.py | 1 - 1 file changed, 1 deletion(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index de1c60c1..e5a63547 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -332,7 +332,6 @@ class update_IRLS(InversionDirective): # Compute change in model objective function and update scaling phim_new = self.reg.eval(self.invProb.curModel) self.reg.gamma = self.phi_m_last / phim_new - print self.reg.gamma # TO DO: Re-scale beta if too much change in misfit self.invProb.beta = self.invProb.beta * self.phi_d_last / self.invProb.phi_d From f92ff1301db92f5dda32125177281064ec59c17f Mon Sep 17 00:00:00 2001 From: D Fournier Date: Thu, 17 Mar 2016 18:45:09 -0700 Subject: [PATCH 18/18] Add reference model in compact term. --- SimPEG/Regularization.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 6c6b4b4b..c20ed973 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -667,7 +667,7 @@ class Sparse(Simple): self.Rs = Utils.speye(self.regmesh.nC) else: - f_m = self.curModel + f_m = self.curModel - self.reg.mref self.rs = self.R(f_m , self.p) #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs)) self.Rs = Utils.sdiag( self.rs )