Changeset 34576

Timestamp:

Oct 25, 2012, 10:20:24 AM (14 years ago)

Author:

eugene

Message:

fix problems with ChisqPos and object flags for relastro

Location:

trunk/Ohana/src/relastro

Files:

• include/relastro.h (modified) (1 diff)
• src/UpdateObjects.c (modified) (7 diffs)
• src/liststats.c (modified) (1 diff)
• src/relastro_objects.c (modified) (1 diff)

trunk/Ohana/src/relastro/include/relastro.h

@@ -271 +271 @@
 void          initstats           PROTO((char *mode));
 int           liststats           PROTO((double *value, double *dvalue, int N, StatType *stats));
+int           liststats_pos       PROTO((double *value, double *dvalue, int N, StatType *stats, int XVERB));
 Catalog      *load_catalogs       PROTO((SkyList *skylist, int *Ncatalog, int subselect, int hostID, char *hostpath));
 int           load_images         PROTO((FITS_DB *db, SkyList *skylist));

trunk/Ohana/src/relastro/src/UpdateObjects.c

@@ -103 +103 @@
       mode = FIT_MODE;
 
+      // XVERB |= (catalog[i].averge[j].objID == 0xc90) && (catalog[i].average[j].catID == 0x2a1e);
+      XVERB |= (catalog[i].average[j].objID == OBJ_ID_SRC) && (catalog[i].average[j].catID == CAT_ID_SRC);
+      XVERB |= (catalog[i].average[j].objID == OBJ_ID_DST) && (catalog[i].average[j].catID == CAT_ID_DST);
+
       // find the basic properties of the detections for this object (Tmin, Tmax, Tmean)
       for (k = 0; k < catalog[i].average[j].Nmeasure; k++) {
+
+        if (XVERB) {
+          char *date = ohana_sec_to_date (measure[k].t);
+          fprintf (stderr, OFF_T_FMT" %f %f %s : 0x%08x : 0x%08x\n",  k, measure[k].dR, measure[k].dD, date, measure[k].dbFlags, measureBig[k].dbFlags);
+          free (date);
+        }
 
         // does the measurement pass the supplied filtering constraints?
@@ -138 +148 @@
         if (isnan(dX[N])) {
           measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
+          if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
           continue;
         }
         if (isnan(dY[N])) {
           measure[k].dbFlags &= ~ID_MEAS_USED_OBJ;
+          if (measureBig) { measureBig[k].dbFlags &= ~ID_MEAS_USED_OBJ; }
           continue;
         }
@@ -207 +219 @@
       coords.crval1 = R[0];
       coords.crval2 = D[0];
-
-      // XVERB |= (catalog[i].averge[j].objID == 0xc90) && (catalog[i].average[j].catID == 0x2a1e);
-      XVERB |= (catalog[i].average[j].objID == OBJ_ID_SRC) && (catalog[i].average[j].catID == CAT_ID_SRC);
-      XVERB |= (catalog[i].average[j].objID == OBJ_ID_DST) && (catalog[i].average[j].catID == CAT_ID_DST);
 
       // to judge the quality of the PM and PAR fits, we need to fit all three models and compare Chisq
@@ -266 +274 @@
       // fit the average model
       if ((mode == FIT_AVERAGE) || (mode == FIT_PM_ONLY) || (mode == FIT_PM_AND_PAR)) {
-        liststats (R, dR, N, &statsR); // WARNING: this function modifies R (do not use after here)
-        liststats (D, dD, N, &statsD); // WARNING: this function modifies D (do not use after here)
+        liststats_pos (R, dR, N, &statsR, XVERB); // WARNING: this function modifies R (do not use after here)
+        liststats_pos (D, dD, N, &statsD, XVERB); // WARNING: this function modifies D (do not use after here)
 
         fitAve.Ro = statsR.mean;
@@ -275 +283 @@
         fitAve.dDo = 3600.0*statsD.sigma;
 
-        fitAve.chisq = 0.5*(statsR.chisq + statsD.chisq);
+        fitAve.chisq = 0.5 * (statsR.chisq + statsD.chisq);
         fitAve.Nfit = N;
 
@@ -358 +366 @@
       catalog[i].average[j].Trange = (Trange * 86400 * 365.26);
       catalog[i].average[j].Npos = fit.Nfit;
-      if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f)\n",
+      if (XVERB) fprintf (stderr, "%f %f -> %f %f (%f,%f) pm=(%f %f) chisq=(%f, %f, %f)\n",
                           catalog[i].average[j].R,
                           catalog[i].average[j].D,
@@ -365 +373 @@
                           3600*(catalog[i].average[j].D - fit.Do),
                           catalog[i].average[j].uR,
-                          catalog[i].average[j].uD);
-
+                          catalog[i].average[j].uD,
+                          fitAve.chisq, fitPM.chisq, fitPAR.chisq);
     }
 

trunk/Ohana/src/relastro/src/liststats.c

@@ -148 +148 @@
 }
 
+int liststats_pos (double *value, double *dvalue, int N, StatType *stats, int XVERB) {
+  
+  int i, ks, ke, Nm;
+  double Mo, dMo, M, dM, X2, dS, *chi;
+
+  ke = ks = dMo = 0;
+
+  stats[0].Nmeas = N;
+  stats[0].mean  = 0;
+  stats[0].sigma = 0;
+  stats[0].error = 0;
+  stats[0].chisq = 0;
+  if (N < 1) return (FALSE);
+
+  dsortpair (value, dvalue, N);
+  if (N % 2) {
+    stats[0].median = value[(int)(0.5*N)];
+  } else {
+    stats[0].median = 0.5*(value[N/2] + value[N/2 - 1]);
+  }
+  stats[0].min    = value[0];
+  stats[0].max    = value[N-1];
+
+  switch (statmode) {
+  case M_MEDIAN:
+    ks = 0;
+    ke = N;
+    Mo = stats[0].median;
+    Nm = N;
+    goto chisq;
+    break;
+  case M_MEAN:
+  case M_WTMEAN:
+    ks = 0;
+    ke = N;
+    break;
+  case M_INNER_MEAN:
+  case M_INNER_WTMEAN:
+  case M_CHI_INNER_MEAN:
+  case M_CHI_INNER_WTMEAN:
+    ks = 0.25*N + 0.50;
+    ke = 0.75*N + 0.25;
+    if (N <= 3) {
+      ks = 0;
+      ke = N;
+    }
+    break;
+  case M_INNER_80_MEAN:
+  case M_INNER_80_WTMEAN:
+  case M_CHI_INNER_80_MEAN:
+  case M_CHI_INNER_80_WTMEAN:
+    ks = 0.1*N + 0.1;
+    ke = 0.8*N + 0.1;
+    if (N <= 10) {
+      ks = 0;
+      ke = N;
+    }
+    break;
+  }    
+
+  if ((statmode == M_CHI_INNER_MEAN) || 
+      (statmode == M_CHI_INNER_WTMEAN) || 
+      (statmode == M_CHI_INNER_80_MEAN) || 
+      (statmode == M_CHI_INNER_80_WTMEAN)) {
+    ALLOCATE (chi, double, N);
+    for (i = 0; i < N; i++) {
+      chi[i] = (value[i] - stats[0].median) / dvalue[i];
+    }
+    dsortthree (chi, value, dvalue, N);
+    free (chi);
+  }
+
+  /* calculating the per-star offset based on the weighted average */
+  M = dM = Nm = 0;
+  if ((statmode == M_WTMEAN) || 
+      (statmode == M_INNER_WTMEAN) || 
+      (statmode == M_INNER_80_WTMEAN) || 
+      (statmode == M_CHI_INNER_WTMEAN) || 
+      (statmode == M_CHI_INNER_80_WTMEAN)) {
+    for (i = ks; i < ke; i++) {
+      M   += value[i] / SQ (dvalue[i]);
+      dM  += 1.0 / SQ (dvalue[i]);
+      Nm  ++;  
+    }   
+    Mo = M / dM;
+    dMo = sqrt (1.0 / dM);
+  } else {
+    for (i = ks; i < ke; i++) {
+      M   += value[i];
+      dM  += SQ (dvalue[i]);
+      Nm  ++;  
+    }   
+    Mo = M / (double) Nm;
+    dMo = sqrt (dM / (double) Nm);
+  }
+
+ chisq:
+  /* find sigma and chisq */
+  X2 = dS = 0;
+  for (i = ks; i < ke; i++) {
+    M  = SQ (value[i] - Mo);
+    dM = SQ (dvalue[i]);
+    X2 += M / dM;
+    dS += M;
+    if (XVERB) fprintf (stderr, "chisq contrib : %f : %f : %f : %f : %f\n", Mo, 3600.0*(value[i] - Mo), 3600.0*dvalue[i], (value[i] - Mo) / dvalue[i], X2);
+  }
+  X2 = X2 / (Nm - 1);
+  dS = sqrt (dS / Nm);
+
+  stats[0].mean  = Mo;
+  stats[0].Nmeas = Nm;
+  stats[0].chisq = X2;
+  stats[0].sigma = dS;
+  stats[0].error = dMo;
+
+  return (TRUE);
+}
+

trunk/Ohana/src/relastro/src/relastro_objects.c

@@ -132 +132 @@
 
     char command[1024];
-    snprintf (command, 1024, "relastro_client -update-objects -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f", 
-              table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax);
+    snprintf (command, 1024, "relastro_client -update-objects -hostID %d -D CATDIR %s -hostdir %s -region %f %f %f %f -statmode %s", 
+              table->hosts[i].hostID, CATDIR, table->hosts[i].pathname, UserPatch.Rmin, UserPatch.Rmax, UserPatch.Dmin, UserPatch.Dmax, STATMODE);
 
     char tmpline[1024];