IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 37444


Ignore:
Timestamp:
Sep 29, 2014, 8:40:51 AM (12 years ago)
Author:
eugene
Message:

fix vstats to better handle combinations of weight and dvalue NULL or not

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/eam_branches/ipp-20140904/Ohana/src/libohana/src/vstats.c

    r37387 r37444  
    6060  int N20 = MAX (  0, 0.2*N);
    6161
    62   if (weight) {
     62  if (weight && dvalue) {
    6363    dsortthree (value, dvalue, weight, N);
    6464  } else {
    6565    if (dvalue) {
    6666      dsortpair (value, dvalue, N);
     67    } else if (weight) {
     68      dsortpair (value, weight, N);
    6769    } else {
    6870      dsort (value, N);
     
    112114  // I'm actually just using chisq to get the correct sorting order
    113115  if ((stats->statmode == VSTATS_CHI_INNER_MEAN) || (stats->statmode == VSTATS_CHI_INNER_WTMEAN)) {
     116    if (!dvalue) myAbort ("invalid combination: Chisq-based mean, no dvalue provided");
    114117    ALLOCATE (chi, double, N);
    115118    for (i = 0; i < N; i++) {
     
    128131  WeightedMean |= (stats->statmode == VSTATS_INNER_WTMEAN);
    129132  WeightedMean |= (stats->statmode == VSTATS_CHI_INNER_WTMEAN);
     133
     134  if (WeightedMean && !dvalue) myAbort ("invalid combination: weighted mean, no dvalue provided");
    130135
    131136  /* calculate stats based on the desired weighting scheme */
     
    149154        M  += value[i] * weight[i];
    150155        W  +=            weight[i];
    151         dM += SQ (weight[i] * dvalue[i]);
     156        dM += dvalue ? SQ (weight[i] * dvalue[i]) : 0.0;
    152157        R  += weight[i];
    153158        Nm += 1.0; 
    154159      }
    155160      Mo  = M / W;
    156       dMo = sqrt (dM) / R;
     161      dMo = dvalue ? sqrt (dM) / R : NAN;
    157162    }
    158163  } else {
     
    171176      for (i = ks; i < ke; i++) {
    172177        M   += value[i];
    173         dM  += SQ (dvalue[i]);
     178        dM  += dvalue ? SQ (dvalue[i]) : 0.0;
    174179        Nm  += 1.0; 
    175180      }
    176181      Mo = M / Nm;
    177       dMo = sqrt (dM) / Nm;
     182      dMo = dvalue ? sqrt (dM) / Nm : NAN;
    178183    }
    179184  }
     
    244249  int N20 = MAX (  0, 0.2*N);
    245250
    246   if (weight) {
     251  if (weight && dvalue) {
    247252    fsortthree (value, dvalue, weight, N);
    248253  } else {
    249254    if (dvalue) {
    250255      fsortpair (value, dvalue, N);
     256    } else if (weight) {
     257      fsortpair (value, weight, N);
    251258    } else {
    252259      fsort (value, N);
     
    296303  // I'm actually just using chisq to get the correct sorting order
    297304  if ((stats->statmode == VSTATS_CHI_INNER_MEAN) || (stats->statmode == VSTATS_CHI_INNER_WTMEAN)) {
     305    if (!dvalue) myAbort ("invalid combination: Chisq-based mean, no dvalue provided");
    298306    ALLOCATE (chi, float, N);
    299307    for (i = 0; i < N; i++) {
     
    312320  WeightedMean |= (stats->statmode == VSTATS_INNER_WTMEAN);
    313321  WeightedMean |= (stats->statmode == VSTATS_CHI_INNER_WTMEAN);
     322
     323  if (WeightedMean && !dvalue) myAbort ("invalid combination: weighted mean, no dvalue provided");
    314324
    315325  /* calculate stats based on the desired weighting scheme */
     
    333343        M  += value[i] * weight[i];
    334344        W  +=            weight[i];
    335         dM += SQ (weight[i] * dvalue[i]);
     345        dM += dvalue ? SQ (weight[i] * dvalue[i]) : 0.0;
    336346        R  += weight[i];
    337347        Nm += 1.0; 
    338348      }
    339349      Mo  = M / W;
    340       dMo = sqrt (dM) / R;
     350      dMo = dvalue ? sqrt (dM) / R : NAN;
    341351    }
    342352  } else {
     
    355365      for (i = ks; i < ke; i++) {
    356366        M   += value[i];
    357         dM  += SQ (dvalue[i]);
     367        dM  += dvalue ? SQ (dvalue[i]) : 0.0;
    358368        Nm  += 1.0; 
    359369      }
    360370      Mo = M / Nm;
    361       dMo = sqrt (dM) / Nm;
     371      dMo = dvalue ? sqrt (dM) / Nm : NAN;
    362372    }
    363373  }
     
    368378  for (i = ks; i < ke; i++) {
    369379    M  = SQ (value[i] - Mo);
    370     dM = SQ (dvalue[i]);
     380    dM = dvalue ? SQ (dvalue[i]) : 1.0;
    371381    X2 += M / dM;
    372382    dS += M;
    373383  }
    374   X2 = X2 / (Nm - 1);
     384  X2 = dvalue ? X2 / (Nm - 1) : NAN;
    375385  dS = sqrt (dS / (Nm - 1));
    376386
Note: See TracChangeset for help on using the changeset viewer.