From 0a0caceaca45da0a542e8f4cc4e9750cebd03409 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 29 Mar 2016 22:49:03 -0700 Subject: [PATCH] Problem.Jvec, Problem.Jtvec, Problem.fields, DataMisfit, survey.dpred take a fields object f (not a solution vector, u) --- SimPEG/DCIP/BaseDC.py | 24 +++++----- SimPEG/DCIP/BaseIP.py | 8 ++-- SimPEG/DataMisfit.py | 46 +++++++++----------- SimPEG/FLOW/Richards/RichardsProblem.py | 58 ++++++++++++------------- SimPEG/Problem.py | 28 ++++++------ SimPEG/Survey.py | 26 +++++------ tests/flow/test_Richards.py | 12 ++--- 7 files changed, 98 insertions(+), 104 deletions(-) diff --git a/SimPEG/DCIP/BaseDC.py b/SimPEG/DCIP/BaseDC.py index e3d353d7..475c621e 100644 --- a/SimPEG/DCIP/BaseDC.py +++ b/SimPEG/DCIP/BaseDC.py @@ -200,11 +200,11 @@ class ProblemDC_CC(Problem.BaseProblem): return F - def Jvec(self, m, v, u=None): + def Jvec(self, m, v, f=None): """ :param numpy.array m: model :param numpy.array v: vector to multiply - :param numpy.array u: fields + :param Fields f: fields :rtype: numpy.array :return: Jv @@ -225,11 +225,10 @@ class ProblemDC_CC(Problem.BaseProblem): # Set current model; clear dependent property $\mathbf{A(m)}$ self.curModel = m sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$ - if u is None: + if f is None: # Run forward simulation if $u$ not provided - u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol'] - else: - u = u[self.survey.srcList, 'phi_sol'] + f = self.fields(self.curModel) + u = f[self.survey.srcList, 'phi_sol'] D = self.mesh.faceDiv G = self.mesh.cellGrad @@ -251,19 +250,18 @@ class ProblemDC_CC(Problem.BaseProblem): if self.Ainv is None: self.Ainv = self.Solver(dA_du, **self.solverOpts) - P = self.survey.getP(self.mesh) + P = self.survey.getP(self.mesh) Jv = - P * mkvc( self.Ainv * dCdm_x_v ) return Jv - def Jtvec(self, m, v, u=None): + def Jtvec(self, m, v, f=None): self.curModel = m sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$ - if u is None: - # Run forward simulation if $u$ not provided - u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol'] - else: - u = u[self.survey.srcList, 'phi_sol'] + if f is None: + # Run forward simulation if $f$ not provided + f = self.fields(self.curModel) + u = f[self.survey.srcList, 'phi_sol'] shp = u.shape P = self.survey.getP(self.mesh) diff --git a/SimPEG/DCIP/BaseIP.py b/SimPEG/DCIP/BaseIP.py index cec0ea2e..a18b5a47 100644 --- a/SimPEG/DCIP/BaseIP.py +++ b/SimPEG/DCIP/BaseIP.py @@ -14,12 +14,12 @@ class SurveyIP(SurveyDC): Survey.BaseSurvey.__init__(self, **kwargs) self._Ps = {} - def dpred(self, m, u=None): + def dpred(self, m, f=None): """ Predicted data. .. math:: - d_\\text{pred} = Pu(m) + d_\\text{pred} = Pf(m) """ return self.prob.forward(m) @@ -143,10 +143,10 @@ class ProblemIP(Problem.BaseProblem): J_x_v = - P * mkvc( self.Ainv * dCdm_x_v ) return -J_x_v - def Jvec(self, m, v, u=None): + def Jvec(self, m, v, f=None): return self.forward(v) - def Jtvec(self, m, v, u=None): + def Jtvec(self, m, v, f=None): self.curModel = m # sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$ diff --git a/SimPEG/DataMisfit.py b/SimPEG/DataMisfit.py index 13574771..1fbf7b2d 100644 --- a/SimPEG/DataMisfit.py +++ b/SimPEG/DataMisfit.py @@ -22,11 +22,11 @@ class BaseDataMisfit(object): Utils.setKwargs(self,**kwargs) @Utils.timeIt - def eval(self, m, u=None): - """eval(m, u=None) + def eval(self, m, f=None): + """eval(m, f=None) :param numpy.array m: geophysical model - :param numpy.array u: fields + :param Fields f: fields :rtype: float :return: data misfit @@ -34,11 +34,11 @@ class BaseDataMisfit(object): raise NotImplementedError('This method should be overwritten.') @Utils.timeIt - def evalDeriv(self, m, u=None): - """evalDeriv(m, u=None) + def evalDeriv(self, m, f=None): + """evalDeriv(m, f=None) :param numpy.array m: geophysical model - :param numpy.array u: fields + :param Fields f: fields :rtype: numpy.array :return: data misfit derivative @@ -47,12 +47,12 @@ class BaseDataMisfit(object): @Utils.timeIt - def eval2Deriv(self, m, v, u=None): - """eval2Deriv(m, v, u=None) + def eval2Deriv(self, m, v, f=None): + """eval2Deriv(m, v, f=None) :param numpy.array m: geophysical model :param numpy.array v: vector to multiply - :param numpy.array u: fields + :param Fields f: fields :rtype: numpy.array :return: data misfit derivative @@ -108,24 +108,20 @@ class l2_DataMisfit(BaseDataMisfit): self._Wd = value @Utils.timeIt - def eval(self, m, u=None): - "eval(m, u=None)" - prob = self.prob - survey = self.survey - R = self.Wd * survey.residual(m, u) + def eval(self, m, f=None): + "eval(m, f=None)" + if f is None: f = self.prob.fields(m) + R = self.Wd * self.survey.residual(m, f) return 0.5*np.vdot(R, R) @Utils.timeIt - def evalDeriv(self, m, u=None): - "evalDeriv(m, u=None)" - prob = self.prob - survey = self.survey - if u is None: u = prob.fields(m) - return prob.Jtvec(m, self.Wd * (self.Wd * survey.residual(m, u)), u) + def evalDeriv(self, m, f=None): + "evalDeriv(m, f=None)" + if f is None: f = self.prob.fields(m) + return self.prob.Jtvec(m, self.Wd * (self.Wd * self.survey.residual(m, f=f)), f=f) @Utils.timeIt - def eval2Deriv(self, m, v, u=None): - "eval2Deriv(m, v, u=None)" - prob = self.prob - if u is None: u = prob.fields(m) - return prob.Jtvec_approx(m, self.Wd * (self.Wd * prob.Jvec_approx(m, v, u=u)), u) + def eval2Deriv(self, m, v, f=None): + "eval2Deriv(m, v, f=None)" + if f is None: f = prob.fields(m) + return self.prob.Jtvec_approx(m, self.Wd * (self.Wd * prob.Jvec_approx(m, v, f=f)), f=f) diff --git a/SimPEG/FLOW/Richards/RichardsProblem.py b/SimPEG/FLOW/Richards/RichardsProblem.py index 4dcabe60..2346f4da 100644 --- a/SimPEG/FLOW/Richards/RichardsProblem.py +++ b/SimPEG/FLOW/Richards/RichardsProblem.py @@ -45,19 +45,19 @@ class RichardsSurvey(Survey.BaseSurvey): @Utils.count @Utils.requires('prob') - def dpred(self, m, u=None): + def dpred(self, m, f=None): """ Create the projected data from a model. - The field, u, (if provided) will be used for the predicted data + The field, f, (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), m) + d_\\text{pred} = P(f(m), m) Where P is a projection of the fields onto the data space. """ - if u is None: u = self.prob.fields(m) - return Utils.mkvc(self.eval(u, m)) + if f is None: f = self.prob.fields(m) + return Utils.mkvc(self.eval(f, m)) @Utils.requires('prob') def eval(self, U, m): @@ -233,16 +233,16 @@ class RichardsProblem(Problem.BaseTimeProblem): return r, J @Utils.timeIt - def Jfull(self, m, u=None): - if u is None: - u = self.fields(m) + def Jfull(self, m, f=None): + if f is None: + f = self.fields(m) - nn = len(u)-1 + nn = len(f)-1 Asubs, Adiags, Bs = range(nn), range(nn), range(nn) for ii in range(nn): dt = self.timeSteps[ii] - bc = self.getBoundaryConditions(ii, u[ii]) - Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, u[ii], u[ii+1], dt, bc) + bc = self.getBoundaryConditions(ii, f[ii]) + Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, f[ii], f[ii+1], dt, bc) Ad = sp.block_diag(Adiags) zRight = Utils.spzeros((len(Asubs)-1)*Asubs[0].shape[0],Adiags[0].shape[1]) zTop = Utils.spzeros(Adiags[0].shape[0], len(Adiags)*Adiags[0].shape[1]) @@ -251,7 +251,7 @@ class RichardsProblem(Problem.BaseTimeProblem): B = np.array(sp.vstack(Bs).todense()) Ainv = self.Solver(A, **self.solverOpts) - P = self.survey.evalDeriv(u, m) + P = self.survey.evalDeriv(f, m) AinvB = Ainv * B z = np.zeros((self.mesh.nC, B.shape[1])) zAinvB = np.vstack((z, AinvB)) @@ -259,41 +259,41 @@ class RichardsProblem(Problem.BaseTimeProblem): return J @Utils.timeIt - def Jvec(self, m, v, u=None): - if u is None: - u = self.fields(m) + def Jvec(self, m, v, f=None): + if f is None: + f = self.fields(m) - JvC = range(len(u)-1) # Cell to hold each row of the long vector. + JvC = range(len(f)-1) # Cell to hold each row of the long vector. # This is done via forward substitution. - bc = self.getBoundaryConditions(0, u[0]) - temp, Adiag, B = self.diagsJacobian(m, u[0], u[1], self.timeSteps[0], bc) + bc = self.getBoundaryConditions(0, f[0]) + temp, Adiag, B = self.diagsJacobian(m, f[0], f[1], self.timeSteps[0], bc) Adiaginv = self.Solver(Adiag, **self.solverOpts) JvC[0] = Adiaginv * (B*v) - for ii in range(1,len(u)-1): - bc = self.getBoundaryConditions(ii, u[ii]) - Asub, Adiag, B = self.diagsJacobian(m, u[ii], u[ii+1], self.timeSteps[ii], bc) + for ii in range(1,len(f)-1): + bc = self.getBoundaryConditions(ii, f[ii]) + Asub, Adiag, B = self.diagsJacobian(m, f[ii], f[ii+1], self.timeSteps[ii], bc) Adiaginv = self.Solver(Adiag, **self.solverOpts) JvC[ii] = Adiaginv * (B*v - Asub*JvC[ii-1]) - P = self.survey.evalDeriv(u, m) + P = self.survey.evalDeriv(f, m) return P * np.concatenate([np.zeros(self.mesh.nC)] + JvC) @Utils.timeIt - def Jtvec(self, m, v, u=None): - if u is None: - u = self.field(m) + def Jtvec(self, m, v, f=None): + if f is None: + f = self.field(m) - P = self.survey.evalDeriv(u, m) + P = self.survey.evalDeriv(f, m) PTv = P.T*v # This is done via backward substitution. minus = 0 BJtv = 0 - for ii in range(len(u)-1,0,-1): - bc = self.getBoundaryConditions(ii-1, u[ii-1]) - Asub, Adiag, B = self.diagsJacobian(m, u[ii-1], u[ii], self.timeSteps[ii-1], bc) + for ii in range(len(f)-1,0,-1): + bc = self.getBoundaryConditions(ii-1, f[ii-1]) + Asub, Adiag, B = self.diagsJacobian(m, f[ii-1], f[ii], self.timeSteps[ii-1], bc) #select the correct part of v vpart = range((ii)*Adiag.shape[0], (ii+1)*Adiag.shape[0]) AdiaginvT = self.Solver(Adiag.T, **self.solverOpts) diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index cd8a4aaa..8ed94e97 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -88,28 +88,28 @@ class BaseProblem(object): return self.survey is not None @Utils.timeIt - def Jvec(self, m, v, u=None): - """Jvec(m, v, u=None) + def Jvec(self, m, v, f=None): + """Jvec(m, v, f=None) Effect of J(m) on a vector v. :param numpy.array m: model :param numpy.array v: vector to multiply - :param numpy.array u: fields + :param Fields f: fields :rtype: numpy.array :return: Jv """ raise NotImplementedError('J is not yet implemented.') @Utils.timeIt - def Jtvec(self, m, v, u=None): - """Jtvec(m, v, u=None) + def Jtvec(self, m, v, f=None): + """Jtvec(m, v, f=None) Effect of transpose of J(m) on a vector v. :param numpy.array m: model :param numpy.array v: vector to multiply - :param numpy.array u: fields + :param Fields f: fields :rtype: numpy.array :return: JTv """ @@ -117,28 +117,28 @@ class BaseProblem(object): @Utils.timeIt - def Jvec_approx(self, m, v, u=None): - """Jvec_approx(m, v, u=None) + def Jvec_approx(self, m, v, f=None): + """Jvec_approx(m, v, f=None) Approximate effect of J(m) on a vector v :param numpy.array m: model :param numpy.array v: vector to multiply - :param numpy.array u: fields + :param Fields f: fields :rtype: numpy.array :return: approxJv """ return self.Jvec(m, v, u) @Utils.timeIt - def Jtvec_approx(self, m, v, u=None): - """Jtvec_approx(m, v, u=None) + def Jtvec_approx(self, m, v, f=None): + """Jtvec_approx(m, v, f=None) Approximate effect of transpose of J(m) on a vector v. :param numpy.array m: model :param numpy.array v: vector to multiply - :param numpy.array u: fields + :param Fields f: fields :rtype: numpy.array :return: JTv """ @@ -224,9 +224,9 @@ class LinearProblem(BaseProblem): def fields(self, m): return self.G.dot(m) - def Jvec(self, m, v, u=None): + def Jvec(self, m, v, f=None): return self.G.dot(v) - def Jtvec(self, m, v, u=None): + def Jtvec(self, m, v, f=None): return self.G.T.dot(v) diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 37024028..9ddd269e 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -295,21 +295,21 @@ class BaseSurvey(object): @Utils.count @Utils.requires('prob') - def dpred(self, m, u=None): - """dpred(m, u=None) + def dpred(self, m, f=None): + """dpred(m, f=None) Create the projected data from a model. - The field, u, (if provided) will be used for the predicted data + The fields, f, (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)) + d_\\text{pred} = P(f(m)) Where P is a projection of the fields onto the data space. """ - if u is None: u = self.prob.fields(m) - return Utils.mkvc(self.eval(u)) + if f is None: f = self.prob.fields(m) + return Utils.mkvc(self.eval(f)) @Utils.count @@ -337,11 +337,11 @@ class BaseSurvey(object): raise NotImplemented('eval is not yet implemented.') @Utils.count - def residual(self, m, u=None): - """residual(m, u=None) + def residual(self, m, f=None): + """residual(m, f=None) :param numpy.array m: geophysical model - :param numpy.array u: fields + :param numpy.array f: fields :rtype: numpy.array :return: data residual @@ -352,14 +352,14 @@ class BaseSurvey(object): \mu_\\text{data} = \mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs} """ - return Utils.mkvc(self.dpred(m, u=u) - self.dobs) + return Utils.mkvc(self.dpred(m, f=f) - self.dobs) @property def isSynthetic(self): "Check if the data is synthetic." return self.mtrue is not None - def makeSyntheticData(self, m, std=0.05, u=None, force=False): + def makeSyntheticData(self, m, std=0.05, f=None, force=False): """ Make synthetic data given a model, and a standard deviation. @@ -372,7 +372,7 @@ class BaseSurvey(object): if getattr(self, 'dobs', None) is not None and not force: raise Exception('Survey already has dobs. You can use force=True to override this exception.') self.mtrue = m - self.dtrue = self.dpred(m, u=u) + self.dtrue = self.dpred(m, f=f) noise = std*abs(self.dtrue)*np.random.randn(*self.dtrue.shape) self.dobs = self.dtrue+noise self.std = self.dobs*0 + std @@ -381,7 +381,7 @@ class BaseSurvey(object): class LinearSurvey(BaseSurvey): def eval(self, u): return u - + @property def nD(self): return self.prob.G.shape[0] diff --git a/tests/flow/test_Richards.py b/tests/flow/test_Richards.py index d63a6210..f67ec71d 100644 --- a/tests/flow/test_Richards.py +++ b/tests/flow/test_Richards.py @@ -116,8 +116,8 @@ class RichardsTests1D(unittest.TestCase): v = np.random.rand(self.survey.nD) z = np.random.rand(self.M.nC) Hs = self.prob.fields(self.Ks) - vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs)) - zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs)) + vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs)) + zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs)) tol = TOL*(10**int(np.log10(np.abs(zJv)))) passed = np.abs(vJz - zJv) < tol print 'Richards Adjoint Test - PressureHead' @@ -188,8 +188,8 @@ class RichardsTests2D(unittest.TestCase): v = np.random.rand(self.survey.nD) z = np.random.rand(self.M.nC) Hs = self.prob.fields(self.Ks) - vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs)) - zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs)) + vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs)) + zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs)) tol = TOL*(10**int(np.log10(np.abs(zJv)))) passed = np.abs(vJz - zJv) < tol print '2D: Richards Adjoint Test - PressureHead' @@ -260,8 +260,8 @@ class RichardsTests3D(unittest.TestCase): v = np.random.rand(self.survey.nD) z = np.random.rand(self.M.nC) Hs = self.prob.fields(self.Ks) - vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs)) - zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs)) + vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs)) + zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs)) tol = TOL*(10**int(np.log10(np.abs(zJv)))) passed = np.abs(vJz - zJv) < tol print '3D: Richards Adjoint Test - PressureHead'