diff --git a/simpegEM/Base.py b/simpegEM/Base.py index 7f938bf6..7f6c9766 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -1,6 +1,14 @@ -from SimPEG import Survey, Problem, Utils, Models, np, sp, Solver as SimpegSolver +from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solver as SimpegSolver from scipy.constants import mu_0 +class EMPropMap(Maps.PropMap): + sigma = Maps.Property("Electrical Conductivity", defaultInvProp = True, propertyLink=('rho',Maps.ReciprocalMap)) + mu = Maps.Property("Inverse Magnetic Permeability", defaultVal = mu_0, propertyLink=('mui',Maps.ReciprocalMap)) + + rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap)) + mui = Maps.Property("Inverse Magnetic Permeability", defaultVal = 1./mu_0, propertyLink=('mu', Maps.ReciprocalMap)) + + class BaseEMProblem(Problem.BaseProblem): def __init__(self, mesh, **kwargs): @@ -9,6 +17,8 @@ class BaseEMProblem(Problem.BaseProblem): surveyPair = Survey.BaseSurvey dataPair = Survey.Data + + PropMap = EMPropMap Solver = SimpegSolver solverOpts = {} @@ -24,94 +34,20 @@ class BaseEMProblem(Problem.BaseProblem): self.__makeASymmetric = True return self.__makeASymmetric - #################################################### - # Phys Props - #################################################### - - # Mu - @property - def mu(self): - if getattr(self, '_mu', None) is None: - # if getattr(self, '_mui', None) is not None: - # self._mu = sel - self._mu = mu_0 - return self._mu - @mu.setter - def mu(self, value): - if getattr(self, '_MfMui', None) is not None: - del self._MfMui - if getattr(self, '_MfMuiI', None) is not None: - del self._MfMuiI - if getattr(self, '_MeMu', None) is not None: - del delf._MeMu - if getattr(self, '_MeMuI', None) is not None: - del self._MeMuI - self._mu = value - - # TODO: hardcoded to assume diagonal mu - @property - def mui(self): - if getattr(self, '_mui', None) is None: - self._mui = 1./mu_0 - return self._mui - @mui.setter - def mui(self, value): - if getattr(self, '_MfMui', None) is not None: - del self._MfMui - if getattr(self, '_MfMuiI', None) is not None: - del self._MfMuiI - if getattr(self, '_MeMu', None) is not None: - del delf._MeMu - if getattr(self, '_MeMuI', None) is not None: - del self._MeMuI - self._mui = value - - # Sigma - # deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmai','_MfSigmaiI'] - @property - def sigma(self): - if getattr(self, '_sigma', None) is None: - self._sigma = self.curModel.transform - return self._sigma - @sigma.setter - def sigma(self, value): - if getattr(self, '_MeSigma', None) is not None: - del self._MeSigma - if getattr(self, '_MeSigmaI', None) is not None: - del self._MeSigmaI - if getattr(self, '_MfSigmai', None) is not None: - del delf._MfSigmai - if getattr(self, '_MfSigmaiI', None) is not None: - del self._MfSigmaiI - self._sigma = value - - # def dsigma_dm(self): - # return self.curModel.transformDeriv - - - # TODO: hardcoded to assume diagonal sigma - @property - def sigmai(self): - if getattr(self, '_sigmai', None) is None: - self._sigmai = 1./self.curModel.transform - return self._sigmai - @sigmai.setter - def sigmai(self, value): - if getattr(self, '_MeSigma', None) is not None: - del self._MeSigma - if getattr(self, '_MeSigmaI', None) is not None: - del self._MeSigmaI - if getattr(self, '_MfSigmai', None) is not None: - del delf._MfSigmai - if getattr(self, '_MfSigmaiI', None) is not None: - del self._MfSigmaiI - self._sigma = value - #################################################### # Mass Matrices #################################################### + @property + def deleteTheseOnModelUpdate(self): + toDelete = [] + if self.mapping.sigmaMap is not None or self.mapping.rhoMap is not None: + toDelete += ['_MeSigma', '_MeSigmaI','_MfRho','_MfRhoI'] + if self.mapping.muMap is not None or self.mapping.muiMap is not None: + toDelete += ['_MeMu', '_MeMuI','_MfMui','_MfMuiI'] + return toDelete + @property def Me(self): if getattr(self, '_Me', None) is None: @@ -129,25 +65,25 @@ class BaseEMProblem(Problem.BaseProblem): @property def MfMui(self): if getattr(self, '_MfMui', None) is None: - self._MfMui = self.mesh.getFaceInnerProduct(self.mui) + self._MfMui = self.mesh.getFaceInnerProduct(self.curModel.mui) return self._MfMui @property def MfMuiI(self): if getattr(self, '_MfMuiI', None) is None: - self._MfMuiI = self.mesh.getFaceInnerProduct(self.mui, invMat=True) + self._MfMuiI = self.mesh.getFaceInnerProduct(self.curModel.mui, invMat=True) return self._MfMuiI @property def MeMu(self): if getattr(self, '_MeMu', None) is None: - self._MeMu = self.mesh.getEdgeInnerProduct(self.mu) + self._MeMu = self.mesh.getEdgeInnerProduct(self.curModel.mu) return self._MeMu @property def MeMuI(self): if getattr(self, '_MeMuI', None) is None: - self._MeMuI = self.mesh.getEdgeInnerProduct(self.mu, invMat=True) + self._MeMuI = self.mesh.getEdgeInnerProduct(self.curModel.mu, invMat=True) return self._MeMuI # ----- Electrical Conductivity ----- # @@ -155,42 +91,43 @@ class BaseEMProblem(Problem.BaseProblem): @property def MeSigma(self): if getattr(self, '_MeSigma', None) is None: - self._MeSigma = self.mesh.getEdgeInnerProduct(self.sigma) + self._MeSigma = self.mesh.getEdgeInnerProduct(self.curModel.sigma) return self._MeSigma - # def dMeSigma_dsigma(self, u): - # return self.mesh.getEdgeInnerProductDeriv(self.sigma)(u) + def MeSigmaDeriv(self, u): + """ + Deriv of MeSigma wrt sigma + """ + return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) + @property def MeSigmaI(self): if getattr(self, '_MeSigmaI', None) is None: - self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.sigma, invMat=True) + self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.curModel.sigma, invMat=True) return self._MeSigmaI - # def dMeSigmaI_dsigma(self,u) + def MeSigmaIDeriv(self, u): + """ + Deriv of MeSigma wrt sigma + """ + return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma, invMat=True)(u) + @property - def MfSigmai(self): - if getattr(self, '_MfSigmai', None) is None: - self._MfSigmai = self.mesh.getFaceInnerProduct(self.sigmai) - return self._MfSigmai + def MfRho(self): + if getattr(self, '_MfRho', None) is None: + self._MfRho = self.mesh.getFaceInnerProduct(self.curModel.rho) + return self._MfRho - # def dMfSigmai_dsigmai(self,u) + def MfRhoDeriv(self,u): + return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) @property - def MfSigmaiI(self): - if getattr(self, '_MfSigmaiI', None) is None: - self._MfSigmaiI = self.mesh.getFaceInnerProduct(self.sigmai, invMat=True) - return self._MfSigmaiI + def MfRhoI(self): + if getattr(self, '_MfRhoI', None) is None: + self._MfRhoI = self.mesh.getFaceInnerProduct(self.curModel.rho, invMat=True) + return self._MfRhoI - # def dMfSigmaiI(self,u) - - - #################################################### - # Fields - #################################################### - - def fields(self, m): - self.curModel = m - F = self.forward(m) - return F + def dMfRhoIDeriv(self,u): + return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True) diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index e83fe0e0..6564183d 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -19,8 +19,8 @@ class BaseFDEMProblem(BaseEMProblem): surveyPair = SurveyFDEM fieldsPair = FieldsFDEM - def forward(self, m): - + def fields(self, m): + self.curModel = m F = self.fieldsPair(self.mesh, self.survey) for freq in self.survey.freqs: @@ -175,9 +175,11 @@ class ProblemFDEM_e(BaseFDEMProblem): def getADeriv(self, freq, u, v, adjoint=False): - sig = self.sigma - dsig_dm = self.curModel.transformDeriv - dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(u) + # sig = self.sigma + # dsig_dm = self.curModel.transformDeriv + # dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(u) + dsig_dm = self.curModel.sigmaDeriv + dMe_dsig = self.MeSigmaDeriv(u) if adjoint: return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) ) @@ -327,14 +329,14 @@ class ProblemFDEM_j(BaseFDEMProblem): """ MeMuI = self.MeMuI - MfSigmai = self.MfSigmai + MfRho = self.MfRho C = self.mesh.edgeCurl iomega = 1j * omega(freq) * sp.eye(self.mesh.nF) - A = C * MeMuI * C.T * MfSigmai + iomega + A = C * MeMuI * C.T * MfRho + iomega if self._makeASymmetric is True: - return MfSigmai.T*A + return MfRho.T*A return A @@ -347,7 +349,7 @@ class ProblemFDEM_j(BaseFDEMProblem): """ MeMuI = self.MeMuI - MfSigmai = self.MfSigmai + MfRho = self.MfRho C = self.mesh.edgeCurl sigi = self.sigmai dsig_dm = self.curModel.transformDeriv @@ -356,11 +358,11 @@ class ProblemFDEM_j(BaseFDEMProblem): if adjoint: if self._makeASymmetric is True: - v = MfSigmai * v + v = MfRho * v return dsig_dm.T * ( dsigi_dsig.T *( dMf_dsigi.T * ( C * ( MeMuI.T * ( C.T * v ) ) ) ) ) if self._makeASymmetric is True: - return MfSigmai.T * ( C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) ) + return MfRho.T * ( C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) ) return C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) @@ -378,8 +380,8 @@ class ProblemFDEM_j(BaseFDEMProblem): RHS = C * (MeMuI * S_m) - 1j * omega(freq) * S_e if self._makeASymmetric is True: - MfSigmai = self.MfSigmai - return MfSigmai.T*RHS + MfRho = self.MfRho + return MfRho.T*RHS return RHS @@ -429,10 +431,10 @@ class ProblemFDEM_h(BaseFDEMProblem): """ MeMu = self.MeMu - MfSigmai = self.MfSigmai + MfRho = self.MfRho C = self.mesh.edgeCurl - return C.T * MfSigmai * C + 1j*omega(freq)*MeMu + return C.T * MfRho * C + 1j*omega(freq)*MeMu def getADeriv(self, freq, u, v, adjoint=False): @@ -458,9 +460,9 @@ class ProblemFDEM_h(BaseFDEMProblem): S_m, S_e = self.getSourceTerm(freq) C = self.mesh.edgeCurl - MfSigmai = self.MfSigmai + MfRho = self.MfRho - RHS = S_m + C.T * ( MfSigmai * S_e ) + RHS = S_m + C.T * ( MfRho * S_e ) return RHS diff --git a/simpegEM/FDEM/FieldsFDEM.py b/simpegEM/FDEM/FieldsFDEM.py index e12752b6..ba4dbdd3 100644 --- a/simpegEM/FDEM/FieldsFDEM.py +++ b/simpegEM/FDEM/FieldsFDEM.py @@ -24,9 +24,9 @@ class FieldsFDEM_e(FieldsFDEM): def _e(self, e_sol, srcList): e = e_sol for i, src in enumerate(srcList): - e_p = src.e_p(self.survey.prob) - if e_p is not None: - e[:,i] += e_p + ePrimary = src.ePrimary(self.survey.prob) + if ePrimary is not None: + e[:,i] += ePrimary return e def _b(self, e_sol, srcList): @@ -38,9 +38,9 @@ class FieldsFDEM_e(FieldsFDEM): if S_m is not None: b[:,i] += 1./(1j*omega(src.freq)) * S_m - b_p = src.b_p(self.survey.prob) - if b_p is not None: - b[:,i] += b_p + bPrimary = src.bPrimary(self.survey.prob) + if bPrimary is not None: + b[:,i] += bPrimary return b @@ -72,9 +72,9 @@ class FieldsFDEM_b(FieldsFDEM): b = b_sol for i, src in enumerate(srcList): - b_p = src.b_p(self.survey.prob) - if b_p is not None: - b[:,i] += b_p + bPrimary = src.bPrimary(self.survey.prob) + if bPrimary is not None: + b[:,i] += bPrimary return b def _e(self, b_sol, srcList): @@ -85,9 +85,9 @@ class FieldsFDEM_b(FieldsFDEM): if S_e is not None: e += -self._MeSigmaI*S_e - e_p = src.e_p(self.survey.prob) - if e_p is not None: - e[:,i] += e_p + ePrimary = src.ePrimary(self.survey.prob) + if ePrimary is not None: + e[:,i] += ePrimary return e @@ -114,23 +114,23 @@ class FieldsFDEM_j(FieldsFDEM): def startup(self): self._edgeCurl = self.survey.prob.mesh.edgeCurl self._MeMuI = self.survey.prob.MeMuI - self._MfSigmai = self.survey.prob.MfSigmai + self._MfRho = self.survey.prob.MfRho self._curModel = self.survey.prob.curModel def _j(self, j_sol, srcList): j = j_sol for i, src in enumerate(srcList): - j_p = src.j_p(self.survey.prob) - if j_p is not None: - j[:,i] += j_p + jPrimary = src.jPrimary(self.survey.prob) + if jPrimary is not None: + j[:,i] += jPrimary return j def _h(self, j_sol, srcList): MeMuI = self._MeMuI C = self._edgeCurl - MfSigmai = self._MfSigmai + MfRho = self._MfRho - h = MeMuI * (C.T * (MfSigmai * j_sol) ) + h = MeMuI * (C.T * (MfRho * j_sol) ) for i, src in enumerate(srcList): h[:,i] *= -1./(1j*omega(src.freq)) @@ -138,9 +138,9 @@ class FieldsFDEM_j(FieldsFDEM): if S_m is not None: h[:,i] += 1./(1j*omega(src.freq)) * MeMuI * S_m - h_p = src.h_p(self.survey.prob) - if h_p is not None: - h[:,i] += h_p + hPrimary = src.hPrimary(self.survey.prob) + if hPrimary is not None: + h[:,i] += hPrimary return h @@ -151,7 +151,7 @@ class FieldsFDEM_j(FieldsFDEM): dsig_dm = self._curModel.transformDeriv dsigi_dsig = -Utils.sdiag(sigi)**2 dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j) - sigi = self._MfSigmai + sigi = self._MfRho S_mDeriv,_ = src.getSourceDeriv(self.survey.prob, v, adjoint) @@ -177,14 +177,14 @@ class FieldsFDEM_h(FieldsFDEM): def startup(self): self._edgeCurl = self.survey.prob.mesh.edgeCurl self._MeMuI = self.survey.prob.MeMuI - self._MfSigmai = self.survey.prob.MfSigmai + self._MfRho = self.survey.prob.MfRho def _h(self, h_sol, srcList): h = h_sol for i, src in enumerate(srcList): - h_p = src.h_p(self.survey.prob) - if h_p is not None: - h[:,i] += h_p + hPrimary = src.hPrimary(self.survey.prob) + if hPrimary is not None: + h[:,i] += hPrimary return h def _j(self, h_sol, srcList): @@ -194,9 +194,9 @@ class FieldsFDEM_h(FieldsFDEM): if S_e is not None: j[:,i] += -S_e - j_p = src.j_p(self.survey.prob) - if j_p is not None: - j[:,i] += j_p + jPrimary = src.jPrimary(self.survey.prob) + if jPrimary is not None: + j[:,i] += jPrimary return j def _jDeriv(self, h_sol, srcList, v, adjoint=False): @@ -215,11 +215,11 @@ class FieldsFDEM_h(FieldsFDEM): # elif fieldType == 'h': # MeMuI = self._MeMuI # C = self.mesh.edgeCurl - # MfSigmai = self._MfSigmai + # MfRho = self._MfRho # if not adjoint: - # h = -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( MfSigmai * j ) ) + # h = -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( MfRho * j ) ) # else: - # h = -(1./(1j*omega(freq))) * MfSigmai.T * ( C * ( MeMuI.T * j ) ) + # h = -(1./(1j*omega(freq))) * MfRho.T * ( C * ( MeMuI.T * j ) ) # return h # raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) @@ -235,7 +235,7 @@ class FieldsFDEM_h(FieldsFDEM): # dsig_dm = self._curModel.transformDeriv # dsigi_dsig = -Utils.sdiag(sigi)**2 # dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j) - # sigi = self._MfSigmai + # sigi = self._MfRho # if not adjoint: # return -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) # else: diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index ff1f50fe..e2df08c9 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -96,51 +96,35 @@ class SrcFDEM(Survey.BaseSrc): rxPair = RxFDEM def eval(self, prob): - S_m = self._getS_m(prob) - S_e = self._getS_e(prob) + S_m = self.S_m(prob) + S_e = self.S_e(prob) return S_m, S_e def evalDeriv(self, prob, v, adjoint=None): - return self._getS_mDeriv(prob,v,adjoint), self._getS_eDeriv(prob,v,adjoint) + return self.S_mDeriv(prob,v,adjoint), self.S_eDeriv(prob,v,adjoint) - def b_p(self,prob): - b_p = self._getb_p(prob) - return b_p + def bPrimary(self,prob): + return None - def h_p(self,prob): - h_p = self._geth_p(prob) - return h_p - - def e_p(self,prob): - e_p = self._gete_p(prob) - return e_p - - def j_p(self,prob): - j_p = self._getj_p(prob) - return j_p - - def _getb_p(self,prob): + def hPrimary(self,prob): return None - def _geth_p(self,prob): + def ePrimary(self,prob): return None - def _gete_p(self,prob): + def jPrimary(self,prob): return None - def _getj_p(self,prob): + def S_m(self,prob): return None - def _getS_m(self,prob): + def S_e(self,prob): return None - def _getS_e(self,prob): + def S_mDeriv(self, prob, v, adjoint = False): return None - def _getS_mDeriv(self, prob, v, adjoint = False): - return None - - def _getS_eDeriv(self, prob, v, adjoint = False): + def S_eDeriv(self, prob, v, adjoint = False): return None @@ -154,12 +138,12 @@ class SrcFDEM_RawVec_e(SrcFDEM): """ def __init__(self, rxList, freq, S_e): - self.S_e = np.array(S_e,dtype=float) + self._S_e = np.array(S_e,dtype=float) self.freq = float(freq) SrcFDEM.__init__(self, rxList) - def _getS_e(self, prob): - return self.S_e + def S_e(self, prob): + return self._S_e class SrcFDEM_RawVec_m(SrcFDEM): @@ -172,12 +156,12 @@ class SrcFDEM_RawVec_m(SrcFDEM): """ def __init__(self, rxList, freq, S_m): - self.S_m = np.array(S_m,dtype=float) + self._S_m = np.array(S_m,dtype=float) self.freq = float(freq) SrcFDEM.__init__(self, rxList) - def _getS_m(self, prob): - return self.S_m + def S_m(self, prob): + return self._S_m class SrcFDEM_RawVec(SrcFDEM): @@ -190,16 +174,16 @@ class SrcFDEM_RawVec(SrcFDEM): :param rxList: receiver list """ def __init__(self, rxList, freq, S_m, S_e): - self.S_m = np.array(S_m,dtype=float) - self.S_e = np.array(S_e,dtype=float) + self._S_m = np.array(S_m,dtype=float) + self._S_e = np.array(S_e,dtype=float) self.freq = float(freq) SrcFDEM.__init__(self, rxList) - def _getS_m(self,prob): - return self.S_m + def S_m(self,prob): + return self._S_m - def _getS_e(self,prob): - return self.S_e + def S_e(self,prob): + return self._S_e class SrcFDEM_MagDipole(SrcFDEM): @@ -212,7 +196,7 @@ class SrcFDEM_MagDipole(SrcFDEM): self.moment = moment SrcFDEM.__init__(self, rxList) - def _getb_p(self,prob): + def bPrimary(self,prob): eqLocs = prob._eqLocs if eqLocs is 'FE': @@ -243,12 +227,12 @@ class SrcFDEM_MagDipole(SrcFDEM): return C*a - def _geth_p(self,prob): - b = self._getb_p(prob) + def hPrimary(self,prob): + b = self.bPrimary(prob) return h_from_b(prob,b) - def _getS_m(self,prob): - b_p = self._getb_p(prob) + def S_m(self,prob): + b_p = self.bPrimary(prob) return -1j*omega(self.freq)*b_p @@ -264,7 +248,7 @@ class SrcFDEM_MagDipole_Bfield(SrcFDEM): self.moment = moment SrcFDEM.__init__(self, rxList) - def _getb_p(self,prob): + def bPrimary(self,prob): eqLocs = prob._eqLocs if eqLocs is 'FE': @@ -295,12 +279,12 @@ class SrcFDEM_MagDipole_Bfield(SrcFDEM): return b - def _geth_p(self,prob): - b = self._getb_p(prob) + def hPrimary(self,prob): + b = self.bPrimary(prob) return h_from_b(prob, b) - def _getS_m(self,prob): - b = self._getb_p(prob) + def S_m(self,prob): + b = self.bPrimary(prob) return -1j*omega(self.freq)*b @@ -313,7 +297,7 @@ class SrcFDEM_CircularLoop(SrcFDEM): self.radius = radius SrcFDEM.__init__(self, rxList) - def _getb_p(self,prob): + def bPrimary(self,prob): eqLocs = prob._eqLocs if eqLocs is 'FE': @@ -343,12 +327,12 @@ class SrcFDEM_CircularLoop(SrcFDEM): return C*a - def _geth_p(self,prob): - b = self._getb_p(prob) + def hPrimary(self,prob): + b = self.bPrimary(prob) return h_from_b - def _getS_m(self, prob): - b = self._getb_p(prob) + def S_m(self, prob): + b = self.bPrimary(prob) return -1j*omega(self.freq)*b diff --git a/simpegEM/Utils/EMUtils.py b/simpegEM/Utils/EMUtils.py index 5cdaf150..916d6ae2 100644 --- a/simpegEM/Utils/EMUtils.py +++ b/simpegEM/Utils/EMUtils.py @@ -16,7 +16,7 @@ def e_from_j(prob,j): if eqLocs is 'FE': MSigmaI = prob.MeSigmaI elif eqLocs is 'EF': - MSigmaI = prob.MfSigmai + MSigmaI = prob.MfRho return MSigmaI*j def j_from_e(prob,e): @@ -24,13 +24,13 @@ def j_from_e(prob,e): if eqLocs is 'FE': MSigma = prob.MeSigma elif eqLocs is 'EF': - MSigma = prob.MfSigma + MSigma = prob.MfRhoI return MSigma*e def b_from_h(prob,h): eqLocs = prob._eqLocs if eqLocs is 'FE': - MMu = prob.MfMu + MMu = prob.MfMuiI elif eqLocs is 'EF': MMu = prob.MeMu return MMu*h