Commit ad49aa52 authored by John Douglas Veitch's avatar John Douglas Veitch
Browse files

Merge branch 'AddBinaryLove' into 'master'

Add binary love EoS model

See merge request !790
parents fdd2a533 6117be18
Pipeline #65417 passed with stages
in 161 minutes and 43 seconds
......@@ -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