Changeset 40400

Timestamp:

Apr 19, 2018, 4:00:59 PM (8 years ago)

Author:

eugene

Message:

updates to avperiodogram

File:

• trunk/Ohana/src/opihi/dvo/avperiodogram.c (modified) (11 diffs)

trunk/Ohana/src/opihi/dvo/avperiodogram.c

@@ -2 +2 @@
 # define MIN_VAR 1e-8
 # define NPERIODS 50
+
+/* basic outline of the code:
+
+ * parse optional arguments
+ * -parallel : launch parallel jobs [below is called for non-parallel or for dvo_client
+
+ * parse optional -where-measure clause
+ * 
+
+ */
 
 typedef struct {
@@ -17 +27 @@
 int avperiodogram (int argc, char **argv) {
   
+# if (0)
+
   int N, next, Nfields;
 
@@ -25 +37 @@
   SkyRegionSelection *selection = NULL;
 
+  // **** parse the optional arguments ****
   if ((N = get_argument (argc, argv, "-h"))) goto help;
   if ((N = get_argument (argc, argv, "--help"))) goto help;
@@ -45 +58 @@
     PARALLEL = TRUE;
   }
+
+  // parse skyregion options.  NOTE: this is stripped off in parallel operation and always
+  // defined for the client via the -skyregion option.  The dvo_client parses this
+  // argument in the main program, before it is passed to the command (like mextract)
+  if ((selection = SetRegionSelection (&argc, argv)) == NULL) {
+    print_error(); 
+    goto escape; 
+  }
+
+  /* load region corresponding to selection above */
+  if ((skylist = SelectRegions (selection)) == NULL) goto escape;
+
+  // **** launch parallel / remote jobs ****
+  // this does all the work of re-packaging the command, calling it on the remote machines, then loading in the results
+  if (PARALLEL && !HOST_ID) {
+    // XXX this is not going to work : the measure args will have been stripped off
+    int status = HostTableParallelOps (skylist, argc, argv, RESULT_FILE, TRUE, 0, VERBOSE);
+
+    dbFreeFields (fields, Nfields);
+    dbFreeStack (stack, Nstack);
+    free (stack);
+    FreeSkyRegionSelection (selection);
+    dvo_catalog_free (&catalog);
+
+    return status;
+  }
+
+  // NOTE: optional arguments below are parsed after the possible remote call so they are also passed to the dvo_client
 
   if ((N = get_argument (argc, argv, "-min-period"))) {
@@ -92 +133 @@
   ALLOCATE (mvalues, dbValue, Nmfields);
 
-  // XXX : I've removed -parallel-local : see avextract for an example
+  // use remote tables, but not dvo_client..
+  int PARALLEL_LOCAL = FALSE;
+  HostTable *table = NULL;
+  if ((N = get_argument (argc, argv, "-parallel-local"))) {
+    remove_argument (N, &argc, argv);
+    PARALLEL_LOCAL = TRUE;
+
+    char *CATDIR = GetCATDIR();
+    if (!CATDIR) {
+      gprint (GP_ERR, "CATDIR is not set\n");
+      return FALSE;
+    }
+    SkyTable *sky = GetSkyTable();
+    if (!sky) {
+      gprint (GP_ERR, "failed to load sky table for database\n");
+      return FALSE;
+    }
+    table = HostTableLoad (CATDIR, sky->hosts);
+    if (!table) {
+      gprint (GP_ERR, "ERROR: failure reading Host Table %s for database %s\n", sky->hosts, CATDIR);
+      return FALSE;
+    }    
+  }
 
   Catalog catalog;
@@ -101 +164 @@
   int Nsecfilt = GetPhotcodeNsecfilt ();
 
-  // parse skyregion options.  NOTE: this is stripped off in parallel operation and always
-  // defined for the client via the -skyregion option.  The dvo_client parses this
-  // argument in the main program, before it is passed to the command (like mextract)
-  if ((selection = SetRegionSelection (&argc, argv)) == NULL) {
-    print_error(); 
-    goto escape; 
-  }
-
   // init locally static variables (time refs)
   dbExtractAveragesInit (); 
 
-  // command-line is of the form: avperiodogram where (field op value)...
+  // command-line is of the form: avperiodogram (output) where (field op value)...
+  if (argc < 2) goto help;
+  char *output = strcreate (argv[1]);
 
   // examine line for 'where' or 'match to'.  'match to' is forbidden
-  int state = dbCmdlineConditions (argc, argv, 1, &next);
+  int state = dbCmdlineConditions (argc, argv, 2, &next);
   if (state == DVO_DB_CMDLINE_ERROR) goto escape;
   if (state == DVO_DB_CMDLINE_IS_MATCH) goto escape; // not allowed for avperiodogram
@@ -140 +197 @@
   ALLOCATE (values, dbValue, Nfields);
 
-  /* load region corresponding to selection above */
-  if ((skylist = SelectRegions (selection)) == NULL) goto escape;
-
-  // this does all the work of re-packaging the command, calling it on the remote machines, then loading in the results
-  if (PARALLEL && !HOST_ID) {
-    int status = HostTableParallelOps (skylist, argc, argv, RESULT_FILE, TRUE, 0, VERBOSE);
-
-    dbFreeFields (fields, Nfields);
-    dbFreeStack (stack, Nstack);
-    free (stack);
-    FreeSkyRegionSelection (selection);
-    dvo_catalog_free (&catalog);
-
-    return status;
-  }
-
   // grab data from all selected sky regions
   struct sigaction *old_sigaction = SetInterrupt();
@@ -180 +221 @@
     }
     dvo_catalog_unlock (&catalog);
+
+    // we save the results in a single FITS file, one extension per object
+    FILE *foutput = NULL;
+    if (RESULT_FILE) {
+      foutput = fopen (RESULT_FILE, "w");
+    } else {
+      foutput = fopen (output, "w");
+    }
+
+    // generate the PHU and write to disk
+    Header header;
+    Matrix matrix;
+    gfits_init_header (&header);
+    gfits_init_matrix (&matrix);
+    header.extend = TRUE;
+    gfits_create_header (&header);
+    gfits_create_matrix (&header, &matrix);
+    gfits_fwrite_header (foutput, &header);
+    gfits_fwrite_matrix (foutput, &matrix);
+    gfits_free_header (&header);
+    gfits_free_matrix (&matrix);
 
     // Scan the catalog for objects which match the WHERE clause
@@ -236 +298 @@
 
       PeriodogramResult *result = periodogram_fm_raw (time, mag, dmag, Npts);
-      // XXX write the periodogram result to a file
+
+      PeriodogramResultSave (result, extname, foutput);
       PeriodogramResultFree (result);
     }
+
+    fclose (foutput);
+    fflush (foutput);
+
     dvo_catalog_free (&catalog);
   }
@@ -261 +328 @@
   return (FALSE);
 
+
  help:
-  gprint (GP_ERR, "USAGE: avperiodogram where (expression) -where-measure (expression)\n");
+
+# endif
+  gprint (GP_ERR, "USAGE: avperiodogram (output) where (expression) -where-measure (expression)\n");
   gprint (GP_ERR, "  the where (expression) is a boolean to limit the objects analysed based on average properties\n");
   gprint (GP_ERR, "  NOTE: the optional -where-measure (expression) MUST come after the average where expression\n");
   return (FALSE);
+
+
 }
 
@@ -397 +469 @@
   return;
 }
+
+/*
+
+// write the period and power vectors to an extension in this output file
+void PeriodogramResultFree (PeriodogramResult *result, char *extname, FILE *foutput) {
+
+  // generate a binary table 
+  Header outheader;
+  FTable outtable;
+  gfits_init_header (&outheader);
+  gfits_init_table  (&outtable);
+  outtable.header = &outheader;
+
+  // XXX add some metadata to the header to identify the object
+  // extname should be related to the object IDs
+  gfits_create_table_header (&outheader, "BINTABLE", extname);
+
+  gfits_define_bintable_column (&outheader, "D", "PERIOD", NULL, NULL, 1.0, 0.0);
+  gfits_define_bintable_column (&outheader, "D", "POWER", NULL, NULL, 1.0, 0.0);
+
+  gfits_create_table (theader, ftable);
+
+  gfits_set_bintable_column (&outheader, &outtable, "PERIOD", result->period, result->Nperiods);
+  gfits_set_bintable_column (&outheader, &outtable, "POWER",  result->power,  result->Nperiods);
+
+  // write the actual table data
+  gfits_fwrite_Theader (foutput, &outheader);
+  gfits_fwrite_table   (foutput, &outtable);
+    
+  gfits_free_header (&outheader);
+  gfits_free_table  (&outtable);
+}
+
+*/