diff --git a/src/BayesWaveLikelihood.c b/src/BayesWaveLikelihood.c
index 773a47ee82fd9f972f89174859841e783a8124b8..4af2c3bddda93585c97a56125d0aaff28f6190e2 100644
--- a/src/BayesWaveLikelihood.c
+++ b/src/BayesWaveLikelihood.c
@@ -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);
diff --git a/src/BayesWaveMCMC.c b/src/BayesWaveMCMC.c
index dd20b8589650b234417f51637c3ab4bb742f20ee..defaf95a23c930bf3903ef4c4fbe46e81ab911f5 100644
--- a/src/BayesWaveMCMC.c
+++ b/src/BayesWaveMCMC.c
@@ -74,6 +74,20 @@ 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
+ 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 +127,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);
}
}
@@ -182,8 +181,6 @@ static void restart_sampler(struct Data *data, struct Chain *chain, struct Prior
if(data->signalFlag)
{
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;
/*
WaveformProject() was split into two functions:
@@ -210,6 +207,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);
}
}
@@ -511,10 +511,10 @@ static void resume_sampler(struct Data *data, struct Chain *chain, struct Model
fprintf(stderr," Exiting to system (1)\n");
exit(1);
}
- parse_signal_parameters(data, model[ic], fptr2, h);
- for(n=0; n<data->Npol; n++) fclose(fptr2[n]);
- free(fptr2);
}
+ parse_signal_parameters(data, model[ic], fptr2, h);
+ for(n=0; n<data->Npol; n++) fclose(fptr2[n]);
+ free(fptr2);
}
// PSD model
@@ -1579,22 +1579,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 +1822,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 +1846,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 +1928,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 +1969,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 +2146,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 +2245,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 +2294,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 +2311,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 +2323,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
@@ -2402,8 +2424,6 @@ 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;
//TODO: Why can't I use a band-limited response for xtrinsic moves?
//Related to computeProjectionCoeffs()?
@@ -2496,8 +2516,6 @@ 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;
//Rescale intrinsic amplitudes and phases with extrinsic parameters
for(i=1; i<ienddim; i++)
diff --git a/src/BayesWavePrior.c b/src/BayesWavePrior.c
index 0ebf2a5c9b4dad5852b18ad8a8a9c0ccb7bc75b5..8421dd5b7d65d0082e7dee4958b2aa9e8853c899 100644
--- a/src/BayesWavePrior.c
+++ b/src/BayesWavePrior.c
@@ -130,14 +130,11 @@ int extrinsic_checkrange(double *p)
if(p[1] < -1.0 || p[1] >= 1.0) return 1;
- //psi ranges from [0,pi/2]
- while(p[2]<0.0 || p[2] >= LAL_PI_2)
+ //psi ranges from [0,pi]
+ 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?
diff --git a/src/BayesWaveProposal.c b/src/BayesWaveProposal.c
index 3422911e2cd07912688a4bdea04d5b7fedb0456d..a2ced926414ea9d19b6ee0dec6f86597bd830b77 100644
--- a/src/BayesWaveProposal.c
+++ b/src/BayesWaveProposal.c
@@ -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;
diff --git a/src/BayesWaveToLALPSD.c b/src/BayesWaveToLALPSD.c
index 3c2a79fef1e80cd5a96de29b42ba4c727fa2eacb..7cc955fef770e493ee63e102f2b2308ae7b413ce 100644
--- a/src/BayesWaveToLALPSD.c
+++ b/src/BayesWaveToLALPSD.c
@@ -34,8 +34,17 @@
int main(int argc, char* argv[])
{
+ if(argc<2)
+ {
+ fprintf(stdout,"USAGE: BayesWaveToLALPSD IFO [runName]\n");
+ return 0;
+ }
+
FILE *waveFile=NULL;
char filename[1024];
+ char ifoname[4];
+
+ sprintf(ifoname,"%s",argv[1]);
int i,n;
@@ -62,8 +71,8 @@ int main(int argc, char* argv[])
{
//sprintf(filename,"BayesWave_PSD_%s_clean_psd_%i.dat.0",argv[1],i);
- if(argc==2)sprintf(filename,"waveforms/%s_clean_psd_%i.dat.0",argv[1],i);
- else sprintf(filename,"waveforms/clean_psd_%i.dat.0",i);
+ if(argc==3)sprintf(filename,"waveforms/%s_clean_psd_%s_%i.dat",argv[2],ifoname,i);
+ else sprintf(filename,"waveforms/clean_psd_%s_%i.dat",ifoname,i);
waveFile=fopen(filename,"r");
if(waveFile==NULL)
@@ -96,14 +105,18 @@ int main(int argc, char* argv[])
fclose(waveFile);
}
+
+ sprintf(filename,"%s_psd_intervals.dat",ifoname);
+ FILE *outfile = fopen(filename,"w");
- if(argc==2)sprintf(filename,"%s_IFO0_asd_median.dat",argv[1]);
- else sprintf(filename,"IFO0_asd_median.dat");
+ sprintf(filename,"%s_psd_median.dat",ifoname);
+ FILE *outfile2 = fopen(filename,"w");
- FILE *outfile = fopen("psd_intervals.dat","w");
- FILE *outfile2 = fopen("psd_median.dat","w");
+ if(argc==3)sprintf(filename,"%s_%s_asd_median.dat",argv[2],ifoname);
+ else sprintf(filename,"%s_asd_median.dat",ifoname);
FILE *outfile3 = fopen(filename,"w");
+
if(outfile3==NULL)printf("Failed to open %s\n",filename);
double median,upper90,lower90,upper50,lower50;
diff --git a/src/BayesWaveWavelet.c b/src/BayesWaveWavelet.c
index 61fdd778c8b1bc327ea6bcf230f481e5736ba21a..ba663125d2da0b1d332a246be672c7cdc28d096e 100644
--- a/src/BayesWaveWavelet.c
+++ b/src/BayesWaveWavelet.c
@@ -22,8 +22,10 @@
#include <math.h>
#include "BayesWave.h"
+#include "BayesLine.h"
#include "BayesWaveMath.h"
#include "BayesWaveWavelet.h"
+#include "BayesWaveIO.h"
/* ********************************************************************************** */