Skip to content

GitLab

Commit 37ea10d6 authored by Marcella Wijngaarden's avatar Marcella Wijngaarden
Browse files
Options

First version seperated extrinsic burnin 

Moved extrinsic burnin helper functions to BayesCBCExtrinsic.c
which uses LAL functions for time delays etc.

The extrinsic burnin is seperated in BayesWaveMCMC. todo: decide
where this function should live permantly.
parent a878d5bd
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 316 additions and 1579 deletions
src/BayesCBC.c
View file @ 37ea10d6
This diff is collapsed.
src/BayesCBC.h
View file @ 37ea10d6
......@@ -22,9 +22,14 @@ struct Net
int *lstart;
int *lstop;
double **location; // The three components, in an Earth-fixed Cartesian coordinate system, of the position vector from the center of the Earth to the detector in meters
double ***response; // The Earth-fixed Cartesian components of the detector's response tensor \f$d^{ab}\f$
// These are only used locally and updated per iteration
double *Fplus; // antenna pattern functions
double *Fcross; // antenna pattern functions
double *dtimes; //time delay between reference location and IFOs
};
// Structure to hold BayesCBC run settings (that were previously hardcoded in QuickCBC header)
......@@ -124,6 +129,33 @@ void print_projected_cbc_waveform(double **SN, double Tobs, double ttrig, double
double z_DL(double DL);
double DL_fit(double z);
// Helper functions for extrinsic burnin
void skymcmc(struct Net *net, int MCX, int *mxc, FILE *chain, double **paramx, double **skyx, double **pallx, int *who, double *heat, double dtx, int nt, double *DD, double **WW, double ***DHc, double ***DHs, double ***HH, double Tobs, gsl_rng * r, struct bayesCBC *rundata);
// Templates
void templates(struct Net *net, double **hwave, RealVector *freq, double *params, int N, struct bayesCBC *rundata);
void geotemplate(double *gwave, RealVector *freq, double *params, int N, struct bayesCBC *rundata);
void skyring(struct Net *net, double *params, double *sky, double *pall, double *SNRsq, RealVector *freq, double **D, double **SN, int N, double Tobs, gsl_rng * r, int NX, int NS, double gmst);
void ringfind(struct Net *net, double *tdelays, double *params, double *SNRsq, gsl_rng * r, double gmst);
double skylike(struct Net *net, double *params, double *D, double *H, double **DHc, double **DHs, double dt, int nt, int flag, struct bayesCBC *rundata);
void skylikesetup(struct Net *net, double **data, double **wave, double *D, double *H, double **DHc, double **DHs, double Tobs, int n, int bn, int nt, int imin, int imax);
void fisherskysetup(struct Net *net, double **wave, double **HH, double Tobs, int n);
double log_likelihood_full(struct Net *net, double **D, double *params, RealVector *freq, double **SN, double *rho, int N, double Tobs, struct bayesCBC *rundata);
void upsample(int n, double T, int *nt, int *bn);
// Conversion routines (between geocenter and detector frames)
void dshifts(struct Net *net, double *sky, double *params, int NX, struct bayesCBC *rundata);
void pmap_all(struct Net *net, double *pall, double *param, double *sky, int NX, double gmst);
void TimeDelays(struct Net *net, double alpha, double sindelta, double *dtimes, double gmst);
void ComputeDetFant(struct Net *net, double psi, double alpha, double sindelta, double *Fplus, double *Fcross, int id, double gmst);
// Temporarily setting these here for computing D in main BayesWave
void tukey(double *data, double alpha, int N);
void tukey_scale(double *s1, double *s2, double alpha, int N);
\ No newline at end of file
src/BayesCBC_internals.h
View file @ 37ea10d6
......@@ -85,7 +85,6 @@ 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);
......@@ -109,8 +108,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 templates(struct Net *net, double **hwave, RealVector *freq, double *params, int N, struct bayesCBC *rundata);
void geotemplate(double *gwave, RealVector *freq, double *params, int N, struct bayesCBC *rundata);
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);
......@@ -119,18 +116,15 @@ void Fisher_All(struct Net *net, double **fish, double *params, RealVector *freq
double log_likelihood(struct Net *net, double **D, double *params, RealVector *freq, double **SN, int N, double Tobs, struct bayesCBC *rundata);
double log_likelihood_max(struct Net *net, double **D, double *params, RealVector *freq, double **SN, int N, double Tobs, double tmin, double tmax, struct bayesCBC *rundata);
double log_likelihood_full(struct Net *net, double **D, double *params, RealVector *freq, double **SN, double *rho, int N, double Tobs, struct bayesCBC *rundata);
double log_likelihood_ext(struct Net *net, double **D, double *params, RealVector *freq, double *hcos, double *hsin, double **SN, double *rho, int N, double Tobs);
void jacobi(double **a, int n, double e[], double **v, int *nrot);
void fisherskysetup(struct Net *net, double **wave, double **HH, double Tobs, int n);
void RingPlot(double *skyx, int d1, int d2, double *theta, double *phi, int M, int NX, double gmst);
void Ring(double *skyx, double *skyy, int d1, int d2, gsl_rng * r, double gmst);
void ringfind(struct Net *net, double *tdelays, double *params, double *SNRsq, gsl_rng * r, double gmst);
double skydensity(double *paramsx, double *paramsy, int ifo1, int ifo2, int iref, int NS, double gmst);
void skymap(double *paramsx, double *paramsy, int ifo1, int ifo2, int iref, double gmst);
void Ring(struct Net *net, double *skyx, double *skyy, int d1, int d2, gsl_rng * r, double gmst);
double skydensity(struct Net *net, double *paramsx, double *paramsy, int ifo1, int ifo2, int iref, int NS, double gmst);
void skymap(struct Net *net, double *paramsx, double *paramsy, int ifo1, int ifo2, int iref, double gmst);
void Ring_all(double *paramx, double *paramy, int d1, int d2, gsl_rng * r, int NX, double gmst);
void uvwz_all(double *u, double *v, double *w, double *z, double *params, int NXs);
......@@ -144,14 +138,9 @@ void F_LHO(double psi, double alpha, double sindelta, double *Fplus, double *Fcr
void F_LLO(double psi, double alpha, double sindelta, double *Fplus, double *Fcross, double gmst);
void F_VIRGO(double psi, double alpha, double sindelta, double *Fplus, double *Fcross, double gmst);
void skyring(struct Net *net, double *params, double *sky, double *pall, double *SNRsq, RealVector *freq, double **D, double **SN, int N, double Tobs, gsl_rng * r, int NX, int NS, double gmst);
void skystart(struct Net *net, double *params, double *sky, double *pall, double *SNRsq, RealVector *freq, double **D, double **SN, int N, double Tobs, gsl_rng * r, struct bayesCBC *rundata);
double skylike(struct Net *net, double *params, double *D, double *H, double **DHc, double **DHs, double dt, int nt, int flag, struct bayesCBC *rundata);
void skylikesetup(struct Net *net, double **data, double **wave, double *D, double *H, double **DHc, double **DHs, double Tobs, int n, int bn, int nt, int imin, int imax);
void upsample(int n, double T, int *nt, int *bn);
void uvwz(double *u, double *v, double *w, double *z, double *params);
void exsolve(double *phiy, double *psiy, double *Ay, double *ey, double uy, double vy, double wy, double zy);
void uvwz_sol(double *uy, double *vy, double *wy, double *zy, double ux, double vx, double wx, double zx, \
......@@ -183,16 +172,12 @@ 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, 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);
void CBC_update(struct Net *net, int *mxc, FILE *chainI, FILE *chainE, double **paramx, double **skyx, double **pallx, int *who, double *heat, double ***history, double ***global, RealVector *freq, double **D, double **SN, int N, double Tobs, gsl_rng *r, struct bayesCBC *rundata);
void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, double **paramx, double **skyx, double **pallx, int *who, double *heat, double ***history, double ***global, RealVector *freq, double **D, double **SN, int N, double Tobs, gsl_rng *r, struct bayesCBC *rundata);
void skymcmc(struct Net *net, int MCX, int *mxc, FILE *chain, double **paramx, double **skyx, double **pallx, int *who, double *heat, double dtx, int nt, double *DD, double **WW, double ***DHc, double ***DHs, double ***HH, double Tobs, gsl_rng * r, struct bayesCBC *rundata);
void updatei(int k, struct Net *net, int lmax, double *logLx, double **paramx, double **paramy, double *min, double *max, double *Fscale, int *who, double *heat, double ***history, double ***global, RealVector *freq, double **D, double **SN, double **ejump, double ***evec, int N, double Tobs, int **cv, int **av, gsl_rng *r, struct bayesCBC *rundata);
void update(int k, struct Net *net, double *logLx, double *logPx, double **rhox, double **paramx, double **paramy, double *min, double *max, int *who, double *heat, double ***history, double ***global, RealVector *freq, double **D, double **SN, double **ejump, double ***evec, int N, double Tobs, int **cv, int **av, gsl_rng *r, int *evolveparams, struct bayesCBC *rundata);
......
src/BayesWaveMCMC.c
View file @ 37ea10d6
......@@ -1255,10 +1255,20 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
// Quick search to find initial signal
fprintf(stdout, " ======= Start CBC_burnin() ======\n");
fprintf(stdout, " Start #1 intrinsic 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, 0, bayescbc);
fclose(bayescbc->chainS);
// Quick search to find initial signal extrinisc parameters
// Currently is not executed with sky fixed flag (assuming sky paramters are given)
// TODO: Allow for extrinsic updates in burnin with sky-fix.
if (bayescbc->intrinsic_only == 0 && data->skyFixFlag == 0)
{
fprintf(stdout, " Start #2 (optional) extrinisic burnin \n");
extrinsic_burnin_cbc(bayescbc, data, model, chain);
}
// Pass updated cbc parameters to model
for(ic=0; ic<chain->NC; ic++)
{
......@@ -2905,6 +2915,241 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
free_double_matrix(intParams,dim);
}
/*
* Initial (optional) MCMC search to find extrinsic model parameters
*/
void extrinsic_burnin_cbc(struct bayesCBC *bayescbc, struct Data *data, struct Model **model, struct Chain *chain)
{
int i, j, k, id;
int nt, bn;
int imin, imax;
double x, dtx;
int mty, mtid, mti;
double *logL, *logLsky, *logLfull;
double **ext_data, ***wave;
double *DD, **WW;
double *rho;
double ***DHc, ***DHs, ***HH;
double *heat;
double Lmax;
FILE *chainFile;
// local pointers
struct Net *net = bayescbc->net;
int N = data->N;
int NX = bayescbc->NX;
int NP = bayescbc->NP;
int NS = bayescbc->NS;
int NC = bayescbc->NC; // number of chains
int NCC = bayescbc->NCC; // number of cold chains
int NH = bayescbc->NH; // length of history
int NQ = bayescbc->NQ; // number of mass ratios in global proposal
int NM = bayescbc->NM; // number of chirp masses in global proposal
double **paramx = bayescbc->paramx; // intrinsic burnin parameters
double **skyx = bayescbc->skyx; // intrinsic burnin parameters
double **pallx = bayescbc->pallx; // intrinsic burnin parameters
double **D = bayescbc->D;
double **SN = bayescbc->SN;
logLsky = double_vector(NC);
logLfull = double_vector(NC);
// local heat for burnin only
heat = double_vector(NC);
// whitened signal in each detector
wave = double_tensor(NC,net->Nifo,N);
// whitened data in each detector
ext_data = double_matrix(net->Nifo,N);
rho = double_vector(net->Nifo);
// convert data to local format (needed for local fourier transforms)
for (id = 0; id < net->Nifo; ++id)
{
ext_data[id][0] = 0.0;
ext_data[id][N/2] = 0.0;
}
for (i = 1; i < N/2; ++i)
{
for (id = 0; id < net->Nifo; ++id)
{
x = 1.0/sqrt(SN[id][i]);
ext_data[id][i] = D[id][2*i]*x;
ext_data[id][N-i] = D[id][2*i+1]*x;
}
}
Lmax = 0.0;
upsample(N, data->Tobs, &nt, &bn);
dtx = data->Tobs/(double)(bn);
imin = (int)((data->fmin*data->Tobs));
imax = (int)((data->fmax*data->Tobs));
DD = double_vector(net->Nifo);
WW = double_matrix(NC,net->Nifo);
DHc = double_tensor(NC,net->Nifo,nt+1);
DHs = double_tensor(NC,net->Nifo,nt+1);
HH = double_tensor(NC,net->Nifo,3);
printf("\nFinding sky location\n");
// // Find sky location and set up whitend signal arrays
#pragma omp parallel for
for(j = 0; j < NC; j++)
{
int mtid, mti;
double mty, Scale;
double *SNRsq;
double **hwave;
SNRsq = double_vector(net->Nifo);
hwave = double_matrix(net->Nifo,N);
// hwave in [0,N/2] real and [N/2,N] imag
templates(net, hwave, bayescbc->freq, paramx[j], N, bayescbc);
// // might want to put a catch here for any chains that didn't reach a decent logL and re-run the rearch phase
// // logLstart[j] = log_likelihood(net, D, paramx[j], freq, SN, N, Tobs, rundata);
mty = 0.0;
for (mtid = 0; mtid < net->Nifo; ++mtid)
{
SNRsq[mtid] = 0.0;
for (mti = 1; mti < N/2; ++mti)
{
// The sky mapping code puts in the amplitude and phase shifts. The signals only differ due to the whitening
SNRsq[mtid] += 4.0*(hwave[0][mti]*hwave[0][mti]+hwave[0][N-mti]*hwave[0][N-mti])/SN[mtid][mti];
}
Scale = 1.0;
if(mtid > 0) Scale = paramx[j][(mtid-1)*3+NX+2]*paramx[j][(mtid-1)*3+NX+2];
SNRsq[mtid] *= Scale;
mty += SNRsq[mtid];
}
printf("%f %f\n", SNRsq[0], SNRsq[1]);
// find a sky location roughly consistent with the time delays
skyring(net, paramx[j], skyx[j], pallx[j], SNRsq, bayescbc->freq, D, SN, N, data->Tobs, bayescbc->rvec[j], NX, NS, bayescbc->gmst);
dshifts(net, skyx[j], paramx[j], NX, bayescbc);
// map to the geocenter params
pmap_all(net, pallx[j], paramx[j], skyx[j], NX, bayescbc->gmst);
// get the geocenter reference template. This does depend on the assumed sky location, hence follows skystart
geotemplate(hwave[0], bayescbc->freq, pallx[j], N, bayescbc);
for (mtid = 0; mtid < net->Nifo; ++mtid)
{
wave[j][mtid][0] = 0.0;
wave[j][mtid][N/2] = 0.0;
for (mti = 1; mti < N/2; ++mti)
{
// The sky mapping code puts in the amplitude and phase shifts. The signals only differ due to the whitening
mty = 1.0/sqrt(SN[mtid][mti]);
wave[j][mtid][mti] = hwave[0][mti]*mty;
wave[j][mtid][N-mti] = hwave[0][N-mti]*mty;
}
}
skylikesetup(net, ext_data, wave[j], DD, WW[j], DHc[j], DHs[j], data->Tobs, N, bn, nt, imin, imax);
fisherskysetup(net, wave[j], HH[j], data->Tobs, N);
logLsky[j] = skylike(net, skyx[j], DD, WW[j], DHc[j], DHs[j], dtx, nt, 0, bayescbc);
logLfull[j] = log_likelihood_full(net, D, pallx[j], bayescbc->freq, SN, rho, N, data->Tobs, bayescbc);
printf("%d logL full %f sky %f\n", j, logLfull[j], logLsky[j]);
free_double_vector(SNRsq);
free_double_matrix(hwave,net->Nifo);
}
// intialize the local heat array
// run cold to force ML
for(i=0; i<NCC; i++) heat[i] = 0.25;
for(i=NCC; i<NC; i++) heat[i] = heat[i-1]*1.25; // spacing can be big here since we start at low temperature
// with 10 hot chains this gets the hottest up to heat > 1.
chainFile = fopen("sky.dat","w");
skymcmc(net, 1000000, bayescbc->mxc, chainFile, paramx, skyx, pallx, bayescbc->who, heat, dtx, nt, DD, WW, DHc, DHs, HH, data->Tobs ,chain->seed, bayescbc);
fclose(chainFile);
for(j = 0; j < NC; j++)
{
pmap_all(net, pallx[j], paramx[j], skyx[j], NX, bayescbc->gmst);
logLfull[j] = log_likelihood_full(net, D, pallx[j], bayescbc->freq, SN, rho, N, data->Tobs, bayescbc);
}
if (bayescbc->debug)
{
for(j = 0; j < NC; j++)
{
printf("\n");
printf("%d logL full %f\n", j, logLfull[j]);
}
printf("\n");
}
x = 0.0;
k = bayescbc->who[0];
for(j = 0; j < NC; j++)
{
if(logLfull[j] > x)
{
x = logLfull[j];
k = j;
}
}
// set threshold that all chains start withing 20% of highest likelihood
x *= 0.8;
for(j = 0; j < NC; j++)
{
// printf("%d %f ", jj, logLfull[jj]);
if(logLfull[j] < x)
{
for(i = 0; i < NP; i++)
{
pallx[j][i] = pallx[k][i];
}
}
logLfull[j] = log_likelihood_full(net, D, pallx[j], bayescbc->freq, SN, rho, N, data->Tobs, bayescbc);
bayescbc->logLx[j] = logLfull[j]; // BW shared array
printf("%f\n", logLfull[j]);
}
free_double_vector(heat);
free_double_vector(rho);
free_double_vector(logLsky);
free_double_vector(logLfull);
free_double_tensor(wave,NC,net->Nifo);
free_double_vector(DD);
free_double_matrix(WW,NC);
free_double_tensor(DHc,NC,net->Nifo);
free_double_tensor(DHs,NC,net->Nifo);
free_double_tensor(HH,NC,net->Nifo);
free_double_matrix(ext_data,net->Nifo);
printf("after skyring and freed 2\n");
}
/* ********************************************************************************** */
/* */
/* MCMC tools */
......
src/BayesWaveMCMC.h
View file @ 37ea10d6
......@@ -31,6 +31,8 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Model **model, struct Chain *chain, gsl_rng *seed, int ic);
void EvolveBayesCBCParameters(struct Data *data, struct Model **model, struct bayesCBC *bayescbc, struct Chain *chain, int ic);
void extrinsic_burnin_cbc(struct bayesCBC *bayescbc, struct Data *data, struct Model **model, struct Chain *chain);
void kickstart_glitch_model(struct Data *data, struct Prior *prior, struct Model *model);
/* ********************************************************************************** */
......
src/BayesWaveModel.c
View file @ 37ea10d6
......@@ -1002,7 +1002,7 @@ void free_bayescbc(struct bayesCBC *bayescbc, struct Data *data, struct Chain *c
void initialize_bayescbc(struct bayesCBC *bayescbc, struct Data *data, struct Chain *chain, struct Model *model)
{
ProcessParamsTable *ppt = NULL;
int i, ifo, ic;
int i, j, ifo, ic;
int imin, imax;
double x, y;
double fac;
......@@ -1067,6 +1067,19 @@ void initialize_bayescbc(struct bayesCBC *bayescbc, struct Data *data, struct Ch
// Transfer Network struct to CBC Net
initialize_net(bayescbc->net, data->Tobs, data->trigtime, NI+3, data->ifos);
bayescbc->net->location = double_matrix(NI,3);
bayescbc->net->response = double_tensor(NI,3, 3);
// Add detector locations
for(ifo=0; ifo<NI; ifo++)
{
for(i=0; i<3; i++)
{
bayescbc->net->location[ifo][i] = (double) data->detector[ifo]->location[i];
for (j=0; j<3; j++) bayescbc->net->response[ifo][i][j] = (double) data->detector[ifo]->response[i][j];
}
}
// detLogL stores the likelihood in each detector (for quick computation of full logL) for current chain
bayescbc->detLogL = double_vector(NI-1);
bayescbc->who = chain->index;
......
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