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

Tweak the framework.

Browse Source
This commit is contained in:
rowanc1
2014-05-18 10:57:06 -07:00
parent 2392b7d6b9
commit 8325919e63
8 changed files with 102 additions and 105 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/DataMisfit.py
+44 -59
View File
@@ -1,20 +1,6 @@
import Utils, Survey, Problem, numpy as np, scipy.sparse as sp, gc
def _splitForward(forward):
assert forward.ispaired, 'The problem and survey must be paired.'
if isinstance(forward, Survey.BaseSurvey):
survey = forward
prob = forward.prob
elif isinstance(forward, Problem.BaseProblem):
prob = forward
survey = forward.survey
else:
raise Exception('The forward simulation must either be a problem or a survey.')
return prob, survey
class BaseDataMisfit(object):
"""BaseDataMisfit
@@ -28,25 +14,16 @@ class BaseDataMisfit(object):
debug = False #: Print debugging information
counter = None #: Set this to a SimPEG.Utils.Counter() if you want to count things
def __init__(self):
pass
def splitForward(self, forward):
"""splitForward(forward)
Split the forward simulation into a problem and a survey
:param Problem,Survey forward: forward simulation
:rtype: Problem,Survey
:return: (prob, survey)
"""
prob, survey = _splitForward(forward)
return prob, survey
def __init__(self, survey, **kwargs):
assert survey.ispaired, 'The survey must be paired to a problem.'
if isinstance(survey, Survey.BaseSurvey):
self.survey = survey
self.prob = survey.prob
Utils.setKwargs(self,**kwargs)
@Utils.timeIt
def dataObj(self, forward, m, u=None):
"""dataObj(forward, m, u=None)
def dataObj(self, m, u=None):
"""dataObj(m, u=None)
:param Problem,Survey forward: forward simulation
:param numpy.array m: geophysical model
@@ -58,8 +35,8 @@ class BaseDataMisfit(object):
raise NotImplementedError('This method should be overwritten.')
@Utils.timeIt
def dataObjDeriv(self, forward, m, u=None):
"""dataObjDeriv(forward, m, u=None)
def dataObjDeriv(self, m, u=None):
"""dataObjDeriv(m, u=None)
:param Problem,Survey forward: forward simulation
:param numpy.array m: geophysical model
@@ -72,8 +49,8 @@ class BaseDataMisfit(object):
@Utils.timeIt
def dataObj2Deriv(self, forward, m, v, u=None):
"""dataObj2Deriv(forward, m, v, u=None)
def dataObj2Deriv(self, m, v, u=None):
"""dataObj2Deriv(m, v, u=None)
:param Problem,Survey forward: forward simulation
:param numpy.array m: geophysical model
@@ -100,7 +77,7 @@ class BaseDataMisfit(object):
return survey.nD
class l2_DataMisfit(object):
class l2_DataMisfit(BaseDataMisfit):
"""
The data misfit with an l_2 norm:
@@ -111,44 +88,52 @@ class l2_DataMisfit(object):
"""
def __init__(self, **kwargs):
pass
def __init__(self, survey, **kwargs):
BaseDataMisfit.__init__(self, survey, **kwargs)
def getWd(self, survey):
@property
def Wd(self):
"""getWd(survey)
Get the data weighting matrix.
The data weighting matrix.
This is based on the norm of the data plus a noise floor.
The default is based on the norm of the data plus a noise floor.
:param Survey survey: geophysical survey
:rtype: scipy.sparse.csr_matrix
:return: Wd
"""
eps = np.linalg.norm(Utils.mkvc(survey.dobs),2)*1e-5
return Utils.sdiag(1/(abs(survey.dobs)*survey.std+eps))
if getattr(self, '_Wd', None) is None:
print 'SimPEG.l2_DataMisfit is creating default weightings for Wd.'
survey = self.survey
eps = np.linalg.norm(Utils.mkvc(survey.dobs),2)*1e-5
self._Wd = Utils.sdiag(1/(abs(survey.dobs)*survey.std+eps))
return self._Wd
@Wd.setter
def Wd(self, value):
self._Wd = value
@Utils.timeIt
def dataObj(self, forward, m, u=None):
"dataObj2Deriv(forward, m, u=None)"
prob, survey = _splitForward(forward)
Wd = self.getWd(survey)
R = Wd * survey.residual(m, u=u)
def dataObj(self, m, u=None):
"dataObj2Deriv(m, u=None)"
prob = self.prob
survey = self.survey
R = self.Wd * survey.residual(m, u=u)
return 0.5*np.vdot(R, R)
@Utils.timeIt
def dataObjDeriv(self, forward, m, u=None):
"dataObj2Deriv(forward, m, u=None)"
prob, survey = _splitForward(forward)
def dataObjDeriv(self, m, u=None):
"dataObj2Deriv(m, u=None)"
prob = self.prob
survey = self.survey
if u is None: u = prob.fields(m)
Wd = self.getWd(survey)
return prob.Jtvec(m, Wd * (Wd * survey.residual(m, u=u)), u=u)
return prob.Jtvec(m, self.Wd * (self.Wd * survey.residual(m, u=u)), u=u)
@Utils.timeIt
def dataObj2Deriv(self, forward, m, v, u=None):
"dataObj2Deriv(forward, m, v, u=None)"
prob, survey = _splitForward(forward)
def dataObj2Deriv(self, m, v, u=None):
"dataObj2Deriv(m, v, u=None)"
prob = self.prob
if u is None: u = prob.fields(m)
Wd = self.getWd(survey)
return prob.Jtvec_approx(m, Wd * (Wd * prob.Jvec_approx(m, v, u=u)), u=u)
return prob.Jtvec_approx(m, self.Wd * (self.Wd * prob.Jvec_approx(m, v, u=u)), u=u)
SimPEG/Directives.py
+3 -3
View File
@@ -27,9 +27,9 @@ class InversionDirective(object):
@property
def dmisfit(self): return self.invProb.dmisfit
@property
def survey(self): return self.invProb.survey
def survey(self): return self.dmisfit.survey
@property
def prob(self): return self.invProb.prob
def prob(self): return self.dmisfit.prob
def initialize(self):
pass
@@ -126,7 +126,7 @@ class BetaEstimate_ByEig(InversionDirective):
u = self.invProb.u_current or self.prob.fields(m)
x0 = np.random.rand(*m.shape)
t = x0.dot(self.dmisfit.dataObj2Deriv(self.prob,m,x0,u=u))
t = x0.dot(self.dmisfit.dataObj2Deriv(m,x0,u=u))
b = x0.dot(self.reg.modelObj2Deriv(m, v=x0))
self.beta0 = self.beta0_ratio*(t/b)
SimPEG/Examples/Linear.py
+2 -2
View File
@@ -51,9 +51,9 @@ def run(N, plotIt=True):
M = prob.mesh
reg = Regularization.Tikhonov(mesh)
dmis = DataMisfit.l2_DataMisfit()
dmis = DataMisfit.l2_DataMisfit(survey)
opt = Optimization.InexactGaussNewton(maxIter=20)
invProb = InvProblem.BaseInvProblem(prob, reg, dmis, opt)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
beta = Directives.BetaSchedule()
betaest = Directives.BetaEstimate_ByEig()
inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest])
SimPEG/InvProblem.py
+14 -8
View File
@@ -1,9 +1,11 @@
import Utils, Survey, Problem, numpy as np, scipy.sparse as sp, gc
from Utils.SolverUtils import *
from DataMisfit import _splitForward
import DataMisfit
import Regularization
class BaseInvProblem(object):
"""BaseInvProblem(forward, reg, **kwargs)"""
"""BaseInvProblem(dmisfit, reg, opt)"""
__metaclass__ = Utils.SimPEGMetaClass
@@ -20,12 +22,16 @@ class BaseInvProblem(object):
m_current = None #: The most current model
def __init__(self, forward, reg, dmisfit, opt, **kwargs):
def __init__(self, dmisfit, reg, opt, **kwargs):
Utils.setKwargs(self, **kwargs)
self.prob, self.survey = _splitForward(forward)
self.reg = reg
assert isinstance(dmisfit, DataMisfit.BaseDataMisfit), 'dmisfit must be a DataMisfit class.'
assert isinstance(reg, Regularization.BaseRegularization), 'reg must be a Regularization class.'
self.dmisfit = dmisfit
self.reg = reg
self.opt = opt
self.prob, self.survey = dmisfit.prob, dmisfit.survey
#TODO: Remove: (and make iteration printers better!)
self.opt.parent = self
@Utils.callHooks('startup')
def startup(self, m0):
@@ -60,7 +66,7 @@ class BaseInvProblem(object):
u = self.prob.fields(m)
self.u_current = u
phi_d = self.dmisfit.dataObj(forward, m, u=u)
phi_d = self.dmisfit.dataObj(m, u=u)
phi_m = self.reg.modelObj(m)
self.dpred = self.survey.dpred(m, u=u) # This is a cheap matrix vector calculation.
@@ -72,7 +78,7 @@ class BaseInvProblem(object):
out = (f,)
if return_g:
phi_dDeriv = self.dmisfit.dataObjDeriv(forward, m, u=u)
phi_dDeriv = self.dmisfit.dataObjDeriv(m, u=u)
phi_mDeriv = self.reg.modelObjDeriv(m)
g = phi_dDeriv + self.beta * phi_mDeriv
@@ -80,7 +86,7 @@ class BaseInvProblem(object):
if return_H:
def H_fun(v):
phi_d2Deriv = self.dmisfit.dataObj2Deriv(forward, m, v, u=u)
phi_d2Deriv = self.dmisfit.dataObj2Deriv(m, v, u=u)
phi_m2Deriv = self.reg.modelObj2Deriv(m, v=v)
return phi_d2Deriv + self.beta * phi_m2Deriv
SimPEG/Inversion.py
+7 -6
View File
@@ -33,7 +33,8 @@ class BaseInversion(object):
self._directiveList = value
self._directiveList.inversion = self
def __init__(self, invProb, **kwargs):
def __init__(self, invProb, directiveList=[], **kwargs):
self.directiveList = directiveList
Utils.setKwargs(self, **kwargs)
self.invProb = invProb
@@ -43,11 +44,11 @@ class BaseInversion(object):
self.stoppers = [StoppingCriteria.iteration]
# # Check if we have inserted printers into the optimization
# if IterationPrinters.phi_d not in self.opt.printers:
# self.opt.printers.insert(1,IterationPrinters.beta)
# self.opt.printers.insert(2,IterationPrinters.phi_d)
# self.opt.printers.insert(3,IterationPrinters.phi_m)
# Check if we have inserted printers into the optimization
if IterationPrinters.phi_d not in self.opt.printers:
self.opt.printers.insert(1,IterationPrinters.beta)
self.opt.printers.insert(2,IterationPrinters.phi_d)
self.opt.printers.insert(3,IterationPrinters.phi_m)
@Utils.timeIt
def run(self, m0):
SimPEG/Maps.py
-23
View File
@@ -312,29 +312,6 @@ class Mesh2Mesh(IdentityMap):
"""
Takes a model on one mesh are translates it to another mesh.
.. plot::
from SimPEG import *
M = Mesh.TensorMesh([100,100])
h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)])
h1 = h1/h1.sum()
M2 = Mesh.TensorMesh([h1,h1])
V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50)
v = Utils.mkvc(V)
modh = Maps.Mesh2Mesh([M,M2])
modH = Maps.Mesh2Mesh([M2,M])
H = modH.transform(v)
h = modh.transform(H)
ax = plt.subplot(131)
M.plotImage(v, ax=ax)
ax.set_title('Fine Mesh (Original)')
ax = plt.subplot(132)
M2.plotImage(H,clim=[0,1],ax=ax)
ax.set_title('Course Mesh')
ax = plt.subplot(133)
M.plotImage(h,clim=[0,1],ax=ax)
ax.set_title('Fine Mesh (Interpolated)')
"""
def __init__(self, meshes, **kwargs):
SimPEG/Optimization.py
+5 -4
View File
@@ -5,6 +5,7 @@ norm = np.linalg.norm
__all__ = ['Minimize', 'Remember', 'SteepestDescent', 'BFGS', 'GaussNewton', 'InexactGaussNewton', 'ProjectedGradient', 'NewtonRoot', 'StoppingCriteria', 'IterationPrinters']
SolverICG = SolverWrapI(sp.linalg.cg, checkAccuracy=False)
class StoppingCriteria(object):
"""docstring for StoppingCriteria"""
@@ -72,9 +73,9 @@ class IterationPrinters(object):
bSet = {"title": "bSet", "value": lambda M: np.sum(M.bindingSet(M.xc)), "width": 8, "format": "%d"}
comment = {"title": "Comment", "value": lambda M: M.comment, "width": 12, "format": "%s"}
beta = {"title": "beta", "value": lambda M: M.parent.objFunc.beta, "width": 10, "format": "%1.2e"}
phi_d = {"title": "phi_d", "value": lambda M: M.parent.objFunc.phi_d, "width": 10, "format": "%1.2e"}
phi_m = {"title": "phi_m", "value": lambda M: M.parent.objFunc.phi_m, "width": 10, "format": "%1.2e"}
beta = {"title": "beta", "value": lambda M: M.parent.beta, "width": 10, "format": "%1.2e"}
phi_d = {"title": "phi_d", "value": lambda M: M.parent.phi_d, "width": 10, "format": "%1.2e"}
phi_m = {"title": "phi_m", "value": lambda M: M.parent.phi_m, "width": 10, "format": "%1.2e"}
class Minimize(object):
@@ -772,7 +773,7 @@ class InexactGaussNewton(BFGS, Minimize, Remember):
@Utils.timeIt
def findSearchDirection(self):
Hinv = SolverCG(self.H, M=self.approxHinv, tol=self.tolCG, maxiter=self.maxIterCG)
Hinv = SolverICG(self.H, M=self.approxHinv, tol=self.tolCG, maxiter=self.maxIterCG)
p = Hinv * (-self.g)
return p
docs/api_Maps.rst
+27
View File
@@ -157,6 +157,33 @@ Vertical 1D Map
Mesh to Mesh Map
----------------
.. plot::
from SimPEG import *
import matplotlib.pyplot as plt
M = Mesh.TensorMesh([100,100])
h1 = Utils.meshTensor([(6,7,-1.5),(6,10),(6,7,1.5)])
h1 = h1/h1.sum()
M2 = Mesh.TensorMesh([h1,h1])
V = Utils.ModelBuilder.randomModel(M.vnC, seed=79, its=50)
v = Utils.mkvc(V)
modh = Maps.Mesh2Mesh([M,M2])
modH = Maps.Mesh2Mesh([M2,M])
H = modH.transform(v)
h = modh.transform(H)
ax = plt.subplot(131)
M.plotImage(v, ax=ax)
ax.set_title('Fine Mesh (Original)')
ax = plt.subplot(132)
M2.plotImage(H,clim=[0,1],ax=ax)
ax.set_title('Course Mesh')
ax = plt.subplot(133)
M.plotImage(h,clim=[0,1],ax=ax)
ax.set_title('Fine Mesh (Interpolated)')
plt.show()
.. autoclass:: SimPEG.Maps.Mesh2Mesh
:members:
:undoc-members: