Commit c1822acd authored by Carl-Johan Haster's avatar Carl-Johan Haster Committed by Adam Mercer

Add binary love EoS model

(cherry picked from commit 6117be18)
parent 835f9c80
......@@ -142,7 +142,7 @@ ppEParams=['ppEalpha','ppElowera','ppEupperA','ppEbeta','ppElowerb','ppEupperB',
tigerParams=['dchi%i'%(i) for i in range(8)] + ['dchi%il'%(i) for i in [5,6] ] + ['dxi%d'%(i+1) for i in range(6)] + ['dalpha%i'%(i+1) for i in range(5)] + ['dbeta%i'%(i+1) for i in range(3)] + ['dsigma%i'%(i+1) for i in range(4)]
bransDickeParams=['omegaBD','ScalarCharge1','ScalarCharge2']
massiveGravitonParams=['lambdaG']
tidalParams=['lambda1','lambda2','lam_tilde','dlam_tilde','lambdat','dlambdat']
tidalParams=['lambda1','lambda2','lam_tilde','dlam_tilde','lambdat','dlambdat','lambdas','bluni']
eosParams=['logp1','gamma1','gamma2','gamma3']
energyParams=['e_rad', 'l_peak']
strongFieldParams=ppEParams+tigerParams+bransDickeParams+massiveGravitonParams+tidalParams+eosParams+energyParams
......@@ -343,6 +343,8 @@ def get_prior(name):
'lambda2': 'uniform',
'lam_tilde' : None,
'dlam_tilde': None,
'lambdas':'uniform',
'bluni':'uniform',
'logp1':None,
'gamma1':None,
'gamma2':None,
......@@ -484,7 +486,9 @@ def plot_label(param):
'h1_optimal_snr':r'$\rho^{opt}_{H1}$',
'l1_optimal_snr':r'$\rho^{opt}_{L1}$',
'v1_optimal_snr':r'$\rho^{opt}_{V1}$',
'matched_filter_snr':r'$\rho^{MF}$'
'matched_filter_snr':r'$\rho^{MF}$',
'lambdas':r'$\Lambda_S$',
'bluni' : r'$BL_{uniform}$'
}
# Handle cases where multiple names have been used
......@@ -5393,7 +5397,7 @@ def find_ndownsample(samples, nDownsample):
"time_mean", "time_maxl","sky_frame","psdscaleflag","logdeltaf","flow","f_ref",
"lal_amporder","lal_pnorder","lal_approximant","tideo","spino","signalmodelflag",
"temperature","nifo","nlocaltemps","ntemps","randomseed","samplerate","segmentlength","segmentstart",
"t0", "phase_maxl", "azimuth", "cosalpha", "lal_amporder"] + logParams + snrParams + splineParams
"t0", "phase_maxl", "azimuth", "cosalpha", "lal_amporder", "bluni"] + logParams + snrParams + splineParams
fixedParams = [p for p in samples.colnames if all(x==samples[p][0] for x in samples[p])]
print("Fixed parameters: "+str(fixedParams))
nonParams.extend(fixedParams)
......
......@@ -4314,3 +4314,102 @@ void LALInferencePrintSplineCalibration(FILE *output, LALInferenceThreadState *t
}
}
}
void LALInferenceBinaryLove(LALInferenceVariables *vars, REAL8 *lambda1, REAL8 *lambda2){
/**
* Implelentation of the parametrization from Chatziioannou, Haster, Zimmerman (2018)
* https://doi.org/10.1103/PhysRevD.97.104036
*/
REAL8 lambdaS = *(REAL8*) LALInferenceGetVariable(vars, "lambdaS");
REAL8 lambdaSm1o5 = pow(lambdaS,-1.0/5.0);
REAL8 lambdaSm2o5 = lambdaSm1o5*lambdaSm1o5;
REAL8 lambdaSm3o5 = lambdaSm2o5*lambdaSm1o5;
REAL8 lambdaSsqrt = sqrt(lambdaS);
REAL8 q=1.0;
if (LALInferenceCheckVariable(vars,"q")) {
q = *(REAL8 *)LALInferenceGetVariable(vars,"q");
}
else{
REAL8 m1 = *(REAL8 *)LALInferenceGetVariable(vars,"mass1");
REAL8 m2 = *(REAL8 *)LALInferenceGetVariable(vars,"mass2");
q = m2/m1;
}
REAL8 q2 = q*q;
/* Eqn 2 , incorporating the dependance on mass ratio */
REAL8 BL_n = 0.743;
REAL8 q_for_Fnofq = pow(q,10.0/(3.0-BL_n));
REAL8 Fnofq = (1.0-q_for_Fnofq)/(1.0+q_for_Fnofq);
/* bXX and cXX coefficients are given in Table I */
REAL8 b11 = -27.7408;
REAL8 b12 = 8.42358;
REAL8 b21 = 122.686;
REAL8 b22 = -19.7551;
REAL8 b31 = -175.496;
REAL8 b32 = 133.708;
REAL8 c11 = -25.5593;
REAL8 c12 = 5.58527;
REAL8 c21 = 92.0337;
REAL8 c22 = 26.8586;
REAL8 c31 = -70.247;
REAL8 c32 = -56.3076i;
/* Eqn 1, giving the lambdaA_fitOnly, not yet accounting for the uncertainty in the fit */
REAL8 numerator = 1.0 + (b11*q*lambdaSm1o5) + (b12*q2*lambdaSm1o5) + (b21*q*lambdaSm2o5) + (b22*q2*lambdaSm2o5) + (b31*q*lambdaSm3o5) + (b32*q2*lambdaSm3o5);
REAL8 denominator = 1.0 + (c11*q*lambdaSm1o5) + (c12*q2*lambdaSm1o5) + (c21*q*lambdaSm2o5) + (c22*q2*lambdaSm2o5) + (c31*q*lambdaSm3o5) + (c32*q2*lambdaSm3o5);
REAL8 lambdaA_fitOnly = Fnofq*lambdaS*numerator/denominator;
/* Eqn 6, correction on fit for lambdaA caused by uncertainty in the mean of the lambdaS residual fit,
* using coefficients mu_1, mu_2 and mu_3 from Table II */
REAL8 lambdaA_lambdaS_meanCorr = (137.1252739/(lambdaS*lambdaS)) - (32.8026613/lambdaS) + 0.5168637;
/* Eqn 8, correction on fit for lambdaA caused by uncertainty in the standard deviation of lambdaS residual fit,
* using coefficients sigma_1, sigma_2 sigma_3 and sigma_4 from Table II */
REAL8 lambdaA_lambdaS_stdCorr = (-0.0000739*lambdaS*lambdaSsqrt) + (0.0103778*lambdaS) + (0.4581717*lambdaSsqrt) - 0.8341913;
/* Eqn 7, correction on fit for lambdaA caused by uncertainty in the mean of the q residual fit,
* using coefficients mu_4 and mu_5 from Table II */
REAL8 lambdaA_q_meanCorr = (-11.2765281*q2) + (14.9499544*q) - 4.6638851;
/* Eqn 9, correction on fit for lambdaA caused by uncertainty in the standard deviation of the q residual fit,
* using coefficients sigma_5, sigma_6 and sigma_7 from Table II */
REAL8 lambdaA_q_stdCorr = (-201.4323962*q2) + (273.9268276*q) - 71.2342246;
/* Eqn 4, averaging the corrections from the mean of the residual fits */
REAL8 lambdaA_meanCorr = (lambdaA_lambdaS_meanCorr + lambdaA_q_meanCorr)/2.0;
/* Eqn 5, averaging the corrections from the standard deviations of the residual fits */
REAL8 lambdaA_stdCorr = sqrt((lambdaA_lambdaS_stdCorr*lambdaA_lambdaS_stdCorr) + (lambdaA_q_stdCorr*lambdaA_q_stdCorr));
/* Draw a correction on the fit from a Gaussian distribution wiht width lambdaA_stdCorr
* this is done by sampling a inverse cdf through a U{0,1} variable called BLuni */
REAL8 BLuni = *(REAL8*) LALInferenceGetVariable(vars, "BLuni");
REAL8 lambdaA_scatter = gsl_cdf_gaussian_Pinv(BLuni, lambdaA_stdCorr);
/* Add the correction of the residual mean and the Gaussian scatter to the lambdaA_fitOnly value */
REAL8 lambdaA = lambdaA_fitOnly + (lambdaA_meanCorr + lambdaA_scatter);
*lambda1 = lambdaS - lambdaA;
*lambda2 = lambdaS + lambdaA;
return;
}
......@@ -86,6 +86,8 @@
#include <gsl/gsl_randist.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_cdf.h>
#include <sys/time.h>
/*LIB imports*/
......@@ -1182,6 +1184,12 @@ SWIGLAL_CLEAR(OWNS_THIS_STRING(const CHAR*, value));
*/
void LALInferenceFprintSplineCalibrationHeader(FILE *output, LALInferenceThreadState *thread);
/**
* Compute Tidal deformabilities following BinaryLove Universal relations
*/
void LALInferenceBinaryLove(LALInferenceVariables *vars, REAL8 *lambda1, REAL8 *lambda2);
/**
* Conversion routines between Equatorial (RA,DEC) and detector-based coordinate systems, where
* new "north pole" points along vector from det0 to det1.
......
......@@ -745,6 +745,7 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
(--4PolyEOS) Enables 4-piece polytropic EOS parmeterization, only with LALSimulation.\n\
(--4SpectralDecomp) Enables 4-coeff. spectral decomposition EOS parmeterization, only with LALSimulation.\n\
(--eos EOS) Fix the neutron star EOS. Use \"--eos help\" for allowed names\n\
(--BinaryLove) Enable the Binary Neutron Star common EoS tidal model, expressed through EOS-insensitive relations\n\
(--spinOrder PNorder) Specify twice the PN order (e.g. 5 <==> 2.5PN) of spin effects to use, only for LALSimulation (default: -1 <==> Use all spin effects).\n\
(--tidalOrder PNorder) Specify twice the PN order (e.g. 10 <==> 5PN) of tidal effects to use, only for LALSimulation (default: -1 <==> Use all tidal effects).\n\
(--numreldata FileName) Location of NR data file for NR waveforms (with NR_hdf5 approx).\n\
......@@ -789,6 +790,8 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
gamma1 gamma1.\n\
gamma2 gamma2.\n\
gamma3 gamma3.\n\
(requires --BinaryLove)\n\
lambdaS Symmetric tidal deformability.\n\
* \n\
----------------------------------------------\n\
--- Prior Ranges -----------------------------\n\
......@@ -877,6 +880,8 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
REAL8 gamma2Max=4.5;
REAL8 gamma3Min=1.1;
REAL8 gamma3Max=4.5;
REAL8 lambdaSMin=0.0;
REAL8 lambdaSMax=5000.0;
gsl_rng *GSLrandom=state->GSLrandom;
REAL8 endtime=0.0, timeParam=0.0;
REAL8 timeMin=endtime-dt,timeMax=endtime+dt;
......@@ -1363,9 +1368,9 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
if((!!LALInferenceGetProcParamVal(commandLine,"--tidalT") + !!LALInferenceGetProcParamVal(commandLine,"--tidal")
+ !!LALInferenceGetProcParamVal(commandLine,"--4PolyEOS") + !!LALInferenceGetProcParamVal(commandLine,"--4SpectralDecomp")
+ !!LALInferenceGetProcParamVal(commandLine,"--eos")) > 1 )
+ !!LALInferenceGetProcParamVal(commandLine,"--eos") + !!LALInferenceGetProcParamVal(commandLine,"--BinaryLove")) > 1 )
{
XLALPrintError("Error: cannot use more than one of --tidalT, --tidal, --4PolyEOS, --4SpectralDecomp and --eos.\n");
XLALPrintError("Error: cannot use more than one of --tidalT, --tidal, --4PolyEOS, --4SpectralDecomp, --eos and --BinaryLove.\n");
XLAL_ERROR_NULL(XLAL_EINVAL);
}
......@@ -1382,9 +1387,7 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
LALInferenceRegisterUniformVariableREAL8(state, model->params, "gamma1", zero, gamma1Min, gamma1Max, LALINFERENCE_PARAM_LINEAR);
LALInferenceRegisterUniformVariableREAL8(state, model->params, "gamma2", zero, gamma2Min, gamma2Max, LALINFERENCE_PARAM_LINEAR);
LALInferenceRegisterUniformVariableREAL8(state, model->params, "gamma3", zero, gamma3Min, gamma3Max, LALINFERENCE_PARAM_LINEAR);
}
else if((ppt=LALInferenceGetProcParamVal(commandLine,"--eos")))
{
} else if((ppt=LALInferenceGetProcParamVal(commandLine,"--eos"))){
LALSimNeutronStarEOS *eos=NULL;
errnum=XLAL_SUCCESS;
XLAL_TRY(eos=XLALSimNeutronStarEOSByName(ppt->value), errnum);
......@@ -1395,8 +1398,12 @@ LALInferenceModel *LALInferenceInitCBCModel(LALInferenceRunState *state) {
if(errnum!=XLAL_SUCCESS)
XLAL_ERROR_NULL(errnum,"%s: %s",__func__,XLALErrorString(errnum));
if(!model->eos_fam) XLAL_ERROR_NULL(XLAL_EINVAL, "Unable to initialise EOS family");
}
// Pull in symmetric tidal deformability (lambdaS) and the uniform variable used to marginlise over the BinaryLove fit uncertainty
} else if((ppt=LALInferenceGetProcParamVal(commandLine,"--BinaryLove"))){
LALInferenceRegisterUniformVariableREAL8(state, model->params, "lambdaS", zero, lambdaSMin, lambdaSMax, LALINFERENCE_PARAM_LINEAR);
LALInferenceRegisterUniformVariableREAL8(state, model->params, "BLuni", 0.5, 0.0, 1.0, LALINFERENCE_PARAM_LINEAR);
}
LALSimInspiralSpinOrder spinO = LAL_SIM_INSPIRAL_SPIN_ORDER_ALL;
ppt=LALInferenceGetProcParamVal(commandLine, "--spinOrder");
if(ppt) {
......
......@@ -1372,6 +1372,21 @@ static REAL8 LALInferenceFusedFreqDomainLogLikelihood(LALInferenceVariables *cur
//loglikelihood = -1.0 * chisquared; // note (again): the log-likelihood is unnormalised!
if(LALInferenceCheckVariable(model->params, "lambdaS")){
REAL8 lambda1 = *(REAL8 *)LALInferenceGetVariable(model->params,"lambda1");
REAL8 lambda2 = *(REAL8 *)LALInferenceGetVariable(model->params,"lambda2");
LALInferenceAddVariable(currentParams,"lambda1",&lambda1,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
LALInferenceAddVariable(currentParams,"lambda2",&lambda2,LALINFERENCE_REAL8_t,LALINFERENCE_PARAM_OUTPUT);
// The BinaryLove model is only physically valid where lambda2 > lambda1 as it assumes m1 > m2
// It also assumes both lambda1 and lambda2 to be positive
// This is an explicit feature of the "raw" model fit, but since this implementation also incorporates
// marginalisation over the fit uncertainty, there can be instances where those assumptions are randomly broken
// for those cases, set logL = -inf.
if( (lambda1 > lambda2) || (lambda1 < 0.0) || (lambda2 < 0.0)){
loglikelihood = -INFINITY;
}
}
return(loglikelihood);
}
......
......@@ -392,6 +392,17 @@ void LALInferenceROQWrapperForXLALSimInspiralChooseFDWaveformSequence(LALInferen
XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
}
/* ==== BINARY_LOVE PARAMETERS ==== */
if(LALInferenceCheckVariable(model->params, "lambdaS")){
REAL8 lambda1=0.;
REAL8 lambda2=0.;
LALInferenceBinaryLove(model->params, &lambda1, &lambda2);
XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
LALInferenceAddVariable(model->params, "lambda1", &lambda1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
LALInferenceAddVariable(model->params, "lambda2", &lambda2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
}
/* Only use GR templates */
/* Fill in the extra parameters for testing GR, if necessary */
for (UINT4 k=0; k<N_extra_params; k++)
......@@ -894,6 +905,17 @@ void LALInferenceTemplateXLALSimInspiralChooseWaveform(LALInferenceModel *model)
}
/* ==== BINARY_LOVE PARAMETERS ==== */
if(LALInferenceCheckVariable(model->params, "lambdaS")){
REAL8 lambda1=0.;
REAL8 lambda2=0.;
LALInferenceBinaryLove(model->params, &lambda1, &lambda2);
XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
LALInferenceAddVariable(model->params, "lambda1", &lambda1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
LALInferenceAddVariable(model->params, "lambda2", &lambda2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
}
/* Only use GR templates */
/* Fill in the extra parameters for testing GR, if necessary */
for (UINT4 k=0; k<N_extra_params; k++)
......@@ -1373,6 +1395,16 @@ void LALInferenceTemplateXLALSimInspiralChooseWaveformPhaseInterpolated(LALInfer
XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
}
/* ==== BINARY_LOVE PARAMETERS ==== */
if(LALInferenceCheckVariable(model->params, "lambdaS")){
REAL8 lambda1=0.;
REAL8 lambda2=0.;
LALInferenceBinaryLove(model->params, &lambda1, &lambda2);
XLALSimInspiralWaveformParamsInsertTidalLambda1(model->LALpars, lambda1);
XLALSimInspiralWaveformParamsInsertTidalLambda2(model->LALpars, lambda2);
LALInferenceAddVariable(model->params, "lambda1", &lambda1, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
LALInferenceAddVariable(model->params, "lambda2", &lambda2, LALINFERENCE_REAL8_t, LALINFERENCE_PARAM_OUTPUT);
}
/* ==== 4-PIECE POLYTROPE EOS PARAMETERS ==== */
REAL8 logp1 = 0.;
......
......@@ -349,6 +349,9 @@ adapt-temps=
#gamma2 gamma2 [1.1,4.5]
#gamma3 gamma3 [1.1,4.5]
#Common EoS tidal model for BNS systems (requires BinaryLove=)
#lambdaS Symmetric tidal deformability, (lambda1+lambda2)/2 [0,5000]
# Settings for allowed component masses in Solar Masses, with default values
#comp-max=30.0
#comp-min=1.0
......
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