mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-03 19:32:01 +08:00
- Meshlesses Identity Map (takes nP instead of a mesh)
- Tikhonov regularization if active cells are used (don't take derivs across interfaces between active cells and not) - testing improvements: test 1D, 2D, 3D on a random tensor mesh , also test that for a constant mref, phi_m(ref) = 0
This commit is contained in:
Lindsey Heagy
@@ -118,6 +118,78 @@ class IdentityMap(object):
|
||||
def __str__(self):
|
||||
return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1])
|
||||
|
||||
|
||||
class IdentityMap_Meshless(IdentityMap):
|
||||
|
||||
def __init__(self, nP=None, **kwargs):
|
||||
IdentityMap.__init__(self, None, **kwargs)
|
||||
self._nP = nP
|
||||
|
||||
@property
|
||||
def nP(self):
|
||||
"""
|
||||
:rtype: int
|
||||
:return: number of parameters in the model
|
||||
"""
|
||||
if self._nP is None:
|
||||
return '*'
|
||||
return self._nP
|
||||
|
||||
@property
|
||||
def shape(self):
|
||||
"""
|
||||
The default shape is (mesh.nC, nP).
|
||||
|
||||
:rtype: (int,int)
|
||||
:return: shape of the operator as a tuple
|
||||
"""
|
||||
if self._nP is None:
|
||||
return ('*', '*')
|
||||
return (self.nP, self.nP)
|
||||
|
||||
|
||||
def _transform(self, m):
|
||||
"""
|
||||
Changes the model into the physical property.
|
||||
|
||||
.. note::
|
||||
|
||||
This can be called by the __mul__ property against a numpy.ndarray.
|
||||
|
||||
:param numpy.array m: model
|
||||
:rtype: numpy.array
|
||||
:return: transformed model
|
||||
|
||||
"""
|
||||
return m
|
||||
|
||||
def inverse(self, D):
|
||||
"""
|
||||
Changes the physical property into the model.
|
||||
|
||||
.. note::
|
||||
|
||||
The *transformInverse* may not be easy to create in general.
|
||||
|
||||
:param numpy.array D: physical property
|
||||
:rtype: numpy.array
|
||||
:return: model
|
||||
|
||||
"""
|
||||
raise NotImplementedError('The transformInverse is not implemented.')
|
||||
|
||||
def deriv(self, m):
|
||||
"""
|
||||
The derivative of the transformation.
|
||||
|
||||
:param numpy.array m: model
|
||||
:rtype: scipy.csr_matrix
|
||||
:return: derivative of transformed model
|
||||
|
||||
"""
|
||||
return sp.identity(self.nP)
|
||||
|
||||
|
||||
class ComboMap(IdentityMap):
|
||||
"""Combination of various maps."""
|
||||
|
||||
|
||||
@@ -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 Maps.IdentityMap(mesh)
|
||||
self.mapping._assertMatchesPair(self.mapPair)
|
||||
self.indActive = indActive
|
||||
|
||||
@property
|
||||
def parent(self):
|
||||
@@ -112,8 +113,6 @@ class BaseRegularization(object):
|
||||
return mD.T * ( self.W.T * ( self.W * ( mD * v) ) )
|
||||
|
||||
|
||||
|
||||
|
||||
class Tikhonov(BaseRegularization):
|
||||
"""**Tikhonov Regularization**
|
||||
|
||||
@@ -205,14 +204,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
|
||||
@@ -221,6 +224,13 @@ class Tikhonov(BaseRegularization):
|
||||
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
|
||||
@@ -229,6 +239,13 @@ class Tikhonov(BaseRegularization):
|
||||
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
|
||||
@@ -237,6 +254,13 @@ class Tikhonov(BaseRegularization):
|
||||
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
|
||||
@@ -244,6 +268,11 @@ class Tikhonov(BaseRegularization):
|
||||
"""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
|
||||
@@ -251,6 +280,11 @@ class Tikhonov(BaseRegularization):
|
||||
"""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
|
||||
@@ -258,6 +292,11 @@ class Tikhonov(BaseRegularization):
|
||||
"""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
|
||||
|
||||
@@ -4,11 +4,17 @@ from SimPEG import *
|
||||
from scipy.sparse.linalg import dsolve
|
||||
import inspect
|
||||
|
||||
TOL = 1e-20
|
||||
|
||||
class RegularizationTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
self.mesh2 = Mesh.TensorMesh([3, 2])
|
||||
hx, hy, hz = np.random.rand(10), np.random.rand(9), np.random.rand(8)
|
||||
hx, hy, hz = hx/hx.sum(), hy/hy.sum(), hz/hz.sum()
|
||||
mesh1 = Mesh.TensorMesh([hx])
|
||||
mesh2 = Mesh.TensorMesh([hx, hy])
|
||||
mesh3 = Mesh.TensorMesh([hx, hy, hz])
|
||||
self.meshlist = [mesh1,mesh2, mesh3]
|
||||
|
||||
def test_regularization(self):
|
||||
for R in dir(Regularization):
|
||||
@@ -16,18 +22,63 @@ class RegularizationTests(unittest.TestCase):
|
||||
if not inspect.isclass(r): continue
|
||||
if not issubclass(r, Regularization.BaseRegularization):
|
||||
continue
|
||||
# if 'Regularization' not in R: continue
|
||||
mapping = r.mapPair(self.mesh2)
|
||||
reg = r(self.mesh2, mapping=mapping)
|
||||
m = np.random.rand(mapping.nP)
|
||||
reg.mref = m[:]*np.mean(m)
|
||||
|
||||
for i, mesh in enumerate(self.meshlist):
|
||||
|
||||
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 '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)
|
||||
|
||||
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 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
|
||||
|
||||
for i, mesh in enumerate(self.meshlist):
|
||||
|
||||
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)
|
||||
|
||||
mapping = Maps.IdentityMap_Meshless(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)
|
||||
|
||||
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 2 Deriv:', R
|
||||
passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False)
|
||||
self.assertTrue(passed)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
Reference in New Issue
Block a user