Commit 446dca3f authored by Tyson Littenberg's avatar Tyson Littenberg
Browse files

return to unpolarized case

parent 65e1990a
......@@ -165,7 +165,7 @@ void recompute_residual(struct Data *data, struct Model **model, struct Chain *c
for(j=1; j<=model_x->signal[0]->size; j++)
{
i = model_x->signal[0]->index[j];
model_x->wavelet(model_x->signal[0]->templates, model_x->signal[0]->intParams[i], N, 1, data->Tobs);
for(ip=0; ip<NP; ip++) model_x->wavelet(model_x->signal[ip]->templates, model_x->signal[ip]->intParams[i], N, 1, data->Tobs);
}
combinePolarizations(data, model_x->signal, model_x->h, model_x->extParams, model_x->Npol);
......
......@@ -74,6 +74,22 @@ static void catch_alarm(UNUSED int sig, UNUSED siginfo_t *siginfo,UNUSED void *c
}
void constrain_hplus_hcross(struct Wavelet **wave, int i)
{
/*
give hx special treatment (wavelets at same TFQ as h+, but independenta A,phi
*/
wave[1]->intParams[i][0] = wave[0]->intParams[i][0]; //t0
wave[1]->intParams[i][1] = wave[0]->intParams[i][1]; //f0
wave[1]->intParams[i][2] = wave[0]->intParams[i][2]; //Q
// printf("A+ = %g, Ax = %g\n",wave[0]->intParams[i][3],wave[1]->intParams[i][3]); //A
// printf("Phi+ = %g, Phix = %g\n",wave[0]->intParams[i][4],wave[1]->intParams[i][4]); //A
if(wave[0]->dimension > 5)
{
wave[1]->intParams[i][5] = wave[0]->intParams[i][5]; //beta (~fdot for chirplets)
}
}
static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior *prior, struct Model **model, char modelname[])
{
int i,j,n;
......@@ -113,26 +129,11 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
if(data->adaptTemperatureFlag)// && data->runPhase!=0)
{
chain->temperature[chain->NC-1] = 1.0e6;
//chain->dT[0] = pow(chain->temperature[chain->NC-1],(1./(double)(NC-1)));
//chain->dT[0] = pow(chain->temperature[chain->NC-1],(1./(double)(NC-1))*(1./(double)(NC-1))*(1./(double)(NC-1)));
chain->dT[0] = pow(chain->temperature[chain->NC-1],(1./(double)(NC-1))*(1./(double)(NC-1)));
//chain->dT[0] = pow(chain->temperature[chain->NC-1],(1./(double)(NC-1)));
for(ic=0; ic<chain->NC-1; ic++)
{
//if(chain->NC==15)
//{
//chain->temperature[ic]=chain->temperature[ic-1]*chain->dT[0];
//TODO: DO better initializing chain temperature ladder (hard-wired for NCmax=15!)
//chain->temperature[ic]=exp((double)(ic*ic)/(double)(((chain->NC-1)*(chain->NC-1))/(log(1.0e6)))); //empirically determined initial guess for temp spacing
//chain->temperature[ic] = 1./(1.0 - 0.07*(double)ic);
//}
//else
//{
//chain->temperature[ic] = pow(chain->dT[0],ic*ic*ic);
chain->temperature[ic] = pow(chain->dT[0],ic*ic);
//chain->temperature[ic] = pow(chain->dT[0],ic);
//}
chain->temperature[ic] = pow(chain->dT[0],ic);
}
}
......@@ -183,7 +184,7 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
{
extrinsic_uniform_proposal(chain->seed,model[ic]->extParams);
//FIXME: Fix ugly psi=0 hack for unpolarized case
if(!data->polarizationFlag) model[ic]->extParams[2] = 0.0;
//if(!data->polarizationFlag) model[ic]->extParams[2] = 0.0;
/*
WaveformProject() was split into two functions:
......@@ -210,6 +211,9 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
if(data->amplitudePriorFlag)
data->signal_amplitude_proposal(signal[n]->intParams[j], model[ic]->SnGeo, chain->seed, data->Tobs, prior->range, prior->sSNRpeak);
//void constrain_hplus_hcross(struct Wavelet **wave, int i)
if(n==1)constrain_hplus_hcross(signal, j);
model[ic]->wavelet(signal[n]->templates, signal[n]->intParams[j], data->N, 1, Tobs);
}
}
......@@ -1579,22 +1583,6 @@ static void add_wavelet(struct Wavelet *wave_x, struct Wavelet *wave_y, int *ii)
*ii = index;
}
void constrain_hplus_hcross(struct Wavelet **wave, int i)
{
/*
give hx special treatment (wavelets at same TFQ as h+, but independenta A,phi
*/
wave[1]->intParams[i][0] = wave[0]->intParams[i][0]; //t0
wave[1]->intParams[i][1] = wave[0]->intParams[i][1]; //f0
wave[1]->intParams[i][2] = wave[0]->intParams[i][2]; //Q
// printf("A+ = %g, Ax = %g\n",wave[0]->intParams[i][3],wave[1]->intParams[i][3]); //A
// printf("Phi+ = %g, Phix = %g\n",wave[0]->intParams[i][4],wave[1]->intParams[i][4]); //A
if(wave[0]->dimension > 5)
{
wave[1]->intParams[i][5] = wave[0]->intParams[i][5]; //beta (~fdot for chirplets)
}
}
/**************************************************************/
void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Model **model, struct Chain *chain, struct TimeFrequencyMap *tf, gsl_rng *seed, int ic)
......@@ -1838,12 +1826,21 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
larrayy = wave_y[n]->index;
q = 0.0;
if(det==-1)q += ClusterRate*TFtoProx*TFprop_density(paramsx[ii], range, tf, data->NI);
else q += ClusterRate*TFtoProx*TFprop_density(paramsx[ii], range, tf, det);
if(det==-1)
{
if(n==0)
{
q += ClusterRate*TFtoProx*TFprop_density(paramsx[ii], range, tf, data->NI);
}
}
else q += ClusterRate*TFtoProx*TFprop_density(paramsx[ii], range, tf, det);
q += ClusterRate*(1.-TFtoProx)*wavelet_proximity_density(paramsx[ii][1], paramsx[ii][0], paramsy, larrayy, wave_y[n]->size, prior);
q += (1.0-ClusterRate)/prior->TFV;
logqx += log(q);
if(n==0 || det > -1)
{
q += ClusterRate*(1.-TFtoProx)*wavelet_proximity_density(paramsx[ii][1], paramsx[ii][0], paramsy, larrayy, wave_y[n]->size, prior);
q += (1.0-ClusterRate)/prior->TFV;
logqx += log(q);
}
logqy += 0.0;
}
}
......@@ -1853,7 +1850,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
//loop over Npol for signal model, single IFO (nmax=nmin+1) for glitch model
for(n=nmin; n<nmax; n++)
{
logqx += -prior->logTFV;
if(n==0 || det >-1) logqx += -prior->logTFV;
logqy += 0.0;
}
}
......@@ -1935,14 +1932,31 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
larrayx = wave_x[n]->index;
larrayy = wave_y[n]->index;
if(det==-1)q += ClusterRate*TFtoProx*TFprop_density(paramsy[ii], range, tf, data->NI);
else q += ClusterRate*TFtoProx*TFprop_density(paramsy[ii], range, tf, det);
q += ClusterRate*(1.-TFtoProx)*wavelet_proximity_density(paramsy[ii][1], paramsy[ii][0], paramsx, larrayx, wave_x[n]->size, prior);
q += (1.0-ClusterRate)/prior->TFV;
logqy += log(q);
logqx += 0.0;
//signal model
if(det==-1)
{
//only evaluate TFQ proposal for h+ polarization (TFQ are not free parameters for hx)
if(n==0)
{
q += ClusterRate*TFtoProx*TFprop_density(paramsy[ii], range, tf, data->NI);
q += ClusterRate*(1.-TFtoProx)*wavelet_proximity_density(paramsy[ii][1], paramsy[ii][0], paramsx, larrayx, wave_x[n]->size, prior);
q += (1.0-ClusterRate)/prior->TFV;
logqy += log(q);
logqx += 0.0;
}
}
//glitch model
else
{
q += ClusterRate*TFtoProx*TFprop_density(paramsy[ii], range, tf, det);
q += ClusterRate*(1.-TFtoProx)*wavelet_proximity_density(paramsy[ii][1], paramsy[ii][0], paramsx, larrayx, wave_x[n]->size, prior);
q += (1.0-ClusterRate)/prior->TFV;
logqy += log(q);
logqx += 0.0;
}
}
}//end TF density proposal
......@@ -1959,7 +1973,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
//loop over Npol for signal model, single IFO (nmax=nmin+1) for glitch model
for(n=nmin; n<nmax; n++)
{
logqy += -prior->logTFV;
if(n==0 || det>-1) logqy += -prior->logTFV;
logqx += 0.0;
}
}//end uniform TF density proposal
......@@ -2136,6 +2150,11 @@ 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;
}
else
......@@ -2230,8 +2249,11 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
// uniform prior on time-frequency volume
else
{
wave_x[n]->logp += -(double)(wave_x[n]->size)*prior->logTFV;
wave_y[n]->logp += -(double)(wave_y[n]->size)*prior->logTFV;
if(n==0 || det>-1)
{
wave_x[n]->logp += -(double)(wave_x[n]->size)*prior->logTFV;
wave_y[n]->logp += -(double)(wave_y[n]->size)*prior->logTFV;
}
}
}//end loop over n
......@@ -2276,8 +2298,11 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
{
if(data->waveletPriorFlag)
{
wave_x[n]->logp += prior->Nwavelet[wave_x[n]->size];
wave_y[n]->logp += prior->Nwavelet[wave_y[n]->size];
if(n==0 || det>-1)
{
wave_x[n]->logp += prior->Nwavelet[wave_x[n]->size];
wave_y[n]->logp += prior->Nwavelet[wave_y[n]->size];
}
}
}
......@@ -2290,6 +2315,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
logH += logqx - logqy; //proposal density
for(n=nmin; n<nmax; n++) logH += wave_y[n]->logp - wave_x[n]->logp; //prior ratio
/*
if(model_x->size != model_y->size)
{
......@@ -2301,7 +2327,7 @@ void EvolveIntrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
}
*/
//if(chain->index[ic]==0 && data->stochasticFlag)printf("logLy=%g logLx=%g Ay=%g Ax=%g detC=%g logH=%g\n",model_y->logL-data->logLc,model_x->logL-data->logLc,model_x->background->logamp,model_y->background->logamp, model_y->background->detCij[data->N/4],logH);
//if(ic==0 && data->stochasticFlag)printf("logLy=%g logLx=%g Ay=%g Ax=%g detC=%g logH=%g\n",model_y->logL-data->logLc,model_x->logL-data->logLc,model_x->background->logamp,model_y->background->logamp, model_y->background->detCij[data->N/4],logH);
}// end prior test condition
else logH = -1.0e60; // Rejection sample at prior boundaries
......@@ -2403,7 +2429,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
//Initialize parameter vectors & likelihood for extrinsic parameter MCMC
for(i=0; i<NE; i++) paramsy[i] = paramsx[i];
//FIXME: Fix ugly psi=0 hack for unpolarized case
if(!data->polarizationFlag) paramsy[2] = paramsx[2] = 0.0;
//if(!data->polarizationFlag) paramsy[2] = paramsx[2] = 0.0;
//TODO: Why can't I use a band-limited response for xtrinsic moves?
//Related to computeProjectionCoeffs()?
......@@ -2497,7 +2523,7 @@ void EvolveExtrinsicParameters(struct Data *data, struct Prior *prior, struct Mo
}
if(!data->extrinsicAmplitudeFlag) paramsy[5]=1.0;
//FIXME: Fix ugly psi=0 hack for unpolarized case
if(!data->polarizationFlag) paramsy[2] = paramsx[2] = 0.0;
//if(!data->polarizationFlag) paramsy[2] = paramsx[2] = 0.0;
//Rescale intrinsic amplitudes and phases with extrinsic parameters
for(i=1; i<ienddim; i++)
......
......@@ -131,13 +131,10 @@ int extrinsic_checkrange(double *p)
//psi ranges from [0,pi/2]
while(p[2]<0.0 || p[2] >= LAL_PI_2)
while(p[2]<0.0 || p[2] >= LAL_PI)
{
if(p[2]<0.0) p[2] += LAL_PI_2;
else p[2] -= LAL_PI_2;
//psi->psi+pi/2 needs phi->phi+pi
p[4] += LAL_PI;
if(p[2]<0.0) p[2] += LAL_PI;
else p[2] -= LAL_PI;
}
//TODO: ellipticity should be periodic?
......
......@@ -216,7 +216,7 @@ void uniform_orientation_proposal(double *y, gsl_rng *seed)
double newPsi,newPhi,newEcc;
newPhi = uniform_draw(seed)*LAL_TWOPI;
newPsi = uniform_draw(seed)*LAL_PI_2;
newPsi = uniform_draw(seed)*LAL_PI;
newEcc = -1. + 2.*uniform_draw(seed);
y[2] = newPsi;
......
......@@ -22,8 +22,10 @@
#include <math.h>
#include "BayesWave.h"
#include "BayesLine.h"
#include "BayesWaveMath.h"
#include "BayesWaveWavelet.h"
#include "BayesWaveIO.h"
/* ********************************************************************************** */
......
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