Index: /branches/sj_branches/sj_SDSSaddstar_branch_r23928_20090420/Ohana/src/photdbc/src/liststats.c
===================================================================
--- /branches/sj_branches/sj_SDSSaddstar_branch_r23928_20090420/Ohana/src/photdbc/src/liststats.c	(revision 24978)
+++ /branches/sj_branches/sj_SDSSaddstar_branch_r23928_20090420/Ohana/src/photdbc/src/liststats.c	(revision 24979)
@@ -1,5 +1,5 @@
 # include "photdbc.h"
 
-enum {M_MEAN, M_MEDIAN, M_WT_MEAN, M_INNER_MEAN, 
+enum {M_MEAN, M_MEDIAN, M_WT_MEAN, M_INNER_MEAN,
       M_INNER_WTMEAN, M_CHI_INNER_MEAN, M_CHI_INNER_WTMEAN};
 
@@ -24,7 +24,16 @@
 
 int liststats (double *value, double *dvalue, int N, StatType *stats) {
-  
+
+  // XXXJester: Need to skip NaNs in all calculations here! See comment below.
+
+  // I'm pretty certain that we need to allow taking into account flag values in liststats, too,
+  // e.g. excluding measurements that were saturated, CR etc. Otherwise, fatally flawed measurements in one
+  // epoch will poison an average for an entire filter, and hence colours for the entire object.
   int i, ks, ke, Nm;
   double Mo, dMo, M, dM, X2, dS, *chi;
+
+# IF 0
+  // My new code from #IF 0 ... #ENDIF below should go here
+# ENDIF
 
   stats[0].Nmeas = N;
@@ -32,4 +41,49 @@
   if (N < 2) return (FALSE);
 
+  // XXXJester: This sort function is already going to be broken by any NaNs that might be in the input
+  // vector, because NaN is not bigger nor smaller than anything. I.e. right after the variable declarations
+  // and before setting stats[0].Nmeas, we need to purge the NaNs. Something like the following pseudocode,
+  // which would need to go between the #IF 0 #ENDIF above.
+
+  #IF 0
+  int N_orig, N_noNaN;
+  double *value_noNaN, dvalue_noNaN;
+  N_noNaN = purgeNaNs(N,value,value_noNaN,dvalue_noNaN);
+  // purgeNaNs (e.g.) allocates memory for d?value_noNaN and returns all the non-NaN entries in them, and
+  // the number of non-NaN values as return value (not sure what the IPP conventions are for where memory
+  // should be allocated and whether I should return the pointer or the number of elements or false/true
+  // for failure/success and all the other things by-reference).
+
+  // purgeNaNs looks like the following:
+  int purgeNaNs(int N, double *value, double* value_noNaN, double *dvalue, double *dvalue_noNaN) {
+    allocate(space for N entries in value_noNaN,dvalue_noNaN);
+    int i_in, i_out, N_noNaN;
+    for (i_in=0,i_out=0; i_in<N; i_in++) {
+      if ( (!isNaN(value[i_in])) && (!isnan(dvalue[i_in]))) {
+        value_noNaN[i_out] = value[i_in];
+        dvalue_noNaN[i_out] = dvalue[i_in];
+        i_out++ ;
+      }
+    }
+    N_noNaN = i_out + 1;
+    reallocate(space for N_noNaN entries in value_noNaN, dvalue_noNaN);
+    return N_noNaN;
+  }
+  // However, it is debatable whether we should exclude values purely because their error is NaN. I could
+  // imagine not excluding NaN errors, but then all operations on entries of dlist would have to check for
+  // NaN individually. Most of the stuff that uses dvalue below only makes sense if dvalue[i] is guaranteed
+  // not to be NaN, too (e.g. the chi-squared things), so my default is to purge NaN errors, too.
+
+  // purgeNaNs() at the same time could purge entries failing a predefined "fatal flag" setting - e.g. in
+  // most contexts it doesn't make sense to average saturated fluxes.
+
+  // next commands would then be
+  N = N_noNaN;
+  free(value);
+  free(dvalue);
+  value = value_noNaN;
+  dvalue = dvalue_noNaN;
+  // to avoid remaining N, value, dvalue everywhere from here on down
+# ENDIF
   dsortpair (value, dvalue, N);
   stats[0].median = value[(int)(0.5*N)];
@@ -61,5 +115,5 @@
     }
     break;
-  }    
+  }
 
   if ((statmode == M_CHI_INNER_MEAN) || (statmode == M_CHI_INNER_WTMEAN)) {
@@ -78,6 +132,6 @@
       M   += value[i] / SQ (dvalue[i]);
       dM  += 1.0 / SQ (dvalue[i]);
-      Nm  ++;  
-    }	
+      Nm  ++;
+    }
     Mo = M / dM;
     dMo = sqrt (1.0 / dM);
@@ -86,6 +140,6 @@
       M   += value[i];
       dM  += SQ (dvalue[i]);
-      Nm  ++;  
-    }	
+      Nm  ++;
+    }
     Mo = M / (double) Nm;
     dMo = sqrt (dM / (double) Nm);
@@ -109,5 +163,6 @@
   stats[0].sigma = dS;
   stats[0].error = dMo;
-
+  // XXXJester: if statmode == M_MEDIAN, dMo is *undefined* - neither set in case M_MEDIAN nor after label
+  // chisq: - this is fixed in the relphot version of photdbc already. Not sure if this particular file is even used any more?
   return (TRUE);
 }
Index: /branches/sj_branches/sj_SDSSaddstar_branch_r23928_20090420/Ohana/src/relphot/src/liststats.c
===================================================================
--- /branches/sj_branches/sj_SDSSaddstar_branch_r23928_20090420/Ohana/src/relphot/src/liststats.c	(revision 24978)
+++ /branches/sj_branches/sj_SDSSaddstar_branch_r23928_20090420/Ohana/src/relphot/src/liststats.c	(revision 24979)
@@ -1,5 +1,5 @@
 # include "relphot.h"
 
-enum {M_MEAN, M_MEDIAN, M_WT_MEAN, M_INNER_MEAN, 
+enum {M_MEAN, M_MEDIAN, M_WT_MEAN, M_INNER_MEAN,
       M_INNER_WTMEAN, M_CHI_INNER_MEAN, M_CHI_INNER_WTMEAN};
 
@@ -24,9 +24,19 @@
 
 int liststats (double *value, double *dvalue, int N, StatType *stats) {
-  
+
+  // XXXJester: Need to skip NaNs in all calculations here! See comment below.
+
+  // I'm pretty certain that we need to allow taking into account flag values in liststats, too,
+  // e.g. excluding measurements that were saturated, CR etc. Otherwise, fatally flawed measurements in one
+  // epoch will poison an average for an entire filter, and hence colours for the entire object.
   int i, ks, ke, Nm;
   double Mo, dMo, M, dM, X2, dS, *chi;
 
+# IF 0
+  // My new code from #IF 0 ... #ENDIF below should go here
+# ENDIF
+
   ke = ks = dMo = 0;
+
 
   stats[0].Nmeas = N;
@@ -37,4 +47,49 @@
   if (N < 1) return (FALSE);
 
+  // XXXJester: This sort function is already going to be broken by any NaNs that might be in the input
+  // vector, because NaN is not bigger nor smaller than anything. I.e. right after the variable declarations
+  // and before setting stats[0].Nmeas, we need to purge the NaNs. Something like the following pseudocode,
+  // which would need to go between the #IF 0 #ENDIF above.
+
+  #IF 0
+  int N_orig, N_noNaN;
+  double *value_noNaN, dvalue_noNaN;
+  N_noNaN = purgeNaNs(N,value,value_noNaN,dvalue_noNaN);
+  // purgeNaNs (e.g.) allocates memory for d?value_noNaN and returns all the non-NaN entries in them, and
+  // the number of non-NaN values as return value (not sure what the IPP conventions are for where memory
+  // should be allocated and whether I should return the pointer or the number of elements or false/true
+  // for failure/success and all the other things by-reference).
+
+  // purgeNaNs looks like the following:
+  int purgeNaNs(int N, double *value, double* value_noNaN, double *dvalue, double *dvalue_noNaN) {
+    allocate(space for N entries in value_noNaN,dvalue_noNaN);
+    int i_in, i_out, N_noNaN;
+    for (i_in=0,i_out=0; i_in<N; i_in++) {
+      if ( (!isNaN(value[i_in])) && (!isnan(dvalue[i_in]))) {
+        value_noNaN[i_out] = value[i_in];
+        dvalue_noNaN[i_out] = dvalue[i_in];
+        i_out++ ;
+      }
+    }
+    N_noNaN = i_out + 1;
+    reallocate(space for N_noNaN entries in value_noNaN, dvalue_noNaN);
+    return N_noNaN;
+  }
+  // However, it is debatable whether we should exclude values purely because their error is NaN. I could
+  // imagine not excluding NaN errors, but then all operations on entries of dlist would have to check for
+  // NaN individually. Most of the stuff that uses dvalue below only makes sense if dvalue[i] is guaranteed
+  // not to be NaN, too (e.g. the chi-squared things), so my default is to purge NaN errors, too.
+
+  // purgeNaNs() at the same time could purge entries failing a predefined "fatal flag" setting - e.g. in
+  // most contexts it doesn't make sense to average saturated fluxes.
+
+  // next commands would then be
+  N = N_noNaN;
+  free(value);
+  free(dvalue);
+  value = value_noNaN;
+  dvalue = dvalue_noNaN;
+  // to avoid remaining N, value, dvalue everywhere from here on down
+# ENDIF
   dsortpair (value, dvalue, N);
   stats[0].median = value[(int)(0.5*N)];
@@ -66,5 +121,5 @@
     }
     break;
-  }    
+  }
 
   if ((statmode == M_CHI_INNER_MEAN) || (statmode == M_CHI_INNER_WTMEAN)) {
@@ -83,6 +138,6 @@
       M   += value[i] / SQ (dvalue[i]);
       dM  += 1.0 / SQ (dvalue[i]);
-      Nm  ++;  
-    }	
+      Nm  ++;
+    }
     Mo = M / dM;
     dMo = sqrt (1.0 / dM);
@@ -91,6 +146,6 @@
       M   += value[i];
       dM  += SQ (dvalue[i]);
-      Nm  ++;  
-    }	
+      Nm  ++;
+    }
     Mo = M / (double) Nm;
     dMo = sqrt (dM / (double) Nm);
@@ -114,5 +169,4 @@
   stats[0].sigma = dS;
   stats[0].error = dMo;
-
   return (TRUE);
 }