Commit 338ce526 authored by Leslie Wade's avatar Leslie Wade

Enforce additional prior check on eos gammas in approx prior draw

parent c45d99d7
......@@ -2363,11 +2363,11 @@ XLALDestroySimNeutronStarEOS(eos);
/* Find lambda1,2(m1,2|eos) for spectral EOS model */
void LALInferenceSDGammasMasses2Lambdas(REAL8 gamma[], REAL8 mass1, REAL8 mass2, REAL8 *lambda1, REAL8 *lambda2, int size){
// If unreasonable gammas, do not find lambdas
// If gammas outside prior, do not find lambdas
if(LALInferenceSDGammaCheck(gamma, 4) == XLAL_FAILURE){
*lambda1= 0.;
*lambda2= 0.;
}
}
// Else calculate lambdas
else{
// Convert to SI
......@@ -2395,7 +2395,7 @@ else{
// Clean up
XLALDestroySimNeutronStarFamily(fam);
XLALDestroySimNeutronStarEOS(eos);
}
}
}
......
......@@ -1275,8 +1275,7 @@ REAL8 LALInferenceDrawApproxPrior(LALInferenceThreadState *thread,
const char *flat_params[] = {"q", "eta", "t0", "azimuth", "cosalpha", "time", "phase", "polarisation",
"rightascension", "costheta_jn", "phi_jl",
"phi12", "a_spin1", "a_spin2", "logp1", "gamma1", "gamma2", "gamma3",
"SDgamma0","SDgamma1","SDgamma2","SDgamma3", NULL};
"phi12", "a_spin1", "a_spin2", "logp1", "gamma1", "gamma2", "gamma3", NULL};
LALInferenceVariables *args = thread->proposalArgs;
......@@ -1294,6 +1293,29 @@ REAL8 LALInferenceDrawApproxPrior(LALInferenceThreadState *thread,
} else {
logBackwardJump = approxLogPrior(currentParams);
if (LALInferenceCheckVariableNonFixed(proposedParams, "SDgamma0")) {
// Find EOS spectral params in prior range
const char *gamma_params[] = {"SDgamma0","SDgamma1","SDgamma2","SDgamma3", NULL};
double gamma[]={*(double *)LALInferenceGetVariable(proposedParams,"SDgamma0"),
*(double *)LALInferenceGetVariable(proposedParams,"SDgamma1"),
*(double *)LALInferenceGetVariable(proposedParams,"SDgamma2"),
*(double *)LALInferenceGetVariable(proposedParams,"SDgamma3")};
// Draw from flat prior until GammaCheck (another prior constraint) is passed
do{
for (i = 0; gamma_params[i] != NULL; i++) {
if (LALInferenceCheckVariableNonFixed(proposedParams, gamma_params[i])) {
REAL8 val = draw_flat(thread, gamma_params[i]);
LALInferenceSetVariable(proposedParams, gamma_params[i], &val);
}
}
// Fill gamma array for GammaCheck
gamma[0]=*(double *)LALInferenceGetVariable(proposedParams,"SDgamma0");
gamma[1]=*(double *)LALInferenceGetVariable(proposedParams,"SDgamma1");
gamma[2]=*(double *)LALInferenceGetVariable(proposedParams,"SDgamma2");
gamma[3]=*(double *)LALInferenceGetVariable(proposedParams,"SDgamma3");
}while(LALInferenceSDGammaCheck(gamma, 4) == XLAL_FAILURE);
}
for (i = 0; flat_params[i] != NULL; i++) {
if (LALInferenceCheckVariableNonFixed(proposedParams, flat_params[i])) {
REAL8 val = draw_flat(thread, flat_params[i]);
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment