adding O2_cleaning tag for glitch subtraction

git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/tags/O2_cleaning@734 c56465c9-8126-4a4f-9d7d-ac845eff4865

src/BayesLine.h

+#include <stdio.h>
+#include <string.h>
+#include <stdlib.h>
+#include <math.h>
+#include <unistd.h>
+#include <getopt.h>
+
+#include <gsl/gsl_sort.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_spline.h>
+#include <gsl/gsl_statistics.h>
+
+//#define SAmin  1.0e-48
+//#define SAmax  1.0e-30
+//#define LQmin  1.0e2  // allowing Q's below ~ 1e2 starts a fight with the spline model as lines get too fat
+//#define LQmax 1.0e8
+//#define LAmin  1.0e-44
+//#define LAmax  1.0e-30
+#define lAwidth 2.3    // 2.3 is one decade
+#define zeta 1.0
+#define kappa_BL 0.8
+#define FSTEP 10.0   // the stencil separation in Hz for the spline model. Going below 2 Hz is dangerous - will fight with line model
+//#define NSTEP 2000
+// as we are getting down to the line width of the broadest lines
+
+
+typedef struct
+{
+  int n;
+  int size;
+
+  int *larray;
+
+  double *Q;
+  double *A;
+  double *f;
+
+}lorentzianParams;
+
+typedef struct
+{
+  int tmax;
+  int ncut;
+  int nmin;
+  int tfull;
+  int sgmts;
+
+  double df;
+  double fny;
+  double Tobs;
+  double fmin;
+  double fmax;
+  double flow;
+  double fgrid;
+  double fstep;
+  double fhigh;
+  double cadence;
+
+}dataParams;
+
+typedef struct
+{
+  int n;
+  double *points;
+  double *data;
+
+}splineParams;
+
+typedef struct
+{
+  double SAmin;
+  double SAmax;
+  double LQmin;
+  double LQmax;
+  double LAmin;
+  double LAmax;
+
+  //double *invsigma; //variances for each frequency bin
+  double *sigma; //variances for each frequency bin
+  double *upper; //variances for each frequency bin
+  double *lower; //variances for each frequency bin
+  double *mean;     //means for each frequency bin
+
+}BayesLinePriors;
+
+struct BayesLineParams
+{
+  dataParams *data;
+  splineParams *spline;
+  splineParams *spline_x;
+  lorentzianParams *lines_x;
+  lorentzianParams *lines_full;
+  BayesLinePriors *priors;
+
+  double *Snf;
+  double *Sna;
+  double *fa;
+  double *freq;
+  double *power;
+  double *spow;
+  double *sfreq;
+  double *Sbase;
+  double *Sline;
+
+  int constantLogLFlag;
+  
+  int nstep;
+  
+  double TwoDeltaT;
+  gsl_rng *r;
+};
+
+void BayesLineFree(struct BayesLineParams *bptr);
+void BayesLineSetup(struct BayesLineParams *bptr, double *freqData, double fmin, double fmax, double deltaT, double Tobs, int nstep);
+
+void BayesLineRJMCMC(struct BayesLineParams *bayesline, double *freqData, double *psd, double *invpsd, double *splinePSD, int N, int cycle, double beta, int priorFlag);
+void BayesLineSearch(struct BayesLineParams *bayesline);
+
+void BayesLineNonMarkovianFit          (struct BayesLineParams *bayesline, int *nj);
+void BayesLineLorentzSplineMCMC        (struct BayesLineParams *bayesline, double heat, int steps, int focus, int priorFlag, double *dan);
+void BayesLineMarkovianSplineOnly      (struct BayesLineParams *bayesline, int nspline, int jj);
+void BayesLineMarkovianFocusedAnalysis (struct BayesLineParams *bayesline);
+
+double loglike_fit_spline(double *respow, double *Snf, int ncut);
+
+double loglike_pm        (double *respow, double *Sn, double *Snx, int ilow, int ihigh);
+double loglike_single    (double *respow, double *Sn, double *Snx, int ilowx, int ihighx, int ilowy, int ihighy);
+
+double sample(double *fprop, double pmax, dataParams *data, gsl_rng *r);
+double lprop(double f, double *fprop, dataParams *data);
+
+void full_spectrum_single(double *Sn, double *Snx, double *Sbasex, double *sfreq, dataParams *data, lorentzianParams *line_x, lorentzianParams *line_y, int ii, int *ilowx, int *ihighx, int *ilowy, int *ihighy);
+void full_spectrum_add_or_subtract(double *Snew, double *Sold, double *Sbase, double *sfreq, dataParams *restrict data, lorentzianParams *restrict lines, int ii, int *ilow, int *ihigh, int flag);
+void full_spectrum_spline(double *Sline, double *Sbase, double *sfreq, dataParams *restrict data, lorentzianParams *restrict lines);
+
+void spectrum_spline(double *Sn, double *Sbase, double *sfreq, dataParams *data, lorentzianParams *restrict lines, splineParams *restrict spline);
+
+void SpecFitSpline    (BayesLinePriors *priors, int zeroLogL, int steps, double *freq, double *power, splineParams *spline, double *Snf, int ncut, gsl_rng *r);
+void LorentzSplineFit (BayesLinePriors *priors, int zeroLogL, int steps, dataParams *data, lorentzianParams *lines, splineParams *spline, double *sfreq, double *spow, gsl_rng *r);
+
+void CubicSplineGSL(int N, double *x, double *y, int Nint, double *xint, double *yint);
+
+void create_dataParams(dataParams *data, double *f, int n);
+
+void create_lorentzianParams(lorentzianParams *lines, int size);
+void copy_lorentzianParams(lorentzianParams *origin, lorentzianParams *copy);
+void destroy_lorentzianParams(lorentzianParams *lines);
+
+void create_splineParams(splineParams *spline, int size);
+void copy_splineParams(splineParams *origin, splineParams *copy);
+void destroy_splineParams(splineParams *spline);
+
+void copy_bayesline_params(struct BayesLineParams *origin, struct BayesLineParams *copy);
+void print_line_model(FILE *fptr, struct BayesLineParams *bayesline);
+void print_spline_model(FILE *fptr, struct BayesLineParams *bayesline);
+void parse_line_model(FILE *fptr, struct BayesLineParams *bayesline);
+void parse_spline_model(FILE *fptr, struct BayesLineParams *bayesline);
+
+
+

src/BayesWave.h

+
+/* ********************************************************************************** */
+/*                                                                                    */
+/*                                  LAL Dependencies                                  */
+/*                                                                                    */
+/* ********************************************************************************** */
+
+#include <lal/LALInspiral.h>
+#include <lal/LALCache.h>
+#include <lal/LALFrStream.h>
+#include <lal/TimeFreqFFT.h>
+#include <lal/LALDetectors.h>
+#include <lal/ResampleTimeSeries.h>
+#include <lal/TimeSeries.h>
+#include <lal/FrequencySeries.h>
+#include <lal/Units.h>
+#include <lal/Date.h>
+#include <lal/StringInput.h>
+#include <lal/VectorOps.h>
+#include <lal/Random.h>
+#include <lal/LALNoiseModels.h>
+#include <lal/XLALError.h>
+#include <lal/GenerateInspiral.h>
+#include <lal/LIGOLwXMLRead.h>
+#include <lal/LIGOLwXMLInspiralRead.h>
+#include <lal/LALConstants.h>
+#include <lal/SeqFactories.h>
+#include <lal/DetectorSite.h>
+#include <lal/GenerateInspiral.h>
+#include <lal/GeneratePPNInspiral.h>
+#include <lal/SimulateCoherentGW.h>
+#include <lal/Inject.h>
+#include <lal/LIGOMetadataTables.h>
+#include <lal/LIGOMetadataUtils.h>
+#include <lal/LIGOMetadataInspiralUtils.h>
+#include <lal/LIGOMetadataRingdownUtils.h>
+#include <lal/LALInspiralBank.h>
+#include <lal/FindChirp.h>
+#include <lal/LALInspiralBank.h>
+#include <lal/GenerateInspiral.h>
+#include <lal/NRWaveInject.h>
+#include <lal/GenerateInspRing.h>
+#include <lal/LALError.h>
+
+#include <lal/LALInference.h>
+#include <lal/LALInferenceReadData.h>
+
+/* ********************************************************************************** */
+/*                                                                                    */
+/*                            Globally defined parameters                             */
+/*                                                                                    */
+/* ********************************************************************************** */
+
+#define RTPI 1.772453850905516       // sqrt(Pi)
+#define RT2PI 2.5066282746310005024  // sqrt(2 Pi)
+#define CLIGHT 2.99792458e8          // Speed of light (m/s)
+
+#define NE 6 //Number of extrinsic parameters
+#define NW 5 //Number of intrinsic parameters
+
+//Tuning parameters for wavelet proximity proposal
+//TODO: These should not be global!
+
+#define UNIFORM 0.2    // fraction of proximity prior that is uniform
+#define SPRD 4         // width of outer part of proximity proposal
+#define HLW 1          // width of hollowed-out region
+
+/* ********************************************************************************** */
+/*                                                                                    */
+/*                                  Data structures                                   */
+/*                                                                                    */
+/* ********************************************************************************** */
+struct Wavelet
+{
+   int size;
+   int smax;
+   int *index;
+
+   //clustering prior
+   int Ncluster;
+   double cosy;
+   double logp;
+
+   double *templates;
+   double **intParams;
+
+   char *name;
+};
+
+struct Model
+{
+  int size;
+
+  double *Sn;
+  double **Snf;
+  double **invSnf;
+  double **SnS;
+  double *SnGeo;
+
+  double logL;
+  double logLnorm;
+  double *detLogL;
+  double *extParams;
+  double **response;
+
+  struct Wavelet *signal;
+  struct Wavelet **glitch;
+
+  struct Network *projection;
+
+  struct FisherMatrix *fisher;
+  struct FisherMatrix *intfisher;
+
+  struct Background *background;
+
+  void (*wavelet_bandwidth)(double *, double *, double *);
+  void (*wavelet)(double *, double *, int, int, double);
+
+};
+
+struct Network
+{
+  double **expPhase;
+  double *deltaT;
+  double *Fplus;
+  double *Fcross;
+  double *dtimes;
+};
+
+struct Chain
+{
+  int mc;
+  int NC;     // Number of parallel chains
+  int burn;   // Number of burn in samples
+  int mod0;   // Counter keeping track of noise model
+  int mod1;   // Counter keeping track of glitch model
+  int mod2;   // Counter keeping track of signal model
+  int mod3;   // Counter keeping track of glitch+signal model
+
+  int cycle;  // # of model updates per BurstRJMCMC iteration
+  int count;  // # of BurstRJMCMC iterations
+  int blcount;// # of BayesLine burnin iterations
+  int zcount; // # of logL samples for thermodynamic integration
+  int mcount; // # of MCMC attempts
+  int scount; // # of RJMCMC attempts
+  int xcount; // # of extrinsic attempts
+  int fcount; // # of intrinsic fisher attempts
+  int pcount; // # of intrinsic phase attempts
+  int ccount; // # of cluster proposal attempts
+  int dcount; // # of density proposal attempts
+  int ucount; // # of uniform proposal attempts
+  int macc;   // # of MCMC successess
+  int sacc;   // # of RJMCMC successes
+  int xacc;   // # of extrinsic successes
+  int facc;   // # of intrinsic fisher successes
+  int pacc;   // # of intrinsic phase successes
+  int cacc;   // # of cluster proposal successes
+  int dacc;   // # of cluster proposal successes
+  int uacc;   // # of uniform proposal successes
+
+  int burnFlag;// flag telling proposals if we are in burn-in phase
+
+  int *index; // keep track of which chain is at which temperature
+  int *ptacc;
+  int *ptprop;
+
+  double *A;
+
+  long seed;
+
+  double modelRate;   //fraction of updates to signal model
+  double rjmcmcRate;  //fraction of trans- vs. inter-dimensional moves
+  double intPropRate; //fraction of prior vs. fisher intrinsic moves
+
+  double *dT;
+  double beta;
+  double tempMin;
+  double tempStep;
+  double *temperature;
+  double *avgLogLikelihood;
+  double *varLogLikelihood;
+
+  int nPoints;// = 3*chain->count/4/chain->cycle;
+  double **logLchain;//=dmatrix(0,chain->NC-1,0,nPoints-1);
+
+
+  FILE *runFile;
+  FILE **intChainFile;
+  FILE **intWaveChainFile;
+  FILE ***intGlitchChainFile;
+  FILE ***lineChainFile;
+  FILE ***splineChainFile;
+
+};
+
+struct Prior
+{
+  int smin;
+  int smax;
+  int *gmin;
+  int *gmax;
+  int omax;
+
+  double Snmin;
+  double Snmax;
+
+  double TFV;
+  double logTFV;
+
+  double *bias;
+  double **range;
+
+  //background prior
+    char bkg_name[1024];
+    double bkg_df;
+    double bkg_dQ;
+    int bkg_nf;
+    int bkg_nQ;
+    double **bkg_dist;
+
+  //cluster prior parameters
+  double alpha;
+  double beta;
+  double gamma;
+
+  //amplitude prior parameters
+  double sSNRpeak;
+
+  //glitch prior parameters
+  double gSNRpeak;
+
+  //dimension prior parameters
+  double Dtau;
+    double *Nwavelet;
+
+  char path[100];
+
+};
+
+struct Data
+{
+  int N;
+  int NI;
+  int tsize;
+  int fsize;
+  int imin;
+  int imax;
+  int Dmax;
+  int Dmin;
+
+  int runPhase;
+  int psdFitFlag;
+  int noiseSimFlag;
+  int bayesLineFlag;
+  int stochasticFlag;
+
+  int fullModelFlag;
+  int noiseModelFlag;
+  int glitchModelFlag;
+  int signalModelFlag;
+  int cleanModelFlag;
+  int cleanOnlyFlag;
+
+  int cleanFlag;
+  int glitchFlag;
+  int signalFlag;
+  int noiseFlag;
+
+  int fixIntrinsicFlag;
+  int fixExtrinsicFlag;
+  int geocenterPSDFlag;
+  int clusterPriorFlag;
+  int waveletPriorFlag;
+  int backgroundPriorFlag;
+  int amplitudePriorFlag;
+  int amplitudeProposalFlag;
+  int signalAmplitudePriorFlag;
+  int extrinsicAmplitudeFlag;
+  int orientationProposalFlag;
+  int clusterProposalFlag;
+  int logClusterProposalFlag;
+  int adaptTemperatureFlag;
+  int splineEvidenceFlag;
+  int constantLogLFlag;
+  int loudGlitchFlag;
+  int gnuplotFlag;
+  int verboseFlag;
+  int checkpointFlag;
+  int resumeFlag;
+  int searchFlag;
+  int rjFlag;
+
+  int nukeFlag; //tell BWP to remove chains/ and checkpoint/ directories
+
+  int srate;
+  double dt;
+  double df;
+  double fmin;
+  double fmax;
+  double Tobs;
+  double Twin;
+  double Qmin;
+  double Qmax;
+  double gmst;
+  double logLc;
+  double TwoDeltaToverN;
+
+  double TFtoProx;
+
+  double **r;
+  double **s;
+  double **rh;
+  double *h;
+  double *ORF;
+
+  double logZglitch;
+  double logZnoise;
+  double logZsignal;
+  double logZfull;
+
+  double varZglitch;
+  double varZnoise;
+  double varZsignal;
+  double varZfull;
+
+
+  void (*signal_amplitude_proposal)(double *, double *, double, long *, double, double **, double);
+  void (*glitch_amplitude_proposal)(double *, double *, double, long *, double, double **, double);
+
+  double (*signal_amplitude_prior)   (double *, double *, double, double);
+  double (*glitch_amplitude_prior)   (double *, double *, double, double);
+
+  double (*intrinsic_likelihood)(int, int, int, struct Model*, struct Model*, struct Prior*, struct Chain*, struct Data*);
+  double (*extrinsic_likelihood)(struct Network*, double *, double **, double *, double *, double **, struct Data*, double, double);
+
+  char runName[100];
+  char outputDirectory[1000];
+  char **channels;
+
+  ProcessParamsTable *commandLine;
+  LALDetector **detector;
+  LIGOTimeGPS epoch;
+};
+
+struct TimeFrequencyMap
+{
+  int N;
+  int nQ;
+  int nt;
+  int nf;
+  double ****pdf;
+  double ****snr;
+  double **max;
+
+};
+
+struct FisherMatrix
+{
+  int N;
+   double *sigma;
+   double *evalue;
+   double **evector;
+   double **matrix;
+   double *aamps;
+   double **count;
+};
+
+struct Background
+{
+  double logamp;
+  double index;
+  double fref;
+  double invfref;
+  double ***Cij;    //actually stores the inverse noise covariance matrix
+  double *detCij;   //frequency-dependent determenant of noise covariance matrix
+  double *spectrum;
+};
+

src/BayesWaveEvidence.h

+
+void average_log_likelihood_via_direct_downsampling(struct Chain *chain, double *ACL, int nPoints);
+void average_log_likelihood_via_recursion_downsampling(int M, struct Chain *chain);
+void average_log_likelihood_via_ACF(struct Chain *chain);
+
+void SplineThermodynamicIntegration(double *temperature, double *avgLogLikelihood, double *varLogLikelihood, int NC, double *logEvidence, double *varLogEvidence);
+void CovarianceThermodynamicIntegration(double *temperature, double **loglikelihood, int M, int NC, char modelname[], double *logEvidence, double *varLogEvidence);
+void ThermodynamicIntegration(double *temperature, double **loglikelihood, int M, int NC, char modelname[], double *logEvidence, double *varLogEvidence);
+void TrapezoidIntegration(struct Chain *chain, int counts, char modelname[], double *logEvidence, double *varLogEvidence);
+
+void LaplaceApproximation(struct Data *data, struct Model *model, struct Chain *chain, struct Prior *prior, double *logEvidence);

src/BayesWaveIO.h

+#ifdef __GNUC__
+#define UNUSED __attribute__ ((unused))
+#else
+#define UNUSED
+#endif
+
+
+/* ********************************************************************************** */
+/*                                                                                    */
+/*                        Data handling and injection routines                        */
+/*                                                                                    */
+/* ********************************************************************************** */
+
+REAL8TimeSeries *readTseries(CHAR *cachefile, CHAR *channel, LIGOTimeGPS start, REAL8 length);
+
+void InjectFromMDC(ProcessParamsTable *commandLine, LALInferenceIFOData *IFOdata, double *SNR);
+void BayesWaveInjection(ProcessParamsTable *commandLine, struct Data *data, struct Chain *chain, struct Prior *prior, double **psd, int *NC);
+void getChannels(ProcessParamsTable *commandLine, char **channel);
+int outputFrameFile(ProcessParamsTable *commandLine, struct Data *data, struct Model *model);
+
+/* ********************************************************************************** */
+/*                                                                                    */
+/*                                   Output Files                                     */
+/*                                                                                    */
+/* ********************************************************************************** */
+
+void write_gnuplot_script_header(FILE *script, char modelname[]);
+void write_gnuplot_script_frame (FILE *script, char modelname[], int signal, int glitch, int cycle, int NI);
+void write_bayesline_gnuplot_script_frame(FILE *script, char modelname[], int cycle);
+void write_bayeswave_gnuplot_script_frame(FILE *script, char modelname[], double Tobs, double fmin, double fmax, int phase, int signal, int glitch, int cycle, int NI);
+void write_evidence_gnuplot_script(FILE *script, char modelname[]);
+void print_bayesline_spectrum(char filename[], double *f, double *power, double *Sbase, double *Sline, int N);
+
+void print_version(FILE *fptr);
+void print_run_stats(FILE *fptr, struct Data *data, struct Chain *chain);
+void print_run_flags(FILE *fptr, struct Data *data, struct Prior *prior);
+void print_chain_files(struct Data *data, struct Chain *chain, struct Model **model, struct BayesLineParams ***bayesline, int ic);
+void print_chain_status(struct Data *data, struct Chain *chain, struct Model **model, int searchFlag);
+void flush_chain_files(struct Data *data, struct Chain *chain, int ic);
+
+void print_model(FILE *fptr, struct Data *data, struct Chain *chain, struct Model *model);
+void print_signal_model(FILE *fptr, struct Model *model);
+void print_glitch_model(FILE *fptr, struct Wavelet *glitch);
+
+void print_time_domain_waveforms(char filename[], double *h, int N, double *Snf, double eta, double Tobs, int imin, int imax, double tmin, double tmax);
+void print_time_domain_hdot(char filename[], double *h, int N, double *Snf, double eta, double Tobs, int imin, int imax, double tmin, double tmax);
+void print_frequency_domain_waveforms(char filename[], double *h, int N, double *Snf, double Tobs, double eta, int imin, int imax);
+void print_frequency_domain_data(char filename[], double *h, int N, double Tobs, int imin, int imax);
+void print_colored_time_domain_waveforms(char filename[], double *h, int N, double *Snf, double eta, double Tobs, int imin, int imax, double tmin, double tmax);
+
+void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *prior, ProcessParamsTable *commandLine);
+void parse_glitch_parameters(struct Data *data, struct Model *model, FILE **paramsfile, double **grec);
+void parse_signal_parameters(struct Data *data, struct Model *model, FILE *paramsfile, double **hrec, double **hrecPlus, double **hrecCross);
+void parse_bayesline_parameters(struct Data *data, struct Model *model, FILE **splinechain, FILE **linechain, double **psd);