diff --git a/src/BayesWave.c b/src/BayesWave.c
index 64f3ab21a741454f217a4f007de6529d4dc54995..8728fe924297ca14a1c6cd499fdd022ff7d68dac 100644
--- a/src/BayesWave.c
+++ b/src/BayesWave.c
@@ -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");
diff --git a/src/BayesWave.h b/src/BayesWave.h
index 41c939fb197423e3daad08e7c34c819eca5b018f..d35ac6c9b26d8343c3a57b72b401f147b8a5c7c9 100644
--- a/src/BayesWave.h
+++ b/src/BayesWave.h
@@ -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;
diff --git a/src/BayesWaveIO.c b/src/BayesWaveIO.c
index c72879514348f1a4e9d9bb547d47e7a571f4a966..9b043fb4362017ed7a1a4ca18dd2693542f7a103 100644
--- a/src/BayesWaveIO.c
+++ b/src/BayesWaveIO.c
@@ -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;
diff --git a/src/BayesWaveMCMC.c b/src/BayesWaveMCMC.c
index 2ebbd4dadf3c84982d730b29e861003065617143..55228f6a7c7e94e454feb0bc43611f8174a4c444 100644
--- a/src/BayesWaveMCMC.c
+++ b/src/BayesWaveMCMC.c
@@ -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++)
diff --git a/src/BayesWaveModel.c b/src/BayesWaveModel.c
index ce24e8ad3c9ac7679d301c45ea47fd922491556e..3ce93ba5d41d2a8937151384ecc8290af837794e 100644
--- a/src/BayesWaveModel.c
+++ b/src/BayesWaveModel.c
@@ -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
diff --git a/src/BayesWavePrior.c b/src/BayesWavePrior.c
index 8421dd5b7d65d0082e7dee4958b2aa9e8853c899..acde28f633832b8e68a16b27f5317aa45ae98ec4 100644
--- a/src/BayesWavePrior.c
+++ b/src/BayesWavePrior.c
@@ -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?
diff --git a/src/BayesWaveProposal.c b/src/BayesWaveProposal.c
index a2ced926414ea9d19b6ee0dec6f86597bd830b77..8b1a09ba627d6e8146bc2a280a583050219e3cb0 100644
--- a/src/BayesWaveProposal.c
+++ b/src/BayesWaveProposal.c
@@ -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 */
diff --git a/src/BayesWaveProposal.h b/src/BayesWaveProposal.h
index 85915d3c088d690066783b62bcb79f412da8f924..94e8e25c3afd0632d1329eb42aa78cf8b01f8dcb 100644
--- a/src/BayesWaveProposal.h
+++ b/src/BayesWaveProposal.h
@@ -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);