Commit 445684d8 authored by Marcella Wijngaarden's avatar Marcella Wijngaarden

Print tidal param chains & fix burnin ring updates

Also removed some printed warnings in the internal waveform
generation functions.
parent 1d5515f1
Pipeline #167112 failed with stages
in 45 minutes and 33 seconds
......@@ -2365,8 +2365,9 @@ void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, dou
#pragma omp parallel for
for(k=0; k < NC; k++)
{
updatei(k, net, lmax, logLx, paramx, paramy, min, max, Fscale, who, heat, history, global, freq, D, SN, ejump, evec, N, Tobs, cv, av, rundata->rvec[k], rundata);
updatei(k, net, lmax, logLx, paramx, paramy, min, max, Fscale, who, heat, history, global, freq, D, SN, ejump, evec, N, Tobs, cv, av, r, rundata);
}
// add to the history file
if(mc%20 == 0)
{
......@@ -2388,7 +2389,7 @@ void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, dou
}
if(mc%10 == 0)
if(mc%100 == 0)
{
q = who[0];
Mchirp = exp(paramx[q][0])/MSUN_SI;
......@@ -3525,7 +3526,6 @@ void PDwave(double *wavef, RealVector *freq, double *params, int N, struct bayes
lambda_type_version = rundata->lambda_type_version;
if (NRTidal_version != NoNRT_V) {
// Compute lambda1 & lambda2
if (lambda_type_version == lambdaTilde) {
lambdaT = params[7];
......@@ -3537,12 +3537,12 @@ void PDwave(double *wavef, RealVector *freq, double *params, int N, struct bayes
extraParams[1] = lambda2; // lambda2
} else {
extraParams[0] = params[7]; // lambda1
extraParams[1] = params[8]; // lambda2
extraParams[1] = params[8]; // lambda2
}
}
ret = IMRPhenomDGenerateh22FDAmpPhase(&ap, freq, phi0, fRef_in, m1_SI, m2_SI, chi1, chi2, distance, extraParams, NRTidal_version);
for (i=1; i< N/2; i++)
{
wavef[2*i] = 0.0;
......
......@@ -122,6 +122,7 @@ void print_projected_cbc_waveform(double **SN, double Tobs, double ttrig, double
double z_DL(double DL);
double DL_fit(double z);
void LambdaTsEta2Lambdas(double lambdaT, double dLambdaT, double eta, double *lambda1, double *lambda2);
// Helper functions for extrinsic burnin
......
......@@ -1222,6 +1222,7 @@ void Ring(struct Net *net, double *skyx, double *skyy, int d1, int d2, gsl_rng *
double calpha, salpha;
double H1mag, L1mag, zmag;
double cmu, smu;
double gmst_rad;
H1 = (double*)malloc(sizeof(double)* 3);
L1 = (double*)malloc(sizeof(double)* 3);
......@@ -1298,11 +1299,20 @@ void Ring(struct Net *net, double *skyx, double *skyy, int d1, int d2, gsl_rng *
// Unet->Nifot vector between the sites
for(i = 0; i <3; i++) zv[i] /= zmag;
//remap gmst back to [0:2pi]
double intpart;
double decpart;
gmst_rad = gmst/LAL_TWOPI;
intpart = (int)( gmst_rad );
decpart = gmst_rad - (double)intpart;
gmst_rad = decpart*LAL_TWOPI;
alpha = skyx[0];
sind = skyx[1];
cosd = sqrt(1.0-sind*sind);
calpha = cos(alpha-gmst);
salpha = -sin(alpha-gmst);
calpha = cos(alpha-gmst_rad);
salpha = -sin(alpha-gmst_rad);
// Unusual convention on spherical angle
// propagation direction
......@@ -1340,7 +1350,7 @@ void Ring(struct Net *net, double *skyx, double *skyy, int d1, int d2, gsl_rng *
skyy[1] = kvr[2];
rot = atan2(kvr[1],kvr[0]);
rot = gmst - rot;
rot = gmst_rad - rot;
if(rot < 0.0) rot += TPI;
skyy[0] = rot;
......
......@@ -171,8 +171,6 @@ void updatei(int k, struct Net *net, int lmax, double *logLx, double **paramx, d
void update_chain(int k, struct Net *net, double *logLx, double *rhox, double *paramx, double *paramy, double *min, double *max, double heat, double **history, double ***global, RealVector *freq, double **D, double **SN, double *ejump, double **evec, int N, double Tobs, int *cv, int *av, gsl_rng *r, int *evolveparams, struct bayesCBC *rundata);
void LambdaTsEta2Lambdas(double lambdaT, double dLambdaT, double eta, double *lambda1, double *lambda2);
void ParseInjectionXMLParameters(char *file_path, char **injection_columns, char **injection_values);
void GetXMLInjectionParameter(char *name, double *value, char **injection_columns, char **injection_values);
\ No newline at end of file
......@@ -1134,6 +1134,22 @@ void print_cbc_model(FILE *fptr, struct bayesCBC *bayescbc, int ic, double Tobs)
fprintf(fptr,"%f %f %f %f ", bayescbc->pallx[ic][bayescbc->NX], bayescbc->pallx[ic][bayescbc->NX+1], bayescbc->pallx[ic][bayescbc->NX+2], bayescbc->pallx[ic][bayescbc->NX+3]);//ra, sindec, psi, cosiota
fprintf(fptr,"%f %f %f %f ", m1, m2, Mchirp, Mtot); //mass parameters in source frame
// Add tidal parameters
double lambda1, lambda2, lambdaT, dLambdaT;
if (bayescbc->NX > 7)
{
if (bayescbc->lambda_type_version == lambdaTilde) {
lambdaT = bayescbc->pallx[ic][7];
dLambdaT = bayescbc->pallx[ic][8];
LambdaTsEta2Lambdas(lambdaT, dLambdaT, eta, &lambda1, &lambda2);
fprintf(fptr,"%f %f %f %f ", lambda1, lambda2, lambdaT, dLambdaT); //lambda parameters
} else {
lambda1 = bayescbc->pallx[ic][7];
lambda2 = bayescbc->pallx[ic][8];
fprintf(fptr,"%f %f ", lambda1, lambda2); //lambda parameters
}
}
fprintf(fptr,"\n");
// The below is for recreating the QuickCBC allchain file units
......@@ -1142,40 +1158,14 @@ void print_cbc_model(FILE *fptr, struct bayesCBC *bayescbc, int ic, double Tobs)
if (bayescbc->verbose==1)
{
// TODO: Add suport below for converting & printing neutron star tidal parameters
// if (NX > 7) {
// if (lambda_type_version == lambdaTilde) {
// lambdaT = bayescbc->pallx[7];
// dLambdaT = bayescbc->pallx[8];
// LambdaTsEta2Lambdas(lambdaT, dLambdaT, eta, &lambda1, &lambda2);
// } else {
// lambda1 = bayescbc->pallx[7];
// lambda2 = bayescbc->pallx[8];
// }
// fprintf(chain,"%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e\t %e\t %e\t %e\t %e\t %e\t %e\t %e\t", mxc[2], rundata->logLx[ic], \
// Mchirp, Mtot, chieff, paramx[4], \
// tc, DL, paramx[NX], paramx[NX+1], \
// paramx[NX+2], paramx[NX+3], z, Mchirp/(1.0+z), Mtot/(1.0+z), m1/(1.0+z), m2/(1.0+z), m1, m2,
// lambda1, lambda2, lambdaT, dLambdaT);
// } else {
// Merger time is printed relative to trigger time, which is at Tobs-2.
// [0] log(Mc) [1] log(Mt) [2] chi1 [3] chi2 [4] phi0 [5] tp [6] log(DL) [7] alpha [8] sindelta [9] psi [10] ciota
// counter, log likelihood, chirp mass (DF), total mass (DF), effective spin, geocenter orbital phase, geocenter arrival time, distance, RA , sine of DEC,
// polarization angle, cos inclination, redshift, chirp mass (SF), total mass (SF), m1 (SF), m2 (SF), m1 (DF), m2 (DF), arrival time at each detector, SNR at each detector
fprintf(bayescbc->chainA,"%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e", bayescbc->mxc[2], bayescbc->logLx[ic], Mchirp, Mtot, chieff, bayescbc->pallx[ic][4], \
fprintf(bayescbc->chainA,"%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e", bayescbc->mxc[2], bayescbc->logLx[ic], Mchirp, Mtot, chieff, bayescbc->pallx[ic][4], \
tc, DL, bayescbc->pallx[ic][bayescbc->NX], bayescbc->pallx[ic][bayescbc->NX+1], \
bayescbc->pallx[ic][bayescbc->NX+2], bayescbc->pallx[ic][bayescbc->NX+3], z, Mchirp/(1.0+z), Mtot/(1.0+z), m1/(1.0+z), m2/(1.0+z), m1, m2, bayescbc->net->Fplus[0], bayescbc->net->Fcross[0]);
bayescbc->pallx[ic][bayescbc->NX+2], bayescbc->pallx[ic][bayescbc->NX+3], z, Mchirp/(1.0+z), Mtot/(1.0+z), m1/(1.0+z), m2/(1.0+z), m1, m2);
if (bayescbc->NX > 7) fprintf(bayescbc->chainA," %e %e ", lambda1, lambda2); //lambda parameters
fprintf(bayescbc->chainA,"\n");
// printf("%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e \n", bayescbc->mxc[2], bayescbc->logLx[ic], Mchirp, Mtot, chieff, bayescbc->pallx[ic][4], \
// tc, DL, bayescbc->pallx[ic][bayescbc->NX], bayescbc->pallx[ic][bayescbc->NX+1], \
// bayescbc->pallx[ic][bayescbc->NX+2], bayescbc->pallx[ic][bayescbc->NX+3], z, Mchirp/(1.0+z), Mtot/(1.0+z), m1/(1.0+z), m2/(1.0+z), m1, m2);
}
bayescbc->mxc[2] += 1;
......
......@@ -2655,8 +2655,8 @@ void EvolveBayesCBCParameters(struct Data *data, struct Model **model, struct ba
// Get CBC waveforms projected on the detectors
projectCBCWaveform(model_x->cbcamphase, N, NI, fmin, data->Tobs, model_x->extParams, model_x->cbctemplate, projection->dtimes, projection->Fplus, projection->Fcross, 1);
// Save 200 waveforms temporarily for checking residuals + template
if (ic==0 && bayescbc->debug == 1 && chain->mc<=200) print_projected_cbc_waveform(bayescbc->SN, data->Tobs, data->trigtime, model_x->cbctemplate, bayescbc->D, N, bayescbc->mxc[0], bayescbc);
// Save first 10 waveforms temporarily for checking residuals + template
if (ic==0 && bayescbc->debug == 1 && chain->mc<=1000) print_projected_cbc_waveform(bayescbc->SN, data->Tobs, data->trigtime, model_x->cbctemplate, bayescbc->D, N, bayescbc->mxc[0], bayescbc);
// Recompute likelihoods of current chain
model_x->logLnorm = 0.0;
......@@ -3060,7 +3060,7 @@ void extrinsic_burnin_cbc(struct bayesCBC *bayescbc, struct Data *data, struct M
if(bayescbc->debug == 1) printf("%f %f\n", SNRsq[0], SNRsq[1]);
// find a sky location roughly consistent with the time delays
skyring(net, paramx[j], skyx[j], pallx[j], SNRsq, bayescbc->freq, D, SN, N, data->Tobs, bayescbc->rvec[j], NX, NS, bayescbc->gmst);
skyring(net, paramx[j], skyx[j], pallx[j], SNRsq, bayescbc->freq, D, SN, N, data->Tobs, chain->seed, NX, NS, bayescbc->gmst);
dshifts(net, skyx[j], paramx[j], NX, bayescbc);
......
......@@ -1196,7 +1196,6 @@ void initialize_bayescbc(struct bayesCBC *bayescbc, struct Data *data, struct Ch
if (bayescbc->NRTidal_version != NoNRT_V)
{
bayescbc->cap = 100.0; // cap on log likelihood proposal (sets SNR^2/2 threshold)
// bayescbc->ladder = 1.05; // temperature ladder
}
// Reference frequency
......
......@@ -153,8 +153,8 @@ int IMRPhenomDGenerateFD(
const double q = (m1 > m2) ? (m1 / m2) : (m2 / m1);
if (q > MAX_ALLOWED_MASS_RATIO)
PRINT_WARNING("Warning: The model is not supported for high mass ratio, see MAX_ALLOWED_MASS_RATIO\n");
// if (q > MAX_ALLOWED_MASS_RATIO)
// PRINT_WARNING("Warning: The model is not supported for high mass ratio, see MAX_ALLOWED_MASS_RATIO\n");
if (chi1 > 1.0 || chi1 < -1.0 || chi2 > 1.0 || chi2 < -1.0)
ERROR(PD_EDOM, "Spins outside the range [-1,1] are not supported\n");
......@@ -238,8 +238,8 @@ int IMRPhenomDGenerateh22FDAmpPhase(
const double q = (m1 > m2) ? (m1 / m2) : (m2 / m1);
if (q > MAX_ALLOWED_MASS_RATIO)
PRINT_WARNING("Warning: The model is not supported for high mass ratio, see MAX_ALLOWED_MASS_RATIO\n");
// if (q > MAX_ALLOWED_MASS_RATIO)
// PRINT_WARNING("Warning: The model is not supported for high mass ratio, see MAX_ALLOWED_MASS_RATIO\n");
if (chi1 > 1.0 || chi1 < -1.0 || chi2 > 1.0 || chi2 < -1.0)
ERROR(PD_EDOM, "Spins outside the range [-1,1] are not supported\n");
......
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