IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 24979


Ignore:
Timestamp:
Aug 3, 2009, 4:07:27 AM (17 years ago)
Author:
Sebastian Jester
Message:

Comments (marked by XXXJester) and some pseudocode (in #IF 0 #ENDIF) to ignore NaNs and potentially do flag-checking in liststats()

Location:
branches/sj_branches/sj_SDSSaddstar_branch_r23928_20090420/Ohana/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • branches/sj_branches/sj_SDSSaddstar_branch_r23928_20090420/Ohana/src/photdbc/src/liststats.c

    r16040 r24979  
    11# include "photdbc.h"
    22
    3 enum {M_MEAN, M_MEDIAN, M_WT_MEAN, M_INNER_MEAN, 
     3enum {M_MEAN, M_MEDIAN, M_WT_MEAN, M_INNER_MEAN,
    44      M_INNER_WTMEAN, M_CHI_INNER_MEAN, M_CHI_INNER_WTMEAN};
    55
    … …  
    2424
    2525int liststats (double *value, double *dvalue, int N, StatType *stats) {
    26  
     26
     27  // XXXJester: Need to skip NaNs in all calculations here! See comment below.
     28
     29  // I'm pretty certain that we need to allow taking into account flag values in liststats, too,
     30  // e.g. excluding measurements that were saturated, CR etc. Otherwise, fatally flawed measurements in one
     31  // epoch will poison an average for an entire filter, and hence colours for the entire object.
    2732  int i, ks, ke, Nm;
    2833  double Mo, dMo, M, dM, X2, dS, *chi;
     34
     35# IF 0
     36  // My new code from #IF 0 ... #ENDIF below should go here
     37# ENDIF
    2938
    3039  stats[0].Nmeas = N;
    … …  
    3241  if (N < 2) return (FALSE);
    3342
     43  // XXXJester: This sort function is already going to be broken by any NaNs that might be in the input
     44  // vector, because NaN is not bigger nor smaller than anything. I.e. right after the variable declarations
     45  // and before setting stats[0].Nmeas, we need to purge the NaNs. Something like the following pseudocode,
     46  // which would need to go between the #IF 0 #ENDIF above.
     47
     48  #IF 0
     49  int N_orig, N_noNaN;
     50  double *value_noNaN, dvalue_noNaN;
     51  N_noNaN = purgeNaNs(N,value,value_noNaN,dvalue_noNaN);
     52  // purgeNaNs (e.g.) allocates memory for d?value_noNaN and returns all the non-NaN entries in them, and
     53  // the number of non-NaN values as return value (not sure what the IPP conventions are for where memory
     54  // should be allocated and whether I should return the pointer or the number of elements or false/true
     55  // for failure/success and all the other things by-reference).
     56
     57  // purgeNaNs looks like the following:
     58  int purgeNaNs(int N, double *value, double* value_noNaN, double *dvalue, double *dvalue_noNaN) {
     59    allocate(space for N entries in value_noNaN,dvalue_noNaN);
     60    int i_in, i_out, N_noNaN;
     61    for (i_in=0,i_out=0; i_in<N; i_in++) {
     62      if ( (!isNaN(value[i_in])) && (!isnan(dvalue[i_in]))) {
     63        value_noNaN[i_out] = value[i_in];
     64        dvalue_noNaN[i_out] = dvalue[i_in];
     65        i_out++ ;
     66      }
     67    }
     68    N_noNaN = i_out + 1;
     69    reallocate(space for N_noNaN entries in value_noNaN, dvalue_noNaN);
     70    return N_noNaN;
     71  }
     72  // However, it is debatable whether we should exclude values purely because their error is NaN. I could
     73  // imagine not excluding NaN errors, but then all operations on entries of dlist would have to check for
     74  // NaN individually. Most of the stuff that uses dvalue below only makes sense if dvalue[i] is guaranteed
     75  // not to be NaN, too (e.g. the chi-squared things), so my default is to purge NaN errors, too.
     76
     77  // purgeNaNs() at the same time could purge entries failing a predefined "fatal flag" setting - e.g. in
     78  // most contexts it doesn't make sense to average saturated fluxes.
     79
     80  // next commands would then be
     81  N = N_noNaN;
     82  free(value);
     83  free(dvalue);
     84  value = value_noNaN;
     85  dvalue = dvalue_noNaN;
     86  // to avoid remaining N, value, dvalue everywhere from here on down
     87# ENDIF
    3488  dsortpair (value, dvalue, N);
    3589  stats[0].median = value[(int)(0.5*N)];
    … …  
    61115    }
    62116    break;
    63   }   
     117  }
    64118
    65119  if ((statmode == M_CHI_INNER_MEAN) || (statmode == M_CHI_INNER_WTMEAN)) {
    … …  
    78132      M   += value[i] / SQ (dvalue[i]);
    79133      dM  += 1.0 / SQ (dvalue[i]);
    80       Nm  ++; 
    81     }   
     134      Nm  ++;
     135    }
    82136    Mo = M / dM;
    83137    dMo = sqrt (1.0 / dM);
    … …  
    86140      M   += value[i];
    87141      dM  += SQ (dvalue[i]);
    88       Nm  ++; 
    89     }   
     142      Nm  ++;
     143    }
    90144    Mo = M / (double) Nm;
    91145    dMo = sqrt (dM / (double) Nm);
    … …  
    109163  stats[0].sigma = dS;
    110164  stats[0].error = dMo;
    111 
     165  // XXXJester: if statmode == M_MEDIAN, dMo is *undefined* - neither set in case M_MEDIAN nor after label
     166  // chisq: - this is fixed in the relphot version of photdbc already. Not sure if this particular file is even used any more?
    112167  return (TRUE);
    113168}

  • branches/sj_branches/sj_SDSSaddstar_branch_r23928_20090420/Ohana/src/relphot/src/liststats.c

    r16040 r24979  
    11# include "relphot.h"
    22
    3 enum {M_MEAN, M_MEDIAN, M_WT_MEAN, M_INNER_MEAN, 
     3enum {M_MEAN, M_MEDIAN, M_WT_MEAN, M_INNER_MEAN,
    44      M_INNER_WTMEAN, M_CHI_INNER_MEAN, M_CHI_INNER_WTMEAN};
    55
    … …  
    2424
    2525int liststats (double *value, double *dvalue, int N, StatType *stats) {
    26  
     26
     27  // XXXJester: Need to skip NaNs in all calculations here! See comment below.
     28
     29  // I'm pretty certain that we need to allow taking into account flag values in liststats, too,
     30  // e.g. excluding measurements that were saturated, CR etc. Otherwise, fatally flawed measurements in one
     31  // epoch will poison an average for an entire filter, and hence colours for the entire object.
    2732  int i, ks, ke, Nm;
    2833  double Mo, dMo, M, dM, X2, dS, *chi;
    2934
     35# IF 0
     36  // My new code from #IF 0 ... #ENDIF below should go here
     37# ENDIF
     38
    3039  ke = ks = dMo = 0;
     40
    3141
    3242  stats[0].Nmeas = N;
    … …  
    3747  if (N < 1) return (FALSE);
    3848
     49  // XXXJester: This sort function is already going to be broken by any NaNs that might be in the input
     50  // vector, because NaN is not bigger nor smaller than anything. I.e. right after the variable declarations
     51  // and before setting stats[0].Nmeas, we need to purge the NaNs. Something like the following pseudocode,
     52  // which would need to go between the #IF 0 #ENDIF above.
     53
     54  #IF 0
     55  int N_orig, N_noNaN;
     56  double *value_noNaN, dvalue_noNaN;
     57  N_noNaN = purgeNaNs(N,value,value_noNaN,dvalue_noNaN);
     58  // purgeNaNs (e.g.) allocates memory for d?value_noNaN and returns all the non-NaN entries in them, and
     59  // the number of non-NaN values as return value (not sure what the IPP conventions are for where memory
     60  // should be allocated and whether I should return the pointer or the number of elements or false/true
     61  // for failure/success and all the other things by-reference).
     62
     63  // purgeNaNs looks like the following:
     64  int purgeNaNs(int N, double *value, double* value_noNaN, double *dvalue, double *dvalue_noNaN) {
     65    allocate(space for N entries in value_noNaN,dvalue_noNaN);
     66    int i_in, i_out, N_noNaN;
     67    for (i_in=0,i_out=0; i_in<N; i_in++) {
     68      if ( (!isNaN(value[i_in])) && (!isnan(dvalue[i_in]))) {
     69        value_noNaN[i_out] = value[i_in];
     70        dvalue_noNaN[i_out] = dvalue[i_in];
     71        i_out++ ;
     72      }
     73    }
     74    N_noNaN = i_out + 1;
     75    reallocate(space for N_noNaN entries in value_noNaN, dvalue_noNaN);
     76    return N_noNaN;
     77  }
     78  // However, it is debatable whether we should exclude values purely because their error is NaN. I could
     79  // imagine not excluding NaN errors, but then all operations on entries of dlist would have to check for
     80  // NaN individually. Most of the stuff that uses dvalue below only makes sense if dvalue[i] is guaranteed
     81  // not to be NaN, too (e.g. the chi-squared things), so my default is to purge NaN errors, too.
     82
     83  // purgeNaNs() at the same time could purge entries failing a predefined "fatal flag" setting - e.g. in
     84  // most contexts it doesn't make sense to average saturated fluxes.
     85
     86  // next commands would then be
     87  N = N_noNaN;
     88  free(value);
     89  free(dvalue);
     90  value = value_noNaN;
     91  dvalue = dvalue_noNaN;
     92  // to avoid remaining N, value, dvalue everywhere from here on down
     93# ENDIF
    3994  dsortpair (value, dvalue, N);
    4095  stats[0].median = value[(int)(0.5*N)];
    … …  
    66121    }
    67122    break;
    68   }   
     123  }
    69124
    70125  if ((statmode == M_CHI_INNER_MEAN) || (statmode == M_CHI_INNER_WTMEAN)) {
    … …  
    83138      M   += value[i] / SQ (dvalue[i]);
    84139      dM  += 1.0 / SQ (dvalue[i]);
    85       Nm  ++; 
    86     }   
     140      Nm  ++;
     141    }
    87142    Mo = M / dM;
    88143    dMo = sqrt (1.0 / dM);
    … …  
    91146      M   += value[i];
    92147      dM  += SQ (dvalue[i]);
    93       Nm  ++; 
    94     }   
     148      Nm  ++;
     149    }
    95150    Mo = M / (double) Nm;
    96151    dMo = sqrt (dM / (double) Nm);
    … …  
    114169  stats[0].sigma = dS;
    115170  stats[0].error = dMo;
    116 
    117171  return (TRUE);
    118172}
Note: See TracChangeset for help on using the changeset viewer.