Commit 6ba0c6b9 authored by Marcella Wijngaarden's avatar Marcella Wijngaarden
Browse files

Widen epsilon prior + add cbc option flat prior ciota + clean output

parent 165800ee
...@@ -1137,9 +1137,7 @@ void burnin_bayescbc(int N, double Tobs, double ttrig, double **D, double **SN, ...@@ -1137,9 +1137,7 @@ void burnin_bayescbc(int N, double Tobs, double ttrig, double **D, double **SN,
free_double_matrix(rho,NC); free_double_matrix(rho,NC);
free_double_vector(logLfull); free_double_vector(logLfull);
free_double_vector(heat); free_double_vector(heat);
free_double_vector(params); free_double_vector(params);
printf("calling log_likelihood_full 8\n");
} }
} }
...@@ -1286,13 +1284,13 @@ void CBC_start(struct Net *net, int *mxc, FILE *chainS, double **paramx, double ...@@ -1286,13 +1284,13 @@ void CBC_start(struct Net *net, int *mxc, FILE *chainS, double **paramx, double
} }
// set threshold that all chains start withing 20% of highest likelihood // set threshold that all chains start withing 20% of highest likelihood
if (rundata->debug) printf(" Chains CBC logL after intrinsic burnin:\n");
x *= 0.8; x *= 0.8;
for(jj = 0; jj < NC; jj++) for(jj = 0; jj < NC; jj++)
{ {
// printf("%d %f ", jj, logLfull[jj]);
if(logLfull[jj] < x) if(logLfull[jj] < x)
{ {
for(i = 0; i < NP; i++) for(i = 0; i < NP; i++)
...@@ -1303,8 +1301,7 @@ void CBC_start(struct Net *net, int *mxc, FILE *chainS, double **paramx, double ...@@ -1303,8 +1301,7 @@ void CBC_start(struct Net *net, int *mxc, FILE *chainS, double **paramx, double
logLfull[jj] = log_likelihood_full(net, D, pallx[jj], freq, SN, rho[jj], N, Tobs, rundata); logLfull[jj] = log_likelihood_full(net, D, pallx[jj], freq, SN, rho[jj], N, Tobs, rundata);
rundata->logLx[jj] = logLfull[jj]; // BW shared array rundata->logLx[jj] = logLfull[jj]; // BW shared array
printf("%f\n", logLfull[jj]); if(rundata->debug) printf(" %d logL full %f\n", jj, logLfull[jj]);
} }
...@@ -2224,7 +2221,7 @@ void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, dou ...@@ -2224,7 +2221,7 @@ void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, dou
// random start drawn from the global // random start drawn from the global
if(rundata->injectionStartFlag==0) if(rundata->injectionStartFlag==0)
{ {
printf("\ngetting global inits for parameters\n"); printf("\n getting global inits for intrinsic parameters burnin\n");
for(k = 0; k < NC; k++) for(k = 0; k < NC; k++)
{ {
qx = globe(global, max, min, Tobs, paramx[k], N, r, lmax, rundata); qx = globe(global, max, min, Tobs, paramx[k], N, r, lmax, rundata);
...@@ -2232,7 +2229,7 @@ void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, dou ...@@ -2232,7 +2229,7 @@ void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, dou
} }
else // start from injection parameters else // start from injection parameters
{ {
printf("\nstarting burnin from injection parameters\n"); printf("\n starting intrinsic parameter burnin from injection parameters\n");
for(k = 0; k < NC; k++) for(k = 0; k < NC; k++)
{ {
qx = globe(global, max, min, Tobs, paramx[k], N, r, lmax, rundata); qx = globe(global, max, min, Tobs, paramx[k], N, r, lmax, rundata);
...@@ -2261,7 +2258,7 @@ void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, dou ...@@ -2261,7 +2258,7 @@ void MCMC_intrinsic(struct Net *net, int lmax, int *mxc, int M, FILE *chain, dou
} }
printf("Intrinsic MCMC\n"); printf(" Intrinsic MCMC\n");
for(j = 1; j <= 4; j++) for(j = 1; j <= 4; j++)
{ {
...@@ -4197,10 +4194,7 @@ void upsample(int n, double Tobs, int *nt, int *bn) ...@@ -4197,10 +4194,7 @@ void upsample(int n, double Tobs, int *nt, int *bn)
// figure out the upscaling, must be a power of 2 // figure out the upscaling, must be a power of 2
k = (int)(pow(2.0,ceil(log(dt/dtres)/log(2.0)))); k = (int)(pow(2.0,ceil(log(dt/dtres)/log(2.0))));
printf("\n upsampling sqrt(k)=%i, k=%i \n\n", (int) ceil(sqrt(k)), k);
// k = 1;
// upsampled // upsampled
*bn = n*k; *bn = n*k;
......
...@@ -127,7 +127,7 @@ void skymcmc(struct Net *net, int MCX, int *mxc, FILE *chain, double **paramx, d ...@@ -127,7 +127,7 @@ void skymcmc(struct Net *net, int MCX, int *mxc, FILE *chain, double **paramx, d
} }
printf("Extrinsic MCMC\n"); printf(" Extrinsic MCMC\n");
ac = 0; ac = 0;
...@@ -390,7 +390,7 @@ void skymcmc(struct Net *net, int MCX, int *mxc, FILE *chain, double **paramx, d ...@@ -390,7 +390,7 @@ void skymcmc(struct Net *net, int MCX, int *mxc, FILE *chain, double **paramx, d
} }
if(mc%10000 == 0) if(mc%10000 == 0 && rundata->debug == 1)
{ {
ic = who[0]; ic = who[0];
DL = exp(pallx[ic][6])/(1.0e6*PC_SI*skyx[ic][4]); DL = exp(pallx[ic][6])/(1.0e6*PC_SI*skyx[ic][4]);
...@@ -457,13 +457,14 @@ void skymcmc(struct Net *net, int MCX, int *mxc, FILE *chain, double **paramx, d ...@@ -457,13 +457,14 @@ void skymcmc(struct Net *net, int MCX, int *mxc, FILE *chain, double **paramx, d
} }
printf("Swap Acceptance = %f\n", (double)sacc/(double)(scount)); if(rundata->debug == 1)
printf("MCMC Acceptance = %f\n", (double)ac/(double)(mcount)); { printf(" Extrinsic burnin ended with:\n");
printf("Ring Acceptance = %f\n", (double)rca/(double)(rc)); printf(" Swap Acceptance = %f\n", (double)sacc/(double)(scount));
printf("Fisher Acceptance = %f\n", (double)fac/(double)(fc)); printf(" MCMC Acceptance = %f\n", (double)ac/(double)(mcount));
printf("Jiggle Acceptance = %f\n", (double)uac/(double)(uc)); printf(" Ring Acceptance = %f\n", (double)rca/(double)(rc));
printf(" Fisher Acceptance = %f\n", (double)fac/(double)(fc));
printf(" Jiggle Acceptance = %f\n", (double)uac/(double)(uc));
}
free_double_matrix(skyy,NC); free_double_matrix(skyy,NC);
free_double_matrix(skyh,NC); free_double_matrix(skyh,NC);
......
...@@ -1304,8 +1304,6 @@ void RJMCMC(struct Data *data, struct Model **model, struct BayesLineParams ***b ...@@ -1304,8 +1304,6 @@ 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); Shf_Geocenter_full(data, model[chain->index[ic]]->projection, model[chain->index[ic]]->Snf, model[chain->index[ic]]->SnGeo, model[chain->index[ic]]->extParams);
} }
printf("model chain %i = index[%i]\n", ic, chain->index[ic]);
// TODO: Intrinsic template here // TODO: Intrinsic template here
// cbcGeoAmpPhase(bayescbc, chain->index[ic], data->N, model[chain->index[ic]]->cbcamphase); // 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); // 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);
...@@ -2872,7 +2870,24 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo ...@@ -2872,7 +2870,24 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if (data->signalFlag) 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; logH = (logLy - logLx)*chain->beta + logpy - logpx + logJ;
// Make prior flat in cos(iota) instead of epsilon
if(data->cbcFlag)
{
double ciotax = ((1.0 - sqrt(1.0 - paramsx[3]*paramsx[3]))/paramsx[3]);
double ciotay = ((1.0 - sqrt(1.0 - paramsy[3]*paramsy[3]))/paramsy[3]);
// Jacobian to make cos(iota) prior flat
// MCMC variable is epsilon, so p(epsilon) = dCiota/depsilon p(ciota)
// = (2(1-ciota**2)**2 / (1+ciota**2)**2)**-1 p(ciota)
double dCioataEx = log(((1.0 + ciotax*ciotax)*(1.0 + ciotax*ciotax))/(2.0*(1.0 - ciotax*ciotax)));
double dCioataEy = log(((1.0 + ciotay*ciotay)*(1.0 + ciotay*ciotay))/(2.0*(1.0 - ciotay*ciotay)));
logH += dCioataEy - dCioataEx;
}
if(paramsy[0]!=paramsy[0]) if(paramsy[0]!=paramsy[0])
{ {
printf("NAN!\n");fflush(stdout); printf("NAN!\n");fflush(stdout);
...@@ -3015,7 +3030,7 @@ void extrinsic_burnin_cbc(struct bayesCBC *bayescbc, struct Data *data, struct M ...@@ -3015,7 +3030,7 @@ void extrinsic_burnin_cbc(struct bayesCBC *bayescbc, struct Data *data, struct M
HH = double_tensor(NC,net->Nifo,3); HH = double_tensor(NC,net->Nifo,3);
printf("\nFinding sky location\n"); if(bayescbc->debug == 1) printf("\nFinding sky location\n");
// // Find sky location and set up whitend signal arrays // // Find sky location and set up whitend signal arrays
#pragma omp parallel for #pragma omp parallel for
...@@ -3050,7 +3065,7 @@ void extrinsic_burnin_cbc(struct bayesCBC *bayescbc, struct Data *data, struct M ...@@ -3050,7 +3065,7 @@ void extrinsic_burnin_cbc(struct bayesCBC *bayescbc, struct Data *data, struct M
mty += SNRsq[mtid]; mty += SNRsq[mtid];
} }
printf("%f %f\n", SNRsq[0], SNRsq[1]); if(bayescbc->debug == 1) printf("%f %f\n", SNRsq[0], SNRsq[1]);
// find a sky location roughly consistent with the time delays // 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); skyring(net, paramx[j], skyx[j], pallx[j], SNRsq, bayescbc->freq, D, SN, N, data->Tobs, bayescbc->rvec[j], NX, NS, bayescbc->gmst);
...@@ -3083,7 +3098,7 @@ void extrinsic_burnin_cbc(struct bayesCBC *bayescbc, struct Data *data, struct M ...@@ -3083,7 +3098,7 @@ void extrinsic_burnin_cbc(struct bayesCBC *bayescbc, struct Data *data, struct M
logLsky[j] = skylike(net, skyx[j], DD, WW[j], DHc[j], DHs[j], dtx, nt, 0, bayescbc); 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); 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]); if(bayescbc->debug == 1) printf("%d logL full %f sky %f\n", j, logLfull[j], logLsky[j]);
free_double_vector(SNRsq); free_double_vector(SNRsq);
free_double_matrix(hwave,net->Nifo); free_double_matrix(hwave,net->Nifo);
...@@ -3106,16 +3121,6 @@ void extrinsic_burnin_cbc(struct bayesCBC *bayescbc, struct Data *data, struct M ...@@ -3106,16 +3121,6 @@ void extrinsic_burnin_cbc(struct bayesCBC *bayescbc, struct Data *data, struct M
logLfull[j] = log_likelihood_full(net, D, pallx[j], bayescbc->freq, SN, rho, N, data->Tobs, bayescbc); 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; x = 0.0;
k = bayescbc->who[0]; k = bayescbc->who[0];
...@@ -3131,6 +3136,8 @@ void extrinsic_burnin_cbc(struct bayesCBC *bayescbc, struct Data *data, struct M ...@@ -3131,6 +3136,8 @@ void extrinsic_burnin_cbc(struct bayesCBC *bayescbc, struct Data *data, struct M
// set threshold that all chains start withing 20% of highest likelihood // set threshold that all chains start withing 20% of highest likelihood
x *= 0.8; x *= 0.8;
if (bayescbc->debug) printf(" Chains CBC logL before starting main MCMC:\n");
for(j = 0; j < NC; j++) for(j = 0; j < NC; j++)
{ {
...@@ -3143,11 +3150,10 @@ void extrinsic_burnin_cbc(struct bayesCBC *bayescbc, struct Data *data, struct M ...@@ -3143,11 +3150,10 @@ void extrinsic_burnin_cbc(struct bayesCBC *bayescbc, struct Data *data, struct M
pallx[j][i] = pallx[k][i]; pallx[j][i] = pallx[k][i];
} }
} }
logLfull[j] = log_likelihood_full(net, D, pallx[j], bayescbc->freq, SN, rho, N, data->Tobs, bayescbc); 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 bayescbc->logLx[j] = logLfull[j]; // BW shared array
printf("%f\n", logLfull[j]); if(bayescbc->debug) printf(" %d logL full %f\n", j, logLfull[j]);
} }
free_double_vector(heat); free_double_vector(heat);
......
...@@ -142,7 +142,7 @@ int extrinsic_checkrange(double *p) ...@@ -142,7 +142,7 @@ int extrinsic_checkrange(double *p)
} }
//TODO: ellipticity should be periodic? //TODO: ellipticity should be periodic?
if(p[3] < -0.99 || p[3] >= 0.99) return 1; if(p[3] <= -1.00 || p[3] >= 1.00) return 1;
//TODO: phi from [0,2pi], or pi? //TODO: phi from [0,2pi], or pi?
while(p[4]<0.0 || p[4] >= LAL_TWOPI) while(p[4]<0.0 || p[4] >= LAL_TWOPI)
......
...@@ -205,10 +205,11 @@ void extrinsic_uniform_proposal(gsl_rng *seed, double *y) ...@@ -205,10 +205,11 @@ void extrinsic_uniform_proposal(gsl_rng *seed, double *y)
y[0] = uniform_draw(seed)*LAL_TWOPI; y[0] = uniform_draw(seed)*LAL_TWOPI;
y[1] = -1.0 + 2.0*uniform_draw(seed); y[1] = -1.0 + 2.0*uniform_draw(seed);
y[2] = uniform_draw(seed)*LAL_PI; y[2] = uniform_draw(seed)*LAL_PI;
y[3] = -0.99 + 1.98*uniform_draw(seed); y[3] = -1.0 + 2.0*uniform_draw(seed);
y[4] = uniform_draw(seed)*LAL_TWOPI; y[4] = uniform_draw(seed)*LAL_TWOPI;
//y[4] = uniform_draw(seed)*LAL_PI; //y[4] = uniform_draw(seed)*LAL_PI;
y[5] = 1.0; y[5] = 1.0;
} }
void uniform_orientation_proposal(double *y, gsl_rng *seed) void uniform_orientation_proposal(double *y, gsl_rng *seed)
...@@ -327,7 +328,7 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa ...@@ -327,7 +328,7 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa
range = LAL_PI; range = LAL_PI;
break; break;
case 3: case 3:
range = 0.99*2; range = 1.0*2;
break; break;
case 4: case 4:
range = LAL_TWOPI; range = LAL_TWOPI;
......
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