IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 3547


Ignore:
Timestamp:
Mar 29, 2005, 12:34:59 PM (21 years ago)
Author:
gusciora
Message:

Fixes bugs 338. Deeals with robust stats, and min->lastDelta.

Location:
trunk/psLib/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/dataManip/psMinimize.c

    r3540 r3547  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.108 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-03-29 19:41:56 $
     11 *  @version $Revision: 1.109 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-03-29 22:34:59 $
    1313 *
    1414 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    786786        //
    787787        if (currChi2 > newChi2) {
    788             min->lastDelta = currChi2 - newChi2;
     788            min->lastDelta = (currChi2 - newChi2)/currChi2;
    789789            min->value = newChi2;
    790790
    … …  
    14061406with an initial guess at the parameters of that function and vector "line"
    14071407of the same size as the parameter vector.  It will minimize the function
    1408 along that vector anr returns the offset along that vector at which the
     1408along that vector and returns the offset along that vector at which the
    14091409minimum is determined.
    14101410 
    … …  
    17131713            psFree(v);
    17141714            min->iter = iterationNumber;
     1715            // XXX: Ensure that currFuncVal is the correct value to use here.
     1716            min->value = currFuncVal;
    17151717            psTrace(".psLib.dataManip.psMinimizePowell", 4,
    17161718                    "---- psMinimizePowell() end (1)(true) ----\n");
    … …  
    17491751        if (fabs(baseFuncVal - currFuncVal) <= min->tol) {
    17501752            psFree(v);
     1753            // XXX: Ensure that currFuncVal is the correct value to use here.
     1754            min->value = currFuncVal;
    17511755            min->iter = iterationNumber;
    17521756            psTrace(".psLib.dataManip.psMinimizePowell", 4,

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

    r3540 r3547  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.122 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-03-29 19:41:56 $
     11 *  @version $Revision: 1.123 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-03-29 22:34:59 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    14891489    if (fabs(binSize) <= FLT_EPSILON) {
    14901490        if (stats->options & PS_STAT_ROBUST_MEAN) {
    1491             stats->robustMean = stats->clippedMean;
     1491            stats->robustMean = tmpStats->clippedMean;
    14921492        }
    14931493        if (stats->options & PS_STAT_ROBUST_MEDIAN) {
    1494             stats->robustMedian = stats->clippedMean;
     1494            stats->robustMedian = tmpStats->clippedMean;
    14951495        }
    14961496        if (stats->options & PS_STAT_ROBUST_MODE) {
    1497             stats->robustMode = stats->clippedMean;
     1497            stats->robustMode = tmpStats->clippedMean;
    14981498        }
    14991499        if (stats->options & PS_STAT_ROBUST_STDEV) {
    … …  
    15011501        }
    15021502        if (stats->options & PS_STAT_ROBUST_QUARTILE) {
    1503             stats->robustUQ = stats->clippedMean;
    1504             stats->robustLQ = stats->clippedMean;
     1503            stats->robustUQ = tmpStats->clippedMean;
     1504            stats->robustLQ = tmpStats->clippedMean;
    15051505        }
    15061506        // XXX: Set these to the number of unmasked data points?
    … …  
    15121512
    15131513    // Determine minimum and maximum values in the data vector.
    1514     if (isnan(stats->min)) {
    1515         if (0 != p_psVectorMin(myVector, maskVector, maskVal, stats)) {
     1514    if (isnan(tmpStats->min)) {
     1515        if (0 != p_psVectorMin(myVector, maskVector, maskVal, tmpStats)) {
    15161516            psLogMsg(__func__, PS_LOG_WARN,
    15171517                     "WARNING: p_psVectorMin(): p_psVectorMin() reported a NAN mean.\n");
    … …  
    15191519        }
    15201520    }
    1521     if (isnan(stats->max)) {
    1522         if (0 != p_psVectorMax(myVector, maskVector, maskVal, stats)) {
     1521    if (isnan(tmpStats->max)) {
     1522        if (0 != p_psVectorMax(myVector, maskVector, maskVal, tmpStats)) {
    15231523            psLogMsg(__func__, PS_LOG_WARN,
    15241524                     "WARNING: p_psVectorMin(): p_psVectorMax() reported a NAN mean.\n");
    … …  
    15311531    // bins, not the binSize.  Also, if we get here, we know that
    15321532    // binSize != 0.0.
    1533     numBins = (psS32)((stats->max - stats->min) / binSize);
    1534     robustHistogram = psHistogramAlloc(stats->min, stats->max, numBins);
     1533    numBins = (psS32)((tmpStats->max - tmpStats->min) / binSize);
     1534    robustHistogram = psHistogramAlloc(tmpStats->min, tmpStats->max, numBins);
    15351535
    15361536    // Populate the histogram array.
    … …  
    16201620    myStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
    16211621
    1622     // Using the above (myMean, myStdev) as initial estimates, we fit a
    1623     // Gaussian to the robustHistogramVector.
    1624     psMinimization *min = psMinimizationAlloc(100, 0.1);
    1625     psVector *myParams = psVectorAlloc(2, PS_TYPE_F32);
    1626     psArray *myCoords = psArrayAlloc(robustHistogramVector->n);
    1627     psVector *y = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
    1628 
    16291622    p_psNormalizeVectorRangeF32(robustHistogramVector, 0.0, 1.0);
    1630     for (i=0;i<robustHistogramVector->n;i++) {
    1631         myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
    1632         ((psVector *) (myCoords->data[i]))->data.F32[0] = (psF32) i;
    1633         y->data.F32[i] = robustHistogramVector->data.F32[i];
    1634     }
    1635 
    1636     myParams->data.F32[0] = myMean;
    1637     myParams->data.F32[1] = myStdev;
    1638     psImage *covar = NULL;
    1639     psMinimizeLMChi2(min,
    1640                      NULL,
    1641                      myParams,
    1642                      NULL,
    1643                      myCoords,
    1644                      y,
    1645                      NULL,
    1646                      (psMinimizeLMChi2Func) psMinimizeLMChi2Gauss1D);
    1647     psFree(covar);
    1648     psFree(min);
    1649     psFree(myParams);
    1650     psFree(myCoords);
    1651     psFree(y);
     1623
     1624    if ((stats->options & PS_STAT_ROBUST_MEAN) ||
     1625            (stats->options & PS_STAT_ROBUST_STDEV)) {
     1626
     1627        // Using the above (myMean, myStdev) as initial estimates, we fit a
     1628        // Gaussian to the robustHistogramVector.
     1629        psImage *covar = NULL;
     1630        psMinimization *min = psMinimizationAlloc(100, 0.1);
     1631        psVector *myParams = psVectorAlloc(2, PS_TYPE_F32);
     1632        psArray *myCoords = psArrayAlloc(robustHistogramVector->n);
     1633        psVector *y = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
     1634
     1635
     1636        myParams->data.F32[0] = myMean;
     1637        myParams->data.F32[1] = myStdev;
     1638
     1639        /* XXX: old.  Remove this.
     1640        for (i=0;i<robustHistogramVector->n;i++) {
     1641            myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
     1642            ((psVector *) (myCoords->data[i]))->data.F32[0] = (psF32) i;
     1643            y->data.F32[i] = robustHistogramVector->data.F32[i];
     1644        }
     1645        */
     1646
     1647        for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) {
     1648            int index = i - modeBinNum + dL;
     1649            myCoords->data[index] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
     1650            ((psVector *) (myCoords->data[index]))->data.F32[0] = PS_BIN_MIDPOINT(robustHistogram, i);
     1651            y->data.F32[index] = robustHistogramVector->data.F32[i];
     1652        }
     1653
     1654        bool rc = psMinimizeLMChi2(min,
     1655                                   NULL,
     1656                                   myParams,
     1657                                   NULL,
     1658                                   myCoords,
     1659                                   y,
     1660                                   NULL,
     1661                                   (psMinimizeLMChi2Func) psMinimizeLMChi2Gauss1D);
     1662        if (rc == false) {
     1663            psLogMsg(__func__, PS_LOG_WARN,
     1664                     "WARNING: failed to minimize with psMinimizeLMChi2().\n");
     1665        }
     1666
     1667        // XXX: Verify this with IfA
     1668        // XXX: The check on the minimization is better than the difference from myMean.
     1669        //      Do they still want this code?
     1670
     1671        if (stats->options & PS_STAT_ROBUST_MEAN) {
     1672            if (fabs((myParams->data.F32[0] - myMean)/myMean) > 0.1) {
     1673                psLogMsg(__func__, PS_LOG_WARN,
     1674                         "WARNING: the fitted Gaussian has more than 10%% error for the mean.\n");
     1675                psLogMsg(__func__, PS_LOG_WARN,
     1676                         "WARNING: Using the calculated mean instead of Gaussian-fitted mean.");
     1677                //                     "WARNING: Using the calculated mean instead of Gaussian-fitted mean.(calc, fit) is (%f, %f)\n",
     1678                //                     myMean, myParams->data.F32[0]);
     1679                stats->robustMean = myMean;
     1680            } else {
     1681                stats->robustMean = myParams->data.F32[0];
     1682            }
     1683        }
     1684
     1685
     1686        if (stats->options & PS_STAT_ROBUST_STDEV) {
     1687            if (fabs((myParams->data.F32[1] - myStdev)/myStdev) > 0.1) {
     1688                psLogMsg(__func__, PS_LOG_WARN,
     1689                         "WARNING: the fitted Gaussian has more than 10%% error for the stdev.\n");
     1690                psLogMsg(__func__, PS_LOG_WARN,
     1691                         "WARNING: Using the calculated stdev instead of Gaussian-fitted stdev.");
     1692                //                     "WARNING: Using the calculated stdev instead of Gaussian-fitted stdev.(calc, fit) is (%f, %f)\n",
     1693                //                     myStdev, myParams->data.F32[1]);
     1694                stats->robustStdev = myStdev;
     1695            } else {
     1696                stats->robustStdev = myParams->data.F32[1];
     1697            }
     1698        }
     1699        psFree(covar);
     1700        psFree(min);
     1701        psFree(myCoords);
     1702        psFree(y);
     1703        psFree(myParams);
     1704    }
     1705
     1706
    16521707    /**************************************************************************
    16531708    Set the appropriate members in the output stats struct.
    16541709    **************************************************************************/
    1655     if (stats->options & PS_STAT_ROBUST_MEAN) {
    1656         if (fabs((myParams->data.F32[0] - myMean)/myMean) > 0.1) {
    1657             psLogMsg(__func__, PS_LOG_WARN,
    1658                      "WARNING: the fitted Gaussian has more than 10%% error for the mean.\n");
    1659             psLogMsg(__func__, PS_LOG_WARN,
    1660                      "WARNING: Using the calculated mean instead of Gaussian-fitted mean.");
    1661             //                     "WARNING: Using the calculated mean instead of Gaussian-fitted mean.(calc, fit) is (%f, %f)\n",
    1662             //                     myMean, myParams->data.F32[0]);
    1663             stats->robustMean = myMean;
    1664         } else {
    1665             stats->robustMean = myParams->data.F32[0];
    1666         }
    1667     }
    16681710
    16691711    if (stats->options & PS_STAT_ROBUST_MODE) {
    … …  
    16711713    }
    16721714
    1673     if (stats->options & PS_STAT_ROBUST_STDEV) {
    1674         if (fabs((myParams->data.F32[1] - myStdev)/myStdev) > 0.1) {
    1675             psLogMsg(__func__, PS_LOG_WARN,
    1676                      "WARNING: the fitted Gaussian has more than 10%% error for the stdev.\n");
    1677             psLogMsg(__func__, PS_LOG_WARN,
    1678                      "WARNING: Using the calculated stdev instead of Gaussian-fitted stdev.");
    1679             //                     "WARNING: Using the calculated stdev instead of Gaussian-fitted stdev.(calc, fit) is (%f, %f)\n",
    1680             //                     myStdev, myParams->data.F32[1]);
    1681             stats->robustStdev = myStdev;
    1682         } else {
    1683             stats->robustStdev = myParams->data.F32[1];
    1684         }
    1685     }
    16861715
    16871716    // To determine the median, we fit a quadratic y=f(x) to the three bins
    … …  
    21742203 *****************************************************************************/
    21752204psStats* psVectorStats(psStats* stats,
    2176                        psVector* in,
    2177                        psVector* errors,
    2178                        psVector* mask,
     2205                       const psVector* in,
     2206                       const psVector* errors,
     2207                       const psVector* mask,
    21792208                       psU32 maskVal)
    21802209{
    … …  
    21972226    inF32 = p_psConvertToF32((psVector *) in);
    21982227    if (inF32 == NULL) {
    2199         inF32 = in;
     2228        inF32 = (psVector *) in;
    22002229        mustFreeVectorIn = 0;
    22012230    }
    22022231    errorsF32 = p_psConvertToF32((psVector *) errors);
    22032232    if (errorsF32 == NULL) {
    2204         errorsF32 = errors;
     2233        errorsF32 = (psVector *) errors;
    22052234        mustFreeVectorErrors = 0;
    22062235    }

  • trunk/psLib/src/dataManip/psStats.h

    r3264 r3547  
    44 *  @ingroup Stats
    55 *
    6  *  This file will hold the definition of the histogram and stats data
     6xd *  This file will hold the definition of the histogram and stats data
    77 *  structures.  It also contains prototypes for procedures which operate
    88 *  on those data structures.
    … …  
    1010 *  @author George Gusciora, MHPCC
    1111 *
    12  *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2005-02-17 19:26:23 $
     12 *  @version $Revision: 1.40 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2005-03-29 22:34:59 $
    1414 *
    1515 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    8787psStats* psVectorStats(
    8888    psStats* stats,    ///< stats structure defines stats to be calculated and how
    89     psVector* in,      ///< Vector to be analysed.
    90     psVector* errors,  ///< Errors.
    91     psVector* mask,    ///< Ignore elements where (maskVector & maskVal) != 0: must be INT or NULL
     89    const psVector* in,      ///< Vector to be analysed.
     90    const psVector* errors,  ///< Errors.
     91    const psVector* mask,    ///< Ignore elements where (maskVector & maskVal) != 0: must be INT or NULL
    9292    psU32 maskVal      ///< Only mask elements with one of these bits set in maskVector
    9393);

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

    r3540 r3547  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.108 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-03-29 19:41:56 $
     11 *  @version $Revision: 1.109 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-03-29 22:34:59 $
    1313 *
    1414 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    786786        //
    787787        if (currChi2 > newChi2) {
    788             min->lastDelta = currChi2 - newChi2;
     788            min->lastDelta = (currChi2 - newChi2)/currChi2;
    789789            min->value = newChi2;
    790790
    … …  
    14061406with an initial guess at the parameters of that function and vector "line"
    14071407of the same size as the parameter vector.  It will minimize the function
    1408 along that vector anr returns the offset along that vector at which the
     1408along that vector and returns the offset along that vector at which the
    14091409minimum is determined.
    14101410 
    … …  
    17131713            psFree(v);
    17141714            min->iter = iterationNumber;
     1715            // XXX: Ensure that currFuncVal is the correct value to use here.
     1716            min->value = currFuncVal;
    17151717            psTrace(".psLib.dataManip.psMinimizePowell", 4,
    17161718                    "---- psMinimizePowell() end (1)(true) ----\n");
    … …  
    17491751        if (fabs(baseFuncVal - currFuncVal) <= min->tol) {
    17501752            psFree(v);
     1753            // XXX: Ensure that currFuncVal is the correct value to use here.
     1754            min->value = currFuncVal;
    17511755            min->iter = iterationNumber;
    17521756            psTrace(".psLib.dataManip.psMinimizePowell", 4,

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

    r3540 r3547  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.122 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-03-29 19:41:56 $
     11 *  @version $Revision: 1.123 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-03-29 22:34:59 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    14891489    if (fabs(binSize) <= FLT_EPSILON) {
    14901490        if (stats->options & PS_STAT_ROBUST_MEAN) {
    1491             stats->robustMean = stats->clippedMean;
     1491            stats->robustMean = tmpStats->clippedMean;
    14921492        }
    14931493        if (stats->options & PS_STAT_ROBUST_MEDIAN) {
    1494             stats->robustMedian = stats->clippedMean;
     1494            stats->robustMedian = tmpStats->clippedMean;
    14951495        }
    14961496        if (stats->options & PS_STAT_ROBUST_MODE) {
    1497             stats->robustMode = stats->clippedMean;
     1497            stats->robustMode = tmpStats->clippedMean;
    14981498        }
    14991499        if (stats->options & PS_STAT_ROBUST_STDEV) {
    … …  
    15011501        }
    15021502        if (stats->options & PS_STAT_ROBUST_QUARTILE) {
    1503             stats->robustUQ = stats->clippedMean;
    1504             stats->robustLQ = stats->clippedMean;
     1503            stats->robustUQ = tmpStats->clippedMean;
     1504            stats->robustLQ = tmpStats->clippedMean;
    15051505        }
    15061506        // XXX: Set these to the number of unmasked data points?
    … …  
    15121512
    15131513    // Determine minimum and maximum values in the data vector.
    1514     if (isnan(stats->min)) {
    1515         if (0 != p_psVectorMin(myVector, maskVector, maskVal, stats)) {
     1514    if (isnan(tmpStats->min)) {
     1515        if (0 != p_psVectorMin(myVector, maskVector, maskVal, tmpStats)) {
    15161516            psLogMsg(__func__, PS_LOG_WARN,
    15171517                     "WARNING: p_psVectorMin(): p_psVectorMin() reported a NAN mean.\n");
    … …  
    15191519        }
    15201520    }
    1521     if (isnan(stats->max)) {
    1522         if (0 != p_psVectorMax(myVector, maskVector, maskVal, stats)) {
     1521    if (isnan(tmpStats->max)) {
     1522        if (0 != p_psVectorMax(myVector, maskVector, maskVal, tmpStats)) {
    15231523            psLogMsg(__func__, PS_LOG_WARN,
    15241524                     "WARNING: p_psVectorMin(): p_psVectorMax() reported a NAN mean.\n");
    … …  
    15311531    // bins, not the binSize.  Also, if we get here, we know that
    15321532    // binSize != 0.0.
    1533     numBins = (psS32)((stats->max - stats->min) / binSize);
    1534     robustHistogram = psHistogramAlloc(stats->min, stats->max, numBins);
     1533    numBins = (psS32)((tmpStats->max - tmpStats->min) / binSize);
     1534    robustHistogram = psHistogramAlloc(tmpStats->min, tmpStats->max, numBins);
    15351535
    15361536    // Populate the histogram array.
    … …  
    16201620    myStdev = PS_SQRT_F32((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1));
    16211621
    1622     // Using the above (myMean, myStdev) as initial estimates, we fit a
    1623     // Gaussian to the robustHistogramVector.
    1624     psMinimization *min = psMinimizationAlloc(100, 0.1);
    1625     psVector *myParams = psVectorAlloc(2, PS_TYPE_F32);
    1626     psArray *myCoords = psArrayAlloc(robustHistogramVector->n);
    1627     psVector *y = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
    1628 
    16291622    p_psNormalizeVectorRangeF32(robustHistogramVector, 0.0, 1.0);
    1630     for (i=0;i<robustHistogramVector->n;i++) {
    1631         myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
    1632         ((psVector *) (myCoords->data[i]))->data.F32[0] = (psF32) i;
    1633         y->data.F32[i] = robustHistogramVector->data.F32[i];
    1634     }
    1635 
    1636     myParams->data.F32[0] = myMean;
    1637     myParams->data.F32[1] = myStdev;
    1638     psImage *covar = NULL;
    1639     psMinimizeLMChi2(min,
    1640                      NULL,
    1641                      myParams,
    1642                      NULL,
    1643                      myCoords,
    1644                      y,
    1645                      NULL,
    1646                      (psMinimizeLMChi2Func) psMinimizeLMChi2Gauss1D);
    1647     psFree(covar);
    1648     psFree(min);
    1649     psFree(myParams);
    1650     psFree(myCoords);
    1651     psFree(y);
     1623
     1624    if ((stats->options & PS_STAT_ROBUST_MEAN) ||
     1625            (stats->options & PS_STAT_ROBUST_STDEV)) {
     1626
     1627        // Using the above (myMean, myStdev) as initial estimates, we fit a
     1628        // Gaussian to the robustHistogramVector.
     1629        psImage *covar = NULL;
     1630        psMinimization *min = psMinimizationAlloc(100, 0.1);
     1631        psVector *myParams = psVectorAlloc(2, PS_TYPE_F32);
     1632        psArray *myCoords = psArrayAlloc(robustHistogramVector->n);
     1633        psVector *y = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32);
     1634
     1635
     1636        myParams->data.F32[0] = myMean;
     1637        myParams->data.F32[1] = myStdev;
     1638
     1639        /* XXX: old.  Remove this.
     1640        for (i=0;i<robustHistogramVector->n;i++) {
     1641            myCoords->data[i] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
     1642            ((psVector *) (myCoords->data[i]))->data.F32[0] = (psF32) i;
     1643            y->data.F32[i] = robustHistogramVector->data.F32[i];
     1644        }
     1645        */
     1646
     1647        for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) {
     1648            int index = i - modeBinNum + dL;
     1649            myCoords->data[index] = (psPtr *) psVectorAlloc(2, PS_TYPE_F32);
     1650            ((psVector *) (myCoords->data[index]))->data.F32[0] = PS_BIN_MIDPOINT(robustHistogram, i);
     1651            y->data.F32[index] = robustHistogramVector->data.F32[i];
     1652        }
     1653
     1654        bool rc = psMinimizeLMChi2(min,
     1655                                   NULL,
     1656                                   myParams,
     1657                                   NULL,
     1658                                   myCoords,
     1659                                   y,
     1660                                   NULL,
     1661                                   (psMinimizeLMChi2Func) psMinimizeLMChi2Gauss1D);
     1662        if (rc == false) {
     1663            psLogMsg(__func__, PS_LOG_WARN,
     1664                     "WARNING: failed to minimize with psMinimizeLMChi2().\n");
     1665        }
     1666
     1667        // XXX: Verify this with IfA
     1668        // XXX: The check on the minimization is better than the difference from myMean.
     1669        //      Do they still want this code?
     1670
     1671        if (stats->options & PS_STAT_ROBUST_MEAN) {
     1672            if (fabs((myParams->data.F32[0] - myMean)/myMean) > 0.1) {
     1673                psLogMsg(__func__, PS_LOG_WARN,
     1674                         "WARNING: the fitted Gaussian has more than 10%% error for the mean.\n");
     1675                psLogMsg(__func__, PS_LOG_WARN,
     1676                         "WARNING: Using the calculated mean instead of Gaussian-fitted mean.");
     1677                //                     "WARNING: Using the calculated mean instead of Gaussian-fitted mean.(calc, fit) is (%f, %f)\n",
     1678                //                     myMean, myParams->data.F32[0]);
     1679                stats->robustMean = myMean;
     1680            } else {
     1681                stats->robustMean = myParams->data.F32[0];
     1682            }
     1683        }
     1684
     1685
     1686        if (stats->options & PS_STAT_ROBUST_STDEV) {
     1687            if (fabs((myParams->data.F32[1] - myStdev)/myStdev) > 0.1) {
     1688                psLogMsg(__func__, PS_LOG_WARN,
     1689                         "WARNING: the fitted Gaussian has more than 10%% error for the stdev.\n");
     1690                psLogMsg(__func__, PS_LOG_WARN,
     1691                         "WARNING: Using the calculated stdev instead of Gaussian-fitted stdev.");
     1692                //                     "WARNING: Using the calculated stdev instead of Gaussian-fitted stdev.(calc, fit) is (%f, %f)\n",
     1693                //                     myStdev, myParams->data.F32[1]);
     1694                stats->robustStdev = myStdev;
     1695            } else {
     1696                stats->robustStdev = myParams->data.F32[1];
     1697            }
     1698        }
     1699        psFree(covar);
     1700        psFree(min);
     1701        psFree(myCoords);
     1702        psFree(y);
     1703        psFree(myParams);
     1704    }
     1705
     1706
    16521707    /**************************************************************************
    16531708    Set the appropriate members in the output stats struct.
    16541709    **************************************************************************/
    1655     if (stats->options & PS_STAT_ROBUST_MEAN) {
    1656         if (fabs((myParams->data.F32[0] - myMean)/myMean) > 0.1) {
    1657             psLogMsg(__func__, PS_LOG_WARN,
    1658                      "WARNING: the fitted Gaussian has more than 10%% error for the mean.\n");
    1659             psLogMsg(__func__, PS_LOG_WARN,
    1660                      "WARNING: Using the calculated mean instead of Gaussian-fitted mean.");
    1661             //                     "WARNING: Using the calculated mean instead of Gaussian-fitted mean.(calc, fit) is (%f, %f)\n",
    1662             //                     myMean, myParams->data.F32[0]);
    1663             stats->robustMean = myMean;
    1664         } else {
    1665             stats->robustMean = myParams->data.F32[0];
    1666         }
    1667     }
    16681710
    16691711    if (stats->options & PS_STAT_ROBUST_MODE) {
    … …  
    16711713    }
    16721714
    1673     if (stats->options & PS_STAT_ROBUST_STDEV) {
    1674         if (fabs((myParams->data.F32[1] - myStdev)/myStdev) > 0.1) {
    1675             psLogMsg(__func__, PS_LOG_WARN,
    1676                      "WARNING: the fitted Gaussian has more than 10%% error for the stdev.\n");
    1677             psLogMsg(__func__, PS_LOG_WARN,
    1678                      "WARNING: Using the calculated stdev instead of Gaussian-fitted stdev.");
    1679             //                     "WARNING: Using the calculated stdev instead of Gaussian-fitted stdev.(calc, fit) is (%f, %f)\n",
    1680             //                     myStdev, myParams->data.F32[1]);
    1681             stats->robustStdev = myStdev;
    1682         } else {
    1683             stats->robustStdev = myParams->data.F32[1];
    1684         }
    1685     }
    16861715
    16871716    // To determine the median, we fit a quadratic y=f(x) to the three bins
    … …  
    21742203 *****************************************************************************/
    21752204psStats* psVectorStats(psStats* stats,
    2176                        psVector* in,
    2177                        psVector* errors,
    2178                        psVector* mask,
     2205                       const psVector* in,
     2206                       const psVector* errors,
     2207                       const psVector* mask,
    21792208                       psU32 maskVal)
    21802209{
    … …  
    21972226    inF32 = p_psConvertToF32((psVector *) in);
    21982227    if (inF32 == NULL) {
    2199         inF32 = in;
     2228        inF32 = (psVector *) in;
    22002229        mustFreeVectorIn = 0;
    22012230    }
    22022231    errorsF32 = p_psConvertToF32((psVector *) errors);
    22032232    if (errorsF32 == NULL) {
    2204         errorsF32 = errors;
     2233        errorsF32 = (psVector *) errors;
    22052234        mustFreeVectorErrors = 0;
    22062235    }

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

    r3264 r3547  
    44 *  @ingroup Stats
    55 *
    6  *  This file will hold the definition of the histogram and stats data
     6xd *  This file will hold the definition of the histogram and stats data
    77 *  structures.  It also contains prototypes for procedures which operate
    88 *  on those data structures.
    … …  
    1010 *  @author George Gusciora, MHPCC
    1111 *
    12  *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2005-02-17 19:26:23 $
     12 *  @version $Revision: 1.40 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2005-03-29 22:34:59 $
    1414 *
    1515 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    8787psStats* psVectorStats(
    8888    psStats* stats,    ///< stats structure defines stats to be calculated and how
    89     psVector* in,      ///< Vector to be analysed.
    90     psVector* errors,  ///< Errors.
    91     psVector* mask,    ///< Ignore elements where (maskVector & maskVal) != 0: must be INT or NULL
     89    const psVector* in,      ///< Vector to be analysed.
     90    const psVector* errors,  ///< Errors.
     91    const psVector* mask,    ///< Ignore elements where (maskVector & maskVal) != 0: must be INT or NULL
    9292    psU32 maskVal      ///< Only mask elements with one of these bits set in maskVector
    9393);
Note: See TracChangeset for help on using the changeset viewer.