Commit 2176551c authored by Marcella Wijngaarden's avatar Marcella Wijngaarden

Remove old BayesCBC init from main BW

parent 284a7ca9
......@@ -894,129 +894,7 @@ int main(int argc, char *argv[])
/* */
/******************************************************************************/
if(data->bayesCBCFlag && LALInferenceGetProcParamVal(runState->commandLine, "--bayesCBC-only"))
{
fprintf(stdout, " ======= Start BayesCBC() ======\n");
// These need to be passed back and forth between BayesWave and the CBC PE code
// The CBC PE code also needs the PSDs and glitch removed residuals D
double x, y;
FILE *BWData; // Temporarily for debugging
// BayesCBC control default parameters
int start_at_injection = 0;
// Set BayesCBC rundata (this replaces the header file settings of QuickCBC)
struct bayesCBC *bayescbc = malloc(sizeof(struct bayesCBC));
initialize_bayescbc(bayescbc, data, chain, model[0]);
// Use XML injection parameter file from BayesWave command line
if(LALInferenceGetProcParamVal(runState->commandLine, "--inj"))
{
sprintf(filename, LALInferenceGetProcParamVal(runState->commandLine,"--inj")->value);
start_at_injection = 1;
// If injecting tidal parameters, read values from command line and setup version + sampler type
if(LALInferenceGetProcParamVal(runState->commandLine, "--inj-lambda1")||LALInferenceGetProcParamVal(runState->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(runState->commandLine,"--inj-lambda1")&&LALInferenceGetProcParamVal(runState->commandLine,"--inj-dLambda2"))
{
bayescbc->lambda_type_version = (Lambda_type) lambda;
bayescbc->init_tidal1 = atof(LALInferenceGetProcParamVal(runState->commandLine, "--inj-lambda1")->value);
bayescbc->init_tidal2 = atof(LALInferenceGetProcParamVal(runState->commandLine, "--inj-lambda2")->value);
}
else if(LALInferenceGetProcParamVal(runState->commandLine,"--inj-lambdaT")&&LALInferenceGetProcParamVal(runState->commandLine,"--inj-dLambdaT"))
{
bayescbc->lambda_type_version = (Lambda_type) lambdaTilde;
bayescbc->init_tidal1 = atof(LALInferenceGetProcParamVal(runState->commandLine, "--inj-lambdaT")->value);
bayescbc->init_tidal2 = atof(LALInferenceGetProcParamVal(runState->commandLine, "--inj-dLambdaT")->value);
}
else
{
fprintf(stdout, "Error: Injection needs two tidal parameters (i.e., lambda1 & lambda2 or lambdaT & dlambdaT.\nEnding.\n");
return 1;
}
fprintf(stdout, " Tidal injection parameters: %f, %f\n", bayescbc->init_tidal1, bayescbc->init_tidal2);
}
}
// Temporarily for debugging
int dump_cbc_data = 1;
if (dump_cbc_data)
{
// output data
BWData = fopen("dataD.dat","w");
for (j = 0; j < N; ++j)
{
fprintf(BWData,"%e ", (double)(j)/Tobs);
for (i = 0; i < bayescbc->net->Nifo; ++i) fprintf(BWData,"%e ", bayescbc->D[i][j]);
fprintf(BWData,"\n");
}
fclose(BWData);
// output psd data
BWData = fopen("dataPSD.dat","w");
for (j = 0; j < N/2; ++j)
{
fprintf(BWData,"%e ", (double)(j)/Tobs);
for (i = 0; i < NI; ++i) fprintf(BWData,"%e ", bayescbc->SN[i][j]);
fprintf(BWData,"\n");
}
fclose(BWData);
}
// Set & initialize seed chains (temporarily, shoud use the version from Chain structure)
gsl_rng_set(chain->seed, 18346443564);
bayescbc->rvec = (gsl_rng **)malloc(sizeof(gsl_rng *) * (bayescbc->NC));
for(i=0; i<= bayescbc->NC; i++)
{
bayescbc->rvec[i] = gsl_rng_alloc(T);
gsl_rng_set(bayescbc->rvec[i] , gsl_rng_uniform(chain->seed));
}
fprintf(stdout, " ======= Start CBC_initialize() ======\n");
bayescbc->chainS = fopen("chains/cbc_searchchain.dat","w");
// Quick search to find signal
burnin_bayescbc(N, Tobs, data->trigtime, bayescbc->D, bayescbc->SN, chain->seed, filename, start_at_injection, bayescbc);
fclose(bayescbc->chainS);
// Temporary flag for if we want to test a seperate run of BayesCBC without
// doing the BayesWave MCMC.
if(LALInferenceGetProcParamVal(runState->commandLine, "--bayesCBC-only"))
{
fprintf(stdout, " ======= Start MCMC_all() ======\n");
bayescbc->chainA = fopen("chains/cbc_allchain.dat","w");
// This whole function can be removed, as it's just intended for testing here to run
// BayesCBC without running full Bayeswave
// MCMC_all_wrapper((1500000/bayescbc->NCC), N, Tobs, data->trigtime, chain->seed, 0, bayescbc);
// Previous version (which updates all chains in parallel): To test internal updates in BayesCBC
MCMC_all(bayescbc->net, bayescbc->mxc, (1500000/bayescbc->NCC), bayescbc->chainA, bayescbc->pallx, bayescbc->who, bayescbc->heat, bayescbc->historyall, bayescbc->global, bayescbc->freq, bayescbc->D, bayescbc->SN, N, Tobs, data->trigtime, chain->seed, 0, bayescbc);
fclose(bayescbc->chainA);
for(i=0 ;i<= bayescbc->NC; i++){
gsl_rng_free(bayescbc->rvec[i]);
}
free(bayescbc->rvec);
fprintf(stdout, " ======= End BayesCBC() ======\n");
return 0;
}
}
// Currently all BayesCBC initialization is moved to RJMCMC
/******************************************************************************/
......
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