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

DC inversion working again. Almost in the new framework..!

Browse Source
This commit is contained in:
rowanc1
2014-01-16 14:29:08 -08:00
parent fa8a5cd7cb
commit 687ef0c20a
6 changed files with 73 additions and 52 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/Data.py
+29 -8
View File
@@ -1,5 +1,5 @@
import Utils
import numpy as np
def requiresProblem(f):
"""
@@ -13,12 +13,13 @@ def requiresProblem(f):
If a problem is not bound an Exception will be raised, and an nice error message printed.
"""
extra = """
This function requires that a problem be bound to the data.
To use data.%s(), SimPEG requires that a problem be bound to the data.
If a problem has not been bound, an Exception will be raised.
To bind a problem to the Data object::
data.setProblem(myProblem)
"""
""" % f.__name__
from functools import wraps
@wraps(f)
def requiresProblemWrapper(self,*args,**kwargs):
@@ -37,11 +38,11 @@ class BaseData(object):
__metaclass__ = Utils.Save.Savable
prob = None #: The geophysical problem that explains this data, use data.setProblem(prob)
std = None #: Estimated Standard Deviations
dobs = None #: Observed data
dtrue = None #: True data, if data is synthetic
mtrue = None #: True model, if data is synthetic
prob = None #: The geophysical problem that explains this data
counter = None #: A SimPEG.Utils.Counter object
@@ -49,19 +50,39 @@ class BaseData(object):
Utils.setKwargs(self, **kwargs)
def setProblem(self, prob):
# Bind these two instances together using pointers
self.prob = prob
prob.data = self
@Utils.count
@requiresProblem
def dpred(self, m, u=None):
"""
Projection matrix.
Create the projected data from a model.
The field, u, (if provided) will be used for the predicted data
instead of recalculating the fields (which may be expensive!).
.. math::
d_\\text{pred} = P(u(m))
Where P is a projection of the fields onto the data space.
"""
if u is None: u = self.prob.field(m)
return self.projectField(u)
@Utils.count
def projectField(self, u):
"""
This function projects the fields onto the data space.
.. math::
d_\\text{pred} = Pu(m)
"""
if u is None: u = self.prob.field(m)
return self.P*u
return u
#TODO: def projectFieldDeriv(self, u): Does this need to be made??!
@Utils.count
def residual(self, m, u=None):
@@ -122,7 +143,7 @@ class BaseData(object):
"""
Source matrix.
"""
return self._RHS
return getattr(self, '_RHS', None)
@RHS.setter
def RHS(self, value):
self._RHS = value
SimPEG/Examples/DC.py
+28 -31
View File
@@ -10,9 +10,10 @@ class DCData(Data.BaseData):
"""
def __init__(self, mesh, model, **kwargs):
problem.BaseProblem.__init__(self, mesh, model)
self.mesh.setCellGradBC('neumann')
P = None #: projection
def __init__(self, **kwargs):
Data.BaseData.__init__(self, **kwargs)
Utils.setKwargs(self, **kwargs)
def reshapeFields(self, u):
@@ -20,18 +21,14 @@ class DCData(Data.BaseData):
u = u.reshape([-1, self.RHS.shape[1]], order='F')
return u
def dpred(self, m, u=None):
def projectField(self, u):
"""
Predicted data.
.. math::
d_\\text{pred} = Pu(m)
"""
if u is None:
u = self.field(m)
u = self.reshapeFields(u)
return Utils.mkvc(self.P*u)
@@ -47,7 +44,7 @@ class DCProblem(Problem.BaseProblem):
dataPair = DCData
def __init__(self, mesh, model, **kwargs):
problem.BaseProblem.__init__(self, mesh, model)
Problem.BaseProblem.__init__(self, mesh, model)
self.mesh.setCellGradBC('neumann')
Utils.setKwargs(self, **kwargs)
@@ -75,7 +72,7 @@ class DCProblem(Problem.BaseProblem):
def field(self, m):
A = self.createMatrix(m)
solve = Solver(A)
phi = solve.solve(self.RHS)
phi = solve.solve(self.data.RHS)
return Utils.mkvc(phi)
def J(self, m, v, u=None):
@@ -103,9 +100,9 @@ class DCProblem(Problem.BaseProblem):
if u is None:
u = self.field(m)
u = self.reshapeFields(u)
u = self.data.reshapeFields(u)
P = self.P
P = self.data.P
D = self.mesh.faceDiv
G = self.mesh.cellGrad
A = self.createMatrix(m)
@@ -128,10 +125,10 @@ class DCProblem(Problem.BaseProblem):
if u is None:
u = self.field(m)
u = self.reshapeFields(u)
v = self.reshapeFields(v)
u = self.data.reshapeFields(u)
v = self.data.reshapeFields(v)
P = self.P
P = self.data.P
D = self.mesh.faceDiv
G = self.mesh.cellGrad
A = self.createMatrix(m)
@@ -186,7 +183,7 @@ if __name__ == '__main__':
# Create the mesh
h1 = np.ones(20)
h2 = np.ones(100)
M = mesh.TensorMesh([h1,h2])
M = Mesh.TensorMesh([h1,h2])
# Create some parameters for the model
sig1 = np.log(1)
@@ -198,7 +195,7 @@ if __name__ == '__main__':
condVals = [sig1, sig2]
mSynth = Utils.ModelBuilder.defineBlockConductivity(p0,p1,M.gridCC,condVals)
plt.colorbar(M.plotImage(mSynth))
plt.show()
# plt.show()
# Set up the projection
nelec = 50
@@ -211,29 +208,29 @@ if __name__ == '__main__':
q, Q, rxmidloc = genTxRxmat(nelec, spacelec, surfloc, elecini, M)
P = Q.T
# Create some data
problem = DCProblem(M)
problem.P = P
problem.RHS = q
data = problem.createSyntheticData(mSynth, std=0.05)
model = Model.LogModel()
prob = DCProblem(M, model)
u = problem.field(mSynth)
u = problem.reshapeFields(u)
# Create some data
data = prob.createSyntheticData(mSynth, std=0.05, P=P, RHS=q)
u = prob.field(mSynth)
u = data.reshapeFields(u)
M.plotImage(u[:,10])
# plt.show()
# Now set up the problem to do some minimization
# problem.dobs = dobs
# problem.std = dobs*0 + 0.05
# Now set up the prob to do some minimization
# prob.dobs = dobs
# prob.std = dobs*0 + 0.05
m0 = M.gridCC[:,0]*0+sig2
opt = inverse.InexactGaussNewton(maxIterLS=20, maxIter=3, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
reg = inverse.Regularization(M)
inv = inverse.Inversion(problem, reg, opt, data, beta0=1e4)
opt = Inverse.InexactGaussNewton(maxIterLS=20, maxIter=3, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
reg = Inverse.Regularization(M)
inv = Inverse.Inversion(prob, reg, opt, data, beta0=1e4)
# Check Derivative
derChk = lambda m: [inv.dataObj(m), inv.dataObjDeriv(m)]
tests.checkDerivative(derChk, mSynth)
# Tests.checkDerivative(derChk, mSynth)
SimPEG/Inverse/Inversion.py
+6 -6
View File
@@ -201,7 +201,7 @@ class BaseInversion(object):
phi_d = self.dataObj(m, u)
phi_m = self.reg.modelObj(m)
self.dpred = self.prob.dpred(m, u=u) # This is a cheap matrix vector calculation.
self.dpred = self.data.dpred(m, u=u) # This is a cheap matrix vector calculation.
self.phi_d = phi_d
self.phi_m = phi_m
@@ -245,7 +245,7 @@ class BaseInversion(object):
u is the field of interest; d_obs is the observed data; and W is the weighting matrix.
"""
# TODO: ensure that this is a data is vector and Wd is a matrix.
R = self.Wd*self.prob.dataResidual(m, self.data, u=u)
R = self.data.residualWeighted(m, u=u)
R = Utils.mkvc(R)
return 0.5*np.vdot(R, R)
@@ -285,9 +285,9 @@ class BaseInversion(object):
if u is None:
u = self.prob.field(m)
R = self.Wd*self.prob.dataResidual(m, self.data, u=u)
R = self.data.residualWeighted(m, u=u)
dmisfit = self.prob.Jt(m, self.Wd * R, u=u)
dmisfit = self.prob.Jt(m, self.data.Wd * R, u=u)
return dmisfit
@@ -330,11 +330,11 @@ class BaseInversion(object):
if u is None:
u = self.prob.field(m)
R = self.Wd*self.prob.dataResidual(m, self.data, u=u)
R = self.data.residualWeighted(m, u=u)
# TODO: abstract to different norms a little cleaner.
# \/ it goes here. in l2 it is the identity.
dmisfit = self.prob.Jt_approx(m, self.Wd * self.Wd * self.prob.J_approx(m, v, u=u), u=u)
dmisfit = self.prob.Jt_approx(m, self.data.Wd * self.data.Wd * self.prob.J_approx(m, v, u=u), u=u)
return dmisfit
SimPEG/Problem.py
+8 -6
View File
@@ -130,7 +130,7 @@ class BaseProblem(object):
"""
pass
def createSyntheticData(self, m, std=0.05, u=None):
def createSyntheticData(self, m, std=0.05, u=None, **geometry_kwargs):
"""
Create synthetic data given a model, and a standard deviation.
@@ -142,11 +142,13 @@ class BaseProblem(object):
Returns the observed data with random Gaussian noise
and Wd which is the same size as data, and can be used to weight the inversion.
"""
dtrue = self.dpred(m,u=u)
noise = std*abs(dtrue)*np.random.randn(*dtrue.shape)
dobs = dtrue+noise
stdev = dobs*0 + std
return self.dataPair(dobs=dobs, std=stdev, dtrue=dtrue, mtrue=m)
data = self.dataPair(mtrue=m, **geometry_kwargs)
data.setProblem(self)
data.dtrue = self.data.dpred(m,u=u)
noise = std*abs(data.dtrue)*np.random.randn(*data.dtrue.shape)
data.dobs = data.dtrue+noise
data.std = data.dobs*0 + std
return data
SimPEG/Utils/ModelBuilder.py
+1 -1
View File
@@ -196,7 +196,7 @@ def randomModel(shape, seed=None, anisotropy=None, its=100, bounds=[0,1]):
if __name__ == '__main__':
from SimPEG.mesh import TensorMesh
from SimPEG.Mesh import TensorMesh
from matplotlib import pyplot as plt
# Define the mesh
SimPEG/Utils/__init__.py
+1
View File
@@ -7,6 +7,7 @@ from ipythonutils import easyAnimate as animate
from Solver import Solver
import Save
import Geophysics
import ModelBuilder
import types
import time