Skip to content

Commit

Permalink
fix to cdf proposal to guard against repeated samples
Browse files Browse the repository at this point in the history
  • Loading branch information
Tyson Littenberg committed Nov 13, 2018
1 parent 5a0719f commit 12c0e35
Show file tree
Hide file tree
Showing 5 changed files with 33 additions and 29 deletions.
2 changes: 1 addition & 1 deletion galactic_binaries/src/GalacticBinary.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -69,7 +69,7 @@ struct Flags
int DMAX; //max dimension of signal model
int zeroNoise;
int fixSky;
int skyPrior;
int galaxyPrior;
int snrPrior;
int emPrior;
int knownSource;
Expand Down
10 changes: 5 additions & 5 deletions galactic_binaries/src/GalacticBinaryIO.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@ void print_usage()
fprintf(stdout," --data : strain data file \n");
fprintf(stdout," --frac-freq : fractional frequency data (phase) \n");
fprintf(stdout," --fix-sky : pin sky params to injection \n");
fprintf(stdout," --sky-prior : use galaxy model for sky prior \n");
fprintf(stdout," --galaxy-prior: use galaxy model for sky prior \n");
fprintf(stdout," --snr-prior : use SNR-based amplitude prior \n");
fprintf(stdout," --em-prior : update prior ranges from other obs \n");
fprintf(stdout," --known-source: injection is VB (draw orientation) \n");
Expand Down Expand Up @@ -129,7 +129,7 @@ void parse(int argc, char **argv, struct Data **data, struct Orbit *orbit, struc
flags->zeroNoise = 0;
flags->confNoise = 0;
flags->fixSky = 0;
flags->skyPrior = 0;
flags->galaxyPrior = 0;
flags->snrPrior = 0;
flags->emPrior = 0;
flags->cheat = 0;
Expand Down Expand Up @@ -212,7 +212,7 @@ void parse(int argc, char **argv, struct Data **data, struct Orbit *orbit, struc
{"conf-noise", no_argument, 0, 0 },
{"frac-freq", no_argument, 0, 0 },
{"fix-sky", no_argument, 0, 0 },
{"sky-prior", no_argument, 0, 0 },
{"galaxy-prior",no_argument, 0, 0 },
{"snr-prior", no_argument, 0, 0 },
{"known-source",no_argument, 0, 0 },
{"f-double-dot",no_argument, 0, 0 },
Expand Down Expand Up @@ -250,7 +250,7 @@ void parse(int argc, char **argv, struct Data **data, struct Orbit *orbit, struc
if(strcmp("zero-noise", long_options[long_index].name) == 0) flags->zeroNoise = 1;
if(strcmp("conf-noise", long_options[long_index].name) == 0) flags->confNoise = 1;
if(strcmp("fix-sky", long_options[long_index].name) == 0) flags->fixSky = 1;
if(strcmp("sky-prior", long_options[long_index].name) == 0) flags->skyPrior = 1;
if(strcmp("galaxy-prior",long_options[long_index].name) == 0) flags->galaxyPrior= 1;
if(strcmp("snr-prior", long_options[long_index].name) == 0) flags->snrPrior = 1;
if(strcmp("prior", long_options[long_index].name) == 0) flags->prior = 1;
if(strcmp("f-double-dot",long_options[long_index].name) == 0) data_ptr->NP = 9;
Expand Down Expand Up @@ -447,7 +447,7 @@ void parse(int argc, char **argv, struct Data **data, struct Orbit *orbit, struc
else fprintf(stdout," Sky parameters are... ENABLED\n");
if(flags->calibration) fprintf(stdout," Calibration is....... ENABLED\n");
else fprintf(stdout," Calibration is....... DISABLED\n");
if(flags->skyPrior) fprintf(stdout," Galaxy prior is ..... ENABLED\n");
if(flags->galaxyPrior) fprintf(stdout," Galaxy prior is ..... ENABLED\n");
else fprintf(stdout," Galaxy prior is ..... DISABLED\n");
if(flags->snrPrior) fprintf(stdout," SNR prior is ........ ENABLED\n");
else fprintf(stdout," SNR prior is ........ DISABLED\n");
Expand Down
40 changes: 21 additions & 19 deletions galactic_binaries/src/GalacticBinaryMCMC.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,7 +122,7 @@ int main(int argc, char *argv[])

/* Initialize priors */
struct Prior *prior = malloc(sizeof(struct Prior));
if(flags->skyPrior) set_galaxy_prior(flags, prior);
if(flags->galaxyPrior) set_galaxy_prior(flags, prior);


/* Initialize data models */
Expand Down Expand Up @@ -217,13 +217,15 @@ int main(int argc, char *argv[])
else flags->burnin=0;

//set annealinging tempurature during burnin
// if(flags->burnin)
// {
// chain->annealing = data[0]->SNR2*pow(data[0]->SNR2,-((double)mcmc+(double)flags->NBURN)/((double)flags->NBURN/(double)10))/40.;
// if(chain->annealing<1.0)chain->annealing=1.0;
// chain->annealing=1.0;
// //printf("annealing=%g\n",chain->annealing);
// }
/*
if(flags->burnin)
{
chain->annealing = data[0]->SNR2*pow(data[0]->SNR2,-((double)mcmc+(double)flags->NBURN)/((double)flags->NBURN/(double)10))/40.;
if(chain->annealing<1.0)chain->annealing=1.0;
chain->annealing=1.0;
//printf("annealing=%g\n",chain->annealing);
}
*/
chain->annealing=1.0;

// (parallel) loop over chains
Expand Down Expand Up @@ -252,7 +254,9 @@ int main(int argc, char *argv[])
if(flags->rj)galactic_binary_rjmcmc(orbit, data_ptr, model_ptr, trial_ptr, chain, flags, prior, proposal[i], ic);

//delayed rejection mode-hopper
//if(model_ptr->Nlive>0 && mcmc<0 && ic<NC/2)galactic_binary_drmc(orbit, data_ptr, model_ptr, trial_ptr, chain, flags, prior, proposal[i], ic);
/*
if(model_ptr->Nlive>0 && mcmc<0 && ic<NC/2)galactic_binary_drmc(orbit, data_ptr, model_ptr, trial_ptr, chain, flags, prior, proposal[i], ic);
*/

//update fisher matrix for each chain
if(mcmc%100==0)
Expand Down Expand Up @@ -589,6 +593,7 @@ void galactic_binary_mcmc(struct Orbit *orbit, struct Data *data, struct Model *
generate_calibration_model(data, model_y);
apply_calibration_model(data, model_y);
}

//get likelihood for y
model_y->logL = gaussian_log_likelihood(orbit, data, model_y);

Expand All @@ -609,15 +614,12 @@ void galactic_binary_mcmc(struct Orbit *orbit, struct Data *data, struct Model *
loga = log(gsl_rng_uniform(chain->r[ic]));
if(logH > loga)
{
if(!strcmp(proposal[nprop]->name,"cdf draw") && ic==0 && model_y->logL - model_x->logL < -20.)
{
printf("cdf logH=%g, logLx=%g, logLy=%g\n",logH,model_x->logL , model_y->logL);
printf(" dlogQ=%g, logQxy=%g, logQyx=%g\n",logQxy - logQyx,logQxy,logQyx);
exit(1);
}



if(!strcmp(proposal[nprop]->name,"cdf draw") && ic==0 && model_y->logL - model_x->logL < -20.)
{
printf("cdf logH=%g, logLx=%g, logLy=%g\n",logH,model_x->logL , model_y->logL);
printf(" dlogQ=%g, logQxy=%g, logQyx=%g\n",logQxy - logQyx,logQxy,logQyx);
//exit(1);
}
proposal[nprop]->accept[ic]++;
copy_model(model_y,model_x);
}
Expand Down Expand Up @@ -661,7 +663,7 @@ void galactic_binary_rjmcmc(struct Orbit *orbit, struct Data *data, struct Model
else
{
draw_from_prior(data, model_y, model_y->source[create], proposal[0], model_y->source[create]->params, chain->r[ic]);
if(flags->skyPrior) draw_from_galaxy_prior(model_y, prior, model_y->source[create]->params, chain->r[ic]);
if(flags->galaxyPrior) draw_from_galaxy_prior(model_y, prior, model_y->source[create]->params, chain->r[ic]);

logQyx = evaluate_prior(flags, data, model_y, prior, model_y->source[create]->params);
}
Expand Down
2 changes: 1 addition & 1 deletion galactic_binaries/src/GalacticBinaryPrior.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -464,7 +464,7 @@ double evaluate_prior(struct Flags *flags, struct Data *data, struct Model *mode
if(params[0]<uniform_prior[0][0] || params[0]>uniform_prior[0][1]) return -INFINITY;
else logP -= log(uniform_prior[0][1]-uniform_prior[0][0]);

if(flags->skyPrior)
if(flags->galaxyPrior)
{
if(params[1]<uniform_prior[1][0] || params[1]>uniform_prior[1][1]) return -INFINITY;

Expand Down
8 changes: 5 additions & 3 deletions galactic_binaries/src/GalacticBinaryProposal.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -513,20 +513,22 @@ double cdf_density(struct Model *model, struct Source *source, struct Proposal *
double logP=0.0;
double *params = source->params;

int i;
int i,j;

for(int n=0; n<NP; n++)
{
if(params[n]<model->prior[n][0] || params[n]>=model->prior[n][1]) return -INFINITY;

//find samples either end of p-value
if(params[n]<cdf[n][0] || params[n]>= cdf[n][N-1])
return -INFINITY;
else
{
i=0;
while(params[n]>cdf[n][i]) i++;
logP += log( (1./N) / (cdf[n][i+1]-cdf[n][i]) );
j=i+1;
while(cdf[n][j]==cdf[n][i]) j++;
logP += log( ((double)(j-i)/N) / (cdf[n][j]-cdf[n][i]) );
}
}

Expand Down

0 comments on commit 12c0e35

Please sign in to comment.