Commit 837d4168 authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

Create O2 tag


git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@609 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent 991bb378
......@@ -920,6 +920,7 @@ void print_help_message(void)
fprintf(stdout," --clusterAlpha distance between wavelets to be considered a cluster (2)\n");
fprintf(stdout," --clusterBeta cluster prior exp(-beta K) (4)\n");
fprintf(stdout," --clusterGamma occam penalty exp(gamma*J) (0)\n");
fprintf(stdout," --updateGeocenterPSD geocenter PSD depends on extrinsic parameters\n");
fprintf(stdout," --backgroundPrior name of 2-column bkg frequency distribution file\n");
fprintf(stdout," --noOrientationProposal disable MCMC proposal for psi/ecc\n");
fprintf(stdout," --uniformAmplitudePrior don't use SNR-dependent amplitude prior\n");
......
......@@ -1401,7 +1401,7 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
if(ppt) data->extrinsicAmplitudeFlag = 1;
else data->extrinsicAmplitudeFlag = 0;
ppt = LALInferenceGetProcParamVal(commandLine, "--fixGeocenterPSD");
ppt = LALInferenceGetProcParamVal(commandLine, "--updateGeocenterPSD");
if(ppt) data->geocenterPSDFlag = 1;
else data->geocenterPSDFlag = 0;
......
......@@ -2254,39 +2254,6 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
//Set extrinsic parameter fisher matrix (numerical derivatives of response to geocenter waveform)
extrinsic_fisher_update(data, model_x);
//=============================================
//TODO: Testing extrinsic proposal (alloc)
//=============================================
/*
char root[128];
double **hx,**hy;
hx = malloc(NI*sizeof(double *));
hy = malloc(NI*sizeof(double *));
for(ifo=0; ifo<NI; ifo++)
{
hx[ifo] = malloc(data->N*sizeof(double));
hy[ifo] = malloc(data->N*sizeof(double));
}
double hihi; //Detector (hinj|hinj)
double hrhr; //Detector (hrec|hrec)
double hrhi; //Detector (hrec|hinj)
double overlap;
double *fyplus, *fycross;
double *fxplus, *fxcross;
double x, y,p0;
fyplus = malloc(NI*sizeof(double));
fycross = malloc(NI*sizeof(double));
fxplus = malloc(NI*sizeof(double));
fxcross = malloc(NI*sizeof(double));
double Fp[2],Fx[2];
double Fp_test[2], Fx_test[2];
*/
//=============================================
glitch = malloc(NI*sizeof(double*));
for(ifo=0; ifo<NI; ifo++) glitch[ifo] = model_x->glitch[ifo]->templates;
......@@ -2294,8 +2261,6 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
for(i=0; i<NE; i++) paramsy[i] = paramsx[i];
//Find minimum and maximum frequencies of signal model and only compute logL over that range
// double fmin = data->fmin;
// double fmax = data->fmax;
double fmin = data->fmax;
double fmax = data->fmin;
double fi,ff;
......@@ -2323,7 +2288,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if(data->amplitudePriorFlag)
{
//model_x->projection has be updated by data->extrinsic_likelihood, so it stores F+ and Fx for params_y.
Shf_Geocenter(data, model_x, SnGeox, paramsx);
Shf_Geocenter_full(data, model_x->projection, model_x->Snf, model_x->SnGeo, model_x->extParams);
for(i=1; i<ienddim; i++)
{
logpx += (data->signal_amplitude_prior(smodel->intParams[smodel->index[i]],SnGeox, 1.0, data->Tobs, prior->sSNRpeak));
......@@ -2381,6 +2346,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
}
}
if(!data->extrinsicAmplitudeFlag) paramsy[5]=1.0;
//else paramsy[5] = 1.0 + 0.1*gasdev2(seed);
//Rescale intrinsic amplitudes and phases with extrinsic parameters
for(i=1; i<ienddim; i++)
......@@ -2397,130 +2363,12 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if(data->amplitudePriorFlag)
{
//model_x->projection has be updated by data->extrinsic_likelihood, so it stores F+ and Fx for params_y.
Shf_Geocenter(data, model_x, SnGeoy, paramsy);
if(data->geocenterPSDFlag)Shf_Geocenter(data, model_x, SnGeoy, paramsy);
for(i=1; i<ienddim; i++) logpy += (data->signal_amplitude_prior(intParams[i],SnGeoy, 1.0,data->Tobs, prior->sSNRpeak));
}
//logpy += dim*log(paramsy[5])/2.;
logH = (logLy - logLx)*chain->beta + logpy - logpx + logJ;
//=============================================
//TODO: Testing extrinsic proposal (print)
//=============================================
/*
if(flag==1)
{
FILE *temp = fopen("temp.dat","a");
overlap=0.0;
for(i=0;i<NE;i++)fprintf(temp,"%lg ",paramsx[i]);
for(i=0;i<NE;i++)fprintf(temp,"%lg ",paramsy[i]);
sprintf(root,"temp_x.dat");
dump_time_domain_waveform(data, model_x, prior, root);
for(ifo=0; ifo<NI; ifo++) for(i=0; i<data->N; i++)hx[ifo][i] = model_x->response[ifo][i];
//paramsx
for(i=0; i< NI; i++)XLALComputeDetAMResponse(&Fp[i], &Fx[i], (const REAL4(*)[3])data->detector[i]->response, paramsx[0], asin(paramsx[1]), paramsx[2], XLALGreenwichMeanSiderealTime(&data->epoch));
// Get detector frame orientation for current location
for(i=0; i< NI; i++)
{
XLALComputeDetAMResponse(&x, &y, (const REAL4(*)[3])data->detector[i]->response, paramsx[0], asin(paramsx[1]), 0.0, XLALGreenwichMeanSiderealTime(&data->epoch));
p0=0.5*atan2(y,x);
XLALComputeDetAMResponse(&x, &y, (const REAL4(*)[3])data->detector[i]->response, paramsx[0], asin(paramsx[1]), p0, XLALGreenwichMeanSiderealTime(&data->epoch));
fxplus[i] = x*cos(2.0*p0);
fxcross[i] = fxplus[i]*tan(2.0*p0);
Fp_test[i] = fxplus[i]*cos(2.0*paramsx[2]) + fxcross[i]*sin(2.0*paramsx[2]);
Fx_test[i] = -fxplus[i]*sin(2.0*paramsx[2]) + fxcross[i]*cos(2.0*paramsx[2]);
}
for(i=0; i<NI; i++)
{
fprintf(temp,"%lg %lg ",fxplus[i],fxcross[i]);
}
for(i=0; i<NI; i++)
{
fprintf(temp,"%lg %lg ",Fp[i],Fx[i]);
}
for(i=0; i<NI; i++)
{
fprintf(temp,"%lg %lg ",Fp_test[i],Fx_test[i]);
}
for(i=0; i<NI; i++)
{
printf("x IFO %i\n",i);
printf(" True F+,Fx = %g,%g\n",Fp[i],Fx[i]);
printf(" Test F+,Fx = %g,%g\n",Fp_test[i],Fx_test[i]);
}
//paramsy
for(i=0; i< NE; i++) paramsx[i] = paramsy[i];
sprintf(root,"temp_y.dat");
dump_time_domain_waveform(data, model_x, prior, root);
for(ifo=0; ifo<NI; ifo++) for(i=0; i<data->N; i++)hy[ifo][i] = model_x->response[ifo][i];
for(i=0; i< NI; i++)XLALComputeDetAMResponse(&Fp[i], &Fx[i], (const REAL4(*)[3])data->detector[i]->response, paramsy[0], asin(paramsy[1]), paramsy[2], XLALGreenwichMeanSiderealTime(&data->epoch));
// Get detector frame orientation for current location
for(i=0; i< NI; i++)
{
XLALComputeDetAMResponse(&x, &y, (const REAL4(*)[3])data->detector[i]->response, paramsy[0], asin(paramsy[1]), 0.0, XLALGreenwichMeanSiderealTime(&data->epoch));
p0=0.5*atan2(y,x);
XLALComputeDetAMResponse(&x, &y, (const REAL4(*)[3])data->detector[i]->response, paramsy[0], asin(paramsy[1]), p0, XLALGreenwichMeanSiderealTime(&data->epoch));
fxplus[i] = x*cos(2.0*p0);
fxcross[i] = fxplus[i]*tan(2.0*p0);
Fp_test[i] = fxplus[i]*cos(2.0*paramsy[2]) + fxcross[i]*sin(2.0*paramsy[2]);
Fx_test[i] = -fxplus[i]*sin(2.0*paramsy[2]) + fxcross[i]*cos(2.0*paramsy[2]);
}
for(i=0; i<NI; i++)
{
fprintf(temp,"%lg %lg ",fxplus[i],fxcross[i]);
}
for(i=0; i<NI; i++)
{
fprintf(temp,"%lg %lg ",Fp[i],Fx[i]);
}
for(i=0; i<NI; i++)
{
fprintf(temp,"%lg %lg ",Fp_test[i],Fx_test[i]);
}
for(i=0; i<NI; i++)
{
printf("y IFO %i\n",i);
printf(" True F+,Fx = %g,%g\n",Fp[i],Fx[i]);
printf(" Test F+,Fx = %g,%g\n",Fp_test[i],Fx_test[i]);
}
for(ifo=0; ifo<NI; ifo++)
{
hihi = fourier_nwip(data->imin, data->N/2, hx[ifo], hx[ifo], model_x->invSnf[ifo]);
hrhr = fourier_nwip(data->imin, data->N/2, hy[ifo], hy[ifo], model_x->invSnf[ifo]);
hrhi = fourier_nwip(data->imin, data->N/2, hx[ifo], hy[ifo], model_x->invSnf[ifo]);
//Store inj-rec overlap of each detector to file
overlap = hrhi/sqrt(hihi*hrhr); // (hrec|hinj)/sqrt((hrec|hrec)(hinj|hinj))
}
fprintf(temp,"%lg\n",overlap);
fflush(temp);
fclose(temp);
//system_pause();
}
*/
//=============================================
alpha = log(ran2(seed));
......@@ -2534,7 +2382,6 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
if(sky==1 && ic==0) acc++;
}
}
sky=0;
......@@ -2543,7 +2390,6 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
}
//Now update the full model with current extrinsic parameters
for(i=0; i<data->N; i++) geo[i] *= model_x->extParams[5];
......@@ -2563,26 +2409,6 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
free_dvector(paramsy,0,NE);
free_dvector(SnGeoy,0,iendN);
free_dmatrix(intParams,1,dim,0,4);
//=============================================
//TODO: Testing extrinsic proposal (free)
//=============================================
/*
for(ifo=0; ifo<NI; ifo++)
{
free(hx[ifo]);
free(hy[ifo]);
}
free(hx);
free(hy);
free(fyplus);
free(fycross);
free(fxplus);
free(fxcross);
*/
//=============================================
}
/* ********************************************************************************** */
......
......@@ -85,7 +85,8 @@ void Shf_Geocenter_full(struct Data *data, struct Network *projection, double **
for(ifo=0; ifo<NI; ifo++)
{
AntennaPattern = projection->Fplus[ifo]*projection->Fplus[ifo] + ecc*ecc*projection->Fcross[ifo]*projection->Fcross[ifo];
AntennaPattern = 1.0;
if(data->geocenterPSDFlag) AntennaPattern = projection->Fplus[ifo]*projection->Fplus[ifo] + ecc*ecc*projection->Fcross[ifo]*projection->Fcross[ifo];
for(i=imin; i<imax; i++) SnGeo[i] += AntennaPattern/(Snf[ifo][i]);
}
......@@ -94,29 +95,6 @@ void Shf_Geocenter_full(struct Data *data, struct Network *projection, double **
for(i=imin; i<imax; i++) SnGeo[i] = 1./SnGeo[i];
}
void Shf_Geocenter_inv(struct Data *data, struct Network *projection, double **invSnf, double *SnGeo, double *params)
{
int ifo,i, halfN=data->N/2;
int NI = data->NI;
int imin = data->imin;
int imax = data->imax;
double AntennaPattern;
double ecc = params[3];
for(i=0; i<halfN; i++) SnGeo[i] = 0.0;
for(ifo=0; ifo<NI; ifo++)
{
AntennaPattern = projection->Fplus[ifo]*projection->Fplus[ifo] + ecc*ecc*projection->Fcross[ifo]*projection->Fcross[ifo];
for(i=imin; i<imax; i++) SnGeo[i] += AntennaPattern*invSnf[ifo][i];
}
for(i=0; i<imin; i++)SnGeo[i]=1.0;
for(i=imax; i<halfN; i++)SnGeo[i]=1.0;
for(i=imin; i<imax; i++) SnGeo[i] = 1./SnGeo[i];
}
void OverlapReductionFunction(struct Data *data)
{
......
......@@ -179,8 +179,8 @@ void extrinsic_uniform_proposal(long *seed, double *y)
y[1] = -1.0 + 2.0*ran2(seed);
y[2] = ran2(seed)*LAL_PI_2;
y[3] = -0.99 + 1.98*ran2(seed);
//y[4] = ran2(seed)*LAL_TWOPI;
y[4] = ran2(seed)*LAL_PI;
y[4] = ran2(seed)*LAL_TWOPI;
//y[4] = ran2(seed)*LAL_PI;
y[5] = 1.0;
}
......@@ -304,7 +304,8 @@ void extrinsic_fisher_information_matrix(struct FisherMatrix *fisher, double *pa
range = LAL_TWOPI;
break;
case 5:
range = 10.0;
//TODO: HACK! testing extrinsic amplitude update
range = 2.0;
break;
}
......@@ -652,6 +653,7 @@ void intrinsic_fisher_proposal(UNUSED struct Model *model, double **range, doubl
free_dvector(dparamsy,0,5);
}
//TODO: TEST FULL INTRINSIC FISHER MATRIX
void intrinsic_fisher_proposal_2(struct Model *model, double **range, double *paramsx, double *paramsy, double *Snf, double Snx, double *logpx, double *logpy, long *seed, double Tobs, int *test)
{
int j,k;
......
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