Commit 1d5515f1 authored by Marcella Wijngaarden's avatar Marcella Wijngaarden

Compute CBC response after burnin so that MCMC can start at signal/glitch updates

- The CBC templates held in the models after the burnin were still zero, now
  they hold the actual burnin templates so that the MCMC can start with
  CBC templates subtracted from glitch/signal residuals

- Also moved the dtimes convention conversion to IN the projection routine
  so that we don't have to repeat it everywhere. The swapTimesFlag tells
  the routine which convention to use (=1 when called from BW, =0 when
  called internally from BayesCBC routines).
parent 6ba0c6b9
......@@ -3570,7 +3570,12 @@ void PDwave(double *wavef, RealVector *freq, double *params, int N, struct bayes
}
void projectCBCWaveform(double *ampphase, int N, int NI, double fmin, double Tobs, double *extParams, double **response, double *dtimes, double *Fplus, double *Fcross)
/*
* Projects the CBC model amplitudate and phase onto the detectors.
* Uses the CBC convention for timedelays, unless swapTimesFlag is used.
* `response` holds the computed response for each detector.
*/
void projectCBCWaveform(double *ampphase, int N, int NI, double fmin, double Tobs, double *extParams, double **response, double *dtimes, double *Fplus, double *Fcross, int swapTimesFlag)
{
int i,id;
double Fs, lambda;
......@@ -3589,8 +3594,10 @@ void projectCBCWaveform(double *ampphase, int N, int NI, double fmin, double Tob
{
Fs = sqrt(Ap*Ap*Fplus[id]*Fplus[id]+Ac*Ac*Fcross[id]*Fcross[id]); // magnitude of response
lambda = atan2(Ac*Fcross[id],Ap*Fplus[id]); // get phase from polarizations
// td = -1*dtimes[id];
td = dtimes[id];
// When calling with BW/LAL timedelays, swap sign convention to CBC internal
if (swapTimesFlag == 1) td = -1*dtimes[id];
else td = dtimes[id];
if(lambda < 0.0) lambda += TPI;
......@@ -3615,6 +3622,12 @@ void projectCBCWaveform(double *ampphase, int N, int NI, double fmin, double Tob
}
}
/*
* Computes the CBC model intrinsic amplitudes and phase at geocenter.
* `gwave` holds the computed amplitude (at 2*i indices) and phase (at
* 2*i+1 indices).
*/
void geowave(double *gwave, RealVector *freq, double *params, int N, struct bayesCBC *rundata)
{
......
......@@ -116,7 +116,8 @@ void set_bayescbc_priors(struct Net *net, struct bayesCBC *rundata);
void BayesCBC_MCMC(int M, int N, double Tobs, double ttrig, gsl_rng *r, int istart, struct bayesCBC *rundata, int ic, double *amphase);
void MCMC_all_wrapper(int M, int N, double Tobs, double ttrig, gsl_rng *r, int istart, struct bayesCBC *rundata);
void projectCBCWaveform(double *AmpPhase, int N, int NI, double fmin, double Tobs, double *extParams, double **response, double *dtimes, double *Fplus, double *Fcross);
void geowave(double *gwave, RealVector *freq, double *params, int N, struct bayesCBC *rundata);
void projectCBCWaveform(double *ampphase, int N, int NI, double fmin, double Tobs, double *extParams, double **response, double *dtimes, double *Fplus, double *Fcross, int swapTimesFlag);
void print_projected_cbc_waveform(double **SN, double Tobs, double ttrig, double **projected, double **data, int N, int iter, struct bayesCBC *rundata);
double z_DL(double DL);
......
......@@ -1466,7 +1466,7 @@ double log_likelihood_full(struct Net *net, double **D, double *params, RealVect
// Compute template with current parameters
geowave(ampphase, freq, params, N, rundata);
projectCBCWaveform(ampphase, N, net->Nifo, fmin, Tobs, sky, twave, net->dtimes, net->Fplus, net->Fcross);
projectCBCWaveform(ampphase, N, net->Nifo, fmin, Tobs, sky, twave, net->dtimes, net->Fplus, net->Fcross, 0);
logL = 0.0;
for(id = 0; id < net->Nifo; id++)
......
......@@ -107,7 +107,6 @@ double fb_nwip(double *a, double *b, int n, int imin, int imax);
void cossin(struct Net *net, double *hcos, double *hsin, RealVector *freq, double *params, int N);
void extrinsictemplates(struct Net *net, double **hwave, RealVector *freq, double *hcos, double *hsin, double *params, int N);
void geowave(double *gwave, RealVector *freq, double *params, int N, struct bayesCBC *rundata);
void fulltemplates(struct Net *net, double **hwave, RealVector *freq, double *params, int N, struct bayesCBC *rundata);
void fullphaseamp(struct Net *net, double **amp, double **phase, RealVector *freq, double *params, int N, struct bayesCBC *rundata);
void Fisher_All(struct Net *net, double **fish, double *params, RealVector *freq, double **SN, int N, double Tobs, struct bayesCBC *rundata);
......
......@@ -887,10 +887,7 @@ double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, doubl
if(data->cbcFlag)
{
// for(i=0; i<NI; i++) for(n=0; n<N; n++) c[i][n] = 0.0;
dtimes = double_vector(NI);
for(ifo=0; ifo<NI; ifo++) dtimes[ifo] = -1*projection->dtimes[ifo];
projectCBCWaveform(amphase, N, NI, fmin, data->Tobs, params, c, dtimes, projection->Fplus, projection->Fcross);
free_double_vector(dtimes);
projectCBCWaveform(amphase, N, NI, fmin, data->Tobs, params, c, projection->dtimes, projection->Fplus, projection->Fcross, 1);
}
//Form up residual
......@@ -952,10 +949,7 @@ double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *
//compute instrument response from cbc waveform at geocenter with extrinsic params
if(data->cbcFlag)
{
dtimes = double_vector(NI);
for(ifo=0; ifo<NI; ifo++) dtimes[ifo] = -1*projection->dtimes[ifo];
projectCBCWaveform(amphase, N, NI, fmin, data->Tobs, params, c, dtimes, projection->Fplus, projection->Fcross);
free_double_vector(dtimes);
projectCBCWaveform(amphase, N, NI, fmin, data->Tobs, params, c, projection->dtimes, projection->Fplus, projection->Fcross, 1);
}
//Form up residual
......@@ -1002,10 +996,7 @@ double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *
// compute instrument response from cbc waveform at geocenter with extrinsic params
if(data->cbcFlag)
{
dtimes = double_vector(NI);
for(ifo=0; ifo<NI; ifo++) dtimes[ifo] = -1*projection->dtimes[ifo];
projectCBCWaveform(amphase, N, NI, fmin, data->Tobs, params, c, dtimes, projection->Fplus, projection->Fcross);
free_double_vector(dtimes);
projectCBCWaveform(amphase, N, NI, fmin, data->Tobs, params, c, projection->dtimes, projection->Fplus, projection->Fcross, 1);
}
......
......@@ -1304,9 +1304,10 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
Shf_Geocenter_full(data, model[chain->index[ic]]->projection, model[chain->index[ic]]->Snf, model[chain->index[ic]]->SnGeo, model[chain->index[ic]]->extParams);
}
// 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);
// Store CBC templates in model after burnin
geowave(model[chain->index[ic]]->cbcamphase, bayescbc->freq, model[chain->index[ic]]->cbcParams, data->N, bayescbc);
projectCBCWaveform(model[chain->index[ic]]->cbcamphase, data->N, data->NI, data->fmin, data->Tobs, model[chain->index[ic]]->extParams, model[chain->index[ic]]->cbctemplate, model[chain->index[ic]]->projection->dtimes, model[chain->index[ic]]->projection->Fplus, model[chain->index[ic]]->projection->Fcross, 1);
}
// Only evolve extrinsic params from now on
......@@ -1357,7 +1358,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
if((data->glitchFlag || data->signalFlag) && !data->fixIntrinsicFlag) EvolveIntrinsicParameters(data, prior, model, chain, tf, chain->seed, ic);
//This function executes [cycle] extrinsic parameter updates, common to all geocenter wavelets.
if((model[chain->index[ic]]->signal[0]->size>0 || data->cbcFlag) && !data->fixExtrinsicFlag && chain->mc > 1) EvolveExtrinsicParameters(data, prior, model, chain, chain->seed, ic);
if((model[chain->index[ic]]->signal[0]->size>0 || data->cbcFlag) && !data->fixExtrinsicFlag) EvolveExtrinsicParameters(data, prior, model, chain, chain->seed, ic);
//update PSD with BayesLine
if(data->bayesLineFlag) EvolveBayesLineParameters(data, model, bayesline, chain, prior, ic);
......@@ -2651,13 +2652,8 @@ void EvolveBayesCBCParameters(struct Data *data, struct Model **model, struct ba
model_x->cbcParams[i] = bayescbc->pallx[icbc][i];
}
dtimes = double_vector(NI);
for(ifo=0; ifo<NI; ifo++) dtimes[ifo] = -1 * projection->dtimes[ifo];
// Get CBC waveforms projected on the detectors
projectCBCWaveform(model_x->cbcamphase, N, NI, fmin, data->Tobs, model_x->extParams, model_x->cbctemplate, dtimes, projection->Fplus, projection->Fcross);
free_double_vector(dtimes);
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);
......@@ -2931,11 +2927,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if (data->cbcFlag)
{
double *dtimes;
dtimes = double_vector(NI);
for(ifo=0; ifo<NI; ifo++) dtimes[ifo] = -1 * model_x->projection->dtimes[ifo];
projectCBCWaveform(model_x->cbcamphase, data->N, NI, data->fmin, data->Tobs, model_x->extParams, model_x->cbctemplate, dtimes, model_x->projection->Fplus, model_x->projection->Fcross);
free_double_vector(dtimes);
projectCBCWaveform(model_x->cbcamphase, data->N, NI, data->fmin, data->Tobs, model_x->extParams, model_x->cbctemplate, model_x->projection->dtimes, model_x->projection->Fplus, model_x->projection->Fcross, 1);
}
......
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