From 783cd405387b15ecb95f46462e2e7bbb41c5fadb Mon Sep 17 00:00:00 2001 From: rowanc1 Date: Fri, 6 Dec 2013 16:54:38 -0800 Subject: [PATCH] Bug fixes in Save, incorporation of data into the inversion object, taken out of problem. --- SimPEG/examples/Linear.py | 30 +++++++++----------------- SimPEG/forward/Problem.py | 43 +++++++++++-------------------------- SimPEG/inverse/Inversion.py | 26 +++++++++++----------- SimPEG/utils/Save.py | 28 ++++++++++++++---------- 4 files changed, 53 insertions(+), 74 deletions(-) diff --git a/SimPEG/examples/Linear.py b/SimPEG/examples/Linear.py index b593edf8..79bec85c 100644 --- a/SimPEG/examples/Linear.py +++ b/SimPEG/examples/Linear.py @@ -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() diff --git a/SimPEG/forward/Problem.py b/SimPEG/forward/Problem.py index 1f3d817c..1412ff5b 100644 --- a/SimPEG/forward/Problem.py +++ b/SimPEG/forward/Problem.py @@ -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) diff --git a/SimPEG/inverse/Inversion.py b/SimPEG/inverse/Inversion.py index 19280e0a..926f719a 100644 --- a/SimPEG/inverse/Inversion.py +++ b/SimPEG/inverse/Inversion.py @@ -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) diff --git a/SimPEG/utils/Save.py b/SimPEG/utils/Save.py index d9b3dd58..70ce1d2e 100644 --- a/SimPEG/utils/Save.py +++ b/SimPEG/utils/Save.py @@ -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)