Maintenance will be performed on git.ligo.org, chat.ligo.org, and docs.ligo.org, starting at approximately 10am CDT Tuesday 20 August 2019. The maintenance is expected to take around an hour and here will be two short periods of downtime, one at the beginning of the maintenance and another at the end.

Commit a92f09e4 authored by James Clark's avatar James Clark

Merge branch 'bayesline_errors' into 'master'

fix bug when trying to use BL on sim. data

See merge request lscsoft/bayeswave!66
parents ace7c56d 82ca34a8
Pipeline #47686 passed with stages
in 1 minute and 37 seconds
......@@ -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