Commit 08845948 authored by Marcella Wijngaarden's avatar Marcella Wijngaarden
Browse files

Cleanup unused code in BayesCBC + headers + inj option

parent b7c87abc
......@@ -27,7 +27,6 @@ along with this program. If not, see <http://www.gnu.org/licenses/>.
#include "BayesCBC.h"
#include "BayesCBC_internals.h"
#include "ConstCBC.h"
#include "Constants.h"
#include "IMRPhenomD.h"
......@@ -1014,15 +1013,6 @@ void burnin_bayescbc(int N, double Tobs, double ttrig, double **D, double **SN,
}
else
{
// TODO: The initialization here is a DRY violation as a lot of this is also done in CBC_start
// Temporary double so that it still works as stand-alone (i.e., when not calling from
// BayesWave), but should update.
// intialize the who, heat, history and counter arrays
// run cold to force ML
for(i=0; i<NCC; i++) rundata->heat[i] = 0.25;
for(i=NCC; i<NC; i++) rundata->heat[i] = rundata->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.
for(j = 0; j < NC; j++)
{
......@@ -1038,10 +1028,7 @@ void burnin_bayescbc(int N, double Tobs, double ttrig, double **D, double **SN,
rundata->mxc[0] = 0;
rundata->mxc[1] = 0;
logL = double_vector(NC);
logLstart = double_vector(NC);
logLsky = double_vector(NC);
// logLfull = double_vector(NC);
logLfull = double_vector(NC);
rho = double_matrix(NC,rundata->net->Nifo);
params = (double*)malloc(sizeof(double)* (NP));
......@@ -1052,72 +1039,70 @@ void burnin_bayescbc(int N, double Tobs, double ttrig, double **D, double **SN,
char filename[1024];
// Read in injection parameters from injection XML file.
// Note: assuming always starting from XML when injecting
if (true) {
printf("Reading XML injection parameters from file: %s\n", xml_file_path);
sprintf(filename, xml_file_path);
ParseInjectionXMLParameters(filename, injection_columns, injection_values);
// Set chirp mass and total mass
GetXMLInjectionParameter("sim_inspiral:mchirp", &params[0], injection_columns, injection_values);
params[0] = log(params[0]*MSUN_SI); // log(Mc)
// Set total mass from component masses & chirp mass
GetXMLInjectionParameter("sim_inspiral:mass1", &m1, injection_columns, injection_values);
GetXMLInjectionParameter("sim_inspiral:mass2", &m2, injection_columns, injection_values);
q = m1/m2; // q = m1/m2
mc = exp(params[0]);
m2 = mc*pow(1.0+q, 0.2)/pow(q,0.6);
m1 = q*m2;
mt = m1+m2;
params[1] = log(mt); // log(Mt)
// Set spins
GetXMLInjectionParameter("sim_inspiral:spin1z", &params[2], injection_columns, injection_values);
GetXMLInjectionParameter("sim_inspiral:spin2z", &params[3], injection_columns, injection_values);
// Set geocenter gravitational wave phase (2phi0)
GetXMLInjectionParameter("sim_inspiral:phi0", &params[4], injection_columns, injection_values);
params[4] = 2*params[4];
// Set peak time geocenter (currently set at 2s before end of segment)
// GetXMLInjectionParameter("sim_inspiral:end_time_gmst", &params[5], injection_columns, injection_values);
params[5] = Tobs-2.0000;
// Set distance
GetXMLInjectionParameter("sim_inspiral:distance", &params[6], injection_columns, injection_values);
params[6] = log(params[6]*1e6*PC_SI);
// Add tidal parameters (these are not in XML injection file for historical reasons)
lambda_type_version = rundata->lambda_type_version;
if (NX > 7) {
params[7] = rundata->init_tidal1;
params[8] = rundata->init_tidal2;
}
printf("Reading XML injection parameters from file: %s\n", xml_file_path);
sprintf(filename, xml_file_path);
ParseInjectionXMLParameters(filename, injection_columns, injection_values);
// Set chirp mass and total mass
GetXMLInjectionParameter("sim_inspiral:mchirp", &params[0], injection_columns, injection_values);
params[0] = log(params[0]*MSUN_SI); // log(Mc)
// Set total mass from component masses & chirp mass
GetXMLInjectionParameter("sim_inspiral:mass1", &m1, injection_columns, injection_values);
GetXMLInjectionParameter("sim_inspiral:mass2", &m2, injection_columns, injection_values);
q = m1/m2; // q = m1/m2
mc = exp(params[0]);
m2 = mc*pow(1.0+q, 0.2)/pow(q,0.6);
m1 = q*m2;
mt = m1+m2;
params[1] = log(mt); // log(Mt)
// Set sky parameter: alpha, RA
GetXMLInjectionParameter("sim_inspiral:longitude", &sky[0], injection_columns, injection_values);
// Set spins
GetXMLInjectionParameter("sim_inspiral:spin1z", &params[2], injection_columns, injection_values);
GetXMLInjectionParameter("sim_inspiral:spin2z", &params[3], injection_columns, injection_values);
// Set sky parameter: sin(delta) = sin(declination)
GetXMLInjectionParameter("sim_inspiral:latitude", &sky[1], injection_columns, injection_values);
sky[1] = sin(sky[1]);
// Set geocenter gravitational wave phase (2phi0)
GetXMLInjectionParameter("sim_inspiral:phi0", &params[4], injection_columns, injection_values);
params[4] = 2*params[4];
// Set polarization phase psi
GetXMLInjectionParameter("sim_inspiral:polarization", &sky[2], injection_columns, injection_values);
// Set peak time geocenter (currently set at 2s before end of segment)
// GetXMLInjectionParameter("sim_inspiral:end_time_gmst", &params[5], injection_columns, injection_values);
params[5] = Tobs-2.0000;
// Set inclination angle orbital plane: cos(iota)
GetXMLInjectionParameter("sim_inspiral:inclination", &sky[3], injection_columns, injection_values);
sky[3] = cos(sky[3]);
// Set distance
GetXMLInjectionParameter("sim_inspiral:distance", &params[6], injection_columns, injection_values);
params[6] = log(params[6]*1e6*PC_SI);
// Add tidal parameters (these are not in XML injection file for historical reasons)
lambda_type_version = rundata->lambda_type_version;
if (NX > 7) {
params[7] = rundata->init_tidal1;
params[8] = rundata->init_tidal2;
}
// Set sky parameter: alpha, RA
GetXMLInjectionParameter("sim_inspiral:longitude", &sky[0], injection_columns, injection_values);
// Set sky parameter: sin(delta) = sin(declination)
GetXMLInjectionParameter("sim_inspiral:latitude", &sky[1], injection_columns, injection_values);
sky[1] = sin(sky[1]);
// Set polarization phase psi
GetXMLInjectionParameter("sim_inspiral:polarization", &sky[2], injection_columns, injection_values);
// Set inclination angle orbital plane: cos(iota)
GetXMLInjectionParameter("sim_inspiral:inclination", &sky[3], injection_columns, injection_values);
sky[3] = cos(sky[3]);
// Initialize [4] scale [5] dphi [6] dt
// Scale changes the overall distance reference at geocenter
// dt shifts the geocenter time
// dphi shifts the geocenter phase
sky[4] = 1.0;
sky[5] = 0.0;
sky[6] = 0.0;
// Initialize [4] scale [5] dphi [6] dt
// Scale changes the overall distance reference at geocenter
// dt shifts the geocenter time
// dphi shifts the geocenter phase
sky[4] = 1.0;
sky[5] = 0.0;
sky[6] = 0.0;
}
printf("\n");
......@@ -1126,7 +1111,7 @@ void burnin_bayescbc(int N, double Tobs, double ttrig, double **D, double **SN,
// Start chains at given parameter values
for (i=0; i<NC; i++) {
// Set full paramater array
// Set full parameter array
for (j = 0; j < NP; ++j) {
// Quasi-intrinsic parameters
......@@ -1142,96 +1127,18 @@ void burnin_bayescbc(int N, double Tobs, double ttrig, double **D, double **SN,
rundata->freq->data[0] = 1.0/Tobs;
for (i=1; i< N/2; i++) rundata->freq->data[i] = (double)(i)/Tobs;
fmin = rundata->fmin;
fmax = rundata->fmax;
imin = (int)((fmin*Tobs));
imax = (int)((fmax*Tobs));
upsample(N, Tobs, &nt, &bn);
dtx = Tobs/(double)(bn);
printf("Time steps = %d time resolution = %f\n", nt, dtx);
double *dtimes;
double alpha, sindelta;
dtimes = double_vector(5);
// Initial find t0 and setting up parameters with respect to detector reference frame
// TODO: Handle when search_t0 is not true.
if (search_t0) {
for (i=0; i<NC; i++) {
logLstart[i] = log_likelihood_full(rundata->net, D, rundata->pallx[i], rundata->freq, SN, rho[i], N, Tobs, rundata);
}
for(k = 0; k < NQ; k++)
{
for(j = 0; j < NM; j++)
{
for(i = -N/2; i < N/2; i++)
{
rundata->global[k][j][i+N/2] = rundata->cap;
}
}
}
// Quick search for t0, phi0, D0
MCMC_all(rundata->net, rundata->mxc, 30000, rundata->chainS, rundata->pallx, rundata->who, rundata->heat, rundata->history, rundata->global, rundata->freq, D, SN, N, Tobs, ttrig, r, 1, rundata);
// Get likelihoods after initial search
for(j = 0; j < NC; j++) logLfull[j] = log_likelihood_full(rundata->net, D, rundata->pallx[j], rundata->freq, SN, rho[j], N, Tobs, rundata);
// Find maximum likelihood chain
x = 0.0;
k = 1;
for(jj = 0; jj < NC; jj++)
{
if(logLfull[jj] > x)
{
x = logLfull[jj];
k = jj;
}
}
// Since we are starting with injection parameters,
// we pick the chain with the highest likelihood
x *= 1.0;
for(jj = 0; jj < NC; jj++)
{
if(logLfull[jj] < x)
{
for(i = 0; i < NP; i++)
{
rundata->pallx[jj][i] = rundata->pallx[k][i];
}
}
logLfull[jj] = log_likelihood_full(rundata->net, D, rundata->pallx[jj], rundata->freq, SN, rho[jj], N, Tobs, rundata);
}
}
chain = fopen("chain_init.dat", "w");
for(k=0; k < NC; k++)
for(jj = 0; jj < NC; jj++)
{
for(j=0; j<NP; j++) {
fprintf(chain,"%e ", rundata->pallx[k][j]);
}
fprintf(chain,"\n");
logLfull[jj] = log_likelihood_full(rundata->net, D, rundata->pallx[jj], rundata->freq, SN, rho[jj], N, Tobs, rundata);
printf("logLfull[%i]: %f\n", jj, logLfull[jj]);
}
fclose(chain);
// now re-set for regular MCMC
for(i=0; i<NCC; i++) rundata->heat[i] = 1.0;
for(i=NCC; i<NC; i++) rundata->heat[i] = rundata->heat[i-1]*rundata->ladder;
free_double_matrix(rho,NC);
free_double_vector(logL);
free_double_vector(logLfull);
free_double_vector(logLstart);
free_double_vector(logLsky);
free_double_vector(heat);
free_double_vector(params);
printf("calling log_likelihood_full 8\n");
}
}
......@@ -2294,28 +2201,42 @@ void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, dou
}
// reference antenna pattern scaling used to convert reference detector frame amplitude to a physical distance
if(lmax == 0)
{
for(k = 0; k < NC; k++) Fscale[k] = Fmag(skyx[k], net->labels[0], rundata->gmst);
}
else
// Only using in lmax mode.
// if(lmax == 0)
// {
// // for(k = 0; k < NC; k++) Fscale[k] = Fmag(skyx[k], net->labels[0], rundata->gmst);
// }
// else
// {
for(k = 0; k < NC; k++)
{
for(k = 0; k < NC; k++)
{
Fscale[k] = 1.0;
}
Fscale[k] = 1.0;
}
// }
//
max[4] = PI; // phi0. Note: Should be TPI in main cbc, and PI for intrinsic param array
// random start drawn from the global
if(lmax == 1)
{
printf("\ngetting global inits for parameters\n");
for(k = 0; k < NC; k++)
{
// random start drawn from the global
if(rundata->injectionStartFlag==0)
{
printf("\ngetting global inits for parameters\n");
for(k = 0; k < NC; k++)
{
qx = globe(global, max, min, Tobs, paramx[k], N, r, lmax, rundata);
}
}
else // start from injection parameters
{
qx = globe(global, max, min, Tobs, paramx[k], N, r, lmax, rundata);
printf("\nstarting burnin from injection parameters\n");
for(k = 0; k < NC; k++)
{
qx = globe(global, max, min, Tobs, paramx[k], N, r, lmax, rundata);
}
}
}
......@@ -2768,554 +2689,107 @@ void updatei(int k, struct Net *net, int lmax, double *logLx, double **paramx,
}
void MCMC_all(struct Net *net, int *mxc, int M, FILE *chain, double **paramx, int *who, double *heat, double ***history, double ***global, RealVector *freq, double **D, double **SN, int N, double Tobs, double ttrig, gsl_rng *r, int istart, struct bayesCBC *rundata)
{
int i, j, q, k, mc, flag;
int id1, id2;
int *m, *evolveparams;
double *logLx, *logPx, logLy, logL, x;
double **paramy;
double *pref;
double *max, *min;
double *href;
double *pml;
double alpha, beta, H;
double Ap, Ac, psi, sindelta;
double Fcross, Fplus, Fp;
double Mchirp, Mtot, M1, M2, ciota, tc;
double eta, dm, m1, m2, chieff, DL;
double pMcMx, pMcMy;
double a, b, c, d;
double sqh;
double Lmax, mch;
double **tp;
double logpx, logpy;
double Jack;
double dA, dTD, dP, z, y;
double leta;
int bin;
int c1, c2, c3, c4;
int a1, a2, a3, a4;
int typ, hold;
int scount, sacc, mcount;
double *Fscale;
double pDx, pDy, qx, qy;
double *rtheta, *rphi;
double **rhox;
int MR;
double mapL, lmapx, lmapy;
double *pmax;
double *dtimes;
double lambda1, lambda2, lambdaT, dLambdaT;
Lambda_type lambda_type_version;
FILE *like;
FILE *ring;
FILE *out;
double ***fish, ***evec;
double **ejump;
int **av, **cv;
char filename[1024];
/**
* Set prior boundaries for bayescbc full parameter array (NP).
'* Boundaries are stored in min and max array.
* This could be moved to BayesWavePrior.c.
*/
void set_bayescbc_priors(struct Net *net, struct bayesCBC *rundata)
{
// If no tidal parameters:
// [0] log(Mc) [1] log(Mt) [2] chi1 [3] chi2 [4] phi0 [5] tp [6] log(DL) [7] alpha [8] sindelta [9] psi [10] ciota
// With tidal parameters:
// [0] log(Mc) [1] log(Mt) [2] chi1 [3] chi2 [4] phi0 [5] tp [6] log(DL) [7] tidal1 [8] tidal2 [9] alpha [10] sindelta [11] psi [12] ciota
int NX = rundata->NX;
int NP = rundata->NP;
int NS = rundata->NS;
int NC = rundata->NC; // number of chains
int NCC = rundata->NCC; // number of cold chains
int NH = rundata->NH; // length of history
int NQ = rundata->NQ; // number of mass ratios in global proposal
int NM = rundata->NM; // number of chirp masses in global proposal
av = int_matrix(4,NC);
cv = int_matrix(4,NC);
// number of points printed for each sky ring
MR = 1000;
rtheta = double_vector(MR+1);
rphi = double_vector(MR+1);
m = int_vector(NC);
tp = double_matrix(net->Nifo,2);
dtimes = double_vector(6);
rhox = double_matrix(NC,net->Nifo);
ejump = double_matrix(NC,NP);
fish = double_tensor(NC,NP,NP);
evec = double_tensor(NC,NP,NP);
pmax = double_vector(NP);
paramy = double_matrix(NC,NP);
logLx = double_vector(NC);
logPx = double_vector(NC);
// Prior boundaries
rundata->max = (double*)malloc(sizeof(double)* (rundata->NP));
rundata->min = (double*)malloc(sizeof(double)* (rundata->NP));
lambda_type_version=rundata->lambda_type_version;
// Determine which parameters to evolve in MCMC
// Currently hardcoded here for testing
evolveparams = int_vector(NP);
for (i=0; i<NP; i++)
{
evolveparams[i] = 1;
if (istart) evolveparams[i] = 0;
}
// Set default boundaries intrinsic parameters
rundata->DLmin = 1.0e6; // Minimum distance in pc
rundata->DLmax = 1.0e10; // Maximum distance in pc
// When calling from BayesWave, we only want to evolve
// the intrinsic part
if (rundata->intrinsic_only)
{
for (i=NX; i<NP; i++)
{
evolveparams[i] = 0;
}
}
rundata->chi1min = -1.0; // Default min, max for chi1 and chi2
rundata->chi2min = -1.0;
rundata->chi1max = 1.0;
rundata->chi2max = 1.0;
// When starting from injection parameters, we only want to
// do a quick search for t0, phi0 (and sky because there is
// a difference between LAL and QuickCBC sky params)
if (istart)
rundata->mcmax = 174.0; // maximum chirp mass in Msun. Mc is at most 0.435 times Mt
rundata->mcmin = 0.23; // minimum chirp mass in Msun
rundata->mtmax = 400.0; // maximum total mass in Msun (should be not more than 2.3 times mcmax)
rundata->mtmin = 1.0; // minimum total mass in Msun
// Set default prior range for BNS
if (rundata->NRTidal_version != (NRTidal_version_type) NoNRT_V)
{
evolveparams[4] = 1;
evolveparams[5] = 1;
evolveparams[6] = 1;
evolveparams[NX] = 1;
evolveparams[NX+1] = 1;
evolveparams[NX+2] = 1;
evolveparams[NX+3] = 1;
}
// Tidal parameter ranges
rundata->lambda1min = 0.0; // Min/max for lambda1 and lambda2
rundata->lambda1max = 1000.0; // Only used if NRTidal_value != NoNRT_V and
rundata->lambda2min = 0.0; // lambda_type_value = lambda
rundata->lambda2max = 1000.0;
// Prior boundaries
max = (double*)malloc(sizeof(double)* (NP));
min = (double*)malloc(sizeof(double)* (NP));
// [0] log(Mc) [1] log(Mt) [2] chi1 [3] chi2 [4] phi0 [5] tp [6] log(DL) [7] alpha [8] sindelta [9] psi [10] ciota
rundata->lambdaTmin = 0.0; // Min/max for lambdaT and dlambdaT
rundata->lambdaTmax = 1000.0; // Only used if NRTidal_value != NoNRT_V and
rundata->dlambdaTmin = -0.1; // lambda_type_value = lambdaTilde
rundata->dlambdaTmax = 0.1;
// Update spin prior range for BNSs
rundata->chi1min = -0.1; // Default min, max for chi1 and chi2
rundata->chi2min = -0.1;
rundata->chi1max = 0.1;
rundata->chi2max = 0.1;
// Update mass range for BNSs
rundata->mcmax = 2.6; // maximum chirp mass in Msun. Mc is at most 0.435 times Mt
rundata->mcmin = 0.23; // minimum chirp mass in Msun
rundata->mtmax = 6.0; // maximum total mass in Msun (should be not more than 2.3 times mcmax)
rundata->mtmin = 1.0; // minimum total mass in Msun
}
max[0] = log(rundata->mcmax*MSUN_SI); // Mc is at most 0.435 times Mt
max[1] = log(rundata->mtmax*MSUN_SI); // Mt
max[2] = rundata->chi1max; // chi1
max[3] = rundata->chi2max; // chi2
max[4] = TPI; // phi0
max[5] = net->tmax; // peak time (Trigger time is Tobs-2)
max[6] = log(rundata->DLmax * PC_SI); //distance
rundata->max[0] = log(rundata->mcmax*MSUN_SI); // Mc is at most 0.435 times Mt
rundata->max[1] = log(rundata->mtmax*MSUN_SI); // Mt
rundata->max[2] = rundata->chi1max; // chi1
rundata->max[3] = rundata->chi2max; // chi2
rundata->max[4] = TPI; // phi0. Note: Should be TPI in main cbc, and PI for intrinsic param array
rundata->max[5] = net->tmax; // peak time (Trigger time is Tobs-2)
rundata->max[6] = log(rundata->DLmax * PC_SI); //distance
if (NX > 7)
{
if (lambda_type_version == lambdaTilde) {
max[7] = rundata->lambdaTmax; // lambdaT
max[8] = rundata->dlambdaTmax; // dlambdaT
if (rundata->lambda_type_version == (Lambda_type) lambdaTilde) {
rundata->max[7] = rundata->lambdaTmax; // lambdaT
rundata->max[8] = rundata->dlambdaTmax; // dlambdaT
} else {
max[7] = rundata->lambda1max; // lambda1
max[8] = rundata->lambda2max; // lambda2
rundata->max[7] = rundata->lambda1max; // lambda1
rundata->max[8] = rundata->lambda2max; // lambda2
}
}
max[NX] = TPI; // alpha
max[NX+1] = 1.0; // sindelta
max[NX+2] = PI; // psi
max[NX+3] = 1.0; // ciota
rundata->max[NX] = TPI; // alpha
rundata->max[NX+1] = 1.0; // sindelta
rundata->max[NX+2] = PI; // psi
rundata->max[NX+3] = 1.0; // ciota
min[0] = log(rundata->mcmin*MSUN_SI); // Mc
min[1] = log(rundata->mtmin*MSUN_SI); // Mt
min[2] = rundata->chi1min; // chi1
min[3] = rundata->chi2min; // chi2
min[4] = 0.0; // phi0
min[5] = net->tmin; // peak time (Trigger time is Tobs-2)
min[6] = log(rundata->DLmin * PC_SI); //distance
rundata->min[0] = log(rundata->mcmin*MSUN_SI); // Mc
rundata->min[1] = log(rundata->mtmin*MSUN_SI); // Mt
rundata->min[2] = rundata->chi1min; // chi1
rundata->min[3] = rundata->chi2min; // chi2
rundata->min[4] = 0.0; // phi0
rundata->min[5] = net->tmin; // peak time (Trigger time is Tobs-2)
rundata->min[6] = log(rundata->DLmin * PC_SI); //distance
if (NX > 7)
{
if (lambda_type_version == lambdaTilde)
if (rundata->lambda_type_version == (Lambda_type) lambdaTilde)
{
min[7] = rundata->lambdaTmin; // lambdaT
min[8] = rundata->dlambdaTmin; // dlambdaT
rundata->min[7] = rundata->lambdaTmin; // lambdaT
rundata->min[8] = rundata->dlambdaTmin; // dlambdaT
} else {
min[7] = rundata->lambda1min; // lambda1
min[8] = rundata->lambda2min; // lambda2
}
}
min[NX] = 0.0; // alpha
min[NX+1] = -1.0; // sindelta
min[NX+2] = 0.0; // psi
min[NX+3] = -1.0; // ciota
//open MP modifications
omp_set_num_threads(NC);
// Initialize the likelihoods
for(k = 1; k <= NC; k++)
{
logLx[k] = 0.0;
if(rundata