Skip to content

Commit

Permalink
Attempt to revert...again
Browse files Browse the repository at this point in the history
  • Loading branch information
@fourndo
fourndo committed Jan 31, 2016
1 parent 54ec718 commit a051d6e
Showing 1 changed file with 54 additions and 15 deletions.
69 changes: 54 additions & 15 deletions SimPEG/Regularization.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,12 +20,13 @@ class BaseRegularization(object):
mesh = None #: A SimPEG.Mesh instance.
mref = None #: Reference model.

def __init__(self, mesh, mapping=None, **kwargs):
def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
Utils.setKwargs(self, **kwargs)
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):
Expand Down Expand Up @@ -112,8 +113,6 @@ def eval2Deriv(self, m, v=None):
return mD.T * ( self.W.T * ( self.W * ( mD * v) ) )




class Tikhonov(BaseRegularization):
"""
"""
Expand All @@ -126,14 +125,18 @@ class Tikhonov(BaseRegularization):
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):
def __init__(self, mesh, mapping=None, indActive = None, **kwargs):
BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs)
self.indActive = indActive

@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.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
return self._Ws

@property
Expand All @@ -142,6 +145,13 @@ def Wx(self):
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

return self._Wx

@property
Expand All @@ -150,6 +160,13 @@ def Wy(self):
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

return self._Wy

@property
Expand All @@ -158,27 +175,49 @@ def Wz(self):
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

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

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

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

return self._Wzz

@property
Expand Down Expand Up @@ -245,9 +284,9 @@ def evalDeriv(self, m):

class Simple(BaseRegularization):
"""
Only for tensor mesh
"""
Only for a 3D 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")
Expand All @@ -258,10 +297,10 @@ class Simple(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, **kwargs):
assert isinstance(mesh, Mesh.TensorMesh)
assert mesh.dim == 3
BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs)



@property
def Ws(self):
"""Regularization matrix Ws"""
Expand Down Expand Up @@ -415,7 +454,7 @@ def Ws(self):

#iteration = self.opt.iter
#dec = (self.coolrate)**iteration
self.Rs = self.R(self.parent.curModel, Utils.speye(self.mesh.nC), self.p, self.eps)
self.Rs = self.R(self.m, 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
Expand All @@ -426,7 +465,7 @@ 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)
self.Rx = self.R(self.m, 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
Expand All @@ -437,7 +476,7 @@ 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)
self.Ry = self.R(self.m, 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
Expand All @@ -448,7 +487,7 @@ 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)
self.Rz = self.R(self.m, 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
Expand All @@ -461,4 +500,4 @@ def R(self, m, G, p, 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
return R

0 comments on commit a051d6e

Please sign in to comment.