wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-30 15:18:22 +08:00
Code Issues Packages Projects Releases Wiki Activity

bug fixes to framework.

Browse Source
This commit is contained in:
Rowan Cockett
2013-10-23 18:19:33 -07:00
parent 9cd6e7ed08
commit 548ba9aa0a
5 changed files with 139 additions and 68 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/forward/DCProblem.py
+35 -32
View File
@@ -3,7 +3,7 @@ from SimPEG.forward import Problem, SyntheticProblem
from SimPEG.tests import checkDerivative
from SimPEG.utils import ModelBuilder, sdiag
import numpy as np
import scipy as sp
import scipy.sparse as sp
import scipy.sparse.linalg as linalg
class DCProblem(Problem):
@@ -104,14 +104,44 @@ class DCProblem(Problem):
return Jtv
def genTxRxmat(nelec, spacelec, surfloc, elecini, mesh):
""" Generate projection matrix (Q) and """
elecend = 0.5+spacelec*(nelec-1)
elecLocR = np.linspace(elecini, elecend, nelec)
elecLocT = elecLocR+1
nrx = nelec-1
ntx = nelec-1
q = np.zeros((mesh.nC, ntx))
Q = np.zeros((mesh.nC, nrx))
for i in range(nrx):
rxind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i]))
rxind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i+1]))
txind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i]))
txind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i+1]))
q[txind1,i] = 1
q[txind2,i] = -1
Q[rxind1,i] = 1
Q[rxind2,i] = -1
Q = sp.csr_matrix(Q)
rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5
return q, Q, rxmidLoc
if __name__ == '__main__':
from SimPEG.regularization import Regularization
from SimPEG import inverse
# Create the mesh
h1 = np.ones(100)
h2 = np.ones(100)
h1 = np.ones(10)
h2 = np.ones(10)
mesh = TensorMesh([h1,h2])
# Create some parameters for the model
@@ -153,14 +183,14 @@ if __name__ == '__main__':
problem = DCProblem(mesh)
problem.P = P
problem.RHS = q
problem.W = Wd
problem.dobs = dobs
problem.std = dobs*0 + 0.05
m0 = mesh.gridCC[:,0]*0+sig1
# print problem.misfit(m0)
# print problem.misfit(mSynth)
opt = inverse.InexactGaussNewton()
opt = inverse.InexactGaussNewton(maxIterLS=20)
reg = Regularization(mesh)
inv = inverse.Inversion(problem, reg, opt)
@@ -181,30 +211,3 @@ if __name__ == '__main__':
def genTxRxmat(nelec, spacelec, surfloc, elecini, mesh):
""" Generate projection matrix (Q) and """
elecend = 0.5+spacelec*(nelec-1)
elecLocR = np.linspace(elecini, elecend, nelec)
elecLocT = elecLocR+1
nrx = nelec-1
ntx = nelec-1
q = np.zeros((mesh.nC, ntx))
Q = np.zeros((mesh.nC, nrx))
for i in range(nrx):
rxind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i]))
rxind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i+1]))
txind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i]))
txind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i+1]))
q[txind1,i] = 1
q[txind2,i] = -1
Q[rxind1,i] = 1
Q[rxind2,i] = -1
Q = sp.csr_matrix(Q)
rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5
return q, Q, rxmidLoc
SimPEG/inverse/Inversion.py
+23 -22
View File
@@ -1,5 +1,6 @@
import numpy as np
from SimPEG.utils import sdiag
import scipy.sparse as sp
from SimPEG.utils import sdiag, mkvc
class Inversion(object):
"""docstring for Inversion"""
@@ -16,7 +17,10 @@ class Inversion(object):
"""
Standard deviation weighting matrix.
"""
return sdiag(1/self.prob.std)
if getattr(self,'_Wd',None) is None:
eps = np.linalg.norm(mkvc(self.prob.dobs),2)*1e-5
self._Wd = 1/(abs(self.prob.dobs)*self.prob.std+eps)
return self._Wd
def run(self, m0):
m = m0
@@ -41,7 +45,7 @@ class Inversion(object):
u = self.prob.field(m)
phi_d = self.dataObj(m, u)
phi_m = self.modelObj(m)
phi_m = self.reg.modelObj(m)
self._phi_d_last = phi_d
self._phi_m_last = phi_m
@@ -50,27 +54,24 @@ class Inversion(object):
out = (f,)
if return_g:
phi_dDeriv = self.dataObjDeriv(m, u)
phi_mDeriv = self.modelObjDeriv(m)
phi_dDeriv = self.dataObjDeriv(m, u=u)
phi_mDeriv = self.reg.modelObjDeriv(m)
g = phi_dDeriv + self._beta * phi_mDeriv
out += (g,)
if return_H:
def H_fun(v):
phi_d2Deriv = self.dataObj2Deriv(m, u, v)
phi_m2Deriv = self.modelObj2Deriv(m)*v
phi_d2Deriv = self.dataObj2Deriv(m, v, u=u)
phi_m2Deriv = self.reg.modelObj2Deriv(m)*v
return phi_d2Deriv + self._beta * phi_m2Deriv
out += (H_fun,)
operator = sp.linalg.LinearOperator( (m.size, m.size), H_fun, dtype=float )
out += (operator,)
return out
def modelObj(self, m, u=None):
self.reg.misfit(m)
def dataObj(self, m, u=None):
"""
:param numpy.array m: geophysical model
@@ -87,7 +88,7 @@ class Inversion(object):
Where P is a projection matrix that brings the field on the full domain to the data measurement locations;
u is the field of interest; d_obs is the observed data; and W is the weighting matrix.
"""
R = self.Wd*self.prob.misfit(u=u)
R = self.Wd*self.prob.misfit(m, u=u)
R = mkvc(R)
return 0.5*R.dot(R)
@@ -123,17 +124,17 @@ class Inversion(object):
"""
if u is None:
u = self.field(m)
u = self.prob.field(m)
R = self.W*(self.dpred(m, u=u) - self.dobs)
R = self.Wd*self.prob.misfit(m, u=u)
dmisfit = 0
for i in range(self.RHS.shape[1]): # Loop over each right hand side
dmisfit += self.Jt(m, self.W[:,i]*R[:,i], u=u[:,i])
for i in range(self.prob.RHS.shape[1]): # Loop over each right hand side
dmisfit += self.prob.Jt(m, self.Wd[:,i]*R[:,i], u=u[:,i])
return dmisfit
def dataObj2Deriv(self, m, u=None):
def dataObj2Deriv(self, m, v, u=None):
"""
:param numpy.array m: geophysical model
:param numpy.array u: fields
@@ -167,12 +168,12 @@ class Inversion(object):
"""
if u is None:
u = self.field(m)
u = self.prob.field(m)
R = self.W*(self.dpred(m, u=u) - self.dobs)
R = self.Wd*self.prob.misfit(m, u=u)
dmisfit = 0
for i in range(self.RHS.shape[1]): # Loop over each right hand side
dmisfit += self.Jt(m, self.W[:,i]*R[:,i], u=u[:,i])
for i in range(self.prob.RHS.shape[1]): # Loop over each right hand side
dmisfit += self.prob.Jt(m, self.Wd[:,i] * self.Wd[:,i] * self.prob.J(m, v, u=u[:,i]), u=u[:,i])
return dmisfit
SimPEG/inverse/Optimize.py
+3 -1
View File
@@ -2,6 +2,7 @@ import numpy as np
import matplotlib.pyplot as plt
from SimPEG.utils import mkvc, sdiag
norm = np.linalg.norm
import scipy.sparse as sp
class Minimize(object):
@@ -149,7 +150,8 @@ class GaussNewton(Minimize):
class InexactGaussNewton(Minimize):
name = 'InexactGaussNewton'
def findSearchDirection(self):
return sparse.linalg.cg(self.H, -self.g, tol=1e-05, maxiter=10)
p, info = sp.linalg.cg(self.H, -self.g, tol=1e-05, maxiter=10)
return p
class SteepestDescent(Minimize):
SimPEG/mesh/DiffOperators.py
+62
View File
@@ -161,6 +161,68 @@ class DiffOperators(object):
_cellGrad = None
cellGrad = property(**cellGrad())
def cellGradx():
doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
def fget(self):
if getattr(self, '_cellGradx', None) is None:
BC = ['neumann', 'neumann']
n = self.n
if(self.dim == 1):
G1 = ddxCellGrad(n[0], BC)
elif(self.dim == 2):
G1 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC))
elif(self.dim == 3):
G1 = kron3(speye(n[2]), speye(n[1]), ddxCellGrad(n[0], BC))
# Compute areas of cell faces & volumes
S = self.r(self.area, 'F','Fx', 'V')
V = self.vol
self._cellGradx = sdiag(S)*G1*sdiag(1/V)
return self._cellGradx
return locals()
cellGradx = property(**cellGradx())
def cellGrady():
doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
def fget(self):
if self.dim < 2:
return None
if getattr(self, '_cellGrady', None) is None:
BC = ['neumann', 'neumann']
n = self.n
if(self.dim == 2):
G2 = sp.kron(speye(n[1]), ddxCellGrad(n[0], BC))
elif(self.dim == 3):
G2 = kron3(speye(n[2]), ddxCellGrad(n[1], BC), speye(n[0]))
# Compute areas of cell faces & volumes
S = self.r(self.area, 'F','Fy', 'V')
V = self.vol
self._cellGrady = sdiag(S)*G2*sdiag(1/V)
return self._cellGrady
return locals()
cellGrady = property(**cellGrady())
def cellGradz():
doc = "Cell centered Gradient in the x dimension. Has neumann boundary conditions."
def fget(self):
if self.dim < 3:
return None
if getattr(self, '_cellGradz', None) is None:
BC = ['neumann', 'neumann']
n = self.n
G3 = kron3(ddxCellGrad(n[2], BC), speye(n[1]), speye(n[0]))
# Compute areas of cell faces & volumes
S = self.r(self.area, 'F','Fz', 'V')
V = self.vol
self._cellGradz = sdiag(S)*G3*sdiag(1/V)
return self._cellGradz
return locals()
cellGradz = property(**cellGradz())
def edgeCurl():
doc = "Construct the 3D curl operator."
SimPEG/regularization/Regularization.py
+16 -13
View File
@@ -1,10 +1,13 @@
from SimPEG.utils import sdiag
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):
@@ -12,25 +15,25 @@ class Regularization(object):
@property
def Wx(self):
if self._Wx is None:
self._Wx = mesh.cellGradx
if getattr(self, '_Wx', None) is None:
self._Wx = self.mesh.cellGradx
return self._Wx
@property
def Wy(self):
if self._Wy is None:
self._Wy = mesh.cellGrady
if getattr(self, '_Wy', None) is None:
self._Wy = self.mesh.cellGrady
return self._Wy
@property
def Wz(self):
if self._Wz is None:
self._Wz = mesh.cellGradz
if getattr(self, '_Wz', None) is None:
self._Wz = self.mesh.cellGradz
return self._Wz
@property
def Ws(self):
if self._Ws is None:
if getattr(self,'_Ws', None) is None:
self._Ws = sdiag(self.mesh.vol)
return self._Ws
@@ -87,12 +90,12 @@ class Regularization(object):
mobjDeriv = self.alpha_s * self.Ws.T * ( self.Ws * mresid)
mobjDeriv += self.alpha_x * self.Wx.T * ( self.Wx * mresid)
mobjDeriv = mobjDeriv + self.alpha_x * self.Wx.T * ( self.Wx * mresid)
if self.mesh.dim > 1:
mobjDeriv += self.alpha_y * self.Wy.T * ( self.Wy * mresid)
mobjDeriv = mobjDeriv + self.alpha_y * self.Wy.T * ( self.Wy * mresid)
if self.mesh.dim > 2:
mobjDeriv += self.alpha_z * self.Wz.T * ( self.Wz * mresid)
mobjDeriv = mobjDeriv + self.alpha_z * self.Wz.T * ( self.Wz * mresid)
return mobjDeriv
@@ -102,12 +105,12 @@ class Regularization(object):
mobj2Deriv = self.alpha_s * self.Ws.T * self.Ws
mobj2Deriv += self.alpha_x * self.Wx.T * self.Wx
mobj2Deriv = mobj2Deriv + self.alpha_x * self.Wx.T * self.Wx
if self.mesh.dim > 1:
mobj2Deriv += self.alpha_y * self.Wy.T * self.Wy
mobj2Deriv = mobj2Deriv + self.alpha_y * self.Wy.T * self.Wy
if self.mesh.dim > 2:
mobj2Deriv += self.alpha_z * self.Wz.T * self.Wz
mobj2Deriv = mobj2Deriv + self.alpha_z * self.Wz.T * self.Wz
return mobj2Deriv