Commit 4305fe92 authored by Meg Millhouse's avatar Meg Millhouse Committed by Katerina Chatziioannou
Browse files

Fixed sky location functionality, and bug fix for extrinsic parameter proposal

parent 77092c68
......@@ -1169,6 +1169,9 @@ void print_help_message(void)
fprintf(stdout," --fixIntrinsicParams hold intrinsic parameters fixed\n");
fprintf(stdout," --fixExtrinsicParams hold extrinsic parameters fixed\n");
fprintf(stdout," --loudGlitchPrior set floor of PSD prior to .01*Sn(f=200Hz)\n");
fprintf(stdout," --fixSky hold sky locaiton fixed\n");
fprintf(stdout," --fixRA right ascension in radians of fixed sky location\n");
fprintf(stdout," --fixDEC declination in radians of fixed sky location\n");
fprintf(stdout,"\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --- Parallel Tempering parameters ----------------------------------------------\n");
......
......@@ -292,6 +292,8 @@ struct Data
int glitchFlag;
int signalFlag;
int noiseFlag;
int skyFixFlag;
int polarizationFlag;
int fixIntrinsicFlag;
......@@ -321,10 +323,10 @@ struct Data
int nukeFlag; //tell BWP to remove chains/ and checkpoint/ directories
int srate;
int chirpletFlag;
int NW;
int srate;
int chirpletFlag;
int NW;
double dt;
double df;
......@@ -337,6 +339,9 @@ struct Data
double gmst;
double logLc;
double TwoDeltaToverN;
double FIX_RA;
double FIX_DEC;
double TFtoProx;
......
......@@ -1405,6 +1405,15 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
ppt = LALInferenceGetProcParamVal(commandLine, "--tempSpacing");
if(ppt) chain->tempStep = (double)atof(ppt->value);
else chain->tempStep = 1.2;
ppt = LALInferenceGetProcParamVal(commandLine, "--fixRA");
if(ppt) data->FIX_RA = (double)atof(ppt->value);
else data->FIX_RA = 0.0;
ppt = LALInferenceGetProcParamVal(commandLine, "--fixDEC");
if(ppt) data->FIX_DEC = (double)atof(ppt->value);
else data->FIX_DEC = 0.0;
ppt = LALInferenceGetProcParamVal(commandLine, "--noAdaptTemperature");
if(ppt) data->adaptTemperatureFlag = 0;
......@@ -1476,6 +1485,10 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
data->fullModelFlag = 0;
}
else data->stochasticFlag = 0;
ppt = LALInferenceGetProcParamVal(commandLine, "--fixSky");
if(ppt) data->skyFixFlag = 1;
else data->skyFixFlag = 0;
ppt = LALInferenceGetProcParamVal(commandLine, "--clusterPrior");
if(ppt) data->clusterPriorFlag = 1;
......
......@@ -180,7 +180,15 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
if(data->signalFlag)
{
extrinsic_uniform_proposal(chain->seed,model[ic]->extParams);
if(data->skyFixFlag == 1)
{
model[ic]->extParams[0] = data->FIX_RA;
model[ic]->extParams[1] = sin(data->FIX_DEC);
}
/*
WaveformProject() was split into two functions:
......@@ -2417,7 +2425,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
intParams = double_matrix(dim,smodel[0]->dimension-1);
//Set extrinsic parameter fisher matrix (numerical derivatives of response to geocenter waveform)
extrinsic_fisher_update(data, model_x);
if(data->skyFixFlag == 0) extrinsic_fisher_update(data, model_x);
glitch = malloc(NI*sizeof(double*));
for(ifo=0; ifo<NI; ifo++) glitch[ifo] = model_x->glitch[ifo]->templates;
......@@ -2493,29 +2501,40 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
Don't bother for hot chains, which should largely be
sampling from the prior
*/
if(uniform_draw(seed)< 0.5 && chain->index[ic]<10 && NI>1)
if(data->skyFixFlag == 0) // standard mode
{
if(ic==0)cnt++;
if(ic==0)sky=1;
sky_ring_proposal(paramsx,paramsy,data,seed);
if(data->orientationProposalFlag && uniform_draw(seed)<0.5)
network_orientation_proposal(paramsx, paramsy, data, &logJ);
if(uniform_draw(seed)< 0.5 && chain->index[ic]<10 && NI>1)
{
if(ic==0)cnt++;
if(ic==0)sky=1;
sky_ring_proposal(paramsx,paramsy,data,seed);
if(data->orientationProposalFlag && uniform_draw(seed)<0.5){
network_orientation_proposal(paramsx, paramsy, data, &logJ);
}
else{
uniform_orientation_proposal(paramsy, seed);
}
}
else
uniform_orientation_proposal(paramsy, seed);
{
/* Uniform Draw */
draw = uniform_draw(seed);
if(draw<0.1)extrinsic_uniform_proposal(seed,paramsy);
/* Fisher Matrix Proposal */
else
{
fisher_matrix_proposal(fisher, paramsx, seed, paramsy);
}
}
if(!data->extrinsicAmplitudeFlag) paramsy[5]=1.0;
}
else
else // Fixed sky proposal
{
/* Uniform Draw */
draw = uniform_draw(seed);
if(draw<0.1)extrinsic_uniform_proposal(seed,paramsy);
sky_fix_proposal(paramsx,paramsy,data,seed);
if(!data->extrinsicAmplitudeFlag) paramsy[5]=1.0;
/* Fisher Matrix Proposal */
else
{
fisher_matrix_proposal(fisher, paramsx, seed, paramsy);
}
}
if(!data->extrinsicAmplitudeFlag) paramsy[5]=1.0;
//Rescale intrinsic amplitudes and phases with extrinsic parameters
for(i=1; i<ienddim; i++)
......
......@@ -742,6 +742,13 @@ void initialize_model(struct Model *model, struct Data *data, struct Prior *prio
//draw extrinsic parameters from prior
extrinsic_uniform_proposal(seed,model->extParams);
if(data->skyFixFlag == 1)
{
model->extParams[0] = data->FIX_RA;
model->extParams[1] = sin(data->FIX_DEC);
}
//set "delta" parameters to not modify signal model
model->extParams[4] = 0.0; //geocenter phase shift
......
......@@ -135,6 +135,8 @@ int extrinsic_checkrange(double *p)
{
if(p[2]<0.0) p[2] += LAL_PI;
else p[2] -= LAL_PI;
}
//TODO: ellipticity should be periodic?
......
......@@ -203,7 +203,7 @@ void extrinsic_uniform_proposal(gsl_rng *seed, double *y)
{
y[0] = uniform_draw(seed)*LAL_TWOPI;
y[1] = -1.0 + 2.0*uniform_draw(seed);
y[2] = uniform_draw(seed)*LAL_PI_2;
y[2] = uniform_draw(seed)*LAL_PI;
y[3] = -0.99 + 1.98*uniform_draw(seed);
y[4] = uniform_draw(seed)*LAL_TWOPI;
//y[4] = uniform_draw(seed)*LAL_PI;
......@@ -323,7 +323,7 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa
range = 2.0;
break;
case 2:
range = LAL_PI_2;
range = LAL_PI;
break;
case 3:
range = 0.99*2;
......@@ -1393,6 +1393,44 @@ void TFprop_signal(struct Data *data, double **range, struct TimeFrequencyMap *t
/* */
/* ********************************************************************************** */
void sky_fix_proposal(double *x, double *y, struct Data *data, gsl_rng *seed)
{
double scale;
double draw;
// The fixed sky location. Assumes that RA and DEC have been provided in radians
y[0] = data->FIX_RA;
y[1] = sin(data->FIX_DEC);
draw = uniform_draw(seed);
if(draw > 0.8)
{
uniform_orientation_proposal(y, seed);
}
else
{
scale = exp(-9.0+7.0*uniform_draw(seed)); // yields scales between ~ 1e-4 and 1e-1
y[2] = x[2] + gaussian_draw(seed)*scale*LAL_PI;
y[3] = x[3] + gaussian_draw(seed)*scale*2.0;
y[4] = x[4] + gaussian_draw(seed)*scale*LAL_TWOPI;
// apply period BC
while(y[2] < 0.0) y[2] += LAL_PI;
while(y[2] > LAL_PI) y[2] -= LAL_PI;
while(y[3] < -1.0) y[3] += 2.0;
while(y[3] > 1.0) y[3] -= 2.0;
while(y[4] < 0.0) y[4] += LAL_TWOPI;
while(y[4] > LAL_TWOPI) y[4] -= LAL_TWOPI;
}
}
void sky_ring_proposal(double *x, double *y, struct Data *data, gsl_rng *seed)
{
......@@ -1773,6 +1811,8 @@ void network_orientation_proposal(double *paramsx, double *paramsy, struct Data
}
/* ********************************************************************************** */
/* */
/* Stochastic background proposals */
......
......@@ -91,6 +91,8 @@ void TFprop_signal(struct Data *data, double **range, struct TimeFrequencyMap *t
void sky_ring_proposal(double *x, double *y, struct Data *data, gsl_rng *seed);
void sky_fix_proposal(double *x, double *y, struct Data *data, gsl_rng *seed);
void uniform_orientation_proposal(double *y, gsl_rng *seed);
void network_orientation_proposal(double *paramsx, double *paramsy, struct Data *data, double *logJ);
......
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