Commit 58e15430 authored by Marcella Wijngaarden's avatar Marcella Wijngaarden

Update code structure BayesCBC setup

The number of CBC parameters is set based on the cmd options (so
data->NX is always initialized).

Move BayesCBC initialization to main BayesWave program and add
some print settings.

Move setup of BayesCBC that needs to be done before all independent
updates (to make sure BayesCBC is up to date about the latest
model parameters, projections, SN etc.) to subroutine.

Add optional tidal version cmd option (to determine how many semi-intrinsic
parameters are needed).
parent 183ccc51
Pipeline #159458 failed with stages
in 1 minute and 45 seconds
......@@ -75,7 +75,8 @@ struct bayesCBC
int NQ; // number of mass ratios in global proposal
int NM; // number of chirp masses in global proposal
int constantLogLFlag; // fix the likelihood to a constant
int constantLogLFlag; // fix the likelihood to a constant
int injectionStartFlag; // start at injection values (instead of burnin)
double *max; // Array that will hold upper boundary prior range
double *min; // Array that will hold lower boundary prior range
......
......@@ -21,8 +21,6 @@ static double maxarg1, maxarg2;
#define SQR(x) ((x) * (x))
#define CUBE(x) ((x) * (x) * (x))
#define FOURTH(x) ((x) * (x) * (x) * (x))
#define STRINGIFY(x) #x
#define TOSTRING(x) STRINGIFY(x)
#define TINY 1.0e-10
#define IM1 2147483563
......
......@@ -614,6 +614,27 @@ int main(int argc, char *argv[])
fprintf(stdout,"\n");
}
/******************************************************************************/
/* */
/* Setup BayesCBC */
/* */
/******************************************************************************/
struct bayesCBC *bayescbc = malloc(sizeof(struct bayesCBC));
if (data->bayesCBCFlag)
{
// Initialize baeyscbc with default values & allocate memory
initialize_bayescbc(bayescbc, data, chain, model[chain->index[0]]);
// Save BayesCBC run settings
sprintf(filename,"%sbayeswave.run",data->runName);
chain->runFile = fopen(filename,"a");
print_cbc_run_stats(chain->runFile, data, bayescbc);
fclose(chain->runFile);
}
/******************************************************************************/
/* */
/* The Full-time MCMC Run (Cleaning phase) */
......@@ -673,7 +694,7 @@ int main(int argc, char *argv[])
//TODO disable checkpointing and adaptation during cleaning phase
int saveAdaptationFlag = data->adaptTemperatureFlag;
data->adaptTemperatureFlag = 0;
RJMCMC(data, model, bayesline, chain, prior, &cleanLogZ, &cleanVarZ);
RJMCMC(data, model, bayesline, bayescbc, chain, prior, &cleanLogZ, &cleanVarZ);
data->adaptTemperatureFlag = saveAdaptationFlag;
/******************************************************************************/
......@@ -887,14 +908,6 @@ int main(int argc, char *argv[])
}
/******************************************************************************/
/* */
/* Now all data/PSDs are in their final form, test calling BayesCBC */
/* */
/******************************************************************************/
// Currently all BayesCBC initialization is moved to RJMCMC
/******************************************************************************/
......@@ -944,7 +957,7 @@ int main(int argc, char *argv[])
printf("characterizing signal model...\n\n");
RJMCMC(data, model, bayesline, chain, prior, &data->logZsignal, &data->varZsignal);
RJMCMC(data, model, bayesline, bayescbc, chain, prior, &data->logZsignal, &data->varZsignal);
//shut off signal model for checkpointing
data->signalModelFlag = 0;
......@@ -986,7 +999,7 @@ int main(int argc, char *argv[])
printf("characterizing glitch model...\n\n");
RJMCMC(data, model, bayesline, chain, prior, &data->logZglitch, &data->varZglitch);
RJMCMC(data, model, bayesline, bayescbc, chain, prior, &data->logZglitch, &data->varZglitch);
//shut off glitch model for checkpointing
data->glitchModelFlag = 0;
......@@ -1022,7 +1035,7 @@ int main(int argc, char *argv[])
printf("characterizing Gaussian noise model...\n\n");
RJMCMC(data, model, bayesline, chain, prior, &data->logZnoise, &data->varZnoise);
RJMCMC(data, model, bayesline, bayescbc, chain, prior, &data->logZnoise, &data->varZnoise);
//shut off noise model for checkpointing
data->noiseModelFlag = 0;
......@@ -1062,7 +1075,7 @@ int main(int argc, char *argv[])
reset_priors(data, prior);
reset_model(data, chain, prior, model);
RJMCMC(data, model, bayesline, chain, prior, &data->logZfull, &data->varZfull);
RJMCMC(data, model, bayesline, bayescbc, chain, prior, &data->logZfull, &data->varZfull);
fprintf(evidence,"full %.12g %lg\n",data->logZfull,data->varZfull);
fflush(evidence);
......@@ -1185,10 +1198,12 @@ void print_help_message(void)
fprintf(stdout," --noiseOnly use noise model only (no signal or glitches)\n");
fprintf(stdout," --signalOnly use signal model only (no glitches)\n");
fprintf(stdout," --glitchOnly use glitch model only (no signal)\n");
fprintf(stdout," --cbcOnly use full CBC model only (incl. extrinsic parameters)\n");
fprintf(stdout," --noPSDfit keep PSD parameters fixed\n");
fprintf(stdout," --bayesLine use BayesLine for PSD model\n");
fprintf(stdout," --stochastic use stochastic background model\n");
fprintf(stdout," --bayesCBC use BayesCBC to model inspiral\n");
fprintf(stdout," --bayesCBC-tidal specify BayesCBC waveform used: NRTidal_V or NoNRT_V (default)) (\n");
fprintf(stdout,"\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --- Priors & Proposals -----------------------------------------------------------\n");
......
......@@ -956,6 +956,33 @@ void print_run_flags(FILE *fptr, struct Data *data, struct Prior *prior)
fprintf(fptr, "\n");
}
void print_cbc_run_stats(FILE *fptr, struct Data *data, struct bayesCBC *bayescbc)
{
fprintf(fptr, " ============ BayesCBC Settings ===========\n");
fprintf(fptr, " nr. semi-intrinsic params (NX) %i\n", bayescbc->NX);
fprintf(fptr, " tidal version ");
if(bayescbc->NRTidal_version == NoNRT_V) fprintf(fptr,"NoNRT_V");
else if(bayescbc->NRTidal_version == NRTidal_V) fprintf(fptr,"NRTidal_V");
else fprintf(fptr, "UNKNOWN");
fprintf(fptr, "\n");
fprintf(fptr, " evolve extrinsic/sky in burnin ");
if(bayescbc->intrinsic_only == 0) fprintf(fptr,"ENABLED");
else fprintf(fptr, "DISABLED");
fprintf(fptr, "\n");
fprintf(fptr, " debug mode (--debug-cbc) ");
if(bayescbc->debug) fprintf(fptr,"ENABLED");
else fprintf(fptr, "DISABLED");
fprintf(fptr, "\n");
fprintf(fptr, " all chain output (--verbose) ");
if(bayescbc->verbose) fprintf(fptr,"ENABLED");
else fprintf(fptr, "DISABLED");
fprintf(fptr,"\n");
}
void print_chain_files(struct Data *data, struct Chain *chain, struct Model **model, struct BayesLineParams ***bayesline, struct bayesCBC *bayescbc, int ic)
{
int n;
......@@ -1565,6 +1592,8 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
data->fullModelFlag = 0;
data->cbcOnlyFlag = 1;
data->bayesCBCFlag = 1;
// TODO: This should also make fixIntrinsicParameters required
}
ppt = LALInferenceGetProcParamVal(commandLine, "--fullOnly");
......@@ -1683,6 +1712,17 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
if(ppt) data->bayesCBCFlag = 1;
else data->bayesCBCFlag = 0;
ppt = LALInferenceGetProcParamVal(commandLine, "--bayesCBC-tidal");
data->NX = 7;
if(ppt)
{
char command[1024];
sprintf(command,"%s",ppt->value);
printf("'%s'\n", command);
if (strstr(command, "NoNRT_V") != NULL) data->NX = 7;
else data->NX = 9;
}
ppt = LALInferenceGetProcParamVal(commandLine, "--0noise");
if(ppt) data->noiseSimFlag = 0;
else data->noiseSimFlag = 1;
......
......@@ -52,6 +52,7 @@ void print_bayesline_spectrum(char filename[], double *f, double *power, double
void print_version(FILE *fptr);
void print_run_stats(FILE *fptr, struct Data *data, struct Chain *chain);
void print_cbc_run_stats(FILE *fptr, struct Data *data, struct bayesCBC *bayescbc);
void print_run_flags(FILE *fptr, struct Data *data, struct Prior *prior);
void print_chain_files(struct Data *data, struct Chain *chain, struct Model **model, struct BayesLineParams ***bayesline, struct bayesCBC *bayescbc, int ic);
void print_chain_status(struct Data *data, struct Chain *chain, struct Model **model, int searchFlag);
......
......@@ -1020,11 +1020,12 @@ static void save_sampler(struct Data *data, struct Chain *chain, struct Model **
/* */
/* ********************************************************************************** */
void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***bayesline, struct Chain *chain, struct Prior *prior, double *logEvidence, double *varLogEvidence)
void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***bayesline, struct bayesCBC *bayescbc, struct Chain *chain, struct Prior *prior, double *logEvidence, double *varLogEvidence)
{
print_run_stats(stdout, data, chain);
print_run_flags(stdout, data, prior);
if (data->bayesCBCFlag) print_cbc_run_stats(stdout, data, bayescbc);
int N = data->N;
int NI = data->NI;
......@@ -1244,119 +1245,30 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
/* */
/******************************************************************************/
// Set BayesCBC rundata (this replaces the header file settings of QuickCBC)
struct bayesCBC *bayescbc = malloc(sizeof(struct bayesCBC));
// BayesCBC model independent burnin
if (data->runPhase && (data->signalFlag||data->glitchFlag||data->fullModelFlag) && data->bayesCBCFlag)
{
double x, y;
printf("Setting up BayesCBC in RJMCMC\n");
// Set BayesCBC rundata (this replaces the header file settings of QuickCBC)
initialize_bayescbc(bayescbc, data, chain, model[chain->index[0]]);
if(LALInferenceGetProcParamVal(data->commandLine, "--debug-cbc")) {
bayescbc->debug = 1;
}
// Use XML injection parameter file from BayesWave command line
int start_at_injection = 0;
if(LALInferenceGetProcParamVal(data->commandLine, "--inj"))
{
sprintf(filename, LALInferenceGetProcParamVal(data->commandLine,"--inj")->value);
start_at_injection = 1;
// If injecting tidal parameters, read values from command line and setup version + sampler type
if(LALInferenceGetProcParamVal(data->commandLine, "--inj-lambda1")||LALInferenceGetProcParamVal(data->commandLine, "--inj-lambdaT"))
{
bayescbc->NRTidal_version = (NRTidal_version_type) NRTidal_V; // TODO: Set from XML
// Increase parameter numbers
bayescbc->NX = bayescbc->NX+2;
bayescbc->NP = bayescbc->NX+4;
bayescbc->cap = 100.0; // cap on log likelihood proposal (sets SNR^2/2 threshold)
bayescbc->ladder = 1.05; // temperature ladder
// Set tidal sampler type (in lambda1/2 directly or in lambdaT/dlambdaT) and injection values
if(LALInferenceGetProcParamVal(data->commandLine,"--inj-lambda1")&&LALInferenceGetProcParamVal(data->commandLine,"--inj-dLambda2"))
{
bayescbc->lambda_type_version = (Lambda_type) lambda;
bayescbc->init_tidal1 = atof(LALInferenceGetProcParamVal(data->commandLine, "--inj-lambda1")->value);
bayescbc->init_tidal2 = atof(LALInferenceGetProcParamVal(data->commandLine, "--inj-lambda2")->value);
}
else if(LALInferenceGetProcParamVal(data->commandLine,"--inj-lambdaT")&&LALInferenceGetProcParamVal(data->commandLine,"--inj-dLambdaT"))
{
bayescbc->lambda_type_version = (Lambda_type) lambdaTilde;
bayescbc->init_tidal1 = atof(LALInferenceGetProcParamVal(data->commandLine, "--inj-lambdaT")->value);
bayescbc->init_tidal2 = atof(LALInferenceGetProcParamVal(data->commandLine, "--inj-dLambdaT")->value);
}
else
{
fprintf(stdout, "Error: Injection needs two tidal parameters (i.e., lambda1 & lambda2 or lambdaT & dlambdaT.\nEnding.\n");
return ;
}
fprintf(stdout, " Tidal injection parameters: %f, %f\n", bayescbc->init_tidal1, bayescbc->init_tidal2);
}
}
data->NX = bayescbc->NX;
// Copy over chain seeds to temporary bayescbc structure
// This is not actually used because we are not running in parallel at the moment
bayescbc->rvec = (gsl_rng **)malloc(sizeof(gsl_rng *) * (chain->NC+1));
for(i=0; i<= chain->NC; i++)
{
// bayescbc->rvec[i] = gsl_rng_alloc(T);
bayescbc->rvec[chain->index[i]] = chain->rvec[chain->index[i]];
}
// Setup antenna patterns and time delays for current extrinsic params
for(ifo = 0; ifo<data->NI; ifo++)
{
// Time offsets for detectors
bayescbc->net->dtimes[ifo] = model[chain->index[0]]->projection->dtimes[ifo];
// Antenna patterns for detectors
bayescbc->net->Fplus[ifo] = model[chain->index[0]]->projection->Fplus[ifo];
bayescbc->net->Fcross[ifo] = model[chain->index[0]]->projection->Fcross[ifo];
}
// Pass extrinsic params from BW to CBC
// Pass the most recent model parameters to BayesCBC structure
for(ic=0; ic<chain->NC; ic++)
{
model[chain->index[i]]->cbcParams = double_vector(bayescbc->NX);
for(i = bayescbc->NX; i < bayescbc->NP; i++)
{
// Convert Bayeswaves' ellipticity to BayesCBC's cos(inclination)
if (i == bayescbc->NX+3)
{
double epsilon = model[chain->index[0]]->extParams[i-bayescbc->NX];
double ciota = (1.0 - sqrt(1.0-(epsilon*epsilon)))/epsilon;
bayescbc->pallx[chain->index[ic]][i] = ciota;
}
else
{
bayescbc->pallx[chain->index[ic]][i] = model[chain->index[0]]->extParams[i-bayescbc->NX];
}
}
setup_bayescbc_model(bayescbc, data, chain, model[chain->index[ic]], ic);
}
// Quick search to find initial signal
fprintf(stdout, " ======= Start CBC_burnin() ======\n");
bayescbc->chainS = fopen("chains/cbc_searchchain.dat","w");
burnin_bayescbc(data->N, data->Tobs, data->trigtime, bayescbc->D, bayescbc->SN, chain->seed, filename, start_at_injection, bayescbc);
burnin_bayescbc(data->N, data->Tobs, data->trigtime, bayescbc->D, bayescbc->SN, chain->seed, filename, 0, bayescbc);
fclose(bayescbc->chainS);
// Pass initial signal parameters to model
// Pass updated cbc parameters to model
for(ic=0; ic<chain->NC; ic++)
{
for(i = 0; i < bayescbc->NX; i++)
{
model[chain->index[ic]]->cbcParams[i] = bayescbc->pallx[chain->index[ic]][i];
}
// TODO: Intrinsic template here
// cbcGeoAmpPhase(bayescbc, chain->index[ic], data->N, model[chain->index[ic]]->cbcamphase);
// waveformProject(data, model[chain->index[ic]]->projection, model[chain->index[ic]]->extParams, model[chain->index[ic]]->cbctemplate, model[chain->index[ic]]->cbcamphase, data->fmin, data->fmax);
}
......@@ -2681,67 +2593,8 @@ void EvolveBayesCBCParameters(struct Data *data, struct Model **model, struct ba
}
}
// Use latest estimates of PSD if using BayesLine
// if (data->bayesLineFlag)
if (true)
{
for(ifo=0; ifo<NI; ifo++)
{
for (i=0; i<data->N/2; ++i)
{
if(i>imin && i<imax)
{
bayescbc->SN[ifo][i] = model_x->Snf[ifo][i];
}
else
{
bayescbc->SN[ifo][i] = 1.0;
}
}
}
}
// Setup antenna patterns and time delays for current extrinsic params
scale = model_x->extParams[5];
for(ifo = 0; ifo<data->NI; ifo++)
{
// Time offsets for detectors
bayescbc->net->dtimes[ifo] = projection->dtimes[ifo];
// Antenna patterns for detectors
bayescbc->net->Fplus[ifo] = projection->Fplus[ifo];
bayescbc->net->Fcross[ifo] = projection->Fcross[ifo];
}
// Pass intrinsic CBC params from BW to CBC
for(i = 0; i < bayescbc->NX; i++)
{
bayescbc->pallx[icbc][i] = model_x->cbcParams[i];
}
// Pass extrinsic params from BW to CBC
for(i = bayescbc->NX; i <= bayescbc->NP; i++)
{
// Convert Bayeswaves' ellipticity to BayesCBC's cos(inclination)
if (i == bayescbc->NX+3)
{
epsilon = model_x->extParams[i-bayescbc->NX];
ciota = (1.0 - sqrt(1.0-(epsilon*epsilon)))/epsilon;
bayescbc->pallx[icbc][i] = ciota;
}
else
{
bayescbc->pallx[icbc][i] = model_x->extParams[i-bayescbc->NX];
}
}
// Set heat
bayescbc->heat[ic] = chain->temperature[ic];
bayescbc->mxc[0] = chain->mc;
bayescbc->who = chain->index;
// if (bayescbc->debug==1) printf("before BayesCBC_MCMC chain %i, %i, model_x->logL = %f, logLnorm = %f, logLx = %f \n", chain->index[ic], ic, model_x->logL, model_x->logLnorm, bayescbc->logLx[icbc]);
// Pass the updated model parameters to BayesCBC structure
setup_bayescbc_model(bayescbc, data, chain, model_x, ic);
// Do updates
BayesCBC_MCMC(M, N, data->Tobs, data->trigtime, chain->seed, istart, bayescbc, ic, model_x->cbcamphase);
......@@ -2786,8 +2639,6 @@ void EvolveBayesCBCParameters(struct Data *data, struct Model **model, struct ba
}
}
// if (bayescbc->debug==1) printf("after BayesCBC_MCMC chain %i, %i, model_x->logL = %f, logLnorm = %f, logLx = %f \n", chain->index[ic], ic, model_x->logL, model_x->logLnorm, bayescbc->logLx[icbc]);
return;
}
......
......@@ -23,7 +23,7 @@
/* */
/* ********************************************************************************** */
void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***bayesline, struct Chain *chain, struct Prior *prior, double *logEvidence, double *varLogEvidence);
void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***bayesline, struct bayesCBC *bayescbc,struct Chain *chain, struct Prior *prior, double *logEvidence, double *varLogEvidence);
void PTMCMC(struct Model **model, struct Chain *chain, gsl_rng *seed, int NC);
void EvolveBayesLineParameters(struct Data *data, struct Model **model, struct BayesLineParams ***bayesline, struct Chain *chain, struct Prior *prior, int ic);
......
......@@ -710,7 +710,7 @@ void initialize_model(struct Model *model, struct Data *data, struct Prior *prio
int NI = data->NI;
int Npol = data->Npol;
int NW = data->NW; //number of intrinsic parameters / frame (5 for wavelet)
int NX = data->NX;
int NX = data->NX; //number of semi-intrinsic CBC parameters (masses etc.)
//gaussian noise model
//model->Sn = double_vector(NI-1);
......@@ -721,7 +721,6 @@ void initialize_model(struct Model *model, struct Data *data, struct Prior *prio
//instrument cbc template model
model->cbctemplate = double_matrix(NI-1,N-1);
// model->cbcgeotemplate = double_vector(N-1);
model->cbcamphase = double_vector(N-1);
//store sum of 1/psd in 1/SnGeo
......@@ -999,13 +998,15 @@ void free_bayescbc(struct bayesCBC *bayescbc, struct Data *data, struct Chain *c
fclose(bayescbc->chainA);
}
void initialize_bayescbc(struct bayesCBC *bayescbc, struct Data *data, struct Chain *chain, struct Model *model)
{
int i, ifo;
ProcessParamsTable *ppt = NULL;
int i, ifo, ic;
int imin, imax;
double x, y;
double fac;
char filename[1024];
char command[1024];
int N = data->N;
int NI = data->NI;
......@@ -1019,8 +1020,7 @@ void initialize_bayescbc(struct bayesCBC *bayescbc, struct Data *data, struct Ch
bayescbc->constantLogLFlag = data->constantLogLFlag;
bayescbc->NX = 7; // number of quasi-intrinsic parameters (Mc, Mt, chi1, chi2, phic, tc, DL)
data->NX = bayescbc->NX;
bayescbc->NX = data->NX; // number of quasi-intrinsic parameters (Mc, Mt, chi1, chi2, phic, tc, DL)
bayescbc->NP = bayescbc->NX+4; // total number of entries in the full parameter arrays NX+4
bayescbc->NS = 7; // number of quasi-extrinsic parameters (alpha, sin(delta), psi, ellipticity, scale, phi0=2*phic, dt)
bayescbc->NC = chain->NC; // number of chains
......@@ -1029,7 +1029,6 @@ void initialize_bayescbc(struct bayesCBC *bayescbc, struct Data *data, struct Ch
bayescbc->NH = 1000; // length of history
bayescbc->NQ = 4; // number of mass ratios in global proposal
bayescbc->NM = 1000; // number of chirp masses in global proposal
bayescbc->NRTidal_version = (NRTidal_version_type) NoNRT_V;
bayescbc->cap = 10.0; // cap on log likelihood proposal (sets SNR^2/2 threshold)
bayescbc->qfac = 2.0; // geometric progression in mass ratio
......@@ -1050,6 +1049,21 @@ void initialize_bayescbc(struct bayesCBC *bayescbc, struct Data *data, struct Ch
for (i = 0; i < 3; i++) bayescbc->mxc[i] = 0;
bayescbc->pallx = double_matrix(bayescbc->NC,bayescbc->NP);
// Parse/set tidal version
ppt = LALInferenceGetProcParamVal(data->commandLine, "--bayesCBC-tidal");
if(ppt)
{
sprintf(command,"%s",ppt->value);
if (strstr(command, "NoNRT_V") != NULL) bayescbc->NRTidal_version = (NRTidal_version_type) NoNRT_V;
else if (strstr(command, "NRTidal_V") != NULL) bayescbc->NRTidal_version = (NRTidal_version_type) NRTidal_V;
else
{
printf(" BayesCBC only supports tidal version NoNRT_V or NRTidal_V at the moment. Exiting! \n");
exit(1);
}
}
else bayescbc->NRTidal_version = (NRTidal_version_type) NoNRT_V;
// Transfer Network struct to CBC Net
initialize_net(bayescbc->net, data->Tobs, data->trigtime, NI+3, data->ifos);
......@@ -1078,6 +1092,7 @@ void initialize_bayescbc(struct bayesCBC *bayescbc, struct Data *data, struct Ch
}
}
// TODO: Can remove this bayesline option after debugging
if (data->bayesLineFlag)
{
// Use Bayesline PSD if available (same estimation method as BayesCBC)
......@@ -1088,10 +1103,10 @@ void initialize_bayescbc(struct bayesCBC *bayescbc, struct Data *data, struct Ch
// Commented out alternative option for PSD that is not smoothed.
// But smooth seems to work fine (and already has the correct scale)!
// sprintf(filename, "PSD_%d_%s.dat", (int)(Tobs), data->ifos[ifo]);
sprintf(filename, "%s_burnin_psd.dat", data->ifos[ifo]);
fprintf(stdout, "Using BayesLine PSD for BayesCBC from file: %s\n", filename);
sprintf(command, "%s_burnin_psd.dat", data->ifos[ifo]);
fprintf(stdout, "Using BayesLine PSD for BayesCBC from file: %s\n", command);
fptr = fopen(filename,"r");
fptr = fopen(command,"r");
bayescbc->SN[ifo][0] = 1.0; // BayesLine PSD starts at i=1 not 0
for (i=1; i<N/2; ++i)
{
......@@ -1134,18 +1149,138 @@ void initialize_bayescbc(struct bayesCBC *bayescbc, struct Data *data, struct Ch
// all chains in the file, otherwise only cold chain)
bayescbc->chainA = fopen("chains/cbc_allchain.dat","w");
bayescbc->debug = 0;
// Debug output mode
if(LALInferenceGetProcParamVal(data->commandLine, "--debug-cbc")) bayescbc->debug = 1;
else bayescbc->debug = 0;
// Verbose mode (for printing all chains and cbc waveforms)
if(data->verboseFlag) bayescbc->verbose = 1;
else bayescbc->verbose = 0;
if(data->verboseFlag)
// Set run setting when NS in system
if (bayescbc->NRTidal_version != NoNRT_V)
{
bayescbc->verbose = 1;
bayescbc->cap = 100.0; // cap on log likelihood proposal (sets SNR^2/2 threshold)
// bayescbc->ladder = 1.05; // temperature ladder
}
else
// Handle injection settings (TODO!)
if(LALInferenceGetProcParamVal(data->commandLine, "--inj"))
{
bayescbc->injectionStartFlag = 1;
sprintf(command, LALInferenceGetProcParamVal(data->commandLine,"--inj")->value);
// If injecting tidal parameters, read values from command line and setup version + sampler type
if(LALInferenceGetProcParamVal(data->commandLine, "--inj-lambda1")||LALInferenceGetProcParamVal(data->commandLine, "--inj-lambdaT"))
{
// Set tidal sampler type (in lambda1/2 directly or in lambdaT/dlambdaT) and injection values
if(LALInferenceGetProcParamVal(data->commandLine,"--inj-lambda1")&&LALInferenceGetProcParamVal(data->commandLine,"--inj-dLambda2"))
{
bayescbc->lambda_type_version = (Lambda_type) lambda;
bayescbc->init_tidal1 = atof(LALInferenceGetProcParamVal(data->commandLine, "--inj-lambda1")->value);
bayescbc->init_tidal2 = atof(LALInferenceGetProcParamVal(data->commandLine, "--inj-lambda2")->value);
}
else if(LALInferenceGetProcParamVal(data->commandLine,"--inj-lambdaT")&&LALInferenceGetProcParamVal(data->commandLine,"--inj-dLambdaT"))
{
bayescbc->lambda_type_version = (Lambda_type) lambdaTilde;
bayescbc->init_tidal1 = atof(LALInferenceGetProcParamVal(data->commandLine, "--inj-lambdaT")->value);
bayescbc->init_tidal2 = atof(LALInferenceGetProcParamVal(data->commandLine, "--inj-dLambdaT")->value);
}
else
{
fprintf(stdout, "Error: Injection needs two tidal parameters (i.e., lambda1 & lambda2 or lambdaT & dlambdaT.\nEnding.\n");
return ;
}
fprintf(stdout, " Tidal injection parameters: %f, %f\n", bayescbc->init_tidal1, bayescbc->init_tidal2);
}
}
else
{
bayescbc->verbose = 0;
bayescbc->injectionStartFlag = 0;
}
// Copy over chain seeds to temporary bayescbc structure
// This is not actually used because we are not running in parallel at the moment
bayescbc->rvec = (gsl_rng **)malloc(sizeof(gsl_rng *) * (chain->NC+1));
for(i=0; i<= chain->NC; i++)
{
// bayescbc->rvec[i] = gsl_rng_alloc(T);
bayescbc->rvec[chain->index[i]] = chain->rvec[chain->index[i]];
}
}
void setup_bayescbc_model(struct bayesCBC *bayescbc, struct Data *data, struct Chain *chain, struct Model *model, int ic)
{
int ifo, i;
double epsilon, ciota;
int icbc = chain->index[ic];
// This routine is called before any independent bayesCBC MCMC iteration
// run (both burnin and EvolveBayesCBCParameters) and sets up the BayesCBC
// structure with the latest model parameters/(run)data
// Setup antenna patterns and time delays for initial extrinsic params
for(ifo = 0; ifo<data->NI; ifo++)
{
// Time offsets for detectors
bayescbc->net->dtimes[ifo] = model->projection->dtimes[ifo];
// Antenna patterns for detectors
bayescbc->net->Fplus[ifo] = model->projection->Fplus[ifo];