Commit 654af9b3 authored by Marcella Wijngaarden's avatar Marcella Wijngaarden

Add extinsic mcmc in CBC burnin + use same template functions as BW

parent f8c37969
This diff is collapsed.
......@@ -85,6 +85,7 @@ double fbegin(double *param);
void pmap(struct Net *net, double *pallx, double *paramx, double *sky, int NX, double gmst);
void pmap_back(struct Net *net, double *pall, double *param, double *sky, int NX, double gmst);
void pmap_all(struct Net *net, double *pall, double *param, double *sky, int NX, double gmst);
void printwave(struct Net *net, int N, RealVector *freq, double **paramx, int *who, double Tobs, double ttrig, int iter, struct bayesCBC *rundata);
void printwaveall(struct Net *net, int N, RealVector *freq, double *paramx, double **SN, double Tobs, double ttrig, int iter, struct bayesCBC *rundata);
......@@ -182,7 +183,7 @@ void SNRvsf(struct Net *net, double **D, double *params, RealVector *freq, doubl
double globeden(double ***global, double *max, double *min, double Tobs, double *params, int N);
double globe(double ***global, double *max, double *min, double Tobs, double *params, int N, gsl_rng *r, int lmax, struct bayesCBC *rundata);
void dshifts(struct Net *net, double *sky, double *params, int NX, double gmst);
void dshifts(struct Net *net, double *sky, double *params, int NX, struct bayesCBC *rundata);
void CBC_start(struct Net *net, int *mxc, FILE *chainI, double **paramx, double **skyx, double **pallx, double *heat, double ***history, double ***global, RealVector *freq, double **D, double **SN, int N, double Tobs, gsl_rng *r, struct bayesCBC *rundata);
......
......@@ -262,21 +262,10 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
}
}
if(data->signalFlag)
if(data->signalFlag || data->cbcModelFlag)
{
// TEMPORARY ADDITION MARCELLA TO SET INITAL EXTRINISC FIXED FOR GW150914
if(data->fixExtrinsicFlag == 1)
{
model[ic]->extParams[2] = 1.5000;
model[ic]->extParams[3] = -0.975609756097561;
printf("IMPORTANT: FIXED INITIAL EXTRINISC PARAMS TO CHOSEN VALUES\n");
}
else
{
extrinsic_uniform_proposal(chain->seed,model[ic]->extParams);
}
extrinsic_uniform_proposal(chain->seed,model[ic]->extParams);
if(data->skyFixFlag == 1)
{
......@@ -298,7 +287,11 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
computeProjectionCoeffs(data, model[ic]->projection, model[ic]->extParams, data->fmin, data->fmax);
Shf_Geocenter_full(data, model[ic]->projection, model[ic]->Snf, model[ic]->SnGeo, model[ic]->extParams);
}
if(data->signalFlag)
{
for(n=0; n<model[ic]->Npol; n++)
{
signal[n]->size = 1;
......@@ -1255,6 +1248,8 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
setup_bayescbc_model(bayescbc, data, chain, model[chain->index[ic]], ic);
}
bayescbc->intrinsic_only = 0;
// Quick search to find initial signal
fprintf(stdout, " ======= Start CBC_burnin() ======\n");
bayescbc->chainS = fopen("chains/cbc_searchchain.dat","w");
......@@ -1268,6 +1263,33 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
{
model[chain->index[ic]]->cbcParams[i] = bayescbc->pallx[chain->index[ic]][i];
}
if (bayescbc->intrinsic_only == 0)
{
printf("updating model with burnin extrinisc params \n");
// Pass extrinsic parameters from BW model to BayesCBC
for(i = bayescbc->NX; i < bayescbc->NP; i++)
{
// Convert Bayeswaves' ellipticity to BayesCBC's cos(inclination)
if (i == bayescbc->NX+3)
{
double ciota = bayescbc->pallx[chain->index[ic]][bayescbc->NX+3];
double epsilon =2*ciota/(1+ciota*ciota);
model[chain->index[ic]]->extParams[i-bayescbc->NX] = epsilon;
}
else
{
model[chain->index[ic]]->extParams[i-bayescbc->NX] = bayescbc->pallx[chain->index[ic]][i];
}
}
// Update projection based on extrinsic params from burnin
computeProjectionCoeffs(data, model[chain->index[ic]]->projection, model[chain->index[ic]]->extParams, data->fmin, data->fmax);
Shf_Geocenter_full(data, model[chain->index[ic]]->projection, model[chain->index[ic]]->Snf, model[chain->index[ic]]->SnGeo, model[chain->index[ic]]->extParams);
}
printf("model chain %i = index[%i]\n", ic, chain->index[ic]);
// 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);
......@@ -2608,6 +2630,7 @@ 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);
// Save waveform temporarily for checking residuals + template
if (ic==0 && bayescbc->debug == 1) 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
......
......@@ -1278,7 +1278,6 @@ void setup_bayescbc_model(struct bayesCBC *bayescbc, struct Data *data, struct C
bayescbc->heat[ic] = chain->temperature[ic];
bayescbc->mxc[0] = chain->mc;
bayescbc->who = chain->index;
}
void initialize_network(struct Network *projection, int N, int NI)
......
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