Skip to content
Snippets Groups Projects

Fix "PSD" options in help message

Merged Fix "PSD" options in help message
katerina.chatziioannou/bayeswave:master into master
Merged Katerina Chatziioannou requested to merge katerina.chatziioannou/bayeswave:master into master
Compare
and
1 file
+ 1209
− 1209
Compare changes
  • Side-by-side
  • Inline
src/BayesWave.c
+ 1209
− 1209
/*
* Copyright (C) 2018 Neil J. Cornish, Tyson B. Littenberg, Margaret Millhouse,Katerina Chatziioannou
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with with program; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
* MA 02111-1307 USA
*/
/*************************** REQUIRED LIBRARIES ***************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <sys/stat.h>
#include <fftw3.h>
#include "BayesWave.h"
#include "BayesLine.h"
#include "BayesWaveIO.h"
#include "BayesWaveMCMC.h"
#include "BayesWaveMath.h"
#include "BayesWavePrior.h"
#include "BayesWaveModel.h"
#include "BayesWaveWavelet.h"
#include "BayesWaveProposal.h"
#include "BayesWaveEvidence.h"
#include "BayesWaveLikelihood.h"
/************* PROTOTYPE DECLARATIONS FOR INTERNAL FUNCTIONS **************/
void print_help_message(void);
/* ============================ MAIN PROGRAM ============================ */
int main(int argc, char *argv[])
{
/* Variable declaration */
int i, j, ic, ifo, imin,imax;
char filename[MAXSTRINGSIZE];
char modelname[MAXSTRINGSIZE];
/* LAL data structure that will contain all of the data from the frame files */
LALInferenceRunState *runState = XLALCalloc(1,sizeof(LALInferenceRunState));
runState->commandLine = LALInferenceParseCommandLine(argc,argv);
/* Print version number and command line to screen */
print_version(stdout);
fprintf(stdout," %s\n",LALInferencePrintCommandLine(runState->commandLine) );
fprintf(stdout,"\n");
/******************************************************************************/
/* */
/* Usage statment if --help */
/* */
/******************************************************************************/
if(LALInferenceGetProcParamVal(runState->commandLine, "--version"))
{
return 0;
}
if(LALInferenceGetProcParamVal(runState->commandLine, "--help") || argc==1)
{
print_help_message();
return 0;
}
/******************************************************************************/
/* */
/* Read/store data */
/* */
/******************************************************************************/
fprintf(stdout, "\n");
fprintf(stdout, " ======= LALInferenceReadData() ======\n");
runState->data = LALInferenceReadData(runState->commandLine);
fprintf(stdout, "\n");
/******************************************************************************/
/* */
/* Unpack runState into BayesWave data structures */
/* */
/******************************************************************************/
LALInferenceIFOData *dataPtr = NULL;
dataPtr = runState->data;
double *SNRinj = NULL;
int NI;
int N = dataPtr->timeData->data->length;
double Tobs = (double)N*dataPtr->timeData->deltaT;
ifo=0;
while(dataPtr!=NULL)
{
dataPtr = dataPtr->next;
ifo++;
}
NI=ifo;
double **freqData = double_matrix(NI-1,N-1);
double **psd = double_matrix(NI-1,N/2-1);
double **timeData = double_matrix(NI-1,N-1);
struct Data *data = malloc(sizeof(struct Data));
struct Chain *chain = malloc(sizeof(struct Chain));
struct Prior *prior = malloc(sizeof(struct Prior));
const gsl_rng_type *T = gsl_rng_default;
chain->seed = gsl_rng_alloc(T);
gsl_rng_env_setup();
parse_command_line(data, chain, prior, runState->commandLine);
if (data->chirpletFlag) {
data->NW = 6;
}
else{
data->NW = 5;
}
// Move to the directory specified by the user, or stay in ./
chdir(data->outputDirectory);
fprintf(stdout,"\nOutputting all data to: %s \n\n", data->outputDirectory);
// Create working directories for output
mode_t process_mask = umask(0);
mkdir("chains", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
mkdir("waveforms", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
if(data->checkpointFlag) mkdir("checkpoint",S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
sprintf(filename,"%ssnr",data->runName);
mkdir(filename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
umask(process_mask);
// burstMDC injection
if(LALInferenceGetProcParamVal(runState->commandLine, "--MDC-cache"))
{
SNRinj = malloc((NI+1)*sizeof(double));
fprintf(stdout, " ==== LALInferenceInjectFromMDC(): started. ====\n");
InjectFromMDC(runState->commandLine, runState->data, SNRinj);
fprintf(stdout, " ==== LALInferenceInjectFromMDC(): finished. ====\n\n");
free(SNRinj);
}
// CBC injection
if(LALInferenceGetProcParamVal(runState->commandLine, "--inj"))
{
fprintf(stdout, " ==== LALInferenceInjectInspiralSignal(): started. ====\n");
LALInferenceInjectInspiralSignal(runState->data,runState->commandLine);
fprintf(stdout, " ==== LALInferenceInjectInspiralSignal(): finished. ====\n\n");
}
//Output Fourier domain data for BayesLine
dataPtr = runState->data;
double fmin = dataPtr->fLow;
double fmax = ceil(dataPtr->fHigh);
//set frequency resolution to be 32 Hz because why not?
double df,dt;
df = 32.0;
dt = 0.5/df;
int tsize = (int)(Tobs/dt);
//Copy Fourier domain data & PSD from LAL data types into simple arrays
ifo=0;
imin = (int)(fmin*Tobs);
imax = (int)(fmax*Tobs);
bool simDataFlag = False;
while(dataPtr!=NULL)
{
//read in time domain data from LALInference data structre
for(i=0; i<N; i++)
{
timeData[ifo][i] = dataPtr->timeData->data->data[i];
//timeData remains uninitialized if simulating data
if(timeData[ifo][i]!=timeData[ifo][i])
{
simDataFlag = True;
timeData[ifo][i] = 0.0;
}
}
//read in frequency domain data from LALInference data structure
for(i=0; i<N/2; i++)
{
// if(simDataFlag)
// {
// dataPtr->freqData->data->data[imin] = dataPtr->freqData->data->data[imin+1];
// dataPtr->oneSidedNoisePowerSpectrum->data->data[imin] = dataPtr->oneSidedNoisePowerSpectrum->data->data[imin+1];
// }
if(i<imin+1)
{
freqData[ifo][2*i] = creal(dataPtr->freqData->data->data[imin+1]);
freqData[ifo][2*i+1] = cimag(dataPtr->freqData->data->data[imin+1]);
psd[ifo][i] = dataPtr->oneSidedNoisePowerSpectrum->data->data[imin+1]*Tobs/2.0;
}
else
{
freqData[ifo][2*i] = creal(dataPtr->freqData->data->data[i]);
freqData[ifo][2*i+1] = cimag(dataPtr->freqData->data->data[i]);
/*
Innerproducts are expecting <n_i^2> instead of Sn(f):
So psd arrays == <n_i^2> = T/2 * Sn(f)
*/
psd[ifo][i] = dataPtr->oneSidedNoisePowerSpectrum->data->data[i]*Tobs/2.0;
}
}
dataPtr = dataPtr->next;
ifo++;
}
if(simDataFlag)
{
fprintf(stdout,"Detected run is using simulated data\n");
fprintf(stdout,"copied d(fmin+1/t) to d(f<=fmin)\n\n");
}
/******************************************************************************/
/* */
/* Setup DATA, CHAIN, PRIOR, MODEL structures */
/* */
/******************************************************************************/
/*
Setup DATA structure
*/
initialize_data(data,freqData,N,tsize,Tobs,NI,fmin,fmax);
//print run summary file
sprintf(filename,"%sbayeswave.run",data->runName);
chain->runFile = fopen(filename,"w");
print_version(chain->runFile);
fprintf(chain->runFile," %s\n\n",LALInferencePrintCommandLine(data->commandLine) );
print_run_stats(chain->runFile, data, chain);
fclose(chain->runFile);
data->detector = malloc(NI*sizeof(LALDetector*));
data->epoch = runState->data->epoch;
ifo=0;
dataPtr = runState->data;
while(dataPtr!=NULL)
{
sprintf(data->ifos[ifo],"%s",dataPtr->name);
data->detector[ifo] = dataPtr->detector;
ifo++;
dataPtr = dataPtr->next;
}
OverlapReductionFunction(data);
/*
Setup CHAIN structure
*/
//give chain buffer in case we need more
int NC=chain->NC;
chain->NC+=20;
allocate_chain(chain,NI,data->Npol);
for(ic=0; ic<chain->NC; ic++) chain->index[ic]=ic;
/*
Setup PRIOR structure
*/
// int Dmax = data->Dmax;
// data->Dmax = 200;
int omax = 10;
if(data->Dmax < omax) omax = data->Dmax;
initialize_priors(data, prior, omax);
data->glitch_amplitude_proposal = draw_glitch_amplitude;
data->glitch_amplitude_prior = glitch_amplitude_prior;
// choose functions for amplitude priors
if(data->signalAmplitudePriorFlag)
{
data->signal_amplitude_proposal = draw_signal_amplitude;
data->signal_amplitude_prior = signal_amplitude_prior;
}
else
{
data->signal_amplitude_proposal = draw_glitch_amplitude;
data->signal_amplitude_prior = glitch_amplitude_prior;
}
/*
Setup MODEL structure
*/
struct Model **model = malloc(chain->NC*sizeof(struct Model*));
for(ic=0; ic<chain->NC; ic++)
{
model[ic] = malloc(sizeof(struct Model));
initialize_model(model[ic], data, prior, psd, chain->seed);
}
// BayesWave internal injection
if(LALInferenceGetProcParamVal(runState->commandLine,"--BW-inject"))
{
fprintf(stdout, " ==== BayesWaveInjection(): started. ====\n");
BayesWaveInjection(runState->commandLine, data, chain, prior, psd, &NC);
fprintf(stdout, " ==== BayesWaveInjection(): finished. ====\n\n");
}
// Print Fourier domain data to file for post-processing
for(ifo=0; ifo<data->NI; ifo++)
{
sprintf(filename,"waveforms/%sfourier_domain_data_%s.dat",data->runName,data->ifos[ifo]);
print_frequency_domain_data(filename, data->s[ifo], data->N, data->Tobs, data->imin, data->imax);
}
/******************************************************************************/
/* */
/* PSD search has to be done once even on restart */
/* Much of the BayesLine setup is done inside of BayesLineSearch() */
/* */
/******************************************************************************/
struct BayesLineParams ***bayesline = malloc(chain->NC*sizeof(struct BayesLineParams **));
if(data->bayesLineFlag)
{
fprintf(stdout,"\n ============ BayesLine ==============\n");
/*
Setup BayesLine structure
*/
for(ic=0; ic<chain->NC; ic++)
{
bayesline[ic] = malloc(data->NI*sizeof(struct BayesLineParams *));
initialize_bayesline(bayesline[ic], data, psd);
}
/*
BayesLine spectral estimation search-phase
*/
double *asd = malloc(N*sizeof(double));
for(ifo=0; ifo<NI; ifo++)
{
//Trick BayesLine into fitting to off source PSD.
for(i=0; i<data->N/2; i++)
{
asd[2*i] = data->s[ifo][2*i];//sqrt(psd[ifo][i]/2.);
asd[2*i+1] = data->s[ifo][2*i+1];//sqrt(psd[ifo][i]/2.);
}
for(i=0; i<(imax-imin); i++)
{
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];
//pad frequencies below fmin, just in case
for(i=0; i<=imin; i++)
{
timeData[ifo][2*i] = data->s[ifo][2*(imin+1)];
timeData[ifo][2*i+1] = data->s[ifo][2*(imin+1)+1];
}
fftw_wrapper(timeData[ifo],data->N,-1);
//normalize time series as expected by BayesLine
for(i=0; i<N; i++) timeData[ifo][i]/=data->Tobs;
fprintf(stdout,"BayesLine search phase for IFO %s\n", data->ifos[ifo]);
BayesLineBurnin(bayesline[0][ifo], timeData[ifo], asd, data->ifos[ifo]);
//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];
}
//Quit lying to BayesLine about the data. Poor BayesLine.
for(i=0; i<(imax-imin); i++)
{
j = i+imin;
bayesline[0][ifo]->power[i] = (data->s[ifo][2*j]*data->s[ifo][2*j]+data->s[ifo][2*j+1]*data->s[ifo][2*j+1]);
}
}
free(asd);
for(ifo=0; ifo<NI; ifo++)
{
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]);
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];
}
}
}
fprintf(stdout,"\n");
}
else
{
/*
Here is an ugly hack:
BayesWavePost expects a psd for its IFO-cache argument.
The bayeswave_pipe script defaults to using the fairdraw output from BayesLine
e.g. --[IFO]-cache interp:[rundir]/[IFO]_fairdraw_asd.dat
which doesn't exist when running w/out BayesLine, e.g. simulated data for review
purposes. Here we print a file with the right name and content so that bayeswave_pipe
doesn't have to decide what to use for the cache argument.
TODO: rethink naming to keep people out of trouble
*/
export_cleaned_data(data, model[0]);
}
//Restore chain size to requested value
chain->NC=NC;
/******************************************************************************/
/* */
/* PSD search & glitch cleaning stages */
/* Remove excess power from data outside of analysis window */
/* Residuals and PSDs are passed on to model selection phase */
/* Skip this step on restart */
/* */
/******************************************************************************/
// Restart run using stored residual and PSD
if(data->checkpointFlag)
{
fprintf(stdout, " ========= Checkpoint flag found ==========\n");
fflush(stdout);
if(data->stochasticFlag)
{
fprintf(stdout,"Stochastic background model not compatible with checkpointing\n");
fprintf(stdout,"Disabling checkpointFlag\n");
fprintf(stdout,"May the Condor gods be with you...\n");
data->checkpointFlag = 0;
}
else
{
FILE *fptr=NULL;
if( (fptr = fopen("checkpoint/state.dat","r")) == NULL )
{
fprintf(stdout,"No checkpoint files found\n");
fprintf(stdout,"Continue with full run\n");
data->resumeFlag = 0; //set resume flag to 0 so run starts cleanly
}
else
{
/****************************************************/
/* */
/* READ IN CURRENT STATE VECTOR */
/* */
/****************************************************/
fscanf(fptr,"%s", modelname);
fscanf(fptr,"%i",&data->cleanModelFlag);
fscanf(fptr,"%i",&data->noiseModelFlag);
fscanf(fptr,"%i",&data->glitchModelFlag);
fscanf(fptr,"%i",&data->signalModelFlag);
fscanf(fptr,"%i",&data->fullModelFlag);
fscanf(fptr,"%i",&data->cleanOnlyFlag);
fscanf(fptr,"%i",&data->loudGlitchFlag);
fclose(fptr);
/****************************************************/
/* */
/* PTMCMC TEMPERATURE LADDER */
/* */
/****************************************************/
//Read in dT, index, temperature, ptprop, ptacc, A, avg&var loglikelihood?
sprintf(filename,"checkpoint/temperature.dat");
if( (fptr = fopen(filename,"r")) == NULL )
{
fprintf(stderr,"Error: Could not checkpoint temperature model\n");
fprintf(stderr," Parameter file %s does not exist\n",filename);
fprintf(stderr," Exiting to system (1)\n");
exit(1);
}
fscanf(fptr,"%i",&chain->mc);
fscanf(fptr,"%i",&NC);
fscanf(fptr,"%i",&chain->zcount);
fclose(fptr);
sprintf(filename,"checkpoint/rng.bin");
if( (fptr = fopen(filename,"rb")) == NULL )
{
fprintf(stderr,"Error: Could not checkpoint the rng\n");
fprintf(stderr," Parameter file %s does not exist\n",filename);
fprintf(stderr," Exiting to system (1)\n");
exit(1);
}
gsl_rng_fread(fptr,chain->seed);
fclose(fptr);
/****************************************************/
/* */
/* READ IN CURRENT EVIDENCE RESULTS */
/* */
/****************************************************/
if( (fptr = fopen("checkpoint/evidence.dat","r")) == NULL )
{
fprintf(stderr,"Error: Could not checkpoint model settings\n");
fprintf(stderr," Parameter file checkpoint/evidence.dat does not exist\n");
fprintf(stderr," Exiting to system (1)\n");
exit(1);
}
fscanf(fptr,"%lg %lg",&data->logZnoise,&data->varZnoise);
fscanf(fptr,"%lg %lg",&data->logZglitch,&data->varZglitch);
fscanf(fptr,"%lg %lg",&data->logZsignal,&data->varZsignal);
fscanf(fptr,"%lg %lg",&data->logZfull,&data->varZfull);
fclose(fptr);
fprintf(stdout,"Checkpoint files found for model %s\n",modelname);
fprintf(stdout,"Resume previous run\n");
data->resumeFlag = 1; //set resume flag to 1 so run picks up at checkpoint
}
}
fprintf(stdout,"\n");
}
/******************************************************************************/
/* */
/* The Full-time MCMC Run (Cleaning phase) */
/* */
/******************************************************************************/
if(data->stochasticFlag)
{
fprintf(stdout,"Stochastic background model not compatible with glitch model\n");
fprintf(stdout,"Disabling cleanModelFlag\n");
fprintf(stdout,"May the glitch gods be with you...\n");
data->cleanModelFlag = 0;
data->noiseModelFlag = 1;
data->glitchModelFlag = 0;
data->signalModelFlag = 0;
data->fullModelFlag = 0;
}
if(data->cleanModelFlag)
{
//lie to the code about the number of chains & iterations
// data->Dmax = 200;
int M = chain->count;
int NC = chain->NC;
int CP = data->clusterPriorFlag;
data->clusterPriorFlag = 0;
chain->count = M;
chain->NC = NC/5;
if(chain->NC < 1)chain->NC=1;
//Initialize proposal ratios for glitch cleaning phase
chain->modelRate = 0.0; //even chance for each channel
if(data->constantLogLFlag)
chain->rjmcmcRate = 0.01;
else
chain->rjmcmcRate = 0.5;
//Shut off RJMCMC if --fixD argument given
if(data->rjFlag==0) chain->rjmcmcRate = 0.0;
data->cleanFlag = 1;
data->glitchFlag = 1;
data->signalFlag = 0;
data->runPhase = 0;
initialize_chain(chain,1);
reset_likelihood(data);
reset_priors(data, prior);
reset_model(data, chain, prior, model);
double cleanLogZ, cleanVarZ;
printf("characterizing PSD+glitch model...\n\n");
//TODO disable checkpointing and adaptation during cleaning phase
int saveAdaptationFlag = data->adaptTemperatureFlag;
data->adaptTemperatureFlag = 0;
RJMCMC(data, model, bayesline, chain, prior, &cleanLogZ, &cleanVarZ);
data->adaptTemperatureFlag = saveAdaptationFlag;
/******************************************************************************/
/* */
/* Remove excess power outside the trigger time, reset dimension priors */
/* */
/******************************************************************************/
data->runPhase = 1;
//Restore the number of chains & iterations
chain->count = M;
chain->NC = NC;
// data->Dmax = Dmax;
data->clusterPriorFlag = CP;
reset_priors(data, prior);
//Permanently regress any glitches outside of the restricted prior range
fprintf(stdout, "\n");
fprintf(stdout, " Removing excess power outside of t=[%g,%g]\n",prior->range[0][0],prior->range[0][1]);
//SNR of signal remaining in data
double rSNR=0.0;
double ifoSNR[NI];
double snrRatio;
for(ifo=0; ifo<NI; ifo++) ifoSNR[ifo]=0.0;
struct Model *m = NULL;
struct Wavelet *g = NULL;
double **r = malloc(NI*sizeof(double*)); //total residual (just for plotting)
for(ifo=0; ifo<NI; ifo++)
{
r[ifo] = malloc(N*sizeof(double));
for(i=0; i<N; i++) r[ifo][i] = data->s[ifo][i];
m = model[chain->index[0]];
g = m->glitch[ifo];
for(i=1; i<=g->size; i++)
{
if(g->intParams[i][0]<prior->range[0][0] || g->intParams[i][0]>prior->range[0][1])
{
//regress from the data & glitch model se we can compute window SNR
m->wavelet(data->s[ifo], g->intParams[i], N, -1, data->Tobs);
m->wavelet(g->templates, g->intParams[i], N, -1, data->Tobs);
}
//regress all wavelets from residual array for plotting
m->wavelet(r[ifo], g->intParams[i], N, -1, data->Tobs);
}
ifoSNR[ifo] += fourier_nwip(data->imin,data->imax,g->templates,g->templates,m->invSnf[ifo]);
rSNR += ifoSNR[ifo];
fprintf(stdout, " Residual SNR in %s = %g\n",data->ifos[ifo],sqrt(ifoSNR[ifo]));
}
fprintf(stdout, " Residual SNR = %g\n",sqrt(rSNR));
/*
BayesWave can not handle very loud, very assymetric SNRs.
loudGlitchFlag is set to 1 if the SNR > 100 and the SNRratio is > 5:1
This hack is only implemented for 2-detector networks.
*/
if(NI==2)
{
snrRatio = ifoSNR[0]/ifoSNR[1];
if(snrRatio < 1.0) snrRatio = 1./snrRatio;
if(sqrt(rSNR)>100. && snrRatio > 5. )
{
fprintf(stdout,"\n");
fprintf(stdout,"************************ WARNING ************************\n");
fprintf(stdout,"Cleaning phase found loud, assymetric residual power \n");
fprintf(stdout,"Model selection cannot handle this trigger \n");
fprintf(stdout,"Assume it is a glitch but proceed with analysis \n");
fprintf(stdout,"*********************************************************\n");
data->loudGlitchFlag = 1;
}
}
// Print residual to file
for(ifo=0; ifo<data->NI; ifo++)
{
sprintf(filename,"waveforms/%sfourier_domain_residual_%s.dat",data->runName,data->ifos[ifo]);
print_frequency_domain_data(filename, r[ifo], data->N, data->Tobs, data->imin, data->imax);
sprintf(filename,"waveforms/%sfourier_domain_clean_%s.dat",data->runName,data->ifos[ifo]);
print_frequency_domain_data(filename, data->s[ifo], data->N, data->Tobs, data->imin, data->imax);
free(r[ifo]);
}
free(r);
/******************************************************************************/
/* */
/* Increase the number of parallel chains if residual SNR is too high */
/* */
/******************************************************************************/
//if(sqrt(rSNR)>0.0) resize_model(data, chain, prior, model, bayesline, psd, chain->NC+10);
if(sqrt(rSNR)>100.0) chain->NC+=10;
if(sqrt(rSNR)>200.0) chain->NC+=10;
/******************************************************************************/
/* */
/* Copy T=0 chain's PSD model across chains */
/* */
/******************************************************************************/
if(data->bayesLineFlag)
{
for(ifo=0; ifo<NI; ifo++)
{
for(i=0; i<data->N/2; i++) psd[ifo][i] = model[chain->index[0]]->Snf[ifo][i];
for(ic=1; ic<chain->NC; ic++)
{
copy_bayesline_params(bayesline[chain->index[0]][ifo], bayesline[chain->index[ic]][ifo]);
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<data->N/2; i++)
{
model[chain->index[ic]]->Snf[ifo][i] = model[chain->index[0]]->Snf[ifo][i];
model[chain->index[ic]]->SnS[ifo][i] = model[chain->index[0]]->SnS[ifo][i];
model[chain->index[ic]]->invSnf[ifo][i] = model[chain->index[0]]->invSnf[ifo][i];
}
}
}
}
//Shut off cleaning flag for checkpointing
data->cleanModelFlag = 0;
//Export Reference PSD and cleaned data for post-processing
export_cleaned_data(data, model[chain->index[0]]);
}
/******************************************************************************/
/* */
/* Here we redo the injections into the original data just to get an updated */
/* estimate of the SNR using the clean PSD. This does not affect the data */
/* used by BayesWave (which lives in the Data structure) because the */
/* injection codes only look at the LALInferenceRunState structure */
/* */
/******************************************************************************/
if(LALInferenceGetProcParamVal(runState->commandLine, "--inj"))
{
ifo=0;
dataPtr = runState->data;
while(dataPtr!=NULL)
{
for(i=0; i<N/2; i++) dataPtr->oneSidedNoisePowerSpectrum->data->data[i] = 2.0/Tobs/model[chain->index[0]]->invSnf[ifo][i];
dataPtr = dataPtr->next;
ifo++;
}
fprintf(stdout, " ==== LALInferenceInjectInspiralSignal(): started. ====\n");
LALInferenceInjectInspiralSignal(runState->data,runState->commandLine);
fprintf(stdout, " ==== LALInferenceInjectInspiralSignal(): finished. ====\n\n");
}
if(LALInferenceGetProcParamVal(runState->commandLine, "--MDC-cache"))
{
fprintf(stdout, " ====== Recompute burstMDC SNR =======\n");
SNRinj = malloc((NI+1)*sizeof(double));
ifo=0;
dataPtr = runState->data;
while(dataPtr!=NULL)
{
// Fill runState PSD with results from BayesLine
for(i=0; i<N/2; i++) dataPtr->oneSidedNoisePowerSpectrum->data->data[i] = 2.0/Tobs/model[chain->index[0]]->invSnf[ifo][i];
dataPtr = dataPtr->next;
ifo++;
}
//Redo the injection to get injected SNR
InjectFromMDC(runState->commandLine, runState->data, SNRinj);
double SNRmin = SNRinj[NI];
for(i=0; i<NI; i++) if(SNRinj[i]<SNRmin) SNRmin=SNRinj[i];
FILE *snrFile = NULL;
sprintf(filename,"%ssnr/snr.dat",data->runName);
snrFile = fopen(filename,"w");
fprintf(snrFile,"%.12g %.12g\n",SNRmin,SNRinj[NI]);
fclose(snrFile);
if(data->cleanModelFlag==0 && SNRinj[NI]>100.)
{
printf("MDC/XML injection with SNR>100\n");
printf(" Cleaning phase disabled\n");
printf(" Increasing number of chains to %i\n",chain->NC+10);
chain->NC+=10;
}
if(data->cleanModelFlag==0 && SNRinj[NI]>200.)
{
printf("MDC/XML injection with SNR>100\n");
printf(" Cleaning phase disabled\n");
printf(" Increasing number of chains to %i\n",chain->NC+10);
chain->NC+=10;
}
free(SNRinj);
}
/******************************************************************************/
/* */
/* Full PT & RJMCMC run for each model under consideration */
/* Evidence automatically calculated using thermodynamic integration */
/* */
/******************************************************************************/
// data->Dmax = Dmax;
FILE *evidence;
sprintf(filename,"%sevidence.dat",data->runName);
evidence = fopen(filename,"w");
/******************************************************************************/
/* */
/* Signal Characterization phase */
/* */
/******************************************************************************/
if(data->signalModelFlag)
{
data->runPhase = 1;
reset_priors(data, prior);
//Initialize proposal ratios for GW search phase
chain->modelRate = 1; //always do signal model updates
if(data->constantLogLFlag)
chain->rjmcmcRate = 0.01;
else
chain->rjmcmcRate = 0.5;
//Shut off RJMCMC if --fixD argument given
if(data->rjFlag==0) chain->rjmcmcRate = 0.0;
data->cleanFlag = 0;
data->glitchFlag = 0;
data->signalFlag = 1;
initialize_chain(chain,1);
reset_likelihood(data);
reset_priors(data, prior);
reset_model(data, chain, prior, model);
printf("characterizing signal model...\n\n");
RJMCMC(data, model, bayesline, chain, prior, &data->logZsignal, &data->varZsignal);
//shut off signal model for checkpointing
data->signalModelFlag = 0;
}
fprintf(evidence,"signal %.12g %lg\n",data->logZsignal,data->varZsignal);
fflush(evidence);
/******************************************************************************/
/* */
/* Glitch Characterization phase */
/* */
/******************************************************************************/
if(data->glitchModelFlag)
{
data->runPhase = 1;
reset_priors(data, prior);
//Initialize proposal ratios for GW search phase
chain->modelRate = 1; //even chance for each channel
if(data->constantLogLFlag)
chain->rjmcmcRate = 0.01;
else
chain->rjmcmcRate = 0.5;
//Shut off RJMCMC if --fixD argument given
if(data->rjFlag==0) chain->rjmcmcRate = 0.0;
data->cleanFlag = 0;
data->glitchFlag = 1;
data->signalFlag = 0;
initialize_chain(chain,1);
reset_likelihood(data);
reset_priors(data, prior);
reset_model(data, chain, prior, model);
printf("characterizing glitch model...\n\n");
RJMCMC(data, model, bayesline, chain, prior, &data->logZglitch, &data->varZglitch);
//shut off glitch model for checkpointing
data->glitchModelFlag = 0;
}
fprintf(evidence,"glitch %.12g %lg\n",data->logZglitch,data->varZglitch);
fflush(evidence);
/******************************************************************************/
/* */
/* Gaussian Noise Characterization phase */
/* */
/******************************************************************************/
if(data->noiseModelFlag)
{
data->runPhase = 1;
reset_priors(data, prior);
//Initialize proposal ratios for GW search phase
chain->rjmcmcRate = 0.0;
//set flags for no wavelets in model
data->cleanFlag = 0;
data->glitchFlag = 0;
data->signalFlag = 0;
initialize_chain(chain,1);
reset_likelihood(data);
reset_priors(data, prior);
reset_model(data, chain, prior, model);
printf("characterizing Gaussian noise model...\n\n");
RJMCMC(data, model, bayesline, chain, prior, &data->logZnoise, &data->varZnoise);
//shut off noise model for checkpointing
data->noiseModelFlag = 0;
}
fprintf(evidence,"noise %.12g %lg\n",data->logZnoise,data->varZnoise);
fflush(evidence);
/******************************************************************************/
/* */
/* The full model RJMCMC Run (Signal+Glitch) Characterization phase */
/* */
/******************************************************************************/
if(data->fullModelFlag)
{
data->runPhase = 1;
reset_priors(data, prior);
//Initialize proposal ratios for GW search phase
chain->modelRate = 0.5; //always do signal model updates during search
if(data->constantLogLFlag)
chain->rjmcmcRate = 0.01;
else
chain->rjmcmcRate = 0.5;
//Shut off RJMCMC if --fixD argument given
if(data->rjFlag==0) chain->rjmcmcRate = 0.0;
data->cleanFlag = 0;
data->glitchFlag = 1;
data->signalFlag = 1;
initialize_chain(chain,1);
reset_likelihood(data);
reset_priors(data, prior);
reset_model(data, chain, prior, model);
RJMCMC(data, model, bayesline, chain, prior, &data->logZfull, &data->varZfull);
fprintf(evidence,"full %.12g %lg\n",data->logZfull,data->varZfull);
fflush(evidence);
//shut off glitch+signal model for checkpointing
data->fullModelFlag = 0;
}
/******************************************************************************/
/* */
/* Output log- and normalized evidence */
/* */
/******************************************************************************/
fclose(evidence);
//output file format for stacked histogram
if(data->loudGlitchFlag) sprintf(filename,"%sevidence_stacked_ignore.dat",data->runName);
else sprintf(filename,"%sevidence_stacked.dat",data->runName);
evidence = fopen(filename,"w");
fprintf(evidence,"%.12g %.12g %.12g %lg %lg %lg\n",data->logZnoise,data->logZglitch,data->logZsignal,data->varZnoise,data->varZglitch,data->varZsignal);
fclose(evidence);
//output evidence files if loudGlitchFlag is triggered
if(data->loudGlitchFlag)
{
sprintf(filename,"%sevidence_stacked.dat",data->runName);
evidence = fopen(filename,"w");
fprintf(evidence,"%.12g %.12g %.12g %lg %lg %lg\n",0.0,10000.0,0.0,0.0,1.0,0.0);
fclose(evidence);
sprintf(filename,"%sevidence.dat",data->runName);
evidence = fopen(filename,"w");
fprintf(evidence,"signal %.12g %.12g\n",0.0,0.0);
fprintf(evidence,"glitch %.12g %.12g\n",10000.0,1.0);
fprintf(evidence,"noise %.12g %.12g\n",0.0,0.0);
fclose(evidence);
sprintf(filename,"%sevidence_ignore.dat",data->runName);
evidence = fopen(filename,"w");
fprintf(evidence,"signal %.12g %.12g\n",data->logZsignal,data->varZsignal);
fprintf(evidence,"glitch %.12g %.12g\n",data->logZglitch,data->varZglitch);
fprintf(evidence,"noise %.12g %.12g\n",data->logZnoise,data->varZnoise);
fclose(evidence);
}
//write gnuplot script to make single-run evidence histograms
sprintf(filename,"%sevidence.gpi",data->runName);
evidence = fopen(filename,"w");
write_evidence_gnuplot_script(evidence, data->runName);
/******************************************************************************/
/* */
/* Free memory, exit cleanly */
/* */
/******************************************************************************/
for(ic=0; ic<chain->NC; ic++)
{
free_model(model[ic], data, prior);
free(model[ic]);
}
free(model);
free_chain(chain,NI,data->Npol);
return 0;
}
void print_help_message(void)
{
fprintf(stdout,"\n ========== BayesWave USAGE: =========\n\n");
fprintf(stdout,"REQUIRED:\n");
fprintf(stdout," --ifo IFO interferometer (H1,L1,V1)\n");
fprintf(stdout," --IFO-flow minimum frequency (Hz)\n");
fprintf(stdout," --IFO-cache full path to cache file\n");
fprintf(stdout," --IFO-channel channel name (IFO:LSC-STRAIN)\n");
fprintf(stdout," --trigtime GPS trigger time\n");
fprintf(stdout," --srate sampling rate (Hz)\n");
fprintf(stdout," --seglen duration of data (s)\n");
fprintf(stdout," --PSDstart GPS start time for PSD estimation\n");
fprintf(stdout," --PSDlength duration of PSD estimation length\n");
fprintf(stdout," --dataseed Required if using simulated noise e.g. --H1-cache LALAdLIGO\n");
fprintf(stdout,"\n");
fprintf(stdout,"OPTIONAL:\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --- Run parameters -------------------------------------------------------------\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --segment-start GPS start time of segment (trigtime + 2 - seglen) \n");
fprintf(stdout," --Niter number of iterations (4000000)\n");
fprintf(stdout," --Nchain number of parallel chains (20)\n");
fprintf(stdout," --Ncycle number of model updates per RJMCMC iteration (100)\n");
fprintf(stdout," --Nburnin number of burn-in iterations (50000)\n");
fprintf(stdout," --maxLogL use maximized likelihood during burnin\n");
fprintf(stdout," --chainseed random number seed for Markov chain (1234)\n");
fprintf(stdout," --runName run name for output files\n");
fprintf(stdout," --outputDir absolute path of where the output should be written\n");
fprintf(stdout," --0noise no noise realization\n");
fprintf(stdout," --prior sample from prior using logL = constant test\n");
fprintf(stdout," --gnuplot output files for gnuplot animations\n");
fprintf(stdout," --verbose output hot chains\n");
fprintf(stdout," --window duration of time window for model characterization runs\n");
fprintf(stdout," --checkpoint enable self-checkpointing\n");
fprintf(stdout," --version print BayesWave version and exit\n");
fprintf(stdout," --help print BayesWave options and exit\n");
fprintf(stdout,"\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --- Model parameters -----------------------------------------------------------\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --fullOnly require signal && glitch model\n");
fprintf(stdout," --noClean skip cleaning phase and go right to reduced window\n");
fprintf(stdout," --noSignal skip signal phase and quit after glitch model\n");
fprintf(stdout," --cleanOnly run bayesline & glitch cleaning phase only\n");
fprintf(stdout," --noiseOnly use noise model only (no signal or glitches)\n");
fprintf(stdout," --signalOnly use signal model only (no glitches)\n");
fprintf(stdout," --glitchOnly use glitch model only (no signal)\n");
fprintf(stdout," --noPSDfit keep PSD parameters fixed\n");
fprintf(stdout," --bayesLine use BayesLine for PSD model\n");
fprintf(stdout," --stochastic use stochastic background model\n");
fprintf(stdout,"\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --- Priors & Proposals -----------------------------------------------------------\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --Dmin minimum number of wavelets total (1)\n");
fprintf(stdout," --Dmax maximum number of wavelets per channel (100)\n");
fprintf(stdout," --fixD fix # of wavelets for signal/IFO. Overrides --Dmin,--Dmax\n");
fprintf(stdout," --Qmin maximum quality factor for wavelets (0.1)\n");
fprintf(stdout," --Qmax maximum quality factor for wavelets (40)\n");
fprintf(stdout," --noWaveletPrior revert to uniform p(D)\n");
fprintf(stdout," --clusterPrior use metric-based clustering prior (BROKEN)\n");
fprintf(stdout," --clusterPath full path to cluster prior normalization files\n");
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," --waveletPrior use empirical distribution on number of wavelets (O1)\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");
fprintf(stdout," --noSignalAmplitudePrior use same SNR-dependent prior for signal & glitch model\n");
fprintf(stdout," --noAmplitudeProposal don't draw from SNR-dependent amplitude prior\n");
fprintf(stdout," --varyExtrinsicAmplitude update wavelet amplitudes in extrinsic update (FIXED?)\n");
fprintf(stdout," --noPolarization do not assume elliptical polarization for signal model\n");
fprintf(stdout," --noClusterProposal disable clustering & TF density proposal \n");
fprintf(stdout," --clusterWeight fractional weight for TF to proximity proposal (0.5)\n");
fprintf(stdout," --ampPriorPeak SNR where amplitude prior peaks (5)\n");
fprintf(stdout," --signalPriorPeak SNR where signal amplitude prior peaks (5)\n");
fprintf(stdout," --dimensionDecayRate e-folding time for number of wavelets (1000)\n");
fprintf(stdout," --fixIntrinsicParams hold intrinsic parameters fixed\n");
fprintf(stdout," --fixExtrinsicParams hold extrinsic parameters fixed\n");
fprintf(stdout," --loudGlitchPrior set floor of PSD prior to .01*Sn(f=200Hz)\n");
fprintf(stdout,"\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --- Parallel Tempering parameters ----------------------------------------------\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --tempMin minimum temperature chain (1)\n");
fprintf(stdout," --noAdaptTemperature disable adjust PT ladder to maintain acc. rate across chains\n");
fprintf(stdout," --tempSpacing set temperature spacing for geometric ladder (1.5)\n");
fprintf(stdout," --noSplineEvidence disable spline thermodynamic integration\n");
fprintf(stdout,"\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --- LALInference injection options ---------------------------------------------\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --inj injfile.xml Injection XML file to use\n");
fprintf(stdout," --event N Event number from Injection XML file to use\n");
fprintf(stdout,"\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --- Burst MDC injection options -------------------------------------------------\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --MDC-channel IFO1-chan, IFO2-chan, etc\n");
fprintf(stdout," --MDC-cache IFO1-mdcframe, IFO2-mdcframe, etc\n");
fprintf(stdout," --MDC-prefactor Rescale injection amplitude (1.0)\n");
fprintf(stdout,"\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --- BayesWave internal injection options ----------------------------------------\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --BW-inject (signal/glitch)\n");
fprintf(stdout," --BW-injName runName that produced the chain file \n");
fprintf(stdout," --BW-path Path to BW chain file for injection (./chains) \n");
fprintf(stdout," --BW-event Which sample from BayesWave chain (200000)\n");
fprintf(stdout,"\n");
fprintf(stdout,"EXAMPLE:\n");
fprintf(stdout," BayesWave --ifo H1 --H1-flow 32 --H1-cache LALSimAdLIGO --H1-channel LALSimAdLIGO --trigtime 900000000.00 --srate 512 --seglen 4 --PSDstart 900000000 --PSDlength 1024 --NCmin 2 --NCmax 2 --dataseed 1234\n");
fprintf(stdout,"\n");
}
/*
* Copyright (C) 2018 Neil J. Cornish, Tyson B. Littenberg, Margaret Millhouse,Katerina Chatziioannou
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with with program; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
* MA 02111-1307 USA
*/
/*************************** REQUIRED LIBRARIES ***************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <sys/stat.h>
#include <fftw3.h>
#include "BayesWave.h"
#include "BayesLine.h"
#include "BayesWaveIO.h"
#include "BayesWaveMCMC.h"
#include "BayesWaveMath.h"
#include "BayesWavePrior.h"
#include "BayesWaveModel.h"
#include "BayesWaveWavelet.h"
#include "BayesWaveProposal.h"
#include "BayesWaveEvidence.h"
#include "BayesWaveLikelihood.h"
/************* PROTOTYPE DECLARATIONS FOR INTERNAL FUNCTIONS **************/
void print_help_message(void);
/* ============================ MAIN PROGRAM ============================ */
int main(int argc, char *argv[])
{
/* Variable declaration */
int i, j, ic, ifo, imin,imax;
char filename[MAXSTRINGSIZE];
char modelname[MAXSTRINGSIZE];
/* LAL data structure that will contain all of the data from the frame files */
LALInferenceRunState *runState = XLALCalloc(1,sizeof(LALInferenceRunState));
runState->commandLine = LALInferenceParseCommandLine(argc,argv);
/* Print version number and command line to screen */
print_version(stdout);
fprintf(stdout," %s\n",LALInferencePrintCommandLine(runState->commandLine) );
fprintf(stdout,"\n");
/******************************************************************************/
/* */
/* Usage statment if --help */
/* */
/******************************************************************************/
if(LALInferenceGetProcParamVal(runState->commandLine, "--version"))
{
return 0;
}
if(LALInferenceGetProcParamVal(runState->commandLine, "--help") || argc==1)
{
print_help_message();
return 0;
}
/******************************************************************************/
/* */
/* Read/store data */
/* */
/******************************************************************************/
fprintf(stdout, "\n");
fprintf(stdout, " ======= LALInferenceReadData() ======\n");
runState->data = LALInferenceReadData(runState->commandLine);
fprintf(stdout, "\n");
/******************************************************************************/
/* */
/* Unpack runState into BayesWave data structures */
/* */
/******************************************************************************/
LALInferenceIFOData *dataPtr = NULL;
dataPtr = runState->data;
double *SNRinj = NULL;
int NI;
int N = dataPtr->timeData->data->length;
double Tobs = (double)N*dataPtr->timeData->deltaT;
ifo=0;
while(dataPtr!=NULL)
{
dataPtr = dataPtr->next;
ifo++;
}
NI=ifo;
double **freqData = double_matrix(NI-1,N-1);
double **psd = double_matrix(NI-1,N/2-1);
double **timeData = double_matrix(NI-1,N-1);
struct Data *data = malloc(sizeof(struct Data));
struct Chain *chain = malloc(sizeof(struct Chain));
struct Prior *prior = malloc(sizeof(struct Prior));
const gsl_rng_type *T = gsl_rng_default;
chain->seed = gsl_rng_alloc(T);
gsl_rng_env_setup();
parse_command_line(data, chain, prior, runState->commandLine);
if (data->chirpletFlag) {
data->NW = 6;
}
else{
data->NW = 5;
}
// Move to the directory specified by the user, or stay in ./
chdir(data->outputDirectory);
fprintf(stdout,"\nOutputting all data to: %s \n\n", data->outputDirectory);
// Create working directories for output
mode_t process_mask = umask(0);
mkdir("chains", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
mkdir("waveforms", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
if(data->checkpointFlag) mkdir("checkpoint",S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
sprintf(filename,"%ssnr",data->runName);
mkdir(filename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
umask(process_mask);
// burstMDC injection
if(LALInferenceGetProcParamVal(runState->commandLine, "--MDC-cache"))
{
SNRinj = malloc((NI+1)*sizeof(double));
fprintf(stdout, " ==== LALInferenceInjectFromMDC(): started. ====\n");
InjectFromMDC(runState->commandLine, runState->data, SNRinj);
fprintf(stdout, " ==== LALInferenceInjectFromMDC(): finished. ====\n\n");
free(SNRinj);
}
// CBC injection
if(LALInferenceGetProcParamVal(runState->commandLine, "--inj"))
{
fprintf(stdout, " ==== LALInferenceInjectInspiralSignal(): started. ====\n");
LALInferenceInjectInspiralSignal(runState->data,runState->commandLine);
fprintf(stdout, " ==== LALInferenceInjectInspiralSignal(): finished. ====\n\n");
}
//Output Fourier domain data for BayesLine
dataPtr = runState->data;
double fmin = dataPtr->fLow;
double fmax = ceil(dataPtr->fHigh);
//set frequency resolution to be 32 Hz because why not?
double df,dt;
df = 32.0;
dt = 0.5/df;
int tsize = (int)(Tobs/dt);
//Copy Fourier domain data & PSD from LAL data types into simple arrays
ifo=0;
imin = (int)(fmin*Tobs);
imax = (int)(fmax*Tobs);
bool simDataFlag = False;
while(dataPtr!=NULL)
{
//read in time domain data from LALInference data structre
for(i=0; i<N; i++)
{
timeData[ifo][i] = dataPtr->timeData->data->data[i];
//timeData remains uninitialized if simulating data
if(timeData[ifo][i]!=timeData[ifo][i])
{
simDataFlag = True;
timeData[ifo][i] = 0.0;
}
}
//read in frequency domain data from LALInference data structure
for(i=0; i<N/2; i++)
{
// if(simDataFlag)
// {
// dataPtr->freqData->data->data[imin] = dataPtr->freqData->data->data[imin+1];
// dataPtr->oneSidedNoisePowerSpectrum->data->data[imin] = dataPtr->oneSidedNoisePowerSpectrum->data->data[imin+1];
// }
if(i<imin+1)
{
freqData[ifo][2*i] = creal(dataPtr->freqData->data->data[imin+1]);
freqData[ifo][2*i+1] = cimag(dataPtr->freqData->data->data[imin+1]);
psd[ifo][i] = dataPtr->oneSidedNoisePowerSpectrum->data->data[imin+1]*Tobs/2.0;
}
else
{
freqData[ifo][2*i] = creal(dataPtr->freqData->data->data[i]);
freqData[ifo][2*i+1] = cimag(dataPtr->freqData->data->data[i]);
/*
Innerproducts are expecting <n_i^2> instead of Sn(f):
So psd arrays == <n_i^2> = T/2 * Sn(f)
*/
psd[ifo][i] = dataPtr->oneSidedNoisePowerSpectrum->data->data[i]*Tobs/2.0;
}
}
dataPtr = dataPtr->next;
ifo++;
}
if(simDataFlag)
{
fprintf(stdout,"Detected run is using simulated data\n");
fprintf(stdout,"copied d(fmin+1/t) to d(f<=fmin)\n\n");
}
/******************************************************************************/
/* */
/* Setup DATA, CHAIN, PRIOR, MODEL structures */
/* */
/******************************************************************************/
/*
Setup DATA structure
*/
initialize_data(data,freqData,N,tsize,Tobs,NI,fmin,fmax);
//print run summary file
sprintf(filename,"%sbayeswave.run",data->runName);
chain->runFile = fopen(filename,"w");
print_version(chain->runFile);
fprintf(chain->runFile," %s\n\n",LALInferencePrintCommandLine(data->commandLine) );
print_run_stats(chain->runFile, data, chain);
fclose(chain->runFile);
data->detector = malloc(NI*sizeof(LALDetector*));
data->epoch = runState->data->epoch;
ifo=0;
dataPtr = runState->data;
while(dataPtr!=NULL)
{
sprintf(data->ifos[ifo],"%s",dataPtr->name);
data->detector[ifo] = dataPtr->detector;
ifo++;
dataPtr = dataPtr->next;
}
OverlapReductionFunction(data);
/*
Setup CHAIN structure
*/
//give chain buffer in case we need more
int NC=chain->NC;
chain->NC+=20;
allocate_chain(chain,NI,data->Npol);
for(ic=0; ic<chain->NC; ic++) chain->index[ic]=ic;
/*
Setup PRIOR structure
*/
// int Dmax = data->Dmax;
// data->Dmax = 200;
int omax = 10;
if(data->Dmax < omax) omax = data->Dmax;
initialize_priors(data, prior, omax);
data->glitch_amplitude_proposal = draw_glitch_amplitude;
data->glitch_amplitude_prior = glitch_amplitude_prior;
// choose functions for amplitude priors
if(data->signalAmplitudePriorFlag)
{
data->signal_amplitude_proposal = draw_signal_amplitude;
data->signal_amplitude_prior = signal_amplitude_prior;
}
else
{
data->signal_amplitude_proposal = draw_glitch_amplitude;
data->signal_amplitude_prior = glitch_amplitude_prior;
}
/*
Setup MODEL structure
*/
struct Model **model = malloc(chain->NC*sizeof(struct Model*));
for(ic=0; ic<chain->NC; ic++)
{
model[ic] = malloc(sizeof(struct Model));
initialize_model(model[ic], data, prior, psd, chain->seed);
}
// BayesWave internal injection
if(LALInferenceGetProcParamVal(runState->commandLine,"--BW-inject"))
{
fprintf(stdout, " ==== BayesWaveInjection(): started. ====\n");
BayesWaveInjection(runState->commandLine, data, chain, prior, psd, &NC);
fprintf(stdout, " ==== BayesWaveInjection(): finished. ====\n\n");
}
// Print Fourier domain data to file for post-processing
for(ifo=0; ifo<data->NI; ifo++)
{
sprintf(filename,"waveforms/%sfourier_domain_data_%s.dat",data->runName,data->ifos[ifo]);
print_frequency_domain_data(filename, data->s[ifo], data->N, data->Tobs, data->imin, data->imax);
}
/******************************************************************************/
/* */
/* PSD search has to be done once even on restart */
/* Much of the BayesLine setup is done inside of BayesLineSearch() */
/* */
/******************************************************************************/
struct BayesLineParams ***bayesline = malloc(chain->NC*sizeof(struct BayesLineParams **));
if(data->bayesLineFlag)
{
fprintf(stdout,"\n ============ BayesLine ==============\n");
/*
Setup BayesLine structure
*/
for(ic=0; ic<chain->NC; ic++)
{
bayesline[ic] = malloc(data->NI*sizeof(struct BayesLineParams *));
initialize_bayesline(bayesline[ic], data, psd);
}
/*
BayesLine spectral estimation search-phase
*/
double *asd = malloc(N*sizeof(double));
for(ifo=0; ifo<NI; ifo++)
{
//Trick BayesLine into fitting to off source PSD.
for(i=0; i<data->N/2; i++)
{
asd[2*i] = data->s[ifo][2*i];//sqrt(psd[ifo][i]/2.);
asd[2*i+1] = data->s[ifo][2*i+1];//sqrt(psd[ifo][i]/2.);
}
for(i=0; i<(imax-imin); i++)
{
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];
//pad frequencies below fmin, just in case
for(i=0; i<=imin; i++)
{
timeData[ifo][2*i] = data->s[ifo][2*(imin+1)];
timeData[ifo][2*i+1] = data->s[ifo][2*(imin+1)+1];
}
fftw_wrapper(timeData[ifo],data->N,-1);
//normalize time series as expected by BayesLine
for(i=0; i<N; i++) timeData[ifo][i]/=data->Tobs;
fprintf(stdout,"BayesLine search phase for IFO %s\n", data->ifos[ifo]);
BayesLineBurnin(bayesline[0][ifo], timeData[ifo], asd, data->ifos[ifo]);
//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];
}
//Quit lying to BayesLine about the data. Poor BayesLine.
for(i=0; i<(imax-imin); i++)
{
j = i+imin;
bayesline[0][ifo]->power[i] = (data->s[ifo][2*j]*data->s[ifo][2*j]+data->s[ifo][2*j+1]*data->s[ifo][2*j+1]);
}
}
free(asd);
for(ifo=0; ifo<NI; ifo++)
{
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]);
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];
}
}
}
fprintf(stdout,"\n");
}
else
{
/*
Here is an ugly hack:
BayesWavePost expects a psd for its IFO-cache argument.
The bayeswave_pipe script defaults to using the fairdraw output from BayesLine
e.g. --[IFO]-cache interp:[rundir]/[IFO]_fairdraw_asd.dat
which doesn't exist when running w/out BayesLine, e.g. simulated data for review
purposes. Here we print a file with the right name and content so that bayeswave_pipe
doesn't have to decide what to use for the cache argument.
TODO: rethink naming to keep people out of trouble
*/
export_cleaned_data(data, model[0]);
}
//Restore chain size to requested value
chain->NC=NC;
/******************************************************************************/
/* */
/* PSD search & glitch cleaning stages */
/* Remove excess power from data outside of analysis window */
/* Residuals and PSDs are passed on to model selection phase */
/* Skip this step on restart */
/* */
/******************************************************************************/
// Restart run using stored residual and PSD
if(data->checkpointFlag)
{
fprintf(stdout, " ========= Checkpoint flag found ==========\n");
fflush(stdout);
if(data->stochasticFlag)
{
fprintf(stdout,"Stochastic background model not compatible with checkpointing\n");
fprintf(stdout,"Disabling checkpointFlag\n");
fprintf(stdout,"May the Condor gods be with you...\n");
data->checkpointFlag = 0;
}
else
{
FILE *fptr=NULL;
if( (fptr = fopen("checkpoint/state.dat","r")) == NULL )
{
fprintf(stdout,"No checkpoint files found\n");
fprintf(stdout,"Continue with full run\n");
data->resumeFlag = 0; //set resume flag to 0 so run starts cleanly
}
else
{
/****************************************************/
/* */
/* READ IN CURRENT STATE VECTOR */
/* */
/****************************************************/
fscanf(fptr,"%s", modelname);
fscanf(fptr,"%i",&data->cleanModelFlag);
fscanf(fptr,"%i",&data->noiseModelFlag);
fscanf(fptr,"%i",&data->glitchModelFlag);
fscanf(fptr,"%i",&data->signalModelFlag);
fscanf(fptr,"%i",&data->fullModelFlag);
fscanf(fptr,"%i",&data->cleanOnlyFlag);
fscanf(fptr,"%i",&data->loudGlitchFlag);
fclose(fptr);
/****************************************************/
/* */
/* PTMCMC TEMPERATURE LADDER */
/* */
/****************************************************/
//Read in dT, index, temperature, ptprop, ptacc, A, avg&var loglikelihood?
sprintf(filename,"checkpoint/temperature.dat");
if( (fptr = fopen(filename,"r")) == NULL )
{
fprintf(stderr,"Error: Could not checkpoint temperature model\n");
fprintf(stderr," Parameter file %s does not exist\n",filename);
fprintf(stderr," Exiting to system (1)\n");
exit(1);
}
fscanf(fptr,"%i",&chain->mc);
fscanf(fptr,"%i",&NC);
fscanf(fptr,"%i",&chain->zcount);
fclose(fptr);
sprintf(filename,"checkpoint/rng.bin");
if( (fptr = fopen(filename,"rb")) == NULL )
{
fprintf(stderr,"Error: Could not checkpoint the rng\n");
fprintf(stderr," Parameter file %s does not exist\n",filename);
fprintf(stderr," Exiting to system (1)\n");
exit(1);
}
gsl_rng_fread(fptr,chain->seed);
fclose(fptr);
/****************************************************/
/* */
/* READ IN CURRENT EVIDENCE RESULTS */
/* */
/****************************************************/
if( (fptr = fopen("checkpoint/evidence.dat","r")) == NULL )
{
fprintf(stderr,"Error: Could not checkpoint model settings\n");
fprintf(stderr," Parameter file checkpoint/evidence.dat does not exist\n");
fprintf(stderr," Exiting to system (1)\n");
exit(1);
}
fscanf(fptr,"%lg %lg",&data->logZnoise,&data->varZnoise);
fscanf(fptr,"%lg %lg",&data->logZglitch,&data->varZglitch);
fscanf(fptr,"%lg %lg",&data->logZsignal,&data->varZsignal);
fscanf(fptr,"%lg %lg",&data->logZfull,&data->varZfull);
fclose(fptr);
fprintf(stdout,"Checkpoint files found for model %s\n",modelname);
fprintf(stdout,"Resume previous run\n");
data->resumeFlag = 1; //set resume flag to 1 so run picks up at checkpoint
}
}
fprintf(stdout,"\n");
}
/******************************************************************************/
/* */
/* The Full-time MCMC Run (Cleaning phase) */
/* */
/******************************************************************************/
if(data->stochasticFlag)
{
fprintf(stdout,"Stochastic background model not compatible with glitch model\n");
fprintf(stdout,"Disabling cleanModelFlag\n");
fprintf(stdout,"May the glitch gods be with you...\n");
data->cleanModelFlag = 0;
data->noiseModelFlag = 1;
data->glitchModelFlag = 0;
data->signalModelFlag = 0;
data->fullModelFlag = 0;
}
if(data->cleanModelFlag)
{
//lie to the code about the number of chains & iterations
// data->Dmax = 200;
int M = chain->count;
int NC = chain->NC;
int CP = data->clusterPriorFlag;
data->clusterPriorFlag = 0;
chain->count = M;
chain->NC = NC/5;
if(chain->NC < 1)chain->NC=1;
//Initialize proposal ratios for glitch cleaning phase
chain->modelRate = 0.0; //even chance for each channel
if(data->constantLogLFlag)
chain->rjmcmcRate = 0.01;
else
chain->rjmcmcRate = 0.5;
//Shut off RJMCMC if --fixD argument given
if(data->rjFlag==0) chain->rjmcmcRate = 0.0;
data->cleanFlag = 1;
data->glitchFlag = 1;
data->signalFlag = 0;
data->runPhase = 0;
initialize_chain(chain,1);
reset_likelihood(data);
reset_priors(data, prior);
reset_model(data, chain, prior, model);
double cleanLogZ, cleanVarZ;
printf("characterizing PSD+glitch model...\n\n");
//TODO disable checkpointing and adaptation during cleaning phase
int saveAdaptationFlag = data->adaptTemperatureFlag;
data->adaptTemperatureFlag = 0;
RJMCMC(data, model, bayesline, chain, prior, &cleanLogZ, &cleanVarZ);
data->adaptTemperatureFlag = saveAdaptationFlag;
/******************************************************************************/
/* */
/* Remove excess power outside the trigger time, reset dimension priors */
/* */
/******************************************************************************/
data->runPhase = 1;
//Restore the number of chains & iterations
chain->count = M;
chain->NC = NC;
// data->Dmax = Dmax;
data->clusterPriorFlag = CP;
reset_priors(data, prior);
//Permanently regress any glitches outside of the restricted prior range
fprintf(stdout, "\n");
fprintf(stdout, " Removing excess power outside of t=[%g,%g]\n",prior->range[0][0],prior->range[0][1]);
//SNR of signal remaining in data
double rSNR=0.0;
double ifoSNR[NI];
double snrRatio;
for(ifo=0; ifo<NI; ifo++) ifoSNR[ifo]=0.0;
struct Model *m = NULL;
struct Wavelet *g = NULL;
double **r = malloc(NI*sizeof(double*)); //total residual (just for plotting)
for(ifo=0; ifo<NI; ifo++)
{
r[ifo] = malloc(N*sizeof(double));
for(i=0; i<N; i++) r[ifo][i] = data->s[ifo][i];
m = model[chain->index[0]];
g = m->glitch[ifo];
for(i=1; i<=g->size; i++)
{
if(g->intParams[i][0]<prior->range[0][0] || g->intParams[i][0]>prior->range[0][1])
{
//regress from the data & glitch model se we can compute window SNR
m->wavelet(data->s[ifo], g->intParams[i], N, -1, data->Tobs);
m->wavelet(g->templates, g->intParams[i], N, -1, data->Tobs);
}
//regress all wavelets from residual array for plotting
m->wavelet(r[ifo], g->intParams[i], N, -1, data->Tobs);
}
ifoSNR[ifo] += fourier_nwip(data->imin,data->imax,g->templates,g->templates,m->invSnf[ifo]);
rSNR += ifoSNR[ifo];
fprintf(stdout, " Residual SNR in %s = %g\n",data->ifos[ifo],sqrt(ifoSNR[ifo]));
}
fprintf(stdout, " Residual SNR = %g\n",sqrt(rSNR));
/*
BayesWave can not handle very loud, very assymetric SNRs.
loudGlitchFlag is set to 1 if the SNR > 100 and the SNRratio is > 5:1
This hack is only implemented for 2-detector networks.
*/
if(NI==2)
{
snrRatio = ifoSNR[0]/ifoSNR[1];
if(snrRatio < 1.0) snrRatio = 1./snrRatio;
if(sqrt(rSNR)>100. && snrRatio > 5. )
{
fprintf(stdout,"\n");
fprintf(stdout,"************************ WARNING ************************\n");
fprintf(stdout,"Cleaning phase found loud, assymetric residual power \n");
fprintf(stdout,"Model selection cannot handle this trigger \n");
fprintf(stdout,"Assume it is a glitch but proceed with analysis \n");
fprintf(stdout,"*********************************************************\n");
data->loudGlitchFlag = 1;
}
}
// Print residual to file
for(ifo=0; ifo<data->NI; ifo++)
{
sprintf(filename,"waveforms/%sfourier_domain_residual_%s.dat",data->runName,data->ifos[ifo]);
print_frequency_domain_data(filename, r[ifo], data->N, data->Tobs, data->imin, data->imax);
sprintf(filename,"waveforms/%sfourier_domain_clean_%s.dat",data->runName,data->ifos[ifo]);
print_frequency_domain_data(filename, data->s[ifo], data->N, data->Tobs, data->imin, data->imax);
free(r[ifo]);
}
free(r);
/******************************************************************************/
/* */
/* Increase the number of parallel chains if residual SNR is too high */
/* */
/******************************************************************************/
//if(sqrt(rSNR)>0.0) resize_model(data, chain, prior, model, bayesline, psd, chain->NC+10);
if(sqrt(rSNR)>100.0) chain->NC+=10;
if(sqrt(rSNR)>200.0) chain->NC+=10;
/******************************************************************************/
/* */
/* Copy T=0 chain's PSD model across chains */
/* */
/******************************************************************************/
if(data->bayesLineFlag)
{
for(ifo=0; ifo<NI; ifo++)
{
for(i=0; i<data->N/2; i++) psd[ifo][i] = model[chain->index[0]]->Snf[ifo][i];
for(ic=1; ic<chain->NC; ic++)
{
copy_bayesline_params(bayesline[chain->index[0]][ifo], bayesline[chain->index[ic]][ifo]);
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<data->N/2; i++)
{
model[chain->index[ic]]->Snf[ifo][i] = model[chain->index[0]]->Snf[ifo][i];
model[chain->index[ic]]->SnS[ifo][i] = model[chain->index[0]]->SnS[ifo][i];
model[chain->index[ic]]->invSnf[ifo][i] = model[chain->index[0]]->invSnf[ifo][i];
}
}
}
}
//Shut off cleaning flag for checkpointing
data->cleanModelFlag = 0;
//Export Reference PSD and cleaned data for post-processing
export_cleaned_data(data, model[chain->index[0]]);
}
/******************************************************************************/
/* */
/* Here we redo the injections into the original data just to get an updated */
/* estimate of the SNR using the clean PSD. This does not affect the data */
/* used by BayesWave (which lives in the Data structure) because the */
/* injection codes only look at the LALInferenceRunState structure */
/* */
/******************************************************************************/
if(LALInferenceGetProcParamVal(runState->commandLine, "--inj"))
{
ifo=0;
dataPtr = runState->data;
while(dataPtr!=NULL)
{
for(i=0; i<N/2; i++) dataPtr->oneSidedNoisePowerSpectrum->data->data[i] = 2.0/Tobs/model[chain->index[0]]->invSnf[ifo][i];
dataPtr = dataPtr->next;
ifo++;
}
fprintf(stdout, " ==== LALInferenceInjectInspiralSignal(): started. ====\n");
LALInferenceInjectInspiralSignal(runState->data,runState->commandLine);
fprintf(stdout, " ==== LALInferenceInjectInspiralSignal(): finished. ====\n\n");
}
if(LALInferenceGetProcParamVal(runState->commandLine, "--MDC-cache"))
{
fprintf(stdout, " ====== Recompute burstMDC SNR =======\n");
SNRinj = malloc((NI+1)*sizeof(double));
ifo=0;
dataPtr = runState->data;
while(dataPtr!=NULL)
{
// Fill runState PSD with results from BayesLine
for(i=0; i<N/2; i++) dataPtr->oneSidedNoisePowerSpectrum->data->data[i] = 2.0/Tobs/model[chain->index[0]]->invSnf[ifo][i];
dataPtr = dataPtr->next;
ifo++;
}
//Redo the injection to get injected SNR
InjectFromMDC(runState->commandLine, runState->data, SNRinj);
double SNRmin = SNRinj[NI];
for(i=0; i<NI; i++) if(SNRinj[i]<SNRmin) SNRmin=SNRinj[i];
FILE *snrFile = NULL;
sprintf(filename,"%ssnr/snr.dat",data->runName);
snrFile = fopen(filename,"w");
fprintf(snrFile,"%.12g %.12g\n",SNRmin,SNRinj[NI]);
fclose(snrFile);
if(data->cleanModelFlag==0 && SNRinj[NI]>100.)
{
printf("MDC/XML injection with SNR>100\n");
printf(" Cleaning phase disabled\n");
printf(" Increasing number of chains to %i\n",chain->NC+10);
chain->NC+=10;
}
if(data->cleanModelFlag==0 && SNRinj[NI]>200.)
{
printf("MDC/XML injection with SNR>100\n");
printf(" Cleaning phase disabled\n");
printf(" Increasing number of chains to %i\n",chain->NC+10);
chain->NC+=10;
}
free(SNRinj);
}
/******************************************************************************/
/* */
/* Full PT & RJMCMC run for each model under consideration */
/* Evidence automatically calculated using thermodynamic integration */
/* */
/******************************************************************************/
// data->Dmax = Dmax;
FILE *evidence;
sprintf(filename,"%sevidence.dat",data->runName);
evidence = fopen(filename,"w");
/******************************************************************************/
/* */
/* Signal Characterization phase */
/* */
/******************************************************************************/
if(data->signalModelFlag)
{
data->runPhase = 1;
reset_priors(data, prior);
//Initialize proposal ratios for GW search phase
chain->modelRate = 1; //always do signal model updates
if(data->constantLogLFlag)
chain->rjmcmcRate = 0.01;
else
chain->rjmcmcRate = 0.5;
//Shut off RJMCMC if --fixD argument given
if(data->rjFlag==0) chain->rjmcmcRate = 0.0;
data->cleanFlag = 0;
data->glitchFlag = 0;
data->signalFlag = 1;
initialize_chain(chain,1);
reset_likelihood(data);
reset_priors(data, prior);
reset_model(data, chain, prior, model);
printf("characterizing signal model...\n\n");
RJMCMC(data, model, bayesline, chain, prior, &data->logZsignal, &data->varZsignal);
//shut off signal model for checkpointing
data->signalModelFlag = 0;
}
fprintf(evidence,"signal %.12g %lg\n",data->logZsignal,data->varZsignal);
fflush(evidence);
/******************************************************************************/
/* */
/* Glitch Characterization phase */
/* */
/******************************************************************************/
if(data->glitchModelFlag)
{
data->runPhase = 1;
reset_priors(data, prior);
//Initialize proposal ratios for GW search phase
chain->modelRate = 1; //even chance for each channel
if(data->constantLogLFlag)
chain->rjmcmcRate = 0.01;
else
chain->rjmcmcRate = 0.5;
//Shut off RJMCMC if --fixD argument given
if(data->rjFlag==0) chain->rjmcmcRate = 0.0;
data->cleanFlag = 0;
data->glitchFlag = 1;
data->signalFlag = 0;
initialize_chain(chain,1);
reset_likelihood(data);
reset_priors(data, prior);
reset_model(data, chain, prior, model);
printf("characterizing glitch model...\n\n");
RJMCMC(data, model, bayesline, chain, prior, &data->logZglitch, &data->varZglitch);
//shut off glitch model for checkpointing
data->glitchModelFlag = 0;
}
fprintf(evidence,"glitch %.12g %lg\n",data->logZglitch,data->varZglitch);
fflush(evidence);
/******************************************************************************/
/* */
/* Gaussian Noise Characterization phase */
/* */
/******************************************************************************/
if(data->noiseModelFlag)
{
data->runPhase = 1;
reset_priors(data, prior);
//Initialize proposal ratios for GW search phase
chain->rjmcmcRate = 0.0;
//set flags for no wavelets in model
data->cleanFlag = 0;
data->glitchFlag = 0;
data->signalFlag = 0;
initialize_chain(chain,1);
reset_likelihood(data);
reset_priors(data, prior);
reset_model(data, chain, prior, model);
printf("characterizing Gaussian noise model...\n\n");
RJMCMC(data, model, bayesline, chain, prior, &data->logZnoise, &data->varZnoise);
//shut off noise model for checkpointing
data->noiseModelFlag = 0;
}
fprintf(evidence,"noise %.12g %lg\n",data->logZnoise,data->varZnoise);
fflush(evidence);
/******************************************************************************/
/* */
/* The full model RJMCMC Run (Signal+Glitch) Characterization phase */
/* */
/******************************************************************************/
if(data->fullModelFlag)
{
data->runPhase = 1;
reset_priors(data, prior);
//Initialize proposal ratios for GW search phase
chain->modelRate = 0.5; //always do signal model updates during search
if(data->constantLogLFlag)
chain->rjmcmcRate = 0.01;
else
chain->rjmcmcRate = 0.5;
//Shut off RJMCMC if --fixD argument given
if(data->rjFlag==0) chain->rjmcmcRate = 0.0;
data->cleanFlag = 0;
data->glitchFlag = 1;
data->signalFlag = 1;
initialize_chain(chain,1);
reset_likelihood(data);
reset_priors(data, prior);
reset_model(data, chain, prior, model);
RJMCMC(data, model, bayesline, chain, prior, &data->logZfull, &data->varZfull);
fprintf(evidence,"full %.12g %lg\n",data->logZfull,data->varZfull);
fflush(evidence);
//shut off glitch+signal model for checkpointing
data->fullModelFlag = 0;
}
/******************************************************************************/
/* */
/* Output log- and normalized evidence */
/* */
/******************************************************************************/
fclose(evidence);
//output file format for stacked histogram
if(data->loudGlitchFlag) sprintf(filename,"%sevidence_stacked_ignore.dat",data->runName);
else sprintf(filename,"%sevidence_stacked.dat",data->runName);
evidence = fopen(filename,"w");
fprintf(evidence,"%.12g %.12g %.12g %lg %lg %lg\n",data->logZnoise,data->logZglitch,data->logZsignal,data->varZnoise,data->varZglitch,data->varZsignal);
fclose(evidence);
//output evidence files if loudGlitchFlag is triggered
if(data->loudGlitchFlag)
{
sprintf(filename,"%sevidence_stacked.dat",data->runName);
evidence = fopen(filename,"w");
fprintf(evidence,"%.12g %.12g %.12g %lg %lg %lg\n",0.0,10000.0,0.0,0.0,1.0,0.0);
fclose(evidence);
sprintf(filename,"%sevidence.dat",data->runName);
evidence = fopen(filename,"w");
fprintf(evidence,"signal %.12g %.12g\n",0.0,0.0);
fprintf(evidence,"glitch %.12g %.12g\n",10000.0,1.0);
fprintf(evidence,"noise %.12g %.12g\n",0.0,0.0);
fclose(evidence);
sprintf(filename,"%sevidence_ignore.dat",data->runName);
evidence = fopen(filename,"w");
fprintf(evidence,"signal %.12g %.12g\n",data->logZsignal,data->varZsignal);
fprintf(evidence,"glitch %.12g %.12g\n",data->logZglitch,data->varZglitch);
fprintf(evidence,"noise %.12g %.12g\n",data->logZnoise,data->varZnoise);
fclose(evidence);
}
//write gnuplot script to make single-run evidence histograms
sprintf(filename,"%sevidence.gpi",data->runName);
evidence = fopen(filename,"w");
write_evidence_gnuplot_script(evidence, data->runName);
/******************************************************************************/
/* */
/* Free memory, exit cleanly */
/* */
/******************************************************************************/
for(ic=0; ic<chain->NC; ic++)
{
free_model(model[ic], data, prior);
free(model[ic]);
}
free(model);
free_chain(chain,NI,data->Npol);
return 0;
}
void print_help_message(void)
{
fprintf(stdout,"\n ========== BayesWave USAGE: =========\n\n");
fprintf(stdout,"REQUIRED:\n");
fprintf(stdout," --ifo IFO interferometer (H1,L1,V1)\n");
fprintf(stdout," --IFO-flow minimum frequency (Hz)\n");
fprintf(stdout," --IFO-cache full path to cache file\n");
fprintf(stdout," --IFO-channel channel name (IFO:LSC-STRAIN)\n");
fprintf(stdout," --trigtime GPS trigger time\n");
fprintf(stdout," --srate sampling rate (Hz)\n");
fprintf(stdout," --seglen duration of data (s)\n");
fprintf(stdout," --psdstart GPS start time for PSD estimation\n");
fprintf(stdout," --psdlength duration of PSD estimation length\n");
fprintf(stdout," --dataseed Required if using simulated noise e.g. --H1-cache LALAdLIGO\n");
fprintf(stdout,"\n");
fprintf(stdout,"OPTIONAL:\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --- Run parameters -------------------------------------------------------------\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --segment-start GPS start time of segment (trigtime + 2 - seglen) \n");
fprintf(stdout," --Niter number of iterations (4000000)\n");
fprintf(stdout," --Nchain number of parallel chains (20)\n");
fprintf(stdout," --Ncycle number of model updates per RJMCMC iteration (100)\n");
fprintf(stdout," --Nburnin number of burn-in iterations (50000)\n");
fprintf(stdout," --maxLogL use maximized likelihood during burnin\n");
fprintf(stdout," --chainseed random number seed for Markov chain (1234)\n");
fprintf(stdout," --runName run name for output files\n");
fprintf(stdout," --outputDir absolute path of where the output should be written\n");
fprintf(stdout," --0noise no noise realization\n");
fprintf(stdout," --prior sample from prior using logL = constant test\n");
fprintf(stdout," --gnuplot output files for gnuplot animations\n");
fprintf(stdout," --verbose output hot chains\n");
fprintf(stdout," --window duration of time window for model characterization runs\n");
fprintf(stdout," --checkpoint enable self-checkpointing\n");
fprintf(stdout," --version print BayesWave version and exit\n");
fprintf(stdout," --help print BayesWave options and exit\n");
fprintf(stdout,"\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --- Model parameters -----------------------------------------------------------\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --fullOnly require signal && glitch model\n");
fprintf(stdout," --noClean skip cleaning phase and go right to reduced window\n");
fprintf(stdout," --noSignal skip signal phase and quit after glitch model\n");
fprintf(stdout," --cleanOnly run bayesline & glitch cleaning phase only\n");
fprintf(stdout," --noiseOnly use noise model only (no signal or glitches)\n");
fprintf(stdout," --signalOnly use signal model only (no glitches)\n");
fprintf(stdout," --glitchOnly use glitch model only (no signal)\n");
fprintf(stdout," --noPSDfit keep PSD parameters fixed\n");
fprintf(stdout," --bayesLine use BayesLine for PSD model\n");
fprintf(stdout," --stochastic use stochastic background model\n");
fprintf(stdout,"\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --- Priors & Proposals -----------------------------------------------------------\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --Dmin minimum number of wavelets total (1)\n");
fprintf(stdout," --Dmax maximum number of wavelets per channel (100)\n");
fprintf(stdout," --fixD fix # of wavelets for signal/IFO. Overrides --Dmin,--Dmax\n");
fprintf(stdout," --Qmin maximum quality factor for wavelets (0.1)\n");
fprintf(stdout," --Qmax maximum quality factor for wavelets (40)\n");
fprintf(stdout," --noWaveletPrior revert to uniform p(D)\n");
fprintf(stdout," --clusterPrior use metric-based clustering prior (BROKEN)\n");
fprintf(stdout," --clusterPath full path to cluster prior normalization files\n");
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," --waveletPrior use empirical distribution on number of wavelets (O1)\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");
fprintf(stdout," --noSignalAmplitudePrior use same SNR-dependent prior for signal & glitch model\n");
fprintf(stdout," --noAmplitudeProposal don't draw from SNR-dependent amplitude prior\n");
fprintf(stdout," --varyExtrinsicAmplitude update wavelet amplitudes in extrinsic update (FIXED?)\n");
fprintf(stdout," --noPolarization do not assume elliptical polarization for signal model\n");
fprintf(stdout," --noClusterProposal disable clustering & TF density proposal \n");
fprintf(stdout," --clusterWeight fractional weight for TF to proximity proposal (0.5)\n");
fprintf(stdout," --ampPriorPeak SNR where amplitude prior peaks (5)\n");
fprintf(stdout," --signalPriorPeak SNR where signal amplitude prior peaks (5)\n");
fprintf(stdout," --dimensionDecayRate e-folding time for number of wavelets (1000)\n");
fprintf(stdout," --fixIntrinsicParams hold intrinsic parameters fixed\n");
fprintf(stdout," --fixExtrinsicParams hold extrinsic parameters fixed\n");
fprintf(stdout," --loudGlitchPrior set floor of PSD prior to .01*Sn(f=200Hz)\n");
fprintf(stdout,"\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --- Parallel Tempering parameters ----------------------------------------------\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --tempMin minimum temperature chain (1)\n");
fprintf(stdout," --noAdaptTemperature disable adjust PT ladder to maintain acc. rate across chains\n");
fprintf(stdout," --tempSpacing set temperature spacing for geometric ladder (1.5)\n");
fprintf(stdout," --noSplineEvidence disable spline thermodynamic integration\n");
fprintf(stdout,"\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --- LALInference injection options ---------------------------------------------\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --inj injfile.xml Injection XML file to use\n");
fprintf(stdout," --event N Event number from Injection XML file to use\n");
fprintf(stdout,"\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --- Burst MDC injection options -------------------------------------------------\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --MDC-channel IFO1-chan, IFO2-chan, etc\n");
fprintf(stdout," --MDC-cache IFO1-mdcframe, IFO2-mdcframe, etc\n");
fprintf(stdout," --MDC-prefactor Rescale injection amplitude (1.0)\n");
fprintf(stdout,"\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --- BayesWave internal injection options ----------------------------------------\n");
fprintf(stdout," ----------------------------------------------------------------------------------\n");
fprintf(stdout," --BW-inject (signal/glitch)\n");
fprintf(stdout," --BW-injName runName that produced the chain file \n");
fprintf(stdout," --BW-path Path to BW chain file for injection (./chains) \n");
fprintf(stdout," --BW-event Which sample from BayesWave chain (200000)\n");
fprintf(stdout,"\n");
fprintf(stdout,"EXAMPLE:\n");
fprintf(stdout," BayesWave --ifo H1 --H1-flow 32 --H1-cache LALSimAdLIGO --H1-channel LALSimAdLIGO --trigtime 900000000.00 --srate 512 --seglen 4 --PSDstart 900000000 --PSDlength 1024 --NCmin 2 --NCmax 2 --dataseed 1234\n");
fprintf(stdout,"\n");
}