IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 23, 2004, 4:46:45 PM (22 years ago)
Author:
gusciora
Message:

...

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/math/psStats.c

    r1292 r1293  
    131131mean of the input vector.
    132132Inputs
     133    myVector
     134    maskVector
     135    maskVal
     136    stats
    133137Returns
    134138    NULL
    … …  
    404408    float rangeMin = 0.0;                       // Exclude data below this
    405409    float rangeMax = 0.0;                       // Exclude date above this
    406     psStats *stats2 = NULL;                     // Temporary stats structure
    407410
    408411
    409412    // Determine if the number of data points exceed a threshold which will
    410413    // cause to generate robust stats, as opposed to exact stats.
    411 
    412     if (myVector->n > stats->sampleLimit) {
    413         psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
    414 
    415         // Calculate the robust quartiles.
    416         stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
    417         p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
    418 
    419         // Store the robust quartiles into the sample quartile members.
    420         stats->sampleMedian = stats2->robustMedian;
    421 
    422         // Free temporary data buffers.
    423         psFree(stats2);
    424 
    425         // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
    426         stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
    427 
    428         return;
    429     }
     414    // XXX: This code is no longer used.  We now calculate the median exactly
     415    // regardless of the vector size.
     416    /*
     417        if (myVector->n > stats->sampleLimit) {
     418            psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
     419     
     420            // Calculate the robust quartiles.
     421            stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
     422            p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
     423     
     424            // Store the robust quartiles into the sample quartile members.
     425            stats->sampleMedian = stats2->robustMedian;
     426     
     427            // Free temporary data buffers.
     428            psFree(stats2);
     429     
     430            // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
     431            stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
     432     
     433            return;
     434        }
     435    */
    430436
    431437    // Determine how many data points fit inside this min/max range
    … …  
    483489
    484490    // Calculate the median exactly.
     491    // XXX: Is this the correct action?
    485492    if (0 == (nValues % 2)) {
    486493        stats->sampleMedian = 0.5 * (sortedVector->data.F32[(nValues/2)-1] +
    … …  
    498505    This routine smoothes the data in the input robustHistogram with a
    499506    Gaussian of width sigma.
     507 
     508    XXX: Must consult with IfA on the proper values for the Gaussian
     509    smoothing.  Currently, the Gaussian coefficients are such that only
     510    the center coefficient is non-zero.  This is because "e" is being
     511    raised to a very large negative number for all points except the
     512    center point.
     513 
     514    XXX: use a static variable for gaussianCoefs[] and compute them once.
    500515 *****************************************************************************/
    501516psVector *p_psVectorsmoothHistGaussian(psHistogram *robustHistogram,
    … …  
    541556    }
    542557
     558    // Perform the actual smoothing.
    543559    for(i=0;i<robustHistogram->nums->n;i++) {
    544560        smooth->data.F32[i] = 0.0;
    … …  
    576592    float rangeMin = 0.0;                       // Exclude data below this
    577593    float rangeMax = 0.0;                       // Exclude date above this
    578 
    579     // Determine if the number of data points exceed a threshold which will
    580     // cause to generate robust stats, as opposed to exact stats.
    581     if (myVector->n > stats->sampleLimit) {
    582         psAbort(__func__, "Robust Statistic Algorithms have not yet been defined or implemented.");
    583         psStats *stats2 = NULL;
    584         // Calculate the robust quartiles.
    585         stats2 = psStatsAlloc(PS_STAT_ROBUST_QUARTILE);
    586         p_psVectorRobustStats(myVector, maskVector, maskVal, stats2);
    587 
    588         // Store the robust quartiles into the sample quartile members.
    589         stats->sampleUQ = stats2->robustUQ;
    590         stats->sampleLQ = stats2->robustLQ;
    591 
    592         // Free temporary data buffers.
    593         psFree(stats2);
    594 
    595         // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure.
    596         stats->options = stats->options | PS_STAT_ROBUST_FOR_SAMPLE;
    597 
    598         return;
    599     }
    600594
    601595    // Determine how many data points fit inside this min/max range
    … …  
    861855
    862856/*****************************************************************************
     857p_psNormalizeVector(myData): this is a private function which normalizes the
     858elements of a vector to a range between 0.0 and 1.0.
    863859 *****************************************************************************/
    864860void p_psNormalizeVector(psVector *myData)
    … …  
    885881
    886882
     883/*****************************************************************************
     884p_psGaussian(myData, myParams): an internal function, used by robust stats,
     885intended to compute the Gaussian with a specified sigma/mean, at the
     886specified data point.
     887 *****************************************************************************/
    887888float p_psGaussian(const psVector *restrict myData,
    888889                   const psVector *restrict myParams)
    … …  
    894895    tmp/= ((float)sqrt(2.0 * M_PI * (stdev * stdev)));
    895896
    896     printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp);
     897    //    printf("p_psGaussian((%.2f), %.2f, %.2f) is %.2f\n", x, mean, stdev, tmp);
    897898    return(tmp);
    898899}
    899900
     901/*****************************************************************************
     902p_psGaussianDeriv(myData, myParams, whichParam): an internal function, which
     903calculates the specified partial derivative of the above Gaussian function.
     904 *****************************************************************************/
    900905float p_psGaussianDeriv(const psVector *restrict myData,
    901906                        const psVector *restrict myParams,
    … …  
    922927
    923928
     929/*****************************************************************************
     930p_psQuadratic(myData, myParams): an internal function, used by robust stats,
     931intended to compute a quadratic, with the specified parameters, at the
     932specified data point.
     933 *****************************************************************************/
    924934float p_psQuadratic(const psVector *restrict myParams,
    925935                    const psVector *restrict myCoords)
    … …  
    935945}
    936946
     947/*****************************************************************************
     948p_psQuadraticDeriv(myData, myParams, whichParam): an internal function, which
     949calculates the specified partial derivative of the above quadratic function.
     950 *****************************************************************************/
    937951float p_psQuadraticDeriv(const psVector *restrict myParams,
    938952                         const psVector *restrict myCoords,
    … …  
    959973[rangeLow, rangeHigh].  It determines the x-value of that polynomial such
    960974that f(x) == midpoint.  This functions uses a binary-search algorithm on the
    961 range and assumes that the polnomial is monotonically increasing within that
    962 range.
     975range and assumes that the polynomial is monotonically increasing or
     976decreasing within that range.
    963977 *****************************************************************************/
    964978float p_ps1DPolyMedian(psPolynomial1D *myPoly,
    … …  
    9981012
    9991013/******************************************************************************
    1000 p_ps1DPolyMedian(myPoly, rangeLow, rangeHigh, midpoint): This routine takes
    1001 as input a 1-D polynomial of arbitrary order (though we are using 2nd-order
    1002 polynomials here) and a range of x-values for which it is defined:
    1003 [rangeLow, rangeHigh].  It determines the x-value of that polynomial such
    1004 that f(x) == midpoint.  This functions uses a binary-search algorithm on the
    1005 range and assumes that the polnomial is monotonically increasing within that
    1006 range.
    1007  *****************************************************************************/
    1008 float p_psFitQuadratic(psHistogram *robustHistogram,
     1014p_psFitQuadratic(robustHistogram. cumulativeSums, binNum, fitFloat): given
     1015the specified histogram, with the specified pre-calculated cumulativeSums,
     1016this routine fits a quadratic f(x) to the 3 bins surrounding bin "binNum",
     1017and then finds the point x such that f(x) == fitFloat.
     1018 
     1019XXX: This function is currently not being used.
     1020 *****************************************************************************/
     1021float p_psFitQuadratic(psHistogram *histogram,
     1022                       psVector *cumulativeSums,
    10091023                       int binNum,
    10101024                       float fitFloat)
    10111025{
    1012     /*
    1013         if ((binNum > 0) &&
    1014             (binNum < (robustHistogram->nums->n+1))) {
    1015             x->data.F64[0] = (double) 0.5 *
    1016                       (robustHistogram->bounds->data.F32[binNum-1] +
    1017                        robustHistogram->bounds->data.F32[binNum]);
    1018             x->data.F64[1] = (double) 0.5 *
    1019                       (robustHistogram->bounds->data.F32[binNum] +
    1020                        robustHistogram->bounds->data.F32[binNum+1]);
    1021             x->data.F64[2] = (double) 0.5 *
    1022                       (robustHistogram->bounds->data.F32[binNum+1] +
    1023                        robustHistogram->bounds->data.F32[binNum+2]);
    1024      
    1025             y->data.F64[0] = cumulativeRobustSumsDl->data.F32[binNum-1];
    1026             y->data.F64[1] = cumulativeRobustSumsDl->data.F32[binNum];
    1027             y->data.F64[2] = cumulativeRobustSumsDl->data.F32[binNum+1];
    1028      
    1029             if (!((y->data.F64[0] <= fitFloat) &&
    1030                  (fitFloat <= y->data.F64[2]))) {
    1031                 psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n");
    1032             }
    1033      
    1034             yErr->data.F64[0] = 1.0;
    1035             yErr->data.F64[1] = 1.0;
    1036             yErr->data.F64[2] = 1.0;
    1037      
    1038             myPoly = psGetArrayPolynomial(myPoly, x, y, yErr);
    1039             return(p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], fitFloat);
    1040          } else {
    1041             return(0.5 * (robustHistogram->bounds->data.F32[binNum+1] +
    1042                           robustHistogram->bounds->data.F32[binNum]));
    1043         }
    1044     */
     1026    psVector *x = psVectorAlloc(3, PS_TYPE_F64);
     1027    psVector *y = psVectorAlloc(3, PS_TYPE_F64);
     1028    psVector *yErr = psVectorAlloc(3, PS_TYPE_F64);
     1029    psPolynomial1D *myPoly = psPolynomial1DAlloc(2);
     1030
     1031    if ((binNum > 0) &&
     1032            (binNum < (histogram->nums->n+1))) {
     1033        x->data.F64[0] = (double) 0.5 *
     1034                         (histogram->bounds->data.F32[binNum-1] +
     1035                          histogram->bounds->data.F32[binNum]);
     1036        x->data.F64[1] = (double) 0.5 *
     1037                         (histogram->bounds->data.F32[binNum] +
     1038                          histogram->bounds->data.F32[binNum+1]);
     1039        x->data.F64[2] = (double) 0.5 *
     1040                         (histogram->bounds->data.F32[binNum+1] +
     1041                          histogram->bounds->data.F32[binNum+2]);
     1042
     1043        y->data.F64[0] = cumulativeSums->data.F32[binNum-1];
     1044        y->data.F64[1] = cumulativeSums->data.F32[binNum];
     1045        y->data.F64[2] = cumulativeSums->data.F32[binNum+1];
     1046
     1047        if (!((y->data.F64[0] <= fitFloat) &&
     1048                (fitFloat <= y->data.F64[2]))) {
     1049            psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n");
     1050        }
     1051
     1052        yErr->data.F64[0] = 1.0;
     1053        yErr->data.F64[1] = 1.0;
     1054        yErr->data.F64[2] = 1.0;
     1055
     1056        myPoly = psGetArrayPolynomial(myPoly, x, y, yErr);
     1057        return(p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2],
     1058                                fitFloat));
     1059    } else {
     1060        return(0.5 * (histogram->bounds->data.F32[binNum+1] +
     1061                      histogram->bounds->data.F32[binNum]));
     1062    }
     1063
     1064    psFree(x);
     1065    psFree(y);
     1066    psFree(yErr);
     1067    psFree(myPoly);
    10451068    return(0.0);
    10461069}
    … …  
    10491072p_psVectorRobustStats(myVector, maskVector, maskVal, stats): this procedure
    10501073calculates a variety of robust stat measures:
    1051 PS_STAT_ROBUST_MEAN
    1052 PS_STAT_ROBUST_MEDIAN
    1053 PS_STAT_ROBUST_MODE
    1054 PS_STAT_ROBUST_STDEV
    1055 PS_STAT_ROBUST_QUARTILE
     1074    PS_STAT_ROBUST_MEAN
     1075    PS_STAT_ROBUST_MEDIAN
     1076    PS_STAT_ROBUST_MODE
     1077    PS_STAT_ROBUST_STDEV
     1078    PS_STAT_ROBUST_QUARTILE
    10561079I have included all that computation in a single function, as opposed to
    1057 breaking it across several functions for one primary reason:
    1058 they all
     1080breaking it across several functions for one primary reason:  they all
    10591081require the same basic initial processing steps (calculate the histogram,
    1060 etc.)
    1061 however there is no currently defined means, in the SDRS, to
    1062 specify this computation as a separate step.  If the above robust stat
    1063 measures were calcualted in separate functiosn, then much of the initial
    1064 processing would be duplicated.
     1082etc.)  however there is no currently defined means, in the SDRS, to specify
     1083this computation as a separate step.  If the above robust stat measures were
     1084calcualted in separate functiosn, then much of the initial processing would
     1085be duplicated.
     1086 
    10651087Inputs
    1066 myVector
    1067 maskVector
    1068 maskVal
    1069 stats
     1088    myVector
     1089    maskVector
     1090    maskVal
     1091    stats
    10701092Returns
    1071 NULL
     1093    NULL
    10721094*****************************************************************************/
    10731095void p_psVectorRobustStats(const psVector *restrict myVector,
    … …  
    12241246        sumN50+= (float) robustHistogram->nums->data.S32[i];
    12251247    }
     1248
    12261249    // XXX: is dL defined as the value at the LQ/UQ, or the bin number?
    12271250    dL = (UQBinNum - LQBinNum) / 4;
    … …  
    12441267                robustHistogramVector->data.F32[i];
    12451268            cumulativeMedian+= robustHistogramVector->data.F32[i];
    1246             sumNfit+= (float) robustHistogramVector->data.S32[i];
     1269            sumNfit+= (float) robustHistogram->nums->data.S32[i];
    12471270        }
    12481271    }
    … …  
    17801803    }
    17811804
     1805    // XXX: Should we abort if (stats->min == stats->max)?
     1806    if (stats->min >= stats->max) {
     1807        psAbort(__func__, "psVectorStats() called with range: %f to %f\n",
     1808                stats->min, stats->max);
     1809    }
     1810
     1811
     1812
    17821813    PS_CHECK_VECTOR_TYPE(in, PS_TYPE_F32);
    17831814    if (mask != NULL) {
Note: See TracChangeset for help on using the changeset viewer.