Changeset 40408

Timestamp:

Apr 26, 2018, 12:37:03 PM (8 years ago)

Author:

eugene

Message:

clean up memory leaks for avperiodogram and related; finish avperiodomatch

Location:

trunk/Ohana/src/opihi

Files:

• dvo/Makefile (modified) (2 diffs)
• dvo/avperiodogram.c (modified) (11 diffs)
• dvo/avperiodomatch.c (added)
• dvo/dvo_host_utils.c (modified) (3 diffs)
• dvo/init.c (modified) (2 diffs)
• dvo/region_list.c (modified) (2 diffs)
• include/dvoshell.h (modified) (1 diff)

trunk/Ohana/src/opihi/dvo/Makefile

@@ -27 +27 @@
 $(SRC)/mySequence.$(ARCH).o             \
 $(SRC)/photcode_ops.$(ARCH).o           \
+$(SRC)/PeriodogramOps.$(ARCH).o         \
 $(SRC)/find_matches.$(ARCH).o           
 
@@ -46 +47 @@
 $(SRC)/avmatch.$(ARCH).o                \
 $(SRC)/avperiodogram.$(ARCH).o          \
+$(SRC)/avperiodomatch.$(ARCH).o         \
 $(SRC)/badimages.$(ARCH).o              \
 $(SRC)/catdir.$(ARCH).o                 \

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

@@ -1 +1 @@
 # include "dvoshell.h"
-# 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 {
-  int Nperiods;
-  opihi_flt *period;
-  opihi_flt *power;
-} PeriodogramResult;
-
-void PeriodogramResultFree (PeriodogramResult *result);
-void PeriodogramResultSave (PeriodogramResult *result, char *extname, FILE *foutput, Average *average);
-PeriodogramResult *periodogram_fm_raw (opihi_flt *time, opihi_flt *flux, opihi_flt *dflux, int Npts);
-
-float minP = 0.2;
-float maxP = 200.0;
 
 int avperiodogram (int argc, char **argv) {
@@ -57 +32 @@
     PARALLEL = TRUE;
   }
+
+  /* load photcode information (this needs to come before the SelectRegions command below to define the catdir */
+  if (!InitPhotcodes ()) goto escape;
+  int Nsecfilt = GetPhotcodeNsecfilt ();
 
   // parse skyregion options.  NOTE: this is stripped off in parallel operation and always
@@ -70 +49 @@
 
   // **** 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
+  // the result files are left behind for the user to access (no automatic re-merge)
   if (PARALLEL && !HOST_ID) {
-    int status = HostTableParallelOps (skylist, argc, argv, NULL, TRUE, 0, VERBOSE);
+    int status = HostTableParallelOps (skylist, argc, argv, NULL, FALSE, 0, VERBOSE);
     SkyListFree (skylist);
     FreeSkyRegionSelection (selection);
@@ -82 +61 @@
   if ((N = get_argument (argc, argv, "-min-period"))) {
     remove_argument (N, &argc, argv);
-    minP = atof (argv[N]);
+    PeriodogramSetOptions (atof (argv[N]), NAN);
     remove_argument (N, &argc, argv);
   }
@@ -88 +67 @@
   if ((N = get_argument (argc, argv, "-max-period"))) {
     remove_argument (N, &argc, argv);
-    maxP = atof (argv[N]);
+    PeriodogramSetOptions (NAN, atof (argv[N]));
     remove_argument (N, &argc, argv);
   }
@@ -179 +158 @@
   dvo_catalog_init (&catalog, TRUE);
 
-  /* load photcode information */
-  if (!InitPhotcodes ()) goto escape;
-  int Nsecfilt = GetPhotcodeNsecfilt ();
-
   // init locally static variables (time refs)
   dbExtractAveragesInit (); 
+  dbExtractMeasuresInit(HOST_ID);
 
   // examine line for 'where' or 'match to'.  'match to' is forbidden
@@ -207 +183 @@
 
   // parse stack elements into fields and scalars as needed
+  Nfields = 0;
+  ALLOCATE (fields, dbField, 1);
   if (!dbCheckStack (stack, Nstack, DVO_TABLE_AVERAGE, &fields, &Nfields)) goto escape;
 
@@ -288 +266 @@
         // The first 3 values[] / fields[] will be time, mag, dmag
         for (off_t n = 0; n < Nmfields; n++) {
-          mvalues[n] = dbExtractMeasures (average, secfilt, measure, NULL, NULL, &mfields[n]);
+          mvalues[n] = dbExtractMeasures (average, secfilt, &measure[k], NULL, NULL, &mfields[n]);
         }
         if (!dbBooleanCond (mstack, Nmstack, mvalues)) continue;
+
+        // do not allow any NAN or Inf values to polute the result
+        if (!isfinite(mvalues[0].Flt)) continue;
+        if (!isfinite(mvalues[1].Flt)) continue;
+        if (!isfinite(mvalues[2].Flt)) continue;
 
         time[Npts] = mvalues[0].Flt;
@@ -299 +282 @@
       }
 
-      PeriodogramResult *result = periodogram_fm_raw (time, mag, dmag, Npts);
+      PeriodogramResult *result = PeriodogramRawFloatingMean (time, mag, dmag, Npts);
+      free (time);
+      free (mag);
+      free (dmag);
+
+      if (!result) continue;
+
+      if (VERBOSE2) gprint (GP_ERR, "periodogram for %d, %d\n", (int) average->catID, (int) average->objID);
 
       // XXX if we have 1e6 objects, we will have output files with sizes of ~80GB
@@ -308 +298 @@
       Nobject ++;
     }
-
-    fclose (foutput);
-    fflush (foutput);
-
     dvo_catalog_free (&catalog);
   }
   ClearInterrupt (old_sigaction);
 
-  if (values) free (values);
+  fclose (foutput);
+  fflush (foutput);
+
+  // free measure stack stuff
+  FREE (mvalues);
+  dbFreeFields(mfields, Nmfields);
+  dbFreeStack(mstack, Nmstack);
+  FREE (mstack);
+
+  FREE (values);
   dbFreeFields (fields, Nfields);
   dbFreeStack (stack, Nstack);
-  if (stack) free (stack);
+  FREE (stack);
+
   SkyListFree (skylist);
   FreeSkyRegionSelection (selection);
+
+  free (output);
   return (TRUE);
 
@@ -341 +339 @@
   return (FALSE);
 }
-
-// low-level periodogram-fm function
-PeriodogramResult *periodogram_fm_raw (opihi_flt *time, opihi_flt *flux, opihi_flt *dflux, int Npts) {
-
-  /* find the max baseline, sum the inverse variances */
-  double minT = time[0];
-  double maxT = time[0];
-  opihi_flt Weight = 0;
-
-  opihi_flt *tv =  time;
-  opihi_flt *df = dflux;
-  opihi_flt *fv =  flux;
-
-  for (int i = 0; i < Npts; i++, tv++, df++) {
-    minT = MIN (minT, *tv);
-    maxT = MAX (maxT, *tv);
-    // skip points with 0.0 error? or add minimum variance?
-    Weight += 1.0 / (SQ(*df) + MIN_VAR); // MIN_VAR : added in quadrature to avoid Inf (XXX make this a user parameter?)
-  }
-
-  double WeightInv = 1.0 / Weight;
-  double dTime = maxT - minT;
-  if (dTime == 0) {
-    gprint (GP_ERR, "ERROR: time range is zero\n");
-    return NULL;
-  }
-
-  int Np = 0;
-  int NP = NPERIODS;
-  ALLOCATE_PTR (period, double, NP);
-  ALLOCATE_PTR (power, double, NP);
-
-  // for testing, we are going to write out the terms explicitly 
-  // for optimization, some of the trig functions can be calculated more efficiently 
-  // by storing cos,sin(w t) and using some recurrence rules
-
-  // maxP = minP * dP^n [n = 0 -- Nperiods-1]
-  // log(maxP) = log(minP) + (Nperiods - 1) * log (dP)
-  // log(dP) = (log(maxP) - log(minP)) / (Nperiods - 1)
-  // double dP = LINEAR ? (maxP - minP)/(NPERIODS - 1) : pow(10.0, ((log10(maxP) - log10(minP))/(NPERIODS - 1)));
-
-  double P = minP;
-  for (Np = 0; (Np < NPERIODS) && (P <= maxP); Np++) {
-    double w = 2*M_PI/P;
-    
-    /* find the period offset tau  */
-    tv = time;
-    df = dflux;
-
-    double cs = 0.0, sn = 0.0, sn2 = 0.0, cs2 = 0.0;
-
-    for (int i = 0; i < Npts; i++, tv++, df++) {
-      opihi_flt wt = WeightInv / (SQ(*df) + MIN_VAR);
-      cs  += wt*cos (*tv*w);
-      sn  += wt*sin (*tv*w);
-      cs2 += wt*cos (*tv*w*2);
-      sn2 += wt*sin (*tv*w*2);
-    }
-    double ytan = sn2 - 2*cs*sn;
-    double xtan = cs2 - (SQ(cs) - SQ(sn));
-    double tau = 0.5*atan2 (ytan, xtan) / w;
-      
-    /* find the power at this period */
-    tv = time;
-    df = dflux;
-    fv =  flux;
-
-    // YY = YY_s - Y_s*Y_s
-    // CC = CC_s - C_s*C_s
-    // SS = SS_s - S_s*S_s
-    // YC = YC_s - Y_s*C_s
-    // YS = YS_s - Y_s*S_s
-    
-    double YY_s = 0.0, CC_s = 0.0, SS_s = 0.0, YC_s = 0.0, YS_s = 0.0;
-    double Y_s = 0.0, C_s = 0.0, S_s = 0.0;
-    for (int i = 0; i < Npts; i++, tv++, fv++, df++) {
-      opihi_flt wt = WeightInv / (SQ(*df) + MIN_VAR);
-      cs = cos (w*(*tv-tau));
-      sn = sin (w*(*tv-tau));
-
-      double wtf = *fv*wt;
-      Y_s += wtf;
-      C_s += wt*cs;
-      S_s += wt*sn;
-
-      YY_s += wtf*(*fv);
-      YC_s += wtf*cs;
-      YS_s += wtf*sn;
-
-      CC_s += wt*cs*cs;
-      SS_s += wt*sn*sn;
-    }
-    double YY = YY_s - SQ(Y_s);
-    double YC = YC_s - Y_s*C_s;
-    double YS = YS_s - Y_s*S_s;
-    double CC = CC_s - SQ(C_s);
-    double SS = SS_s - SQ(S_s);
-
-    double Pa = SQ(YC) / CC;
-    double Pb = SQ(YS) / SS;
-    double Po = (Pa + Pb) / YY;
-
-    power[Np] = Po;
-    period[Np] = P;
-    if (Np >= NP) {
-      NP += 100;
-      REALLOCATE (power, opihi_flt, NP);
-      REALLOCATE (period, opihi_flt, NP);
-    }
-
-    float dP = 0.4 * SQ(P) / dTime;
-    P += dP;
-  }
-
-  ALLOCATE_PTR (result, PeriodogramResult, 1);
-  result->Nperiods = Np;
-  result->period = period;
-  result->power = power;
-  
-  return result;
-}
-
-void PeriodogramResultFree (PeriodogramResult *result) {
-  if (!result) return;
-
-  FREE (result->period);
-  FREE (result->power);
-  free (result);
-  return;
-}
-
-// write the period and power vectors to an extension in this output file
-void PeriodogramResultSave (PeriodogramResult *result, char *extname, FILE *foutput, Average *average) {
-
-  // generate a binary table 
-  Header outheader;
-  FTable outtable;
-  gfits_init_header (&outheader);
-  gfits_init_table  (&outtable);
-  outtable.header = &outheader;
-
-  gfits_create_table_header (&outheader, "BINTABLE", extname);
-
-  // add some metadata derived from average:
-  gfits_modify (&outheader, "RA_OBJ",  "%lf", 1, average->R);
-  gfits_modify (&outheader, "DEC_OBJ", "%lf", 1, average->D);
-  gfits_modify (&outheader, "OBJ_ID",  "%d", 1, average->objID);
-  gfits_modify (&outheader, "CAT_ID",  "%d", 1, average->catID);
-
-  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 (&outheader, &outtable);
-
-  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);
-}

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

@@ -212 +212 @@
         continue;
       }
-      free (table->hosts[i].results);
-      table->hosts[i].results = NULL;
+      // free (table->hosts[i].results);
+      // table->hosts[i].results = NULL;
       set_int_variable (name, 1); // result file has been read
 
@@ -225 +225 @@
         FreeVectorArray (invec, Ninvec);
       }
+    } else {
+      // free (table->hosts[i].results);
+      // table->hosts[i].results = NULL;
     }
   }
@@ -242 +245 @@
   free (vec);
 
-  free (table);
+  FreeHostTable (table);
   return TRUE;
 }

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

@@ -4 +4 @@
 int avmatch         PROTO((int, char **));
 int avperiodogram   PROTO((int, char **));
+int avperiodomatch  PROTO((int, char **));
 int badimages       PROTO((int, char **));
 int calextract      PROTO((int, char **));
@@ -68 +69 @@
   {1, "avextract",   avextract,    "extract average data values"},
   {1, "avmatch",     avmatch,      "extract average data values matched to RA,DEC points"},
-  {1, "avperiodogram", avperiodogram, "perform periodogram on objects based on restrictions"},
+  {1, "avperiodogram",  avperiodogram, "perform periodogram on objects based on restrictions"},
+  {1, "avperiodomatch", avperiodomatch, "perform periodogram on objects based on ra,dec list"},
   {1, "badimages",   badimages,    "look for images with anomalous astrometry"},
 //  {1, "calextract",  calextract,   "extract photometry calibration"},

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

@@ -181 +181 @@
   SkyList *skylist;
 
+  // the list of regions comes directly from a file
+  if (selection->list != NULL) {
+    skylist = SkyListLoadFile (selection->list);
+    return (skylist);
+  }
+
+  // all other options require sky to be set
+  if (!sky) {
+    gprint (GP_ERR, "CATDIR not set\n");
+    return NULL;
+  }
+
   /* determine region-file names */
   if (selection->name != NULL) {
@@ -186 +198 @@
     return (skylist);
   } 
-
-  if (selection->list != NULL) {
-    skylist = SkyListLoadFile (selection->list);
-    return (skylist);
-  }
 
   if (selection->useDisplay) {

trunk/Ohana/src/opihi/include/dvoshell.h

@@ -41 +41 @@
 char *HOSTDIR;
 char *RESULT_FILE;
+
+typedef struct {
+  int Nmeasure;
+  int Nperiods;
+  opihi_flt *period;
+  opihi_flt *power;
+  double P50;
+  double R50;
+  double R90;
+} PeriodogramResult;
+
+/*** Periodogram functions (avperiodogram and avperiodomatch) ***/
+void PeriodogramResultFree (PeriodogramResult *result);
+void PeriodogramResultSave (PeriodogramResult *result, char *extname, FILE *foutput, Average *average);
+void PeriodogramSetOptions (float minimumPeriod, float maximumPeriod);
+PeriodogramResult *PeriodogramRawFloatingMean (opihi_flt *time, opihi_flt *flux, opihi_flt *dflux, int Npts);
 
 /*** dvo prototypes ***/