Commit 3795d6bc authored by Marcella Wijngaarden's avatar Marcella Wijngaarden

CBC output control options + seperate heat burnin + fix max/min

parent 80b5cad3
This diff is collapsed.
......@@ -42,6 +42,9 @@ struct bayesCBC
// (temporarily as needed when running multitheaded)
// Note: moved to BW Chain structure!
int debug; // if > 0, prints debug statements in BayesCBCs
int verbose; // set by --verbose flag
// Maybe want to move this here too?
struct Net *net;
double **D;
......@@ -53,7 +56,6 @@ struct bayesCBC
double ***history;
double **SN;
int *mxc;
double ***historyall;
RealVector *freq;
FILE *chainA; // Chains file for main MCMC
......
......@@ -185,7 +185,7 @@ double globe(double ***global, double *max, double *min, double Tobs, double *pa
void dshifts(struct Net *net, double *sky, double *params, int NX, double gmst);
void CBC_start(struct Net *net, int *mxc, FILE *chainI, 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 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);
......
......@@ -1075,51 +1075,53 @@ void print_cbc_model(FILE *fptr, struct bayesCBC *bayescbc, int ic, double Tobs)
// The below is for recreating the QuickCBC allchain file units
// Used for testing (but perhaps not necessary after?)
double Mchirp = exp(bayescbc->pallx[ic][0])/LAL_MSUN_SI;
double Mtot = exp(bayescbc->pallx[ic][1])/LAL_MSUN_SI;
double eta = pow((Mchirp/Mtot), (5.0/3.0));
double dm = sqrt(1.0-4.0*eta);
double m1 = Mtot*(1.0+dm)/2.0;
double m2 = Mtot*(1.0-dm)/2.0;
double chieff = (m1*bayescbc->pallx[ic][2]+m2*bayescbc->pallx[ic][3])/Mtot;
double DL = exp(bayescbc->pallx[ic][6])/(1.0e6*LAL_PC_SI); // includes conversion to Mpc
double z = z_DL(DL);
// TODO: Add suport below for converting & printing neutron star tidal parameters
// if (NX > 7) {
// if (lambda_type_version == lambdaTilde) {
// lambdaT = bayescbc->pallx[7];
// dLambdaT = bayescbc->pallx[8];
// LambdaTsEta2Lambdas(lambdaT, dLambdaT, eta, &lambda1, &lambda2);
// } else {
// lambda1 = bayescbc->pallx[7];
// lambda2 = bayescbc->pallx[8];
// }
// fprintf(chain,"%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e\t %e\t %e\t %e\t %e\t %e\t %e\t %e\t", mxc[2], rundata->logLx[ic], \
// Mchirp, Mtot, chieff, paramx[4], \
// tc, DL, paramx[NX], paramx[NX+1], \
// paramx[NX+2], paramx[NX+3], z, Mchirp/(1.0+z), Mtot/(1.0+z), m1/(1.0+z), m2/(1.0+z), m1, m2,
// lambda1, lambda2, lambdaT, dLambdaT);
// } else {
// Merger time is printed relative to trigger time, which is at Tobs-2.
double tc =bayescbc->pallx[ic][5]-(Tobs-2.0);
// [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
// counter, log likelihood, chirp mass (DF), total mass (DF), effective spin, geocenter orbital phase, geocenter arrival time, distance, RA , sine of DEC,
// polarization angle, cos inclination, redshift, chirp mass (SF), total mass (SF), m1 (SF), m2 (SF), m1 (DF), m2 (DF), arrival time at each detector, SNR at each detector
fprintf(bayescbc->chainA,"%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e ", bayescbc->mxc[2], bayescbc->logLx[ic], Mchirp, Mtot, chieff, bayescbc->pallx[ic][4], \
tc, DL, bayescbc->pallx[ic][bayescbc->NX], bayescbc->pallx[ic][bayescbc->NX+1], \
bayescbc->pallx[ic][bayescbc->NX+2], bayescbc->pallx[ic][bayescbc->NX+3], z, Mchirp/(1.0+z), Mtot/(1.0+z), m1/(1.0+z), m2/(1.0+z), m1, m2);
fprintf(bayescbc->chainA,"\n");
// fprintf(stdout,"printing RA/DEC\n");
// fprintf(stdout,"%f %f \n", bayescbc->pallx[ic][bayescbc->NX], bayescbc->pallx[ic][bayescbc->NX+1]);
printf("%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e \n", bayescbc->mxc[2], bayescbc->logLx[ic], Mchirp, Mtot, chieff, bayescbc->pallx[ic][4], \
tc, DL, bayescbc->pallx[ic][bayescbc->NX], bayescbc->pallx[ic][bayescbc->NX+1], \
bayescbc->pallx[ic][bayescbc->NX+2], bayescbc->pallx[ic][bayescbc->NX+3], z, Mchirp/(1.0+z), Mtot/(1.0+z), m1/(1.0+z), m2/(1.0+z), m1, m2);
if (bayescbc->verbose==1)
{
double Mchirp = exp(bayescbc->pallx[ic][0])/LAL_MSUN_SI;
double Mtot = exp(bayescbc->pallx[ic][1])/LAL_MSUN_SI;
double eta = pow((Mchirp/Mtot), (5.0/3.0));
double dm = sqrt(1.0-4.0*eta);
double m1 = Mtot*(1.0+dm)/2.0;
double m2 = Mtot*(1.0-dm)/2.0;
double chieff = (m1*bayescbc->pallx[ic][2]+m2*bayescbc->pallx[ic][3])/Mtot;
double DL = exp(bayescbc->pallx[ic][6])/(1.0e6*LAL_PC_SI); // includes conversion to Mpc
double z = z_DL(DL);
// TODO: Add suport below for converting & printing neutron star tidal parameters
// if (NX > 7) {
// if (lambda_type_version == lambdaTilde) {
// lambdaT = bayescbc->pallx[7];
// dLambdaT = bayescbc->pallx[8];
// LambdaTsEta2Lambdas(lambdaT, dLambdaT, eta, &lambda1, &lambda2);
// } else {
// lambda1 = bayescbc->pallx[7];
// lambda2 = bayescbc->pallx[8];
// }
// fprintf(chain,"%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e\t %e\t %e\t %e\t %e\t %e\t %e\t %e\t", mxc[2], rundata->logLx[ic], \
// Mchirp, Mtot, chieff, paramx[4], \
// tc, DL, paramx[NX], paramx[NX+1], \
// paramx[NX+2], paramx[NX+3], z, Mchirp/(1.0+z), Mtot/(1.0+z), m1/(1.0+z), m2/(1.0+z), m1, m2,
// lambda1, lambda2, lambdaT, dLambdaT);
// } else {
// Merger time is printed relative to trigger time, which is at Tobs-2.
double tc =bayescbc->pallx[ic][5]-(Tobs-2.0);
// [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
// counter, log likelihood, chirp mass (DF), total mass (DF), effective spin, geocenter orbital phase, geocenter arrival time, distance, RA , sine of DEC,
// polarization angle, cos inclination, redshift, chirp mass (SF), total mass (SF), m1 (SF), m2 (SF), m1 (DF), m2 (DF), arrival time at each detector, SNR at each detector
fprintf(bayescbc->chainA,"%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e ", bayescbc->mxc[2], bayescbc->logLx[ic], Mchirp, Mtot, chieff, bayescbc->pallx[ic][4], \
tc, DL, bayescbc->pallx[ic][bayescbc->NX], bayescbc->pallx[ic][bayescbc->NX+1], \
bayescbc->pallx[ic][bayescbc->NX+2], bayescbc->pallx[ic][bayescbc->NX+3], z, Mchirp/(1.0+z), Mtot/(1.0+z), m1/(1.0+z), m2/(1.0+z), m1, m2);
fprintf(bayescbc->chainA,"\n");
// printf("%d %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e %e \n", bayescbc->mxc[2], bayescbc->logLx[ic], Mchirp, Mtot, chieff, bayescbc->pallx[ic][4], \
// tc, DL, bayescbc->pallx[ic][bayescbc->NX], bayescbc->pallx[ic][bayescbc->NX+1], \
// bayescbc->pallx[ic][bayescbc->NX+2], bayescbc->pallx[ic][bayescbc->NX+3], z, Mchirp/(1.0+z), Mtot/(1.0+z), m1/(1.0+z), m2/(1.0+z), m1, m2);
}
bayescbc->mxc[2] += 1;
......
......@@ -883,7 +883,8 @@ double EvaluateExtrinsicMarkovianLogLikelihood(struct Network *projection, doubl
//compute instrument response from cbc waveform at geocenter with extrinsic params
if(data->bayesCBCFlag)
{
{
for(i=0; i<NI; i++) for(n=0; n<N; n++) c[i][n] = 0.0;
waveformProject(data, projection, params, c, geocbc, fmin, fmax);
}
......
......@@ -1255,13 +1255,17 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b
// Set BayesCBC rundata (this replaces the header file settings of QuickCBC)
initialize_bayescbc(bayescbc, data, chain, model[chain->index[0]]);
if(LALInferenceGetProcParamVal(data->commandLine, "--debug-cbc")) {
bayescbc->debug = 1;
}
// Use XML injection parameter file from BayesWave command line
int start_at_injection = 0;
if(LALInferenceGetProcParamVal(data->commandLine, "--inj"))
{
sprintf(filename, LALInferenceGetProcParamVal(data->commandLine,"--inj")->value);
start_at_injection = 1;
// If injecting tidal parameters, read values from command line and setup version + sampler type
if(LALInferenceGetProcParamVal(data->commandLine, "--inj-lambda1")||LALInferenceGetProcParamVal(data->commandLine, "--inj-lambdaT"))
{
......@@ -2660,8 +2664,6 @@ void EvolveBayesCBCParameters(struct Data *data, struct Model **model, struct ba
imin = fmin*data->Tobs;
imax = fmax*data->Tobs;
if (ic==0) printf("chain-cycle: %i \n", chain->cycle);
// Setup residuals for BayesCBC (data substracted with ex/intrinsic
// model template and glitch template)
for(ifo=0; ifo<NI; ifo++)
......@@ -2776,7 +2778,7 @@ void EvolveBayesCBCParameters(struct Data *data, struct Model **model, struct ba
}
}
// printf("after BayesCBC_MCMC chain %i, %i, model_x->logL = %f, logLnorm = %f, logLx = %f \n", chain->index[ic], ic, model_x->logL, model_x->logLnorm, bayescbc->logLx[icbc]);
if (bayescbc->debug==1) printf("after BayesCBC_MCMC chain %i, %i, model_x->logL = %f, logLnorm = %f, logLx = %f \n", chain->index[ic], ic, model_x->logL, model_x->logLnorm, bayescbc->logLx[icbc]);
return;
}
......
......@@ -1036,7 +1036,6 @@ void initialize_bayescbc(struct bayesCBC *bayescbc, struct Data *data, struct Ch
bayescbc->twidth = 0.1; // width of prior on merger time (used when given non-integer merger time)
bayescbc->intrinsic_only = 1; // evolve only semi-intrinsic parameters (NX), not sky
// False for now during testing, but should be true from BW MCMC
bayescbc->logLx = double_vector(bayescbc->NC);
bayescbc->net = malloc(sizeof(struct Net));
bayescbc->global = double_tensor(bayescbc->NQ,bayescbc->NM,N+1);
......@@ -1045,20 +1044,17 @@ void initialize_bayescbc(struct bayesCBC *bayescbc, struct Data *data, struct Ch
bayescbc->heat = double_vector(bayescbc->NC);
bayescbc->paramx = double_matrix(bayescbc->NC,bayescbc->NX+3*NI);
bayescbc->history = double_tensor(bayescbc->NC,bayescbc->NH+1,bayescbc->NX);
bayescbc->historyall = double_tensor(bayescbc->NC,bayescbc->NH+1,bayescbc->NX);
bayescbc->freq = CreateRealVector((N/2));
bayescbc->mxc = int_vector(3);
for (i = 0; i < 3; i++) bayescbc->mxc[i] = 0;
bayescbc->pallx = double_matrix(bayescbc->NC,bayescbc->NP);
// Hold CBC model for each chain and each detector
// bayescbc->template = double_tensor(bayescbc->NC,bayescbc->NI,N);
// Transfer Network struct to CBC Net
initialize_net(bayescbc->net, data->Tobs, data->trigtime, NI+3, data->ifos);
// 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;
// Set prior range BayesCBC
set_bayescbc_priors(bayescbc->net, bayescbc);
......@@ -1076,7 +1072,7 @@ void initialize_bayescbc(struct bayesCBC *bayescbc, struct Data *data, struct Ch
{
for(i=0; i<N; i++)
{
//copy over current residual
//copy over data to start with (used in burnin)
bayescbc->D[ifo][i] = data->s[ifo][i];
}
}
......@@ -1117,7 +1113,6 @@ void initialize_bayescbc(struct bayesCBC *bayescbc, struct Data *data, struct Ch
fprintf(stdout, "Using LAL PSD for BayesCBC:\n");
// Use LAL PSD
// Compute SN from psd arrays == <n_i^2> = T/2 * Sn(f)
// fac = 1.0;//2.0/Tobs;
for(ifo=0; ifo<NI; ifo++)
{
for (i=0; i<data->N/2; ++i)
......@@ -1134,7 +1129,20 @@ void initialize_bayescbc(struct bayesCBC *bayescbc, struct Data *data, struct Ch
}
}
bayescbc->chainA = fopen("chains/cbc_allchain.dat","w");
// CBC chains with parameters in useful units (if --verbose prints
// all chains in the file, otherwise only cold chain)
bayescbc->chainA = fopen("chains/cbc_allchain.dat","w");
bayescbc->debug = 0;
if(data->verboseFlag)
{
bayescbc->verbose = 1;
}
else
{
bayescbc->verbose = 0;
}
}
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