wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-10 20:19:36 +08:00
Code Issues Packages Projects Releases Wiki Activity

Create sensitivity re-weighting directive

Browse Source 
Adapt Map/polymap for actInd (topography)
This commit is contained in:
D Fournier
2016-04-06 07:22:18 -07:00
parent df620b42bd
commit ba977206ea
2 changed files with 54 additions and 22 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Directives.py
+25 -6
View File
@@ -356,10 +356,29 @@ class update_lin_PreCond(InversionDirective):
PC = Utils.sdiag(diagA**-1.)
self.opt.approxHinv = PC
print 'Updated pre-cond'
class update_Wj(InversionDirective):
"""
Create approx-sensitivity base weighting
"""
k = None
#==============================================================================
# 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))
#==============================================================================
def endIter(self):
if self.opt.iter == 2:
m = self.invProb.curModel
if self.k is None:
self.k = int(self.survey.nD/10)
def JtJv(v):
Jv = self.prob.Jvec(m, v)
return self.prob.Jtvec(m,Jv)
JtJdiag = Utils.diagEst(JtJv,6,k=100)
JtJdiag = JtJdiag
JtJdiag = JtJdiag / max(JtJdiag)
self.reg.wght = JtJdiag
SimPEG/Maps.py
+29 -16
View File
@@ -759,15 +759,27 @@ class PolyMap(IdentityMap):
m = [\sigma_1, \sigma_2, c]
"""
def __init__(self, mesh, order, logSigma=True, normal='X'):
"""
def __init__(self, mesh, order, logSigma=True, normal='X', actInd = None):
IdentityMap.__init__(self, mesh)
self.logSigma = logSigma
self.order = order
self.normal = normal
slope = 1e4
self.normal = normal
self.actInd = actInd
if getattr(self, 'actInd', None) is None:
self.actInd = range(self.mesh.nC)
self.nC = self.mesh.nC
else:
self.nC = len(self.actInd)
slope = 1e4
@property
def shape(self):
return (self.nC, self.nP)
@property
def nP(self):
if np.isscalar(self.order):
@@ -785,8 +797,8 @@ class PolyMap(IdentityMap):
sig1, sig2 = np.exp(sig1), np.exp(sig2)
#2D
if self.mesh.dim == 2:
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
X = self.mesh.gridCC[self.actInd,0]
Y = self.mesh.gridCC[self.actInd,1]
if self.normal =='X':
f = polynomial.polyval(Y, c) - X
elif self.normal =='Y':
@@ -795,9 +807,9 @@ class PolyMap(IdentityMap):
raise(Exception("Input for normal = X or Y or Z"))
#3D
elif self.mesh.dim == 3:
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
Z = self.mesh.gridCC[:,2]
X = self.mesh.gridCC[self.actInd,0]
Y = self.mesh.gridCC[self.actInd,1]
Z = self.mesh.gridCC[self.actInd,2]
if self.normal =='X':
f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
elif self.normal =='Y':
@@ -806,6 +818,7 @@ class PolyMap(IdentityMap):
f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z
else:
raise(Exception("Input for normal = X or Y or Z"))
else:
raise(Exception("Only supports 2D"))
@@ -819,8 +832,8 @@ class PolyMap(IdentityMap):
sig1, sig2 = np.exp(sig1), np.exp(sig2)
#2D
if self.mesh.dim == 2:
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
X = self.mesh.gridCC[self.actInd,0]
Y = self.mesh.gridCC[self.actInd,1]
if self.normal =='X':
f = polynomial.polyval(Y, c) - X
@@ -832,9 +845,9 @@ class PolyMap(IdentityMap):
raise(Exception("Input for normal = X or Y or Z"))
#3D
elif self.mesh.dim == 3:
X = self.mesh.gridCC[:,0]
Y = self.mesh.gridCC[:,1]
Z = self.mesh.gridCC[:,2]
X = self.mesh.gridCC[self.actInd,0]
Y = self.mesh.gridCC[self.actInd,1]
Z = self.mesh.gridCC[self.actInd,2]
if self.normal =='X':
f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X