mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-17 11:32:59 +08:00
Removed Examples dir, and put it up a level (outside SimPEG path) renamed to tutorials for future tutorial development.
This commit is contained in:
rowanc1
@@ -1,249 +0,0 @@
|
||||
from SimPEG import *
|
||||
|
||||
|
||||
|
||||
class DCData(Data.BaseData):
|
||||
"""
|
||||
**DCData**
|
||||
|
||||
Geophysical DC resistivity data.
|
||||
|
||||
"""
|
||||
|
||||
P = None #: projection
|
||||
|
||||
def __init__(self, **kwargs):
|
||||
Data.BaseData.__init__(self, **kwargs)
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
|
||||
def reshapeFields(self, u):
|
||||
if len(u.shape) == 1:
|
||||
u = u.reshape([-1, self.RHS.shape[1]], order='F')
|
||||
return u
|
||||
|
||||
def projectField(self, u):
|
||||
"""
|
||||
Predicted data.
|
||||
|
||||
.. math::
|
||||
d_\\text{pred} = Pu(m)
|
||||
"""
|
||||
u = self.reshapeFields(u)
|
||||
return Utils.mkvc(self.P*u)
|
||||
|
||||
|
||||
|
||||
class DCProblem(Problem.BaseProblem):
|
||||
"""
|
||||
**DCProblem**
|
||||
|
||||
Geophysical DC resistivity problem.
|
||||
|
||||
"""
|
||||
|
||||
dataPair = DCData
|
||||
|
||||
def __init__(self, mesh, model, **kwargs):
|
||||
Problem.BaseProblem.__init__(self, mesh, model)
|
||||
self.mesh.setCellGradBC('neumann')
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
|
||||
|
||||
def createMatrix(self, m):
|
||||
"""
|
||||
Makes the matrix A(m) for the DC resistivity problem.
|
||||
|
||||
:param numpy.array m: model
|
||||
:rtype: scipy.csc_matrix
|
||||
:return: A(m)
|
||||
|
||||
.. math::
|
||||
c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0
|
||||
|
||||
Where M() is the mass matrix and mT is the model transform.
|
||||
"""
|
||||
D = self.mesh.faceDiv
|
||||
G = self.mesh.cellGrad
|
||||
sigma = self.model.transform(m)
|
||||
Msig = self.mesh.getFaceMass(sigma)
|
||||
A = D*Msig*G
|
||||
return A.tocsc()
|
||||
|
||||
def fields(self, m):
|
||||
A = self.createMatrix(m)
|
||||
solve = Solver(A)
|
||||
phi = solve.solve(self.data.RHS)
|
||||
return Utils.mkvc(phi)
|
||||
|
||||
def J(self, m, v, u=None):
|
||||
"""
|
||||
:param numpy.array m: model
|
||||
:param numpy.array v: vector to multiply
|
||||
:param numpy.array u: fields
|
||||
:rtype: numpy.array
|
||||
:return: Jv
|
||||
|
||||
.. math::
|
||||
c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0
|
||||
|
||||
\\nabla_u (A(m)u - q) = A(m)
|
||||
|
||||
\\nabla_m (A(m)u - q) = G\\text{sdiag}(Du)\\nabla_m(M(mT(m)))
|
||||
|
||||
Where M() is the mass matrix and mT is the model transform.
|
||||
|
||||
.. math::
|
||||
J = - P \left( \\nabla_u c(m, u) \\right)^{-1} \\nabla_m c(m, u)
|
||||
|
||||
J(v) = - P ( A(m)^{-1} ( G\\text{sdiag}(Du)\\nabla_m(M(mT(m))) v ) )
|
||||
"""
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
|
||||
u = self.data.reshapeFields(u)
|
||||
|
||||
P = self.data.P
|
||||
D = self.mesh.faceDiv
|
||||
G = self.mesh.cellGrad
|
||||
A = self.createMatrix(m)
|
||||
Av_dm = self.mesh.getFaceMassDeriv()
|
||||
mT_dm = self.model.transformDeriv(m)
|
||||
|
||||
dCdu = A
|
||||
|
||||
dCdm = np.empty_like(u)
|
||||
for i, ui in enumerate(u.T): # loop over each column
|
||||
dCdm[:, i] = D * ( Utils.sdiag( G * ui ) * ( Av_dm * ( mT_dm * v ) ) )
|
||||
|
||||
solve = Solver(dCdu)
|
||||
Jv = - P * solve.solve(dCdm)
|
||||
return Utils.mkvc(Jv)
|
||||
|
||||
def Jt(self, m, v, u=None):
|
||||
"""Takes data, turns it into a model..ish"""
|
||||
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
|
||||
u = self.data.reshapeFields(u)
|
||||
v = self.data.reshapeFields(v)
|
||||
|
||||
P = self.data.P
|
||||
D = self.mesh.faceDiv
|
||||
G = self.mesh.cellGrad
|
||||
A = self.createMatrix(m)
|
||||
Av_dm = self.mesh.getFaceMassDeriv()
|
||||
mT_dm = self.model.transformDeriv(m)
|
||||
|
||||
dCdu = A.T
|
||||
solve = Solver(dCdu)
|
||||
|
||||
w = solve.solve(P.T*v)
|
||||
|
||||
Jtv = 0
|
||||
for i, ui in enumerate(u.T): # loop over each column
|
||||
Jtv += Utils.sdiag( G * ui ) * ( D.T * w[:,i] )
|
||||
|
||||
Jtv = - mT_dm.T * ( Av_dm.T * Jtv )
|
||||
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__':
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
# Create the mesh
|
||||
h1 = np.ones(20)
|
||||
h2 = np.ones(100)
|
||||
M = Mesh.TensorMesh([h1,h2])
|
||||
|
||||
# Create some parameters for the model
|
||||
sig1 = np.log(1)
|
||||
sig2 = np.log(0.01)
|
||||
|
||||
# Create a synthetic model from a block in a half-space
|
||||
p0 = [5, 10]
|
||||
p1 = [15, 50]
|
||||
condVals = [sig1, sig2]
|
||||
mSynth = Utils.ModelBuilder.defineBlockConductivity(M.gridCC,p0,p1,condVals)
|
||||
plt.colorbar(M.plotImage(mSynth))
|
||||
# plt.show()
|
||||
|
||||
# Set up the projection
|
||||
nelec = 50
|
||||
spacelec = 2
|
||||
surfloc = 0.5
|
||||
elecini = 0.5
|
||||
elecend = 0.5+spacelec*(nelec-1)
|
||||
elecLocR = np.linspace(elecini, elecend, nelec)
|
||||
rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5
|
||||
q, Q, rxmidloc = genTxRxmat(nelec, spacelec, surfloc, elecini, M)
|
||||
P = Q.T
|
||||
|
||||
model = Model.LogModel(M)
|
||||
prob = DCProblem(M, model)
|
||||
|
||||
# Create some data
|
||||
data = prob.createSyntheticData(mSynth, std=0.05, P=P, RHS=q)
|
||||
|
||||
u = prob.fields(mSynth)
|
||||
u = data.reshapeFields(u)
|
||||
M.plotImage(u[:,10])
|
||||
plt.show()
|
||||
|
||||
# Now set up the prob to do some minimization
|
||||
# prob.dobs = dobs
|
||||
# prob.std = dobs*0 + 0.05
|
||||
m0 = M.gridCC[:,0]*0+sig2
|
||||
|
||||
reg = Regularization.Tikhonov(model)
|
||||
objFunc = ObjFunction.BaseObjFunction(data, reg)
|
||||
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=3, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
|
||||
inv = Inversion.BaseInversion(objFunc, opt)
|
||||
|
||||
# Check Derivative
|
||||
derChk = lambda m: [objFunc.dataObj(m), objFunc.dataObjDeriv(m)]
|
||||
# Tests.checkDerivative(derChk, mSynth)
|
||||
|
||||
print objFunc.dataObj(m0)
|
||||
print objFunc.dataObj(mSynth)
|
||||
|
||||
m = inv.run(m0)
|
||||
|
||||
plt.colorbar(M.plotImage(m))
|
||||
print m
|
||||
plt.show()
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -1,2 +0,0 @@
|
||||
import DC
|
||||
import Linear
|
||||
@@ -1,85 +0,0 @@
|
||||
# import numpy as np
|
||||
# import unittest
|
||||
# from SimPEG.mesh import TensorMesh
|
||||
# from SimPEG.Utils import ModelBuilder, sdiag
|
||||
# from SimPEG.forward import Problem
|
||||
# from SimPEG.examples.DC import *
|
||||
# from TestUtils import checkDerivative
|
||||
# from scipy.sparse.linalg import dsolve
|
||||
# from SimPEG import inverse
|
||||
|
||||
|
||||
# class DCProblemTests(unittest.TestCase):
|
||||
|
||||
# def setUp(self):
|
||||
# # Create the mesh
|
||||
# h1 = np.ones(20)
|
||||
# h2 = np.ones(20)
|
||||
# mesh = TensorMesh([h1,h2])
|
||||
|
||||
# # Create some parameters for the model
|
||||
# sig1 = 1
|
||||
# sig2 = 0.01
|
||||
|
||||
# # Create a synthetic model from a block in a half-space
|
||||
# p0 = [2, 2]
|
||||
# p1 = [5, 5]
|
||||
# condVals = [sig1, sig2]
|
||||
# mSynth = ModelBuilder.defineBlockConductivity(p0,p1,mesh.gridCC,condVals)
|
||||
|
||||
# # Set up the projection
|
||||
# nelec = 10
|
||||
# spacelec = 2
|
||||
# surfloc = 0.5
|
||||
# elecini = 0.5
|
||||
# elecend = 0.5+spacelec*(nelec-1)
|
||||
# elecLocR = np.linspace(elecini, elecend, nelec)
|
||||
# rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5
|
||||
# q, Q, rxmidloc = genTxRxmat(nelec, spacelec, surfloc, elecini, mesh)
|
||||
# P = Q.T
|
||||
|
||||
# # Create some data
|
||||
|
||||
# problem = DCProblem(mesh)
|
||||
# problem.P = P
|
||||
# problem.RHS = q
|
||||
# data = problem.createSyntheticData(mSynth, std=0.05)
|
||||
|
||||
# # Now set up the problem to do some minimization
|
||||
# opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
|
||||
# reg = inverse.Regularization(mesh)
|
||||
# inv = inverse.Inversion(problem, reg, opt, data, beta0=1e4)
|
||||
|
||||
# self.inv = inv
|
||||
# self.reg = reg
|
||||
# self.p = problem
|
||||
# self.mesh = mesh
|
||||
# self.m0 = mSynth
|
||||
# self.data = data
|
||||
|
||||
# def test_misfit(self):
|
||||
# derChk = lambda m: [self.p.dpred(m), lambda mx: self.p.J(self.m0, mx)]
|
||||
# passed = checkDerivative(derChk, self.m0, plotIt=False)
|
||||
# self.assertTrue(passed)
|
||||
|
||||
# def test_adjoint(self):
|
||||
# # Adjoint Test
|
||||
# u = np.random.rand(self.mesh.nC*self.p.RHS.shape[1])
|
||||
# v = np.random.rand(self.mesh.nC)
|
||||
# w = np.random.rand(self.data.dobs.shape[0])
|
||||
# wtJv = w.dot(self.p.J(self.m0, v, u=u))
|
||||
# vtJtw = v.dot(self.p.Jt(self.m0, w, u=u))
|
||||
# passed = (wtJv - vtJtw) < 1e-10
|
||||
# self.assertTrue(passed)
|
||||
|
||||
# def test_dataObj(self):
|
||||
# derChk = lambda m: [self.inv.dataObj(m), self.inv.dataObjDeriv(m)]
|
||||
# checkDerivative(derChk, self.m0, plotIt=False)
|
||||
|
||||
# def test_modelObj(self):
|
||||
# derChk = lambda m: [self.reg.modelObj(m), self.reg.modelObjDeriv(m)]
|
||||
# checkDerivative(derChk, self.m0, plotIt=False)
|
||||
|
||||
|
||||
# if __name__ == '__main__':
|
||||
# unittest.main()
|
||||
@@ -11,7 +11,6 @@ import ObjFunction
|
||||
import Optimization
|
||||
import Inversion
|
||||
import Parameters
|
||||
import Examples
|
||||
import Tests
|
||||
|
||||
|
||||
|
||||
Reference in New Issue
Block a user