wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-28 14:12:43 +08:00
Code Issues Packages Projects Releases Wiki Activity

Bug fixes in Save, incorporation of data into the inversion object, taken out of problem.

Browse Source
This commit is contained in:
rowanc1
2013-12-06 16:54:38 -08:00
parent 09938ed7df
commit 783cd40538
4 changed files with 53 additions and 74 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/examples/Linear.py
+10 -20
View File
@@ -34,35 +34,27 @@ def example(N):
for i in range(nk):
G[i,:] = g(i)
m_true = np.zeros(M.nC)
m_true[M.vectorCCx > 0.3] = 1.
m_true[M.vectorCCx > 0.45] = -0.5
m_true[M.vectorCCx > 0.6] = 0
d_true = G.dot(m_true)
noise = 0.1 * np.random.rand(d_true.size)
d_obs = d_true + noise
mtrue = np.zeros(M.nC)
mtrue[M.vectorCCx > 0.3] = 1.
mtrue[M.vectorCCx > 0.45] = -0.5
mtrue[M.vectorCCx > 0.6] = 0
prob = LinearProblem(M)
prob.G = G
prob.dobs = d_obs
prob.std = np.ones_like(d_obs)*0.1
data = prob.createSyntheticData(mtrue, std=0.01)
return prob, m_true
return prob, data
if __name__ == '__main__':
prob, m_true = example(100)
prob, data = example(100)
M = prob.mesh
reg = regularization.Regularization(M)
opt = inverse.InexactGaussNewton(maxIter=20)
inv = inverse.Inversion(prob,reg,opt,beta0=1e-4)
m0 = np.zeros_like(m_true)
inv = inverse.Inversion(prob,reg,opt,data)
m0 = np.zeros_like(data.mtrue)
mrec = inv.run(m0)
@@ -72,9 +64,7 @@ if __name__ == '__main__':
plt.figure(2)
plt.plot(M.vectorCCx, m_true, 'b-')
plt.plot(M.vectorCCx, data.mtrue, 'b-')
plt.plot(M.vectorCCx, mrec, 'r-')
plt.show()
SimPEG/forward/Problem.py
+12 -31
View File
@@ -1,4 +1,4 @@
from SimPEG import utils, np, sp
from SimPEG import utils, data, np, sp
norm = np.linalg.norm
@@ -37,9 +37,11 @@ class Problem(object):
__metaclass__ = utils.Save.Savable
counter = None
counter = None #: A SimPEG.utils.Counter object
def __init__(self, mesh):
def __init__(self, mesh, *args, **kwargs):
utils.setKwargs(self, **kwargs)
self.mesh = mesh
@property
@@ -65,26 +67,6 @@ class Problem(object):
def P(self, value):
self._P = value
@property
def std(self):
"""
Estimated Standard Deviations.
"""
return self._std
@std.setter
def std(self, value):
self._std = value
@property
def dobs(self):
"""
Observed data.
"""
return self._dobs
@dobs.setter
def dobs(self, value):
self._dobs = value
@utils.count
def dpred(self, m, u=None):
"""
@@ -98,7 +80,7 @@ class Problem(object):
return self.P*u
@utils.count
def dataResidual(self, m, u=None):
def dataResidual(self, m, data, u=None):
"""
:param numpy.array m: geophysical model
:param numpy.array u: fields
@@ -115,7 +97,7 @@ class Problem(object):
u is the field of interest; d_obs is the observed data.
"""
return self.dpred(m, u=u) - self.dobs
return self.dpred(m, u=u) - data.dobs
@utils.timeIt
def J(self, m, v, u=None):
@@ -238,12 +220,11 @@ class Problem(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.
"""
dobs = self.dpred(m,u=u)
noise = std*abs(dobs)*np.random.randn(*dobs.shape)
dobs = dobs+noise
eps = np.linalg.norm(utils.mkvc(dobs),2)*1e-5
Wd = 1/(abs(dobs)*std+eps)
return dobs, Wd
dtrue = self.dpred(m,u=u)
noise = std*abs(dtrue)*np.random.randn(*dtrue.shape)
dobs = dtrue+noise
stdev = dobs*0 + std
return data.SimPEGData(self, dobs=dobs, std=stdev, dtrue=dtrue, mtrue=m)
SimPEG/inverse/Inversion.py
+14 -12
View File
@@ -5,7 +5,8 @@ from BetaSchedule import Cooling
from SimPEG.inverse import IterationPrinters, StoppingCriteria
class BaseInversion(object):
"""docstring for BaseInversion"""
"""BaseInversion(prob, reg, opt, data, **kwargs)
"""
__metaclass__ = utils.Save.Savable
@@ -20,11 +21,12 @@ class BaseInversion(object):
beta0 = None #: The initial Beta (regularization parameter)
beta0_ratio = 0.1 #: When beta0 is set to None, estimateBeta0 is used with this ratio
def __init__(self, prob, reg, opt, **kwargs):
def __init__(self, prob, reg, opt, data, **kwargs):
utils.setKwargs(self, **kwargs)
self.prob = prob
self.reg = reg
self.opt = opt
self.data = data
self.opt.parent = self
self.stoppers = [StoppingCriteria.iteration]
@@ -46,8 +48,8 @@ class BaseInversion(object):
Standard deviation weighting matrix.
"""
if getattr(self,'_Wd',None) is None:
eps = np.linalg.norm(utils.mkvc(self.prob.dobs),2)*1e-5
self._Wd = 1/(abs(self.prob.dobs)*self.prob.std+eps)
eps = np.linalg.norm(utils.mkvc(self.data.dobs),2)*1e-5
self._Wd = 1/(abs(self.data.dobs)*self.data.std+eps)
return self._Wd
@Wd.setter
def Wd(self, value):
@@ -63,7 +65,7 @@ class BaseInversion(object):
Note that we do not set the target if it is None, but we return the default value.
"""
if getattr(self, '_phi_d_target', None) is None:
return self.prob.dobs.size #
return self.data.dobs.size #
return self._phi_d_target
@phi_d_target.setter
@@ -256,7 +258,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, u=u)
R = self.Wd*self.prob.dataResidual(m, self.data, u=u)
R = utils.mkvc(R)
return 0.5*np.vdot(R, R)
@@ -296,7 +298,7 @@ class BaseInversion(object):
if u is None:
u = self.prob.field(m)
R = self.Wd*self.prob.dataResidual(m, u=u)
R = self.Wd*self.prob.dataResidual(m, self.data, u=u)
dmisfit = self.prob.Jt(m, self.Wd * R, u=u)
@@ -341,7 +343,7 @@ class BaseInversion(object):
if u is None:
u = self.prob.field(m)
R = self.Wd*self.prob.dataResidual(m, u=u)
R = self.Wd*self.prob.dataResidual(m, self.data, u=u)
# TODO: abstract to different norms a little cleaner.
# \/ it goes here. in l2 it is the identity.
@@ -360,8 +362,8 @@ class Inversion(Cooling, Remember, BaseInversion):
maxIter = 10
name = "SimPEG Inversion"
def __init__(self, prob, reg, opt, **kwargs):
BaseInversion.__init__(self, prob, reg, opt, **kwargs)
def __init__(self, prob, reg, opt, data, **kwargs):
BaseInversion.__init__(self, prob, reg, opt, data, **kwargs)
self.stoppers.append(StoppingCriteria.phi_d_target_Inversion)
@@ -377,8 +379,8 @@ class TimeSteppingInversion(Remember, BaseInversion):
maxIter = 1
name = "Time-Stepping SimPEG Inversion"
def __init__(self, prob, reg, opt, **kwargs):
BaseInversion.__init__(self, prob, reg, opt, **kwargs)
def __init__(self, prob, reg, opt, data, **kwargs):
BaseInversion.__init__(self, prob, reg, opt, data, **kwargs)
self.stoppers.append(StoppingCriteria.phi_d_target_Inversion)
SimPEG/utils/Save.py
+17 -11
View File
@@ -50,14 +50,15 @@ class SimPEGTable:
node = self.inversions.addGroup('%d'%self.inversions.numChildren)
saveSavable(invObj,node.addGroup('rebuild'))
results = node.addGroup('results')
preIteration(results)
invObj._invNode = results
self.f.flush()
invObj.hook(_startup_hdf5_inv, overwrite=True)
# At the start of every iteration we will create a inversion iteration node.
def _doStartIteration_hdf5_inv(invObj):
invNodeIt = invObj._invNode.addGroup('%d'%(invObj._iter+1))
preIteration(invNodeIt)
invObj._invNodeIt = invNodeIt
invObj._invNodeIt = invObj._invNode.addGroup('%d'%(invObj._iter+1))
preIteration(invObj._invNodeIt)
invObj.hook(_doStartIteration_hdf5_inv, overwrite=True)
def _doEndIteration_hdf5_inv(invObj):
@@ -68,6 +69,7 @@ class SimPEGTable:
# Delete all iterates that did not finish.
def _finish_hdf5_inv(invObj):
postIteration(invObj._invNode)
for it in invObj._invNode:
if not it.attrs['complete']:
del self.f[it.path]
@@ -76,10 +78,8 @@ class SimPEGTable:
invObj.hook(_finish_hdf5_inv, overwrite=True)
def _doStartIteration_hdf5_opt(optObj):
optNodeIt = optObj.parent._invNode.addGroup('%d.%d'%(optObj.parent._iter, optObj._iter))
preIteration(optNodeIt)
optObj._optNodeIt = optNodeIt
optObj._optNodeIt = optObj.parent._invNode.addGroup('%d.%d'%(optObj.parent._iter, optObj._iter))
preIteration(optObj._optNodeIt)
invObj.opt.hook(_doStartIteration_hdf5_opt, overwrite=True)
def _doEndIteration_hdf5_opt(optObj, xt):
@@ -332,10 +332,16 @@ def loadSavable(node, pointers=None):
cls = get(node, '__class__')
if cls in SAVEABLES:
out = SAVEABLES[cls](*ARGS, **KWARGS)
out._savable = node
pointers.append(out) # Because this is recursive.
return out
try:
out = SAVEABLES[cls](*ARGS, **KWARGS)
out._savable = node
pointers.append(out) # Because this is recursive.
return out
except Exception, e:
print 'Warning: %s Class could not be initiated.' % cls
print 'ARGS: ', ARGS
print 'KWARGS: ', KWARGS
return (cls, ARGS, KWARGS, node)
else:
print 'Warning: %s Class not found in SimPEG.utils.Save.SAVABLES' % cls
return (cls, ARGS, KWARGS, node)