Commit 91aafad6 authored by Marcella Wijngaarden's avatar Marcella Wijngaarden

Merge branch 'BW+QuickCBC' into 'BW+QuickCBC'

Flags

See merge request marcella.wijngaarden/bayeswave!5
parents 654af9b3 df4c413b
......@@ -1012,7 +1012,7 @@ void print_chain_files(struct Data *data, struct Chain *chain, struct Model **mo
}
//print bayescbc parameters
if(data->cbcModelFlag)
if(data->cbcFlag)
{
print_cbc_model(chain->cbcChainFile[ic], bayescbc, chain->index[ic], data->Tobs);
}
......@@ -1045,7 +1045,7 @@ void flush_chain_files(struct Data *data, struct Chain *chain, int ic)
}
//flush cbc file
if(data->cbcModelFlag) fflush(chain->cbcChainFile[ic]);
if(data->cbcFlag) fflush(chain->cbcChainFile[ic]);
}
void print_model(FILE *fptr, struct Data *data, struct Chain *chain, struct Model *model)
......
......@@ -197,7 +197,7 @@ void recompute_residual(struct Data *data, struct Model **model, struct Chain *c
data->r[ifo][i] = data->s[ifo][i];
if(data->glitchFlag) data->r[ifo][i] -= model_x->glitch[ifo]->templates[i];
if(data->signalFlag) data->r[ifo][i] -= model_x->response[ifo][i];
if(data->cbcModelFlag) data->r[ifo][i] -= model_x->cbctemplate[ifo][i];
if(data->cbcFlag) data->r[ifo][i] -= model_x->cbctemplate[ifo][i];
}
}
......@@ -390,7 +390,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
}
}
if(data->cbcModelFlag)
if(data->cbcFlag)
{
//subtract current CBC model from residual
for(i=0; i<N; i++)
......@@ -678,7 +678,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
r[ifo][j] -= ry[ifo][j];
r[ifo][k] -= ry[ifo][k];
}
if(data->cbcModelFlag)
if(data->cbcFlag)
{
r[ifo][j] -= my->cbctemplate[ifo][j];
r[ifo][k] -= my->cbctemplate[ifo][k];
......@@ -704,7 +704,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
r[det][j] = s[det][j] - ry[det][j] - gy[det][j];
r[det][k] = s[det][k] - ry[det][k] - gy[det][k];
if (data->cbcModelFlag)
if (data->cbcFlag)
{
r[det][j] -= my->cbctemplate[det][j];
r[det][k] -= my->cbctemplate[det][k];
......@@ -760,7 +760,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
r[ifo][j] -= ry[ifo][j];
r[ifo][k] -= ry[ifo][k];
}
if(data->cbcModelFlag)
if(data->cbcFlag)
{
r[ifo][j] -= my->cbctemplate[ifo][j];
r[ifo][k] -= my->cbctemplate[ifo][k];
......@@ -788,7 +788,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
r[det][j] = s[det][j] - ry[det][j] - gy[det][j];
r[det][k] = s[det][k] - ry[det][k] - gy[det][k];
if (data->cbcModelFlag)
if (data->cbcFlag)
{
r[det][j] -= my->cbctemplate[det][j];
r[det][k] -= my->cbctemplate[det][k];
......@@ -815,7 +815,7 @@ double EvaluateMarkovianLogLikelihood(int typ, int ii, int det, struct Model *mx
{
r[ifo][i] -= ry[ifo][i];
}
if(data->cbcModelFlag)
if(data->cbcFlag)
{
r[ifo][i] -= my->cbctemplate[ifo][i];
}
......@@ -883,7 +883,7 @@ double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, doubl
//compute instrument response from cbc waveform at geocenter with extrinsic params
if(data->cbcModelFlag)
if(data->cbcFlag)
{
// for(i=0; i<NI; i++) for(n=0; n<N; n++) c[i][n] = 0.0;
projectCBCWaveform(amphase, N, NI, fmin, data->Tobs, params, c, projection->dtimes, projection->Fplus, projection->Fcross);
......@@ -897,7 +897,7 @@ double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, doubl
r[ifo][n] = d[ifo][n];
if(data->signalFlag) r[ifo][n] -= h[ifo][n];
if(data->glitchFlag) r[ifo][n] -= g[ifo][n];
if(data->cbcModelFlag) r[ifo][n] -= c[ifo][n];
if(data->cbcFlag) r[ifo][n] -= c[ifo][n];
}
logL += loglike(imin, imax, r[ifo], invSnf[ifo]);
}
......@@ -945,7 +945,7 @@ double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *
waveformProject(data, projection, params, h, t, data->fmin, data->fmax);
//compute instrument response from cbc waveform at geocenter with extrinsic params
if(data->cbcModelFlag)
if(data->cbcFlag)
{
projectCBCWaveform(amphase, N, NI, fmin, data->Tobs, params, c, projection->dtimes, projection->Fplus, projection->Fcross);
}
......@@ -958,7 +958,7 @@ double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *
r[ifo][n] = d[ifo][n];
if(data->signalFlag) r[ifo][n] -= h[ifo][n];
if(data->glitchFlag) r[ifo][n] -= g[ifo][n];
if(data->cbcModelFlag) r[ifo][n] -= c[ifo][n];
if(data->cbcFlag) r[ifo][n] -= c[ifo][n];
}
}
......@@ -992,7 +992,7 @@ double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *
waveformProject(data, projection, params, h, t, data->fmin, data->fmax);
// compute instrument response from cbc waveform at geocenter with extrinsic params
if(data->cbcModelFlag)
if(data->cbcFlag)
{
projectCBCWaveform(amphase, N, NI, fmin, data->Tobs, params, c, projection->dtimes, projection->Fplus, projection->Fcross);
}
......@@ -1005,7 +1005,7 @@ double EvaluateExtrinsicSearchLogLikelihood(struct Network *projection, double *
r[ifo][n] = d[ifo][n];
if(data->signalFlag) r[ifo][n] -= h[ifo][n];
if(data->glitchFlag) r[ifo][n] -= g[ifo][n];
if(data->cbcModelFlag) r[ifo][n] -= c[ifo][n];
if(data->cbcFlag) r[ifo][n] -= c[ifo][n];
}
}
logL = 0.0;
......
......@@ -339,7 +339,7 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
data->r[ifo][i] = data->s[ifo][i];
if(data->signalFlag) data->r[ifo][i] -= model[ic]->response[ifo][i];
if(data->glitchFlag) data->r[ifo][i] -= glitch[ifo]->templates[i];
if(data->cbcModelFlag) data->r[ifo][i] -= model[ic]->cbctemplate[ifo][i];
if(data->cbcFlag) data->r[ifo][i] -= model[ic]->cbctemplate[ifo][i];
}
if(!data->constantLogLFlag)
{
......@@ -400,7 +400,7 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
}
//Parameters for reproducing CBC model
if(data->cbcModelFlag)
if(data->cbcFlag)
{
sprintf(filename,"chains/cbc_model.dat.%i",ic);
chain->cbcChainFile[ic] = fopen(filename,"w");
......@@ -1018,7 +1018,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
print_run_stats(stdout, data, chain);
print_run_flags(stdout, data, prior);
if (data->cbcModelFlag) print_cbc_run_stats(stdout, data, bayescbc);
if (data->cbcFlag) print_cbc_run_stats(stdout, data, bayescbc);
int N = data->N;
int NI = data->NI;
......@@ -1239,7 +1239,7 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
/******************************************************************************/
// BayesCBC model independent burnin
if (data->runPhase && data->cbcModelFlag)
if (data->runPhase && data->cbcFlag)
{
// Pass the most recent model parameters to BayesCBC structure
......@@ -1343,13 +1343,13 @@ 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->cbcModelFlag) && !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 && chain->mc > 1) EvolveExtrinsicParameters(data, prior, model, chain, chain->seed, ic);
//update PSD with BayesLine
if(data->bayesLineFlag) EvolveBayesLineParameters(data, model, bayesline, chain, prior, ic);
// Evolve BayesCBC
if(data->cbcModelFlag) EvolveBayesCBCParameters(data, model, bayescbc, chain, ic);
if(data->cbcFlag) EvolveBayesCBCParameters(data, model, bayescbc, chain, ic);
}
......@@ -1666,7 +1666,7 @@ void EvolveBayesLineParameters(struct Data *data, struct Model **model, struct B
data->r[ifo][i] = data->s[ifo][i];
if(data->signalFlag) data->r[ifo][i] -= model_x->response[ifo][i];
if(data->glitchFlag) data->r[ifo][i] -= model_x->glitch[ifo]->templates[i];
if(data->cbcModelFlag) data->r[ifo][i] -= model_x->cbctemplate[ifo][i];
if(data->cbcFlag) data->r[ifo][i] -= model_x->cbctemplate[ifo][i];
}
//re-run Markovian, full spectrum, full model part of BayesLine
......@@ -1687,7 +1687,7 @@ void EvolveBayesLineParameters(struct Data *data, struct Model **model, struct B
data->r[ifo][i] = data->s[ifo][i];
if(data->signalFlag) data->r[ifo][i] -= model_x->response[ifo][i];
if(data->glitchFlag) data->r[ifo][i] -= model_x->glitch[ifo]->templates[i];
if(data->cbcModelFlag) data->r[ifo][i] -= model_x->cbctemplate[ifo][i];
if(data->cbcFlag) data->r[ifo][i] -= model_x->cbctemplate[ifo][i];
}
model_x->detLogL[ifo] = loglike(data->imin, data->imax, data->r[ifo], model_x->invSnf[ifo]);
......@@ -1857,7 +1857,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
acc=1;
if(data->signalFlag) copy_int_model(model_x, model_y, N, NI,-1);
if(data->glitchFlag) for(ifo=0; ifo<NI; ifo++) copy_int_model(model_x, model_y, N, NI, ifo);
if(data->cbcModelFlag) copy_cbc_model(model_x, model_y, N, NI, data->NX);
if(data->cbcFlag) copy_cbc_model(model_x, model_y, N, NI, data->NX);
for(mc=1; mc<=M; mc++)
{
......@@ -2750,7 +2750,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
Shf_Geocenter_full(data, model_x->projection, model_x->Snf, model_x->SnGeo, model_x->extParams);
for(i=1; i<ienddim; i++)
{
if (!data->cbcModelFlag) logpx += (data->signal_amplitude_prior(smodel[0]->intParams[smodel[0]->index[i]],SnGeox, data->Tobs, prior->sSNRpeak));
if (data->signalFlag) logpx += (data->signal_amplitude_prior(smodel[0]->intParams[smodel[0]->index[i]],SnGeox, data->Tobs, prior->sSNRpeak));
}
}
//logpx += dim*log(paramsx[5]);
......@@ -2823,7 +2823,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
for(i=1; i<ienddim; i++)
{
intParams[i][3] = smodel[0]->intParams[smodel[0]->index[i]][3]*paramsy[5];//amplitude
if (!data->cbcModelFlag) test += checkrange(intParams[i],prior->range, smodel[0]->dimension);
if (data->signalFlag) test += checkrange(intParams[i],prior->range, smodel[0]->dimension);
}
//check extrinsic parameters
......@@ -2838,7 +2838,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
{
//model_x->projection has be updated by data->extrinsic_likelihood, so it stores F+ and Fx for params_y.
Shf_Geocenter_full(data, model_x->projection, model_x->Snf, SnGeoy, paramsy);
if (!data->cbcModelFlag) for(i=1; i<ienddim; i++) logpy += (data->signal_amplitude_prior(intParams[i],SnGeoy, data->Tobs, prior->sSNRpeak));
if (data->signalFlag) for(i=1; i<ienddim; i++) logpy += (data->signal_amplitude_prior(intParams[i],SnGeoy, data->Tobs, prior->sSNRpeak));
}
logH = (logLy - logLx)*chain->beta + logpy - logpx + logJ;
......
......@@ -843,7 +843,7 @@ void initialize_model(struct Model *model, struct Data *data, struct Prior *prio
}
//Initial cbc template must be zero
if(data->cbcModelFlag)
if(data->cbcFlag)
{
for (ifo=0; ifo<NI; ifo++)
{
......
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