Maintenance will be performed on git.ligo.org and chat.ligo.org on the 25th April 2019 between 9-10pm CDT, towards the end of the maintenance window there will be a short (5-10 minute) downtime.

Commit 82ca34a8 authored by Tyson Littenberg's avatar Tyson Littenberg

fix bug when trying to use BL on sim. data

parent ace7c56d
......@@ -9,6 +9,8 @@
#include <string.h>
#include <math.h>
#include <sys/stat.h>
#include <fftw3.h>
#include "BayesWave.h"
#include "BayesLine.h"
......@@ -178,7 +180,7 @@ int main(int argc, char *argv[])
//read in frequency domain data from LALInference data structure
for(i=0; i<N/2; i++)
{
if(i<imin)
if(i<imin)
{
freqData[ifo][2*i] = 0.0;
freqData[ifo][2*i+1] = 0.0;
......@@ -339,26 +341,37 @@ int main(int argc, char *argv[])
j = i+imin;
bayesline[0][ifo]->power[i] = (asd[2*j]*asd[2*j]+asd[2*j+1]*asd[2*j+1]);
}
/*
When simulating data, LALInferenceReadData does not fill the time series.
BayesLineBurnin requires the time-domain data.
So, instead of trying to figure out if we are simulating data or not,
we always iFFT the freqData and feed that to BayesLine.
*/
//copy frequency domain data into timeData array for inplace iFFT
for(i=0; i<N; i++) timeData[ifo][i] = data->s[ifo][i];
fftw_wrapper(timeData[ifo],data->N,-1);
fprintf(stdout,"BayesLine search phase for IFO %s\n", data->ifos[ifo]);
BayesLineBurnin(bayesline[0][ifo], timeData[ifo], asd);
BayesLineBurnin(bayesline[0][ifo], timeData[ifo], asd);
//Copy BayesLine model into BayesWave model
for(i=0; i<data->N/2; i++)
{
if(i>=imin && i<imax)
{
model[0]->Snf[ifo][i] = bayesline[0][ifo]->Snf[i-imin];
model[0]->SnS[ifo][i] = bayesline[0][ifo]->Sbase[i-imin];
}
else
{
model[0]->Snf[ifo][i] = 1.0;
model[0]->SnS[ifo][i] = 1.0;
}
model[0]->invSnf[ifo][i] = 1./model[0]->Snf[ifo][i];
}
for(i=0; i<data->N/2; i++)
{
if(i>=imin && i<imax)
{
model[0]->Snf[ifo][i] = bayesline[0][ifo]->Snf[i-imin];
model[0]->SnS[ifo][i] = bayesline[0][ifo]->Sbase[i-imin];
}
else
{
model[0]->Snf[ifo][i] = 1.0;
model[0]->SnS[ifo][i] = 1.0;
}
model[0]->invSnf[ifo][i] = 1./model[0]->Snf[ifo][i];
}
//Quit lying to BayesLine about the data. Poor BayesLine.
for(i=0; i<(imax-imin); i++)
{
......@@ -373,18 +386,18 @@ int main(int argc, char *argv[])
{
fprintf(stdout,"BayesLine clone PSD to parallel chains for IFO %s\n", data->ifos[ifo]);
for(ic=1; ic<chain->NC; ic++)
{
copy_bayesline_params(bayesline[0][ifo], bayesline[ic][ifo]);
{
copy_bayesline_params(bayesline[0][ifo], bayesline[ic][ifo]);
for(i=0; i<data->N/2; i++)
{
model[ic]->Snf[ifo][i] = model[0]->Snf[ifo][i];
model[ic]->SnS[ifo][i] = model[0]->SnS[ifo][i];
model[ic]->invSnf[ifo][i] = model[0]->invSnf[ifo][i];
}
for(i=0; i< bayesline[chain->index[0]][ifo]->data->ncut; i++)
{
bayesline[chain->index[ic]][ifo]->priors->lower[i] = bayesline[chain->index[0]][ifo]->priors->lower[i];
bayesline[chain->index[ic]][ifo]->priors->upper[i] = bayesline[chain->index[0]][ifo]->priors->upper[i];
}
for(i=0; i< bayesline[chain->index[0]][ifo]->data->ncut; i++)
{
bayesline[chain->index[ic]][ifo]->priors->lower[i] = bayesline[chain->index[0]][ifo]->priors->lower[i];
bayesline[chain->index[ic]][ifo]->priors->upper[i] = bayesline[chain->index[0]][ifo]->priors->upper[i];
}
}
}
......
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <fftw3.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
......@@ -211,3 +212,82 @@ double uniform_draw(gsl_rng *seed)
{
return gsl_rng_uniform(seed);
}
/* ********************************************************************************** */
/* */
/* Fourier Routines */
/* */
/* ********************************************************************************** */
void fftw_wrapper(double *data, int N, int flag)
{
int n;
double *timeData = (double *)malloc(N*sizeof(double));
fftw_complex *freqData = (fftw_complex *) fftw_malloc(sizeof(fftw_complex)*N);
//setup FFTW plans
fftw_plan reverse = fftw_plan_dft_c2r_1d(N, freqData, timeData, FFTW_MEASURE);
fftw_plan forward = fftw_plan_dft_r2c_1d(N, timeData, freqData, FFTW_MEASURE);
switch (flag)
{
//Reverse transform
case -1:
//fill freqData with contents of data
for(n=0; n<N/2; n++)
{
freqData[n][0] = data[2*n];
freqData[n][1] = data[2*n+1];
}
//get DC and Nyquist where FFTW wants them
freqData[N/2][0] = freqData[0][1];
freqData[N/2][1] = freqData[0][1] = 0.0;
//The FFT
fftw_execute(reverse);
//Copy output back into data
for(n=0; n<N; n++) data[n] = timeData[n];
break;
//Forward transform
case 1:
//fill timeData with contents of data
for(n=0; n<N; n++) timeData[n] = data[n];
//fill timeData with contents of data
for(n=0; n<N; n++) timeData[n] = data[n];
//The FFT
fftw_execute(forward);
//Copy output back into data
for(n=0; n<N/2; n++)
{
data[2*n] = freqData[n][0];
data[2*n+1] = freqData[n][1];
}
//get DC and Nyquist where FFTW wants them
data[1] = freqData[N/2][1];
break;
default:
fprintf(stdout,"Error: unsupported type in fftw_wrapper()\n");
exit(1);
break;
}
fftw_destroy_plan(reverse);
fftw_destroy_plan(forward);
free(timeData);
fftw_free(freqData);
}
......@@ -35,5 +35,10 @@ double gaussian_draw(gsl_rng *seed);
double uniform_draw(gsl_rng *seed);
void dfour1(double data[], unsigned long nn, int isign);
void drealft(double data[], unsigned long n, int isign);
/* ********************************************************************************** */
/* */
/* Fourier Routines */
/* */
/* ********************************************************************************** */
void fftw_wrapper(double *data, int N, int flag);
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