diff --git a/SimPEG/Examples/Inversion_IRLS.py b/SimPEG/Examples/Inversion_IRLS.py index 06ef4be4..c925cce2 100644 --- a/SimPEG/Examples/Inversion_IRLS.py +++ b/SimPEG/Examples/Inversion_IRLS.py @@ -84,12 +84,12 @@ def run(N=200, plotIt=True): #============================================================================== #reg.recModel = mrec - reg.wght = np.ones(mesh.nC) + # reg.cell_weight = np.ones(mesh.nC) reg.mref = np.zeros(mesh.nC) reg.eps_p = 5e-2 reg.eps_q = 1e-2 reg.norms = [0., 0., 2., 2.] - reg.wght = wr + reg.cell_weight = wr opt = Optimization.ProjectedGNCG(maxIter=10 ,lower=-2.,upper=2., maxIterLS = 20, maxIterCG= 20, tolCG = 1e-3) invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.) diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 16bb4b9f..5947ed03 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -410,11 +410,11 @@ class Simple(BaseRegularization): Simple regularization that does not include length scales in the derivatives. """ - mrefInSmooth = False #: SMOOTH and SMOOTH_MOD_DIF options + mrefInSmooth = False #: include mref in the smoothness? alpha_s = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "Smallness weight") - alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction") - alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction") - alpha_z = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction") + alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction") + alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction") + alpha_z = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction") cell_weights = 1. def __init__(self, mesh, mapping=None, indActive=None, **kwargs): @@ -454,6 +454,8 @@ class Simple(BaseRegularization): @property def Wsmooth(self): """Full smoothness regularization matrix W""" + print 'wtf why are we using Wsmooth' + raise NotImplementedError if getattr(self, '_Wsmooth', None) is None: wlist = (self.Wx,) if self.regmesh.dim > 1: @@ -482,6 +484,13 @@ class Simple(BaseRegularization): r = self.Wsmall * ( self.mapping * (m - self.mref) ) return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) ) + @Utils.timeIt + def _evalSmall2Deriv(self, m, v = None): + rDeriv = self.Wsmall * ( self.mapping.deriv(m - self.mref) ) + if v is not None: + return rDeriv.T * (rDeriv * v) + return rDeriv.T * rDeriv + @Utils.timeIt def _evalSmoothx(self, m): if self.mrefInSmooth == True: @@ -524,6 +533,17 @@ class Simple(BaseRegularization): r = self.Wx * ( self.mapping * m ) return r.T * ( self.Wx * self.mapping.deriv(m) ) + @Utils.timeIt + def _evalSmoothx2Deriv(self, m, v=None): + if self.mrefInSmooth == True: + rDeriv = self.Wx * ( self.mapping.deriv( m - self.mref ) ) + elif self.mrefInSmooth == False: + rDeriv = self.Wx * ( self.mapping.deriv(m) ) + + if v is not None: + return rDeriv.T * ( rDeriv * v ) + return rDeriv.T * rDeriv + @Utils.timeIt def _evalSmoothyDeriv(self, m): if self.mrefInSmooth == True: @@ -533,6 +553,17 @@ class Simple(BaseRegularization): r = self.Wy * ( self.mapping * m ) return r.T * ( self.Wy * self.mapping.deriv(m) ) + @Utils.timeIt + def _evalSmoothy2Deriv(self, m, v=None): + if self.mrefInSmooth == True: + rDeriv = self.Wy * ( self.mapping.deriv( m - self.mref ) ) + elif self.mrefInSmooth == False: + rDeriv = self.Wy * ( self.mapping.deriv(m) ) + + if v is not None: + return rDeriv.T * ( rDeriv * v ) + return rDeriv.T * rDeriv + @Utils.timeIt def _evalSmoothzDeriv(self, m): if self.mrefInSmooth == True: @@ -542,6 +573,17 @@ class Simple(BaseRegularization): r = self.Wz * ( self.mapping * m ) return r.T * ( self.Wz * self.mapping.deriv(m) ) + @Utils.timeIt + def _evalSmoothz2Deriv(self, m, v=None): + if self.mrefInSmooth == True: + rDeriv = self.Wz * ( self.mapping.deriv( m - self.mref ) ) + elif self.mrefInSmooth == False: + rDeriv = self.Wz * ( self.mapping.deriv(m) ) + + if v is not None: + return rDeriv.T * ( rDeriv * v ) + return rDeriv.T * rDeriv + @Utils.timeIt def _evalSmoothDeriv(self, m): deriv = self._evalSmoothxDeriv(m) @@ -551,6 +593,15 @@ class Simple(BaseRegularization): deriv += self._evalSmoothzDeriv(m) return deriv + @Utils.timeIt + def _evalSmooth2Deriv(self, m, v=None): + deriv = self._evalSmoothx2Deriv(m, v) + if self.regmesh.dim > 1: + deriv += self._evalSmoothy2Deriv(m, v) + if self.regmesh.dim > 2: + deriv += self._evalSmoothz2Deriv(m, v) + return deriv + @Utils.timeIt def eval(self, m): @@ -574,6 +625,10 @@ class Simple(BaseRegularization): """ return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + @Utils.timeIt + def eval2Deriv(self, m, v=None): + return self._evalSmall2Deriv(m, v) + self._evalSmooth2Deriv(m, v) + class Tikhonov(Simple): @@ -740,7 +795,37 @@ class Tikhonov(Simple): return r.T * ( self.Wzz * self.mapping.deriv(m) ) @Utils.timeIt - def _evalSmooth2Deriv(self, m): + def _evalSmoothxx2Deriv(self, m, v=None): + if self.mrefInSmooth == True: + rDeriv = self.Wxx * ( self.mapping.deriv( m - self.mref ) ) + elif self.mrefInSmooth == False: + rDeriv = self.Wxx * self.mapping.deriv(m) + if v is not None: + return rDeriv.T * (rDeriv * v) + return rDeriv.T * rDeriv + + @Utils.timeIt + def _evalSmoothyy2Deriv(self, m, v=None): + if self.mrefInSmooth == True: + rDeriv = self.Wyy * ( self.mapping.deriv( m - self.mref ) ) + elif self.mrefInSmooth == False: + rDeriv = self.Wyy * self.mapping.deriv(m) + if v is not None: + return rDeriv.T * (rDeriv * v) + return rDeriv.T * rDeriv + + @Utils.timeIt + def _evalSmoothzz2Deriv(self, m, v=None): + if self.mrefInSmooth == True: + rDeriv = self.Wzz * ( self.mapping.deriv( m - self.mref ) ) + elif self.mrefInSmooth == False: + rDeriv = self.Wzz * self.mapping.deriv(m) + if v is not None: + return rDeriv.T * (rDeriv * v) + return rDeriv.T * rDeriv + + @Utils.timeIt + def _evalSmoothDeriv2(self, m): deriv = self._evalSmoothxxDeriv(m) if self.regmesh.dim > 1: deriv += self._evalSmoothyyDeriv(m) @@ -748,6 +833,15 @@ class Tikhonov(Simple): deriv += self._evalSmoothzzDeriv(m) return deriv + @Utils.timeIt + def _evalSmooth2Deriv2(self, m, v=None): + deriv = self._evalSmoothxx2Deriv(m, v) + if self.regmesh.dim > 1: + deriv += self._evalSmoothyy2Deriv(m, v) + if self.regmesh.dim > 2: + deriv += self._evalSmoothzz2Deriv(m, v) + return deriv + @Utils.timeIt def eval(self, m): @@ -769,7 +863,24 @@ class Tikhonov(Simple): R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} """ - return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + self._evalSmooth2Deriv(m) + return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m) + self._evalSmoothDeriv2(m) + + def eval2Deriv(self, m): + """ + The regularization is: + + .. math:: + + R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})} + + So the derivative is straight forward: + + .. math:: + + R(m) = \mathbf{W^\\top W (m-m_\\text{ref})} + + """ + return self._evalSmall2Deriv(m) + self._evalSmooth2Deriv(m) + self._evalSmooth2Deriv2(m)