Adding --fairdraw-index as a command line option to BayesWaveCleanFrame

src/BayesWaveCleanFrame.c

@@ -79,6 +79,8 @@ struct CleanFrameData
   double fr_seglen;   //frame segment length
   double fr_start;    //frame start time
 
+  int fairdraw_index;    //index of fair draw (note, this is 0-indexed.
+                         // The point must be from the second half of the chain file)
   char bw_model[MAXSTRINGSIZE]; //bayeswave glitch model
 
   double bw_seglen;   //bayeswave segment length
@@ -261,7 +263,20 @@ int main(int argc, char *argv[]) {
 
   char filename[128];
   FILE *tempfile;
-  
+
+  // when running with fairdraw, we need to pick a sample at random
+  //pick a sample at "random" e.g. final line
+  int fair_draw = N_burnin - 1;
+  if (data->fairdraw_index >= 0)
+  {
+    if ((data->fairdraw_index - N_burnin) < 0)
+    {
+      fprintf(stderr, "ERROR: fairdraw index %i must be in second half of samples (there are %i samples)\n", data->fairdraw_index, N_chain_samples);
+      exit(1);
+    }
+    fair_draw = data->fairdraw_index - N_burnin;
+  }
+
   //parse posterior samples and generate waveforms
   printf("Generating glitch model posterior:\n");
   for(int n=0; n<N_burnin; n++)
@@ -269,19 +284,26 @@ int main(int argc, char *argv[]) {
     for(i=0; i<NBW; i++) glitch_model[i] = 0.0;
     
     fscanf(paramsfile, "%i", &glitch_model_size);
+
+    if (n == fair_draw && !(data->medianFlag)){
+      if (glitch_model_size == 0){
+        fprintf(stderr, "WARNING! Fair draw requested at index %i, but the glitch is 0 dimensional there\n", fair_draw + N_burnin);
+      }
+    }
     
     if(glitch_model_size>0)
     {
       for(int d=1; d<=glitch_model_size; d++)
       {
         for(i=0; i<NP; i++) fscanf(paramsfile,"%lg", &glitch_model_params[i]);
-        
-        SineGaussianTime(glitch_model, glitch_model_params, NBW, 1, data->bw_seglen);
+        if ((data->medianFlag) || (n == fair_draw))
+        {
+          SineGaussianTime(glitch_model, glitch_model_params, NBW, 1, data->bw_seglen);
+        }
       }
-      
+
       //copy glitch_model to master glitch array for safe keeping
       for(i=0; i<NBW; i++) glitch_model_posterior[i][n] = glitch_model[i];
-
     }
     
     printProgress((double)(n+1)/(double)N_burnin);
@@ -300,8 +322,6 @@ int main(int argc, char *argv[]) {
   }
   else
   {
-    //pick a sample at random
-    int fair_draw = N_burnin-1;//0;//N_burnin/2;
     for(i=0; i<NBW; i++)
       glitch_model[i] = glitch_model_posterior[i][fair_draw];
   }
@@ -449,6 +469,7 @@ void parse_command_line_args(int argc, char **argv, struct CleanFrameData *data)
 
   data->medianFlag  = 0;
   data->verboseFlag = 0;
+  data->fairdraw_index = -1; // default to -1 meaning we take the last sample
 
   static struct option long_options[] =
   {
@@ -464,6 +485,7 @@ void parse_command_line_args(int argc, char **argv, struct CleanFrameData *data)
     {"trigtime",      required_argument, 0, 0},
     {"clean-suffix",  required_argument, 0, 0},
     {"outdir",        required_argument, 0, 0},
+    {"fairdraw-index",required_argument, 0, 0},
     {"median",        no_argument, 0, 0},
     {"verbose",       no_argument, 0, 0},
     {"help",          no_argument, 0,'h'},
@@ -565,6 +587,15 @@ void parse_command_line_args(int argc, char **argv, struct CleanFrameData *data)
         {
           data->medianFlag = 1;
         }
+        if(strcmp("fairdraw-index",long_options[long_index].name) == 0)
+        {
+          data->fairdraw_index = atoi(optarg);
+          if (data->fairdraw_index < 0)
+          {
+            fprintf(stderr, "ERROR: When supplied, fairdraw-index must be a positive integer, given instead %i\n", data->fairdraw_index);
+            exit(1);
+          }
+        }
         if(strcmp("help", long_options[long_index].name) == 0)
         {
           fprintf(stdout,"BayesWaveCleanFrame:\n");
@@ -637,6 +668,7 @@ void print_usage()
   fprintf(stdout,"--outdir        Output directory. Default is current directory ('.')\n");
   fprintf(stdout,"--verbose       print intermediate data products to disk\n");
   fprintf(stdout,"--median        use median glitch reconstruction\n");
+  fprintf(stdout,"--fairdraw-index Index of fairdraw in glitch-model file. Note that the index must be in the second half of the file\n");
   fprintf(stdout,"                can only be used on a completed run\n");
   fprintf(stdout,"--help          print help output\n");
   fprintf(stdout,"\n");