From 6fcd826673bcbfc866541ea1a1a1141856a2f3b2 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Wed, 20 Jan 2016 14:23:42 -0800 Subject: [PATCH 1/6] Start branch for regularization Add LinearSurvey Add LinearProblem --- SimPEG/Problem.py | 15 +++++++++++++++ SimPEG/Survey.py | 9 ++++++++- 2 files changed, 23 insertions(+), 1 deletion(-) diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index 607cff5b..4d257fd8 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -213,5 +213,20 @@ class BaseTimeProblem(BaseProblem): if hasattr(self, '_timeMesh'): del self._timeMesh +class LinearProblem(BaseProblem): + + surveyPair = Survey.LinearSurvey + def __init__(self, mesh, G, **kwargs): + Problem.BaseProblem.__init__(self, mesh, **kwargs) + self.G = G + + def fields(self, m): + return self.G.dot(m) + + def Jvec(self, m, v, u=None): + return self.G.dot(v) + + def Jtvec(self, m, v, u=None): + return self.G.T.dot(v) diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 88355df1..98a44ee3 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -1,6 +1,5 @@ import Utils, numpy as np, scipy.sparse as sp, uuid - class BaseRx(object): """SimPEG Receiver Object""" @@ -374,3 +373,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[1] From 85b55139e898475bf72aa512e1b327a6766570d0 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Thu, 28 Jan 2016 18:36:34 -0800 Subject: [PATCH 2/6] Implement simple regularization Modified the Optimization.ProjectedGNCG to allow active cells back in. Fix problem regarding the Directive.TargetMisfit --> Survey.Linear had wrong nD value --- SimPEG/Optimization.py | 16 +++- SimPEG/Regularization.py | 160 +++++++++++++++++++++++++++++++++++++++ SimPEG/Survey.py | 2 +- 3 files changed, 176 insertions(+), 2 deletions(-) diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py index 4f2cb062..20a74d0d 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -989,5 +989,19 @@ class ProjectedGNCG(BFGS, Minimize, Remember): if np.logical_or(norm(resid)/normResid0 <= self.tolCG, cgiter == self.maxIterCG): cgFlag = 1 # End CG Iterations + + # Take a gradient step on the active cells if exist + if temp != self.xc.size: + + rhs_a = (Active) * -self.g + + dm_i = max( abs( delx ) ) + dm_a = max( abs(rhs_a) ) + + delx = delx + rhs_a * dm_i / dm_a /10. + + # Only keep gradients going in the right direction on the active set + indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0)) + delx[indx] = 0. - return delx + return delx \ No newline at end of file diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 6c5b5f97..42c5b359 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -243,3 +243,163 @@ 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") + alpha_xx = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction") + alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction") + alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction") + + def __init__(self, mesh, mapping=None, **kwargs): + BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs) + + @property + def Gx(self): + + n = self.mesh.vnC + gx = Utils.ddx(n[0]-1) + gx_square = sp.vstack((gx,gx[-1,:]*-1), format="csr") + + self._Gx = Utils.kron3(Utils.speye(n[2]), Utils.speye(n[1]), gx_square) + + return self._Gx + + @property + def Gy(self): + + n = self.mesh.vnC + gy = Utils.ddx(n[1]-1) + gy_square = sp.vstack((gy,gy[-1,:]*-1), format="csr") + + self._Gy = Utils.kron3(Utils.speye(n[2]), gy_square, Utils.speye(n[0])) + + return self._Gy + + @property + def Gz(self): + + n = self.mesh.vnC + gz = Utils.ddx(n[2]-1) + gz_square = sp.vstack((gz,gz[-1,:]*-1), format="csr") + + self._Gz = Utils.kron3( gz_square , Utils.speye(n[1]), Utils.speye(n[0])) + + return self._Gz + + @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.Gx + 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.Gy + 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.Gz + 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 + 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 + 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 + return self._Wzz + + @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 + + @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 diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 98a44ee3..f8bdec47 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -380,4 +380,4 @@ class LinearSurvey(BaseSurvey): @property def nD(self): - return self.prob.G.shape[1] + return self.prob.G.shape[0] From 54ec7187cb667b285f55f65e4920682e31c6a3ba Mon Sep 17 00:00:00 2001 From: D Fournier Date: Fri, 29 Jan 2016 00:50:59 -0800 Subject: [PATCH 3/6] Add unitCellGrad to DiffOperators Implement SparseRegularization and test on mag problem... works! --- SimPEG/Mesh/DiffOperators.py | 53 +++++++++++++- SimPEG/Regularization.py | 129 +++++++++++++++++++++++++---------- 2 files changed, 146 insertions(+), 36 deletions(-) diff --git a/SimPEG/Mesh/DiffOperators.py b/SimPEG/Mesh/DiffOperators.py index 00bf0259..c5754c93 100644 --- a/SimPEG/Mesh/DiffOperators.py +++ b/SimPEG/Mesh/DiffOperators.py @@ -565,7 +565,58 @@ 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 --------------------- @property diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 42c5b359..80888306 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -260,38 +260,7 @@ class Simple(BaseRegularization): def __init__(self, mesh, mapping=None, **kwargs): BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs) - @property - def Gx(self): - - n = self.mesh.vnC - gx = Utils.ddx(n[0]-1) - gx_square = sp.vstack((gx,gx[-1,:]*-1), format="csr") - - self._Gx = Utils.kron3(Utils.speye(n[2]), Utils.speye(n[1]), gx_square) - - return self._Gx - - @property - def Gy(self): - - n = self.mesh.vnC - gy = Utils.ddx(n[1]-1) - gy_square = sp.vstack((gy,gy[-1,:]*-1), format="csr") - - self._Gy = Utils.kron3(Utils.speye(n[2]), gy_square, Utils.speye(n[0])) - - return self._Gy - - @property - def Gz(self): - - n = self.mesh.vnC - gz = Utils.ddx(n[2]-1) - gz_square = sp.vstack((gz,gz[-1,:]*-1), format="csr") - - self._Gz = Utils.kron3( gz_square , Utils.speye(n[1]), Utils.speye(n[0])) - - return self._Gz + @property def Ws(self): @@ -304,21 +273,21 @@ class Simple(BaseRegularization): 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.Gx + 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.Gy + 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.Gz + self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z)**0.5)*self.mesh.unitCellGradz return self._Wz @property @@ -403,3 +372,93 @@ class Simple(BaseRegularization): r = self.W * ( self.mapping * (m - self.mref) ) out = mD.T * ( self.W.T * r ) return out + +class SparseRegularization(Simple): + + eps0 = None + eps = 1e-4 + coolrate = 2. + m = None + phim_before = None + 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""" + + #iteration = self.opt.iter + #dec = (self.coolrate)**iteration + self.Rs = self.R(self.parent.curModel, Utils.speye(self.mesh.nC), self.p, self.eps) + + if getattr(self,'_Ws', None) is None: + self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)*self.Rs + return self._Ws + + @property + def Wx(self): + """Regularization matrix Wx""" + #iteration = self.opt.iter + #dec = (self.coolrate)**iteration + self.Rx = self.R(self.parent.curModel, self.mesh.unitCellGradx, self.qx, self.eps) + + if getattr(self, '_Wx', None) is None: + self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x)**0.5)*self.Rx*self.mesh.unitCellGradx + return self._Wx + + @property + def Wy(self): + """Regularization matrix Wy""" + #iteration = self.opt.iter + #dec = (self.coolrate)**iteration + self.Ry = self.R(self.parent.curModel, self.mesh.unitCellGrady, self.qy, self.eps) + + if getattr(self, '_Wy', None) is None: + self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y)**0.5)*self.Ry*self.mesh.unitCellGrady + return self._Wy + + @property + def Wz(self): + """Regularization matrix Wz""" + #iteration = self.opt.iter + #dec = (self.rate)**iteration + self.Rz = self.R(self.parent.curModel, self.mesh.unitCellGradz, self.qz, self.eps) + + if getattr(self, '_Wz', None) is None: + self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z)**0.5)*self.Rz*self.mesh.unitCellGradz + return self._Wz + + + def R(self, m, G, p, dec): + + #self.eps = self.eps0*dec + # Scaling to assure equal contribution of all norms + eta = (self.eps**(1-p/2))**0.5 + R = Utils.sdiag( eta / ((self.mapping *((G*m)**2+self.eps**2)**(1-p/2) ) )**0.5 ) + return R \ No newline at end of file From 254fd1c02935110fea1993b97d34b58c97c56bb8 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Sun, 31 Jan 2016 10:50:06 -0800 Subject: [PATCH 4/6] Improved sparse regularization. Issue with spyder debugger since last merge with dev... --- SimPEG/Regularization.py | 44 ++++++++++++++++++++++++++-------------- 1 file changed, 29 insertions(+), 15 deletions(-) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 605aaf2e..9aad9c32 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -451,10 +451,12 @@ class SparseRegularization(Simple): @property def Ws(self): """Regularization matrix Ws""" - - #iteration = self.opt.iter - #dec = (self.coolrate)**iteration - self.Rs = self.R(self.parent.curModel, Utils.speye(self.mesh.nC), self.p, self.eps) + 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) if getattr(self,'_Ws', None) is None: self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)*self.Rs @@ -463,9 +465,13 @@ class SparseRegularization(Simple): @property def Wx(self): """Regularization matrix Wx""" - #iteration = self.opt.iter - #dec = (self.coolrate)**iteration - self.Rx = self.R(self.parent.curModel, self.mesh.unitCellGradx, self.qx, self.eps) + + 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) if getattr(self, '_Wx', None) is None: self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x)**0.5)*self.Rx*self.mesh.unitCellGradx @@ -474,9 +480,13 @@ class SparseRegularization(Simple): @property def Wy(self): """Regularization matrix Wy""" - #iteration = self.opt.iter - #dec = (self.coolrate)**iteration - self.Ry = self.R(self.parent.curModel, self.mesh.unitCellGrady, self.qy, self.eps) + + 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) if getattr(self, '_Wy', None) is None: self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y)**0.5)*self.Ry*self.mesh.unitCellGrady @@ -485,19 +495,23 @@ class SparseRegularization(Simple): @property def Wz(self): """Regularization matrix Wz""" - #iteration = self.opt.iter - #dec = (self.rate)**iteration - self.Rz = self.R(self.parent.curModel, self.mesh.unitCellGradz, self.qz, self.eps) + + 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) if getattr(self, '_Wz', None) is None: self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z)**0.5)*self.Rz*self.mesh.unitCellGradz return self._Wz - def R(self, m, G, p, dec): + def R(self, f_m , p, dec): #self.eps = self.eps0*dec # Scaling to assure equal contribution of all norms eta = (self.eps**(1-p/2))**0.5 - R = Utils.sdiag( eta / ((self.mapping *((G*m)**2+self.eps**2)**(1-p/2) ) )**0.5 ) + R = Utils.sdiag( eta / (((f_m**2+self.eps**2)**(1-p/2) ) )**0.5 ) return R From 04d977f861486c76a797cc7235f8b80c7c4cd3f0 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Sun, 31 Jan 2016 15:31:20 -0800 Subject: [PATCH 5/6] Implement and test Directive for sparse norm... need clean up. --- SimPEG/Directives.py | 55 +++++++++++++++++++++++++++++++++------- SimPEG/Regularization.py | 45 +++++++++++++++++--------------- 2 files changed, 70 insertions(+), 30 deletions(-) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 46576df5..cab6c02f 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -239,14 +239,51 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration): -# class UpdateReferenceModel(Parameter): +class update_IRLS(InversionDirective): -# mref0 = None + 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.) -# def nextIter(self): -# mref = getattr(self, 'm_prev', None) -# if mref is None: -# if self.debug: print 'UpdateReferenceModel is using mref0' -# mref = self.mref0 -# self.m_prev = self.invProb.m_current -# return mref + 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)) +#============================================================================== diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 9aad9c32..cdcfd59a 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -414,11 +414,10 @@ class Simple(BaseRegularization): class SparseRegularization(Simple): - eps0 = None - eps = 1e-4 - coolrate = 2. + eps = 1e-1 + m = None - phim_before = None + gamma = 1. p = 0. qx = 2. qy = 2. @@ -456,10 +455,12 @@ class SparseRegularization(Simple): else: f_m = self.m - self.Rs = self.R(f_m , self.p, self.eps) - - if getattr(self,'_Ws', None) is None: - self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)*self.Rs + 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 @@ -471,10 +472,11 @@ class SparseRegularization(Simple): else: f_m = self.mesh.unitCellGradx * self.m - self.Rx = self.R( f_m , self.qx, self.eps) - + 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)**0.5)*self.Rx*self.mesh.unitCellGradx + self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x*self.gamma)**0.5)*self.Rx*self.mesh.unitCellGradx return self._Wx @property @@ -486,10 +488,11 @@ class SparseRegularization(Simple): else: f_m = self.mesh.unitCellGrady * self.m - self.Ry = self.R( f_m , self.qy, self.eps) - + 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)**0.5)*self.Ry*self.mesh.unitCellGrady + self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y*self.gamma)**0.5)*self.Ry*self.mesh.unitCellGrady return self._Wy @property @@ -501,17 +504,17 @@ class SparseRegularization(Simple): else: f_m = self.mesh.unitCellGradz * self.m - self.Rz = self.R( f_m , self.qz, self.eps) + 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)**0.5)*self.Rz*self.mesh.unitCellGradz + 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): - #self.eps = self.eps0*dec - # Scaling to assure equal contribution of all norms - eta = (self.eps**(1-p/2))**0.5 - R = Utils.sdiag( eta / (((f_m**2+self.eps**2)**(1-p/2) ) )**0.5 ) - return R + eta = (self.eps**(1-p/2.))**0.5 + r = eta / (f_m**2.+self.eps**2.)**((1-p/2.)/2.) + + return r From b16b1b7526bc4bb9f7f2f8cc497c17447cf664c4 Mon Sep 17 00:00:00 2001 From: D Fournier Date: Mon, 1 Feb 2016 21:02:11 -0800 Subject: [PATCH 6/6] Add ModelBuilder sphere model. Add example for DC pseudo section -> Requires a pull request in SimpegDC for dependancies. --- .../Examples/DC_PseudoSection_Simulation.py | 179 ++++++++++++++++++ SimPEG/Utils/ModelBuilder.py | 38 ++++ 2 files changed, 217 insertions(+) create mode 100644 SimPEG/Examples/DC_PseudoSection_Simulation.py diff --git a/SimPEG/Examples/DC_PseudoSection_Simulation.py b/SimPEG/Examples/DC_PseudoSection_Simulation.py new file mode 100644 index 00000000..c3517360 --- /dev/null +++ b/SimPEG/Examples/DC_PseudoSection_Simulation.py @@ -0,0 +1,179 @@ +from SimPEG import * +import simpegDCIP as DC +import scipy.interpolate as interpolation +import matplotlib.pyplot as plt +import time +import re + +def run(loc=np.c_[[-50.,0.,-50.],[50.,0.,-50.]], sig=np.r_[1e-2,1e-1,1e-3], radi=np.r_[25.,25.], param = np.r_[30.,30.,5], stype = 'dpdp', plotIt=True): + """ + DC Forward Simulation + + Forward model conductive spheres in a half-space and plot a pseudo-section + + Created on Mon Feb 01 19:28:06 2016 + + @fourndo + """ + + # First we need to create a mesh and a model. + + # This is our mesh + dx = 5. + + hxind = [(dx,15,-1.3), (dx, 75), (dx,15,1.3)] + hyind = [(dx,15,-1.3), (dx, 10), (dx,15,1.3)] + hzind = [(dx,15,-1.3),(dx, 15)] + + mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCN') + + + # Set background conductivity + model = np.ones(mesh.nC) * sig[0] + + # First anomaly + ind = Utils.ModelBuilder.getIndicesSphere(loc[:,0],radi[0],mesh.gridCC) + model[ind] = sig[1] + + # Second anomaly + ind = Utils.ModelBuilder.getIndicesSphere(loc[:,1],radi[1],mesh.gridCC) + model[ind] = sig[2] + + # Get index of the center + indy = int(mesh.nCy/2) + + + # Plot the model for reference + # Define core mesh extent + xlim = 200 + zlim = 125 + + # Specify the survey type: "pdp" | "dpdp" + + + # Then specify the end points of the survey. Let's keep it simple for now and survey above the anomalies, top of the mesh + ends = [(-175,0),(175,0)] + ends = np.c_[np.asarray(ends),np.ones(2).T*mesh.vectorNz[-1]] + + # Snap the endpoints to the grid. Easier to create 2D section. + indx = Utils.closestPoints(mesh, ends ) + locs = np.c_[mesh.gridCC[indx,0],mesh.gridCC[indx,1],np.ones(2).T*mesh.vectorNz[-1]] + + # We will handle the geometry of the survey for you and create all the combination of tx-rx along line + [Tx, Rx] = DC.gen_DCIPsurvey(locs, mesh, stype, param[0], param[1], param[2]) + + # Define some global geometry + dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) ) + dl_x = ( Tx[-1][0,1] - Tx[0][0,0] ) / dl_len + dl_y = ( Tx[-1][1,1] - Tx[0][1,0] ) / dl_len + azm = np.arctan(dl_y/dl_x) + + #Set boundary conditions + mesh.setCellGradBC('neumann') + + # Define the differential operators needed for the DC problem + Div = mesh.faceDiv + Grad = mesh.cellGrad + Msig = Utils.sdiag(1./(mesh.aveF2CC.T*(1./model))) + + A = Div*Msig*Grad + + # Change one corner to deal with nullspace + A[0,0] = 1 + A = sp.csc_matrix(A) + + # We will solve the system iteratively, so a pre-conditioner is helpful + # This is simply a Jacobi preconditioner (inverse of the main diagonal) + dA = A.diagonal() + P = sp.spdiags(1/dA,0,A.shape[0],A.shape[0]) + + # Now we can solve the system for all the transmitters + # We want to store the data + data = [] + + # There is probably a more elegant way to do this, but we can just for-loop through the transmitters + for ii in range(len(Tx)): + + start_time = time.time() # Let's time the calculations + + #print("Transmitter %i / %i\r" % (ii+1,len(Tx))) + + # Select dipole locations for receiver + rxloc_M = np.asarray(Rx[ii][:,0:3]) + rxloc_N = np.asarray(Rx[ii][:,3:]) + + + # For usual cases "dpdp" or "gradient" + if not re.match(stype,'pdp'): + inds = Utils.closestPoints(mesh, np.asarray(Tx[ii]).T ) + RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1,1] / mesh.vol[inds] ) + + else: + + # Create an "inifinity" pole + tx = np.squeeze(Tx[ii][:,0:1]) + tinf = tx + np.array([dl_x,dl_y,0])*dl_len*2 + inds = Utils.closestPoints(mesh, np.c_[tx,tinf].T) + RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1] / mesh.vol[inds] ) + + + # Iterative Solve + Ainvb = sp.linalg.bicgstab(P*A,P*RHS, tol=1e-5) + + # We now have the potential everywhere + phi = mkvc(Ainvb[0]) + + # Solve for phi on pole locations + P1 = mesh.getInterpolationMat(rxloc_M, 'CC') + P2 = mesh.getInterpolationMat(rxloc_N, 'CC') + + # Compute the potential difference + dtemp = (P1*phi - P2*phi)*np.pi + + data.append( dtemp ) + print '\rTransmitter {0} of {1} -> Time:{2} sec'.format(ii,len(Tx),time.time()- start_time), + + print 'Transmitter {0} of {1}'.format(ii,len(Tx)) + print 'Forward completed' + + + # Let's just convert the 3D format into 2D (distance along line) and plot + [Tx2d, Rx2d] = DC.convertObs_DC3D_to_2D(Tx,Rx) + + + # Here is an example for the first tx-rx array + if plotIt: + fig = plt.figure() + ax = plt.subplot(2,1,1, aspect='equal') + mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y', ind = indy,grid=True) + ax.set_title('E-W section at '+str(mesh.vectorCCy[indy])+' m') + plt.gca().set_aspect('equal', adjustable='box') + + plt.scatter(Tx[0][0,:],Tx[0][2,:],s=40,c='g', marker='v') + plt.scatter(Rx[0][:,0::3],Rx[0][:,2::3],s=40,c='y') + plt.xlim([-xlim,xlim]) + plt.ylim([-zlim,mesh.vectorNz[-1]+dx]) + + + ax = plt.subplot(2,1,2, aspect='equal') + + # Plot the location of the spheres for reference + circle1=plt.Circle((loc[0,0]-Tx[0][0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3) + circle2=plt.Circle((loc[0,1]-Tx[0][0,0],loc[2,1]),radi[1],color='k',fill=False, lw=3) + ax.add_artist(circle1) + ax.add_artist(circle2) + + # Add the speudo section + DC.plot_pseudoSection(Tx2d,Rx2d,data,mesh.vectorNz[-1],stype) + + plt.scatter(Tx2d[0][:],Tx[0][2,:],s=40,c='g', marker='v') + plt.scatter(Rx2d[0][:],Rx[0][:,2::3],s=40,c='y') + plt.plot(np.r_[Tx2d[0][0],Rx2d[-1][-1,-1]],np.ones(2)*mesh.vectorNz[-1], color='k') + plt.ylim([-zlim,mesh.vectorNz[-1]+dx]) + + plt.show() + + return fig, ax + +if __name__ == '__main__': + run() \ No newline at end of file diff --git a/SimPEG/Utils/ModelBuilder.py b/SimPEG/Utils/ModelBuilder.py index 52d13820..435856f7 100644 --- a/SimPEG/Utils/ModelBuilder.py +++ b/SimPEG/Utils/ModelBuilder.py @@ -118,6 +118,44 @@ def defineElipse(ccMesh, center=[0,0,0], anisotropy=[1,1,1], slope=10., theta=0. D = np.sqrt(np.sum(G**2,axis=1)) return -np.arctan((D-1)*slope)*(2./np.pi)/2.+0.5 +def getIndicesSphere(center,radius,ccMesh): + """ + Creates a vector containing the sphere indices in the cell centers mesh. + Returns a tuple + + The sphere is defined by the points + + p0, describe the position of the center of the cell + + r, describe the radius of the sphere. + + ccMesh represents the cell-centered mesh + + The points p0 must live in the the same dimensional space as the mesh. + + """ + + # Validation: mesh and point (p0) live in the same dimensional space + dimMesh = np.size(ccMesh[0,:]) + assert len(center) == dimMesh, "Dimension mismatch. len(p0) != dimMesh" + + if dimMesh == 1: + # Define the reference points + + ind = np.abs(center[0] - ccMesh[:,0]) < radius + + elif dimMesh == 2: + # Define the reference points + + ind = np.sqrt( ( center[0] - ccMesh[:,0] )**2 + ( center[1] - ccMesh[:,1] )**2 ) < radius + + elif dimMesh == 3: + # Define the points + ind = np.sqrt( ( center[0] - ccMesh[:,0] )**2 + ( center[1] - ccMesh[:,1] )**2 + ( center[2] - ccMesh[:,2] )**2 ) < radius + + # Return a tuple + return ind + def defineTwoLayers(ccMesh,depth,vals=[0,1]): """ Define a two layered model. Depth of the first layer must be specified.