From d8d8915f942ef1f5430ba57dd726a00cb747c4ba Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 31 Mar 2016 09:28:48 -0700 Subject: [PATCH] f for fields in data misfit, directives etc. Previously, f was used in the InvProblem to be the function value for the objective function --> this has been renamed to phi --- SimPEG/DataMisfit.py | 4 ++-- SimPEG/Directives.py | 4 ++-- SimPEG/EM/FDEM/SurveyFDEM.py | 16 ++++++++-------- SimPEG/EM/TDEM/BaseTDEM.py | 22 +++++++++++----------- SimPEG/InvProblem.py | 26 +++++++++++++------------- SimPEG/MT/BaseMT.py | 26 +++++++++++++------------- SimPEG/MT/SurveyMT.py | 6 +++--- SimPEG/Problem.py | 4 ++-- SimPEG/Survey.py | 14 +++++++------- 9 files changed, 61 insertions(+), 61 deletions(-) diff --git a/SimPEG/DataMisfit.py b/SimPEG/DataMisfit.py index 1fbf7b2d..53728c4e 100644 --- a/SimPEG/DataMisfit.py +++ b/SimPEG/DataMisfit.py @@ -123,5 +123,5 @@ class l2_DataMisfit(BaseDataMisfit): @Utils.timeIt 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) + if f is None: f = self.prob.fields(m) + return self.prob.Jtvec_approx(m, self.Wd * (self.Wd * self.prob.Jvec_approx(m, v, f=f)), f=f) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 2ed27c20..3d726712 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -123,10 +123,10 @@ class BetaEstimate_ByEig(InversionDirective): if self.debug: print 'Calculating the beta0 parameter.' m = self.invProb.curModel - u = self.invProb.getFields(m, store=True, deleteWarmstart=False) + f = self.invProb.getFields(m, store=True, deleteWarmstart=False) x0 = np.random.rand(*m.shape) - t = x0.dot(self.dmisfit.eval2Deriv(m,x0,u=u)) + t = x0.dot(self.dmisfit.eval2Deriv(m,x0,f=f)) b = x0.dot(self.reg.eval2Deriv(m, v=x0)) self.beta0 = self.beta0_ratio*(t/b) diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index 163f6914..1552a12c 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -67,7 +67,7 @@ class Rx(SimPEG.Survey.BaseRx): """Grid Location projection (e.g. Ex Fy ...)""" return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1] - def eval(self, src, mesh, u): + def eval(self, src, mesh, f): """ Project fields to recievers to get data. @@ -80,27 +80,27 @@ class Rx(SimPEG.Survey.BaseRx): # projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0]) # projGLoc += self.knownRxTypes[self.rxType][1] - P = self.getP(mesh, self.projGLoc(u)) - u_part_complex = u[src, self.projField] + P = self.getP(mesh, self.projGLoc(f)) + f_part_complex = f[src, self.projField] # get the real or imag component real_or_imag = self.projComp - u_part = getattr(u_part_complex, real_or_imag) + f_part = getattr(f_part_complex, real_or_imag) - return P*u_part + return P*f_part - def evalDeriv(self, src, mesh, u, v, adjoint=False): + def evalDeriv(self, src, mesh, f, v, adjoint=False): """ Derivative of projected fields with respect to the inversion model times a vector. :param Source src: FDEM source :param Mesh mesh: mesh used - :param Fields u: fields object + :param Fields f: fields object :param numpy.ndarray v: vector to multiply :rtype: numpy.ndarray :return: fields projected to recievers """ - P = self.getP(mesh, self.projGLoc(u)) + P = self.getP(mesh, self.projGLoc(f)) if not adjoint: Pv_complex = P * v diff --git a/SimPEG/EM/TDEM/BaseTDEM.py b/SimPEG/EM/TDEM/BaseTDEM.py index 0da22072..15fc19e3 100644 --- a/SimPEG/EM/TDEM/BaseTDEM.py +++ b/SimPEG/EM/TDEM/BaseTDEM.py @@ -108,11 +108,11 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): Ainv.clean() return F - def Jvec(self, m, v, u=None): + def Jvec(self, m, v, f=None): """ :param numpy.array m: Conductivity model :param numpy.ndarray v: vector (model object) - :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m + :param simpegEM.TDEM.FieldsTDEM f: Fields resulting from m :rtype: numpy.ndarray :return: w (data object) @@ -125,15 +125,15 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): """ if self.verbose: print '%s\nCalculating J(v)\n%s'%('*'*50,'*'*50) self.curModel = m - if u is None: - u = self.fields(m) - p = self.Gvec(m, v, u) + if f is None: + f = self.fields(m) + p = self.Gvec(m, v, f) y = self.solveAh(m, p) - Jv = self.survey.evalDeriv(u, v=y) + Jv = self.survey.evalDeriv(f, v=y) if self.verbose: print '%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50) return - mkvc(Jv) - def Jtvec(self, m, v, u=None): + def Jtvec(self, m, v, f=None): """ :param numpy.array m: Conductivity model :param numpy.ndarray,SimPEG.Survey.Data v: vector (data object) @@ -150,15 +150,15 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): """ if self.verbose: print '%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50) self.curModel = m - if u is None: - u = self.fields(m) + if f is None: + f = self.fields(m) if not isinstance(v, self.dataPair): v = self.dataPair(self.survey, v) - p = self.survey.evalDeriv(u, v=v, adjoint=True) + p = self.survey.evalDeriv(f, v=v, adjoint=True) y = self.solveAht(m, p) - w = self.Gtvec(m, y, u) + w = self.Gtvec(m, y, f) if self.verbose: print '%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50) return - mkvc(w) diff --git a/SimPEG/InvProblem.py b/SimPEG/InvProblem.py index 0296bf4b..fd6c48c3 100644 --- a/SimPEG/InvProblem.py +++ b/SimPEG/InvProblem.py @@ -82,23 +82,23 @@ class BaseInvProblem(object): self._warmstart = value def getFields(self, m, store=False, deleteWarmstart=True): - u = None + f = None for mtest, u_ofmtest in self.warmstart: if m is mtest: - u = u_ofmtest + f = u_ofmtest if self.debug: print 'InvProb is Warm Starting!' break - if u is None: - u = self.prob.fields(m) + if f is None: + f = self.prob.fields(m) if deleteWarmstart: self.warmstart = [] if store: - self.warmstart += [(m,u)] + self.warmstart += [(m,f)] - return u + return f @Utils.timeIt def evalFunction(self, m, return_g=True, return_H=True): @@ -109,21 +109,21 @@ class BaseInvProblem(object): gc.collect() # Store fields if doing a line-search - u = self.getFields(m, store=(return_g==False and return_H==False)) + f = self.getFields(m, store=(return_g==False and return_H==False)) - phi_d = self.dmisfit.eval(m, u=u) + phi_d = self.dmisfit.eval(m, f=f) phi_m = self.reg.eval(m) - self.dpred = self.survey.dpred(m, u=u) # This is a cheap matrix vector calculation. + self.dpred = self.survey.dpred(m, f=f) # This is a cheap matrix vector calculation. self.phi_d, self.phi_d_last = phi_d, self.phi_d self.phi_m, self.phi_m_last = phi_m, self.phi_m - f = phi_d + self.beta * phi_m + phi = phi_d + self.beta * phi_m - out = (f,) + out = (phi,) if return_g: - phi_dDeriv = self.dmisfit.evalDeriv(m, u=u) + phi_dDeriv = self.dmisfit.evalDeriv(m, f=f) phi_mDeriv = self.reg.evalDeriv(m) g = phi_dDeriv + self.beta * phi_mDeriv @@ -131,7 +131,7 @@ class BaseInvProblem(object): if return_H: def H_fun(v): - phi_d2Deriv = self.dmisfit.eval2Deriv(m, v, u=u) + phi_d2Deriv = self.dmisfit.eval2Deriv(m, v, f=f) phi_m2Deriv = self.reg.eval2Deriv(m, v=v) return phi_d2Deriv + self.beta * phi_m2Deriv diff --git a/SimPEG/MT/BaseMT.py b/SimPEG/MT/BaseMT.py index 36389430..c201dfb0 100644 --- a/SimPEG/MT/BaseMT.py +++ b/SimPEG/MT/BaseMT.py @@ -27,7 +27,7 @@ class BaseMTProblem(BaseFDEMProblem): # Might need to add more stuff here. ## NEED to clean up the Jvec and Jtvec to use Zero and Identities for None components. - def Jvec(self, m, v, u=None): + def Jvec(self, m, v, f=None): """ Function to calculate the data sensitivities dD/dm times a vector. @@ -39,8 +39,8 @@ class BaseMTProblem(BaseFDEMProblem): """ # Calculate the fields - if u is None: - u = self.fields(m) + if f is None: + f= self.fields(m) # Set current model self.curModel = m # Initiate the Jv object @@ -56,9 +56,9 @@ class BaseMTProblem(BaseFDEMProblem): # We need fDeriv_m = df/du*du/dm + df/dm # Construct du/dm, it requires a solve # NOTE: need to account for the 2 polarizations in the derivatives. - u_src = u[src,:] + f_src = f[src,:] # dA_dm and dRHS_dm should be of size nE,2, so that we can multiply by dA_duI. The 2 columns are each of the polarizations. - dA_dm = self.getADeriv_m(freq, u_src, v) # Size: nE,2 (u_px,u_py) in the columns. + dA_dm = self.getADeriv_m(freq, f_src, v) # Size: nE,2 (u_px,u_py) in the columns. dRHS_dm = self.getRHSDeriv_m(freq, v) # Size: nE,2 (u_px,u_py) in the columns. if dRHS_dm is None: du_dm = dA_duI * ( -dA_dm ) @@ -68,13 +68,13 @@ class BaseMTProblem(BaseFDEMProblem): for rx in src.rxList: # Get the projection derivative # v should be of size 2*nE (for 2 polarizations) - PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, u, t) # wrt u, we don't have have PDeriv wrt m + PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, f, t) # wrt u, we don't have have PDeriv wrt m Jv[src, rx] = PDeriv_u(mkvc(du_dm)) dA_duI.clean() # Return the vectorized sensitivities return mkvc(Jv) - def Jtvec(self, m, v, u=None): + def Jtvec(self, m, v, f=None): """ Function to calculate the transpose of the data sensitivities (dD/dm)^T times a vector. @@ -85,8 +85,8 @@ class BaseMTProblem(BaseFDEMProblem): :return: Data sensitivities wrt m """ - if u is None: - u = self.fields(m) + if f is None: + f = self.fields(m) self.curModel = m @@ -103,15 +103,15 @@ class BaseMTProblem(BaseFDEMProblem): for src in self.survey.getSrcByFreq(freq): ftype = self._fieldType + 'Solution' - u_src = u[src, :] + f_src = f[src, :] for rx in src.rxList: # Get the adjoint evalDeriv # PTv needs to be nE, - PTv = rx.evalDeriv(src, self.mesh, u, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m + PTv = rx.evalDeriv(src, self.mesh, f, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m # Get the dA_duIT = ATinv * PTv - dA_dmT = self.getADeriv_m(freq, u_src, mkvc(dA_duIT), adjoint=True) + dA_dmT = self.getADeriv_m(freq, f_src, mkvc(dA_duIT), adjoint=True) dRHS_dmT = self.getRHSDeriv_m(freq, mkvc(dA_duIT), adjoint=True) # Make du_dmT if dRHS_dmT is None: @@ -129,4 +129,4 @@ class BaseMTProblem(BaseFDEMProblem): raise Exception('Must be real or imag') # Clean the factorization, clear memory. ATinv.clean() - return Jtv \ No newline at end of file + return Jtv diff --git a/SimPEG/MT/SurveyMT.py b/SimPEG/MT/SurveyMT.py index 4e4a8688..0ec91a0e 100644 --- a/SimPEG/MT/SurveyMT.py +++ b/SimPEG/MT/SurveyMT.py @@ -427,15 +427,15 @@ class Survey(SimPEGsurvey.BaseSurvey): assert freq in self._freqDict, "The requested frequency is not in this survey." return self._freqDict[freq] - def eval(self, u): + def eval(self, f): data = Data(self) for src in self.srcList: sys.stdout.flush() for rx in src.rxList: - data[src, rx] = rx.eval(src, self.mesh, u) + data[src, rx] = rx.eval(src, self.mesh, f) return data - def evalDeriv(self, u): + def evalDeriv(self, f): raise Exception('Use Transmitters to project fields deriv.') ################# diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index 8ed94e97..f8520c5c 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -128,7 +128,7 @@ class BaseProblem(object): :rtype: numpy.array :return: approxJv """ - return self.Jvec(m, v, u) + return self.Jvec(m, v, f) @Utils.timeIt def Jtvec_approx(self, m, v, f=None): @@ -142,7 +142,7 @@ class BaseProblem(object): :rtype: numpy.array :return: JTv """ - return self.Jtvec(m, v, u) + return self.Jtvec(m, v, f) def fields(self, m): """ diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 9ddd269e..fbc88276 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -313,20 +313,20 @@ class BaseSurvey(object): @Utils.count - def eval(self, u): - """eval(u) + def eval(self, f): + """eval(f) This function projects the fields onto the data space. .. math:: - d_\\text{pred} = \mathbf{P} u(m) + d_\\text{pred} = \mathbf{P} f(m) """ raise NotImplemented('eval is not yet implemented.') @Utils.count - def evalDeriv(self, u): - """evalDeriv(u) + def evalDeriv(self, f): + """evalDeriv(f) This function s the derivative of projects the fields onto the data space. @@ -379,8 +379,8 @@ class BaseSurvey(object): return self.dobs class LinearSurvey(BaseSurvey): - def eval(self, u): - return u + def eval(self, f): + return f @property def nD(self):