some changes from the polarization branch that are needed for relaxed ellipticity runs

src/BayesWave.c

@@ -1150,6 +1150,7 @@ void print_help_message(void)
   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,"  --noGlitch             skip glitch phase\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");

src/BayesWave.h

@@ -94,73 +94,77 @@
 /* ********************************************************************************** */
 struct Wavelet
 {
-   int size;
-   int smax;
-   int *index;
-   int dimension;
+   int size; // number of wavelets
+   int smax; // prior upper bound on number of wavelets
+   int *index; // keeps track of which of the smax possible wavelets are active
+   int dimension; // number of intrinsic parameters of each wavelet (5 of 6)
 
    //clustering prior
-   int Ncluster;
-   double cosy;
-   double logp;
+   int Ncluster; // DOESN'T WORK
+   double cosy; // DOESN'T WORK
+   double logp; // DOESN'T WORK
 
-   double *templates;
-   double **intParams;
+   double *templates; // fourier amplitudes evaluation of the wavelet (size N)
+   double **intParams; // parameters of each wavelet (smax times dimension)
 
-   char *name;
+   char *name; // not used, meant to be "wavelet" or "chirplet"
 };
 
 struct Model
 {
-  int size;
-  int Npol;
-
-  //double *Sn;
-  double **Snf;
-  double **invSnf;
-  double **SnS;
-  double *SnGeo;
-
-  double logL;
-  double logLnorm;
-  double *detLogL;
-  double *extParams;
-  double **h;
-  double **response;
+  int size; // total number of wavelets in all IFOs (glitch model) or polarization modes (signal model)
+  int N; // number of data points (srate times Tobs)
+  int Npol; // total number of polarization states we consider
+
+  double **Snf; // PSD
+  double **invSnf; // one over PSD at each frequency
+  double **SnS; // spline part of the PSD
+  double *SnGeo; // "geocenter PSD"
+
+  double logL; // log of the likelihood
+  double logLnorm; // log of the likelihood normalized (with the Sn part)
+  double *detLogL; // log of the likelihood in each detector
+  double *extParams; // extrinsic paramerers
+  double **h; // reference waveform at one interferometer (NpolS times N)
+  double **response; // response in each detector (NI times N)
     
-    int chirpletFlag;
-    int NW;
+  int chirpletFlag; // whether to use chirplets
+  int NW; // number of intrinsic wavelet parameters
 
-  struct Wavelet **signal;
-  struct Wavelet **glitch;
+  struct Wavelet **signal; // structure that holds the signal wavelets (see appropriate structure definition)
+  struct Wavelet **glitch; // structure that holds the glitch wavelets (see appropriate structure definition)
 
-  struct Network *projection;
+  struct Network *projection; // structure that holds the network projection (see appropriate structure definition)
 
-  struct FisherMatrix *fisher;
-  struct FisherMatrix *intfisher;
+  struct FisherMatrix *fisher; // structure that holds the fisher matrix (see appropriate structure definition)
+  struct FisherMatrix *intfisher; // structure that holds the fisher matrix of the intrinsic parameters (see appropriate structure definition)
 
-  struct Background *background;
+  struct Background *background; // structure that holds the stochastic background model (see appropriate structure definition)
 
-  void (*wavelet_bandwidth)(double *, double *, double *);
-  void (*wavelet)(double *, double *, int, int, double);
+  void (*wavelet_bandwidth)(double *, double *, double *); // function that points to the type of basis function used (wavelet or chirplet)
+  void (*wavelet)(double *, double *, int, int, double); // function that points to the type of basis function used (wavelet or chirplet)
 
 };
 
 struct Network
 {
-  double **expPhase;
-  double *deltaT;
-  double *Fplus;
-  double *Fcross;
-  double *dtimes;
+  double **expPhase; // dummy parameters for the recursion relations
+  double *deltaT; // dummy parameters for the recursion relations
+  double *Fplus; // antenna pattern functions
+  double *Fcross; // antenna pattern functions
+  double *Fv1; // antenna pattern functions
+  double *Fv2; // antenna pattern functions
+  double *Fb; // antenna pattern functions
+  double *Fl; // antenna pattern functions
+  double *dtimes; //time delay between reference location and IFOs
 };
 
 struct Chain
 {
-  int mc;
+  int mc;     // keeps track of total number of steps in the chain
   int NC;     // Number of parallel chains
   int burn;   // Number of burn in samples
-  int mod0;   // Counter keeping track of noise model
+  int mod0;   // Counter keeping track of noise model  REMOVE
   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
@@ -169,7 +173,7 @@ struct Chain
   int count;  // # of BurstRJMCMC iterations
   int zcount; // # of logL samples for thermodynamic integration
   int mcount; // # of MCMC attempts
-  int scount; // # of RJMCMC attempts
+  int scount; // # of RJMCMC attempts per polarization mode
   int xcount; // # of extrinsic attempts
   int fcount; // # of intrinsic fisher attempts
   int pcount; // # of intrinsic phase attempts
@@ -177,7 +181,7 @@ struct Chain
   int dcount; // # of density proposal attempts
   int ucount; // # of uniform proposal attempts
   int macc;   // # of MCMC successess
-  int sacc;   // # of RJMCMC successes
+  int sacc;   // # of RJMCMC successes per polarization mode
   int xacc;   // # of extrinsic successes
   int facc;   // # of intrinsic fisher successes
   int pacc;   // # of intrinsic phase successes
@@ -188,27 +192,27 @@ struct Chain
   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;
+  int *ptacc; // number of PTMCMC steps accepted
+  int *ptprop; // number of PTMCMC proposals
 
-  double *A;
+  double *A; // adapt temperature
 
-  gsl_rng *seed;
+  gsl_rng *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;
+  double *dT; // temperature difference
+  double beta; // one over the temperature
+  double tempMin; // minimum temperature
+  double tempStep; // initial chain spacing
+  double *temperature; // array that holds the chain temperatures
+  double *avgLogLikelihood; // average log L per chain
+  double *varLogLikelihood; // variance of logL per chain
 
-  int nPoints;// = 3*chain->count/4/chain->cycle;
-  double **logLchain;//=dmatrix(0,chain->NC-1,0,nPoints-1);
+  int nPoints; // how many likelihood samples to use to compute avgLogLikelihood and varLogLikelihood: 3*chain->count/4/chain->cycle
+  double **logLchain;// likelihood samples used for thermodynamic integration
 
 
   FILE *runFile;
@@ -222,22 +226,22 @@ struct Chain
 
 struct Prior
 {
-  int smin;
-  int smax;
-  int *gmin;
-  int *gmax;
-  int omax;
+  int smin; // minimum number of wavelets
+  int smax; // maximum number of wavelets
+  int *gmin; // minimum number of glitch wavelets (NOT USED)
+  int *gmax; // maximum number of glitch wavelets (NOT USED)
+  int omax; //REMOVE
 
-  double Snmin;
-  double Snmax;
+  double Snmin; // minimum spline amplitude
+  double Snmax; // maximum spline amplitude
 
-  double TFV;
-  double logTFV;
+  double TFV; // time-frequency volume
+  double logTFV; // log of time-frequency volume
 
-  double *bias;
-  double **range;
+  double *bias; // DOESN'T WORK
+  double **range; // prior range boundaries
 
-  //background prior
+  //exploratory background prior parameters, they are not currently used, should be revisited
     char bkg_name[1024];
     double bkg_df;
     double bkg_dQ;
@@ -246,19 +250,19 @@ struct Prior
     double **bkg_dist;
 
   //cluster prior parameters
-  double alpha;
-  double beta;
-  double gamma;
+  double alpha; // DOESN'T WORK
+  double beta; // DOESN'T WORK
+  double gamma; // DOESN'T WORK
 
   //amplitude prior parameters
-  double sSNRpeak;
+  double sSNRpeak; // peak of the amplitude prior per signal wavelet
 
   //glitch prior parameters
-  double gSNRpeak;
+  double gSNRpeak; // peak of the amplitude prior per glitch wavelet
 
   //dimension prior parameters
-  double Dtau;
-    double *Nwavelet;
+  double Dtau; // REMOVE and from the cmd line output
+  double *Nwavelet; //prior on the number of wavelets (smax)
 
   char path[100];
 
@@ -266,104 +270,106 @@ struct Prior
 
 struct Data
 {
-  int N;
-  int NI;
-  int Npol;
-  int tsize;
-  int fsize;
-  int imin;
-  int imax;
-  int Dmax;
-  int Dmin;
-
-  int runPhase;
-  int noiseSimFlag;
-  int bayesLineFlag;
-  int stochasticFlag;
-
-  int fullModelFlag;
+  int N; // number of data points
+  int NI; // number of IFOs
+  int Npol; // total number of polarization states we consider
+  int tsize; // time bins in the TF proposal
+  int fsize; // frequency bins in the TF proposal
+  int imin; // minimum frequency bin for inner product
+  int imax; // maximum frequency bin for inner product
+  int Dmax; // maximum number of wavelets
+  int Dmin; // minimum number of wavelets
+
+  int runPhase; // (0: cleaning, 1: not cleaning phase)
+  int noiseSimFlag; // REMOVE
+  int bayesLineFlag; // use BayesLine
+  int stochasticFlag; // whether you use the stochastic model
+
+  // type of run modes available (uses a combination of the models below)
+  int fullModelFlag; 
   int noiseModelFlag;
   int glitchModelFlag;
   int signalModelFlag;
   int cleanModelFlag;
   int cleanOnlyFlag;
 
+  // models available
   int cleanFlag;
   int glitchFlag;
   int signalFlag;
   int noiseFlag;
     
-  int skyFixFlag;
-
-  int polarizationFlag;
-  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 flatBayesLinePriorFlag;
-  int loudGlitchFlag;
-  int gnuplotFlag;
-  int verboseFlag;
-  int checkpointFlag;
-  int resumeFlag;
-  int searchFlag;
+  int skyFixFlag; //fix the sky location to the given values
+
+  int polarizationFlag; // use the elliptical polarization constraint
+  int fixIntrinsicFlag; // fix intrinsic parameters
+  int fixExtrinsicFlag; // fix extrinsic parameters
+  int geocenterPSDFlag; // use the geocenter PSD fo the amplitude prior
+  int clusterPriorFlag; // DOESN'T WORK
+  int waveletPriorFlag; // use the prior on the number of wavelets
+  int backgroundPriorFlag; // use the backgroung prior on the parameters of the glitch wavelets
+  int amplitudePriorFlag; // use the SNR-based prior on the wavelet amplitude
+  int amplitudeProposalFlag; // draw from the SNR-based prior on the wavelet amplitude as a proposal
+  int signalAmplitudePriorFlag; // use the SNR-based prior on the wavelet amplitude for the signal model only
+  int extrinsicAmplitudeFlag; // DOESN'T WORK, revisit
+  int orientationProposalFlag; // proposal that updates psi and ecc
+  int clusterProposalFlag; // propose wavelets close to existing wavelets (depends on SNR)
+  int logClusterProposalFlag; // propose wavelets close to existing wavelets (depends on log SNR)
+  int adaptTemperatureFlag; // adapt the PTMCMC temperature
+  int splineEvidenceFlag; // use splines and MCMC for the thermodynamic integration
+  int constantLogLFlag; // fix the likelihood to a constant
+  int flatBayesLinePriorFlag; // options for when running on loud glitches, need to be revisited
+  int loudGlitchFlag; // options for when running on loud glitches, need to be revisited
+  int gnuplotFlag; // print 200 samples (waveforms and other diagnostics)
+  int verboseFlag; // print every 100 iterations at stdout
+  int checkpointFlag; // use periodic checkpointing
+  int resumeFlag; // resume after checkpointing
+  int searchFlag; //depricated, maximized logL mode
   int waveletStartFlag;
-  int rjFlag;
+  int rjFlag; // whether you do RJ
 
   int nukeFlag; //tell BWP to remove chains/ and checkpoint/ directories
 
-  int srate;
+  int srate; // sampling rate
   
-  int chirpletFlag;
-  int NW;
+  int chirpletFlag; //use chirplets as a basis function 
+  int NW; // number of intrinsic wavelet parameters
     
-  double dt;
-  double df;
-  double fmin;
-  double fmax;
-  double Tobs;
-  double Twin;
-  double Qmin;
-  double Qmax;
-  double gmst;
-  double logLc;
+  double dt; // time resolution
+  double df; // frequency resolution
+  double fmin; // minimum frequency for inner products
+  double fmax; // maximum frequency for inner products
+  double Tobs; // segment length
+  double Twin; // window length
+  double Qmin; // minimum wavelet quality factor
+  double Qmax; // maximum wavelet quality factor
+  double gmst; // Greenwich time
+  double logLc; // constant factor that can be applied to the likelihood
   double seglen;
   double trigtime;
   double starttime;
-  double TwoDeltaToverN;
+  double TwoDeltaToverN; // what it says
   
-  double FIX_RA;
-  double FIX_DEC;
+  double FIX_RA; // fixed right ascension value
+  double FIX_DEC; //fixed declination value
 
-  double TFtoProx;
+  double TFtoProx; // ratio of TF to proximity proposals
 
-  double **r;
-  double **s;
-  double **rh;
-  double *h;
-  double *ORF;
+  double **r; // residual
+  double **s; // data h(t)
+  double **rh; // depricated
+  double *h; // model
+  double *ORF; // overlap reduction function
 
-  double logZglitch;
-  double logZnoise;
-  double logZsignal;
-  double logZfull;
+  double logZglitch; // glitch evidence
+  double logZnoise; // signal evidence
+  double logZsignal; // noise evidence
+  double logZfull; // full model evidence
 
-  double varZglitch;
-  double varZnoise;
-  double varZsignal;
-  double varZfull;
+  double varZglitch; // variation of the glitch evidence
+  double varZnoise; // variation of the signal evidence
+  double varZsignal; // variation of the noise evidence
+  double varZfull; // variation of the full model evidence
 
 
   void (*signal_amplitude_proposal)(double *, double *, gsl_rng *, double, double **, double);
@@ -379,6 +385,8 @@ struct Data
   char outputDirectory[1000];
   char **channels;
   char **ifos;
+  char **PolArray;
+  char PolModes[10];
 
   ProcessParamsTable *commandLine;
   LALDetector **detector;
@@ -387,35 +395,35 @@ struct Data
 
 struct TimeFrequencyMap
 {
-  int N;
-  int nQ;
-  int nt;
-  int nf;
-  double ****pdf;
-  double ****snr;
-  double **max;
+  int N; // number of samples
+  int nQ; // number of Q bins
+  int nt; // number of t bins
+  int nf; // number of f bins
+  double ****pdf; // TF density
+  double ****snr; // placeholder
+  double **max;// max density 
 
 };
 
 struct FisherMatrix
 {
-  int N;
-   double *sigma;
-   double *evalue;
-   double **evector;
-   double **matrix;
-   double *aamps;
-   double **count;
+  int N; // number of parameters
+   double *sigma; // diagonal elements
+   double *evalue; // eigenvalues
+   double **evector; // eigenvectors
+   double **matrix; // 2d NxN fisher matrix
+   double *aamps; // jump sizes
+   double **count; // deprecated
 };
 
 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;
+  double logamp; // amplitude of the background
+  double index; // spectral index
+  double fref; // reference frequency
+  double invfref; // one over the above
+  double ***Cij;    // actually stores the inverse noise covariance matrix
+  double *detCij;   // frequency-dependent determinant of noise covariance matrix
+  double *spectrum; // background spectrum per frequency
 };
 

src/BayesWaveIO.c

@@ -1428,6 +1428,9 @@ void parse_command_line(struct Data *data, struct Chain *chain, struct Prior *pr
   ppt = LALInferenceGetProcParamVal(commandLine, "--noSignal");
   if(ppt) data->signalModelFlag = 0;
 
+  ppt = LALInferenceGetProcParamVal(commandLine, "--noGlitch");
+  if(ppt) data->glitchModelFlag = 0;
+
   ppt = LALInferenceGetProcParamVal(commandLine, "--cleanOnly");
   if(ppt)
   {

src/BayesWaveWavelet.c

@@ -615,6 +615,8 @@ void computeProjectionCoeffs(struct Data *data, struct Network *projection, doub
     double toff = 0.0;
    double trigArgument=0.0;
 
+   if (!data->polarizationFlag) psi=0; //the polarization angle is only used when running with an elliptical polarization, otherwise its effect should be captured by the phase ratio of the plus and cross wavelets
+
    for(ifo = 0; ifo<NI; ifo++)
    {
       //Get time offsets for detectors