Commit 7f5444f2 authored by Tyson Littenberg's avatar Tyson Littenberg

Merge branch 'Fisher' into 'master'

Restore proposal density in Fisher jumps

See merge request !177
parents f9fef287 635ba4b5
Pipeline #140336 passed with stages
in 2 minutes and 17 seconds
......@@ -1694,6 +1694,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
int acc;
double alpha;
double q,logqx,logqy;
double logqxAp, logqyAp, logqxtfQ, logqytfQ;//proposal density for the Fisher proposal, need to split it for the relaxed polarization case
double logH;
//indicies for looping over the number of states being update (1 IFO for glitch model, Npol for signal model)
......@@ -1837,7 +1838,6 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
rj = 0;
test = 0;
/*******************************************************/
/* */
/* TRANS DIMENSION RJMCMC MOVE */
......@@ -2216,6 +2216,9 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
{
//give hx the same phase shift as h+
wave_y[1]->intParams[ii][4] = wave_x[1]->intParams[ii][4] + (wave_y[0]->intParams[ii][4]-wave_x[0]->intParams[ii][4]);
while(wave_y[1]->intParams[ii][4] > LAL_TWOPI) wave_y[1]->intParams[ii][4] -=LAL_TWOPI;
while(wave_y[1]->intParams[ii][4] < LAL_TWOPI) wave_y[1]->intParams[ii][4] +=LAL_TWOPI;
constrain_hplus_hcross(wave_y, ii);
}
......@@ -2244,17 +2247,28 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if(det==-1)
{
intrinsic_fisher_proposal(model_x, range, paramsx[ii], paramsy[ii], model_x->SnGeo, &logqx, &logqy, seed, data->Tobs, &test, wave_y[n]->dimension, tfqFlag);
//rely on (mostly) symmetry of proposal
// -proposal got tricky for h+,hx case because
// t,f,Q are fixed to h+
logqx=logqy=0.0;
intrinsic_fisher_proposal(model_x, range, paramsx[ii], paramsy[ii], model_x->SnGeo, &logqxAp, &logqyAp, &logqxtfQ, &logqytfQ, seed, data->Tobs, &test, wave_y[n]->dimension, tfqFlag);
// t, f, Q are shared between hplus and hcross, so they should be counted only once
if (n==0) {
logqx += logqxAp;
logqx += logqxtfQ;
logqy += logqyAp;
logqy += logqytfQ;
}
else {
logqx += logqxAp;
logqy += logqyAp;
}
}
else
{
intrinsic_fisher_proposal(model_x, range, paramsx[ii], paramsy[ii], model_x->Snf[det], &logqx, &logqy, seed, data->Tobs, &test, wave_y[n]->dimension, tfqFlag);
intrinsic_fisher_proposal(model_x, range, paramsx[ii], paramsy[ii], model_x->Snf[det], &logqxAp, &logqyAp, &logqxtfQ, &logqytfQ, seed, data->Tobs, &test, wave_y[n]->dimension, tfqFlag);
logqx += logqxAp;
logqx += logqxtfQ;
logqy += logqyAp;
logqy += logqytfQ;
}
}
if(det==-1 && model_y->Npol > 1) constrain_hplus_hcross(wave_y, ii);
......
......@@ -629,12 +629,14 @@ void intrinsic_fisher_matrix(double *params, struct FisherMatrix *fisher, int NW
}
void intrinsic_fisher_proposal(UNUSED struct Model *model, double **range, double *paramsx, double *paramsy, double *Snf, double *logpx, double *logpy, gsl_rng *seed, double Tobs, int *test, int NW, int TFQFLAG)
void intrinsic_fisher_proposal(UNUSED struct Model *model, double **range, double *paramsx, double *paramsy, double *Snf, double *logpxAp, double *logpyAp, double *logpxtfQ, double *logpytfQ, gsl_rng *seed, double Tobs, int *test, int NW, int TFQFLAG)
{
int k;
double x,y;
double qyx = *logpy;
double qxy = *logpx;
double qyxAp = 0.;
double qxyAp = 0.;
double qyxtfQ = 0.;
double qxytfQ = 0.;
double *dparamsx = double_vector(NW-1);
double *dparamsy = double_vector(NW-1);
double scale = 1./sqrt((double)NW);//Number of intrinsic parameters
......@@ -659,36 +661,62 @@ void intrinsic_fisher_proposal(UNUSED struct Model *model, double **range, doubl
return;
}
//compute probabiltiy of such a proposal (could just product the draws from gaussian_draw() above
//compute probability of such a proposal (could just product the draws from gaussian_draw() above
//TODO: Does that break chirplets and non elliptical polarization?
// A, Phi
y = 1.0;
for(k=0; k<NW; k++)
for(k=3; k<NW; k++)
{
x = (paramsx[k]-paramsy[k])/(scale*dparamsx[k]);
y *= scale*dparamsx[k];
qyxAp += -0.5*x*x;
}
qyxAp -= log(y);
// t, f, Q
y = 1.0;
for(k=0; k<NW-2; k++)
{
x = (paramsx[k]-paramsy[k])/(scale*dparamsx[k]);
y *= scale*dparamsx[k];
qyx += -0.5*x*x;
qyxtfQ += -0.5*x*x;
}
qyx -= log(y);
qyxtfQ -= log(y);
//need fisher at new location to compute reverse probability
intrinsic_fisher_update(paramsy,dparamsy,Snf,Tobs, NW, TFQFLAG);
//compute probabiltiy of reverse proposal
// A, Phi
y = 1.0;
for(k=0; k<NW; k++)
for(k=3; k<NW; k++)
{
x = (paramsx[k]-paramsy[k])/(scale*dparamsy[k]);
y *= scale*dparamsy[k];
qxyAp+= -0.5*x*x;
}
qxyAp -= log(y);
// t, f, Q
y = 1.0;
for(k=0; k<NW-2; k++)
{
x = (paramsx[k]-paramsy[k])/(scale*dparamsy[k]);
y *= scale*dparamsy[k];
qxy += -0.5*x*x;
qxytfQ+= -0.5*x*x;
}
qxy -= log(y);
qxytfQ -= log(y);
// map phase onto [0,2pi]
while(paramsy[4] > LAL_TWOPI) paramsy[4] -=LAL_TWOPI;
while(paramsy[4] < 0.0) paramsy[4] +=LAL_TWOPI;
//update incoming pointers and return
*logpy += qyx;
*logpx += qxy;
*logpyAp = qyxAp;
*logpxAp = qxyAp;
*logpytfQ = qyxtfQ;
*logpxtfQ = qxytfQ;
free_double_vector(dparamsx);
free_double_vector(dparamsy);
......
......@@ -53,7 +53,7 @@ void extrinsic_fisher_update(struct Data *data, struct Model *model);
void intrinsic_fisher_update(double *params, double *dparams, double *Snf, double Tobs, int NW, int TFQFLAG);
void intrinsic_fisher_proposal(UNUSED struct Model *model, double **range, double *paramsx, double *paramsy, double *Snf, double *logpx, double *logpy, gsl_rng *seed, double Tobs, int *test, int NW, int TFQFLAG);
void intrinsic_fisher_proposal(UNUSED struct Model *model, double **range, double *paramsx, double *paramsy, double *Snf, double *logpxAp, double *logpyAp, double *logpxtfQ, double *logpytfQ, gsl_rng *seed, double Tobs, int *test, int NW, int TFQFLAG);
/* ********************************************************************************** */
/* */
......
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