mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-28 11:53:55 +08:00
Rowan Cockett
2782a6fcfb
e.g. output: ‘’’ Counters: InexactGaussNewton.doEndIteration : 6 InexactGaussNewton.printIter : 7 InexactGaussNewton.scaleSearchDirection : 6 Times: mean sum InexactGaussNewton.findSearchDirection : 1.55e-02, 9.29e-02, 6x InexactGaussNewton.minimize : 1.10e-01, 1.10e-01, 1x InexactGaussNewton.modifySearchDirection: 2.89e-04, 1.73e-03, 6x InexactGaussNewton.projection : 3.69e-06, 1.11e-04, 30x InexactGaussNewton.stoppingCriteria : 1.16e-04, 1.51e-03, 13x Inversion.dataObj : 6.60e-05, 8.58e-04, 13x Inversion.dataObj2Deriv : 1.03e-04, 6.20e-03, 60x Inversion.dataObjDeriv : 5.06e-05, 3.54e-04, 7x Inversion.evalFunction : 7.75e-04, 1.01e-02, 13x Inversion.run : 1.10e-01, 1.10e-01, 1x Regularization.modelObj : 3.56e-04, 4.63e-03, 13x Regularization.modelObj2Deriv : 1.29e-03, 7.76e-02, 60x Regularization.modelObjDeriv : 5.01e-04, 3.51e-03, 7x ‘’’
122 lines
3.3 KiB
Python
122 lines
3.3 KiB
Python
from SimPEG.utils import sdiag, count, timeIt
|
|
import numpy as np
|
|
|
|
class Regularization(object):
|
|
"""docstring for Regularization"""
|
|
|
|
@property
|
|
def mref(self):
|
|
if getattr(self, '_mref', None) is None:
|
|
self._mref = np.zeros(self.mesh.nC);
|
|
return self._mref
|
|
@mref.setter
|
|
def mref(self, value):
|
|
self._mref = value
|
|
|
|
@property
|
|
def Ws(self):
|
|
if getattr(self,'_Ws', None) is None:
|
|
self._Ws = sdiag(self.mesh.vol)
|
|
return self._Ws
|
|
|
|
@property
|
|
def Wx(self):
|
|
if getattr(self, '_Wx', None) is None:
|
|
a = self.mesh.r(self.mesh.area,'F','Fx','V')
|
|
self._Wx = sdiag(a)*self.mesh.cellGradx
|
|
return self._Wx
|
|
|
|
@property
|
|
def Wy(self):
|
|
if getattr(self, '_Wy', None) is None:
|
|
a = self.mesh.r(self.mesh.area,'F','Fy','V')
|
|
self._Wy = sdiag(a)*self.mesh.cellGrady
|
|
return self._Wy
|
|
|
|
@property
|
|
def Wz(self):
|
|
if getattr(self, '_Wz', None) is None:
|
|
a = self.mesh.r(self.mesh.area,'F','Fz','V')
|
|
self._Wz = sdiag(a)*self.mesh.cellGradz
|
|
return self._Wz
|
|
|
|
alpha_s = 1e-6
|
|
alpha_x = 1
|
|
alpha_y = 1
|
|
alpha_z = 1
|
|
|
|
counter = None
|
|
|
|
def __init__(self, mesh):
|
|
self.mesh = mesh
|
|
|
|
def pnorm(self, r):
|
|
return 0.5*r.dot(r)
|
|
|
|
@timeIt
|
|
def modelObj(self, m):
|
|
mresid = m - self.mref
|
|
|
|
mobj = self.alpha_s * self.pnorm( self.Ws * mresid )
|
|
|
|
mobj += self.alpha_x * self.pnorm( self.Wx * mresid )
|
|
|
|
if self.mesh.dim > 1:
|
|
mobj += self.alpha_y * self.pnorm( self.Wy * mresid )
|
|
if self.mesh.dim > 2:
|
|
mobj += self.alpha_z * self.pnorm( self.Wz * mresid )
|
|
|
|
return mobj
|
|
|
|
@timeIt
|
|
def modelObjDeriv(self, m):
|
|
"""
|
|
|
|
In 1D:
|
|
|
|
.. math::
|
|
|
|
m_{\\text{obj}} = {1 \over 2}\\alpha_s \left\| W_s (m- m_{\\text{ref}})\\right\|^2_2
|
|
+ {1 \over 2}\\alpha_x \left\| W_x (m- m_{\\text{ref}})\\right\|^2_2
|
|
|
|
\\frac{ \partial m_{\\text{obj}} }{\partial m} =
|
|
\\alpha_s W_s^{\\top} W_s (m - m_{\\text{ref}}) +
|
|
\\alpha_x W_x^{\\top} W_x (m - m_{\\text{ref}})
|
|
|
|
|
|
\\frac{ \partial^2 m_{\\text{obj}} }{\partial m^2} =
|
|
\\alpha_s W_s^{\\top} W_s +
|
|
\\alpha_x W_x^{\\top} W_x
|
|
|
|
"""
|
|
|
|
mresid = m - self.mref
|
|
|
|
mobjDeriv = self.alpha_s * self.Ws.T * ( self.Ws * mresid)
|
|
|
|
mobjDeriv = mobjDeriv + self.alpha_x * self.Wx.T * ( self.Wx * mresid)
|
|
|
|
if self.mesh.dim > 1:
|
|
mobjDeriv = mobjDeriv + self.alpha_y * self.Wy.T * ( self.Wy * mresid)
|
|
if self.mesh.dim > 2:
|
|
mobjDeriv = mobjDeriv + self.alpha_z * self.Wz.T * ( self.Wz * mresid)
|
|
|
|
return mobjDeriv
|
|
|
|
|
|
@timeIt
|
|
def modelObj2Deriv(self, m):
|
|
mresid = m - self.mref
|
|
|
|
mobj2Deriv = self.alpha_s * self.Ws.T * self.Ws
|
|
|
|
mobj2Deriv = mobj2Deriv + self.alpha_x * self.Wx.T * self.Wx
|
|
|
|
if self.mesh.dim > 1:
|
|
mobj2Deriv = mobj2Deriv + self.alpha_y * self.Wy.T * self.Wy
|
|
if self.mesh.dim > 2:
|
|
mobj2Deriv = mobj2Deriv + self.alpha_z * self.Wz.T * self.Wz
|
|
|
|
return mobj2Deriv
|
|
|