Index: /trunk/psLib/src/dataManip/psConstants.h =================================================================== --- /trunk/psLib/src/dataManip/psConstants.h (revision 2738) +++ /trunk/psLib/src/dataManip/psConstants.h (revision 2739) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.44 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-16 21:42:39 $ + * @version $Revision: 1.45 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-16 23:04:17 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -320,4 +320,7 @@ NAME = psPolynomial1DAlloc(ORDER, TYPE); \ p_psMemSetPersistent(NAME, true); \ + p_psMemSetPersistent(NAME->coeff, true); \ + p_psMemSetPersistent(NAME->coeffErr, true); \ + p_psMemSetPersistent(NAME->mask, true); \ } \ Index: /trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- /trunk/psLib/src/dataManip/psMinimize.c (revision 2738) +++ /trunk/psLib/src/dataManip/psMinimize.c (revision 2739) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.93 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-07 22:28:29 $ + * @version $Revision: 1.94 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-16 23:04:17 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -585,8 +585,7 @@ } + // XXX: Generate code that modifies covar matrix if not-NULL. if (covar != NULL) { PS_IMAGE_CHECK_SIZE(covar, params->n, params->n, NULL); - } else { - covar = psImageAlloc(params->n, params->n, PS_TYPE_F64); } @@ -923,4 +922,5 @@ A = psImageAlloc(polyOrder, polyOrder, PS_TYPE_F64); ALUD = psImageAlloc(polyOrder, polyOrder, PS_TYPE_F64); + B = psVectorAlloc(polyOrder, PS_TYPE_F64); coeffs = psVectorAlloc(polyOrder, PS_TYPE_F64); Index: /trunk/psLib/src/dataManip/psStats.c =================================================================== --- /trunk/psLib/src/dataManip/psStats.c (revision 2738) +++ /trunk/psLib/src/dataManip/psStats.c (revision 2739) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.102 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-16 03:01:40 $ + * @version $Revision: 1.103 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-16 23:04:18 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -1118,4 +1118,7 @@ XXX: After you fit the polynomial, solve for X analytically. + +XXX: the vectors do not have to be the same length. Must insert the proper +tests to ensure that binNum is within acceptable ranges for both vectors. *****************************************************************************/ float fitQuadraticSearchForYThenReturnX(psVector *xVec, @@ -1126,13 +1129,19 @@ PS_VECTOR_CHECK_NULL(xVec, NAN); PS_VECTOR_CHECK_NULL(yVec, NAN); - PS_VECTOR_CHECK_TYPE(xVec, PS_TYPE_F64, NAN); - PS_VECTOR_CHECK_TYPE(yVec, PS_TYPE_F64, NAN); - PS_VECTOR_CHECK_SIZE_EQUAL(xVec, yVec, NAN); + PS_VECTOR_CHECK_TYPE(xVec, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_TYPE(yVec, PS_TYPE_F32, NAN); + // PS_VECTOR_CHECK_SIZE_EQUAL(xVec, yVec, NAN); PS_INT_CHECK_RANGE(binNum, 0, (xVec->n - 1), NAN); - - PS_VECTOR_DECLARE_ALLOC_STATIC(x, 3, PS_TYPE_F64); - PS_VECTOR_DECLARE_ALLOC_STATIC(y, 3, PS_TYPE_F64); - PS_VECTOR_DECLARE_ALLOC_STATIC(yErr, 3, PS_TYPE_F64); - PS_POLY_1D_DECLARE_ALLOC_STATIC(myPoly, 2, PS_POLYNOMIAL_ORD); + PS_INT_CHECK_RANGE(binNum, 0, (yVec->n - 1), NAN); + + // PS_VECTOR_DECLARE_ALLOC_STATIC(x, 3, PS_TYPE_F64); + // PS_VECTOR_DECLARE_ALLOC_STATIC(y, 3, PS_TYPE_F64); + // PS_VECTOR_DECLARE_ALLOC_STATIC(yErr, 3, PS_TYPE_F64); + // PS_POLY_1D_DECLARE_ALLOC_STATIC(myPoly, 2, PS_POLYNOMIAL_ORD); + psVector *x = psVectorAlloc(3, PS_TYPE_F64); + psVector *y = psVectorAlloc(3, PS_TYPE_F64); + psVector *yErr = psVectorAlloc(3, PS_TYPE_F64); + psPolynomial1D *myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); + float tmpFloat; @@ -1162,4 +1171,8 @@ false, PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLYNOMIAL_1D_FIT); + psFree(myPoly); + psFree(x); + psFree(y); + psFree(yErr); return(NAN); } @@ -1171,4 +1184,8 @@ false, PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLY_MEDIAN); + psFree(myPoly); + psFree(x); + psFree(y); + psFree(yErr); return(NAN); } @@ -1192,4 +1209,8 @@ } + psFree(myPoly); + psFree(x); + psFree(y); + psFree(yErr); return(tmpFloat); } @@ -1394,5 +1415,4 @@ myStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); - // Using the above (myMean, myStdev) as initial estimates, we fit a // Gaussian to the robustHistogramVector. @@ -1411,4 +1431,5 @@ myParams->data.F32[0] = myMean; myParams->data.F32[1] = myStdev; + psImage *covar = NULL; psMinimizeLMChi2(min, NULL, @@ -1419,4 +1440,5 @@ NULL, (psMinimizeLMChi2Func) psMinimizeLMChi2Gauss1D); + psFree(covar); psFree(min); psFree(myParams); @@ -1429,7 +1451,7 @@ if (fabs((myParams->data.F32[0] - myMean)/myMean) > 0.1) { psLogMsg(__func__, PS_LOG_WARN, - "WARNING: the fitted Gaussian has more than 10% error for the mean.\n"); + "WARNING: the fitted Gaussian has more than 10%% error for the mean.\n"); psLogMsg(__func__, PS_LOG_WARN, - "WARNING: Using the calculated mean (calc, fit) is (%f, %f)\n", + "WARNING: Using the calculated mean instead of Gaussian-fitted mean.(calc, fit) is (%f, %f)\n", myMean, myParams->data.F32[0]); stats->robustMean = myMean; @@ -1446,8 +1468,8 @@ if (fabs((myParams->data.F32[1] - myStdev)/myStdev) > 0.1) { psLogMsg(__func__, PS_LOG_WARN, - "WARNING: the fitted Gaussian has more than 10% error for the mean.\n"); + "WARNING: the fitted Gaussian has more than 10%% error for the stdev.\n"); psLogMsg(__func__, PS_LOG_WARN, - "WARNING: Using the calculated mean (calc, fit) is (%f, %f)\n", - myMean, myParams->data.F32[0]); + "WARNING: Using the calculated stdev instead of Gaussian-fitted stdev.(calc, fit) is (%f, %f)\n", + myStdev, myParams->data.F32[1]); stats->robustStdev = myStdev; } else { @@ -1463,8 +1485,12 @@ // Take a psVector. Fit a polynomial y = f(x) to the 3 data elements // surrounding medianBinNum, then find the x-value corresponding y = sumRobust/2.0. - stats->robustMedian = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds, - cumulativeRobustSums, - medianBinNum, - sumRobust/2.0); + + //XXX: robustHistogram->bounds and cumulativeRobustSums are of different lengths. + //Determine id that is okay. + stats->robustMedian = fitQuadraticSearchForYThenReturnX( + robustHistogram->bounds, + cumulativeRobustSums, + medianBinNum, + sumRobust/2.0); } @@ -1476,12 +1502,18 @@ countFloat = cumulativeRobustSums->data.F32[robustHistogramVector->n - 1]; - stats->robustLQ = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds, - cumulativeRobustSums, - LQBinNum, - countFloat/4.0); - stats->robustUQ = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds, - cumulativeRobustSums, - UQBinNum, - 3.0 * countFloat/4.0); + //XXX: robustHistogram->bounds and cumulativeRobustSums are of different lengths. + //Determine id that is okay. + stats->robustLQ = fitQuadraticSearchForYThenReturnX( + robustHistogram->bounds, + cumulativeRobustSums, + LQBinNum, + countFloat/4.0); + //XXX: robustHistogram->bounds and cumulativeRobustSums are of different lengths. + //Determine id that is okay. + stats->robustUQ = fitQuadraticSearchForYThenReturnX( + robustHistogram->bounds, + cumulativeRobustSums, + UQBinNum, + 3.0 * countFloat/4.0); } // XXX: I think sumNfit == sumN50 here. Index: /trunk/psLib/src/math/psConstants.h =================================================================== --- /trunk/psLib/src/math/psConstants.h (revision 2738) +++ /trunk/psLib/src/math/psConstants.h (revision 2739) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.44 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-16 21:42:39 $ + * @version $Revision: 1.45 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-16 23:04:17 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -320,4 +320,7 @@ NAME = psPolynomial1DAlloc(ORDER, TYPE); \ p_psMemSetPersistent(NAME, true); \ + p_psMemSetPersistent(NAME->coeff, true); \ + p_psMemSetPersistent(NAME->coeffErr, true); \ + p_psMemSetPersistent(NAME->mask, true); \ } \ Index: /trunk/psLib/src/math/psMinimize.c =================================================================== --- /trunk/psLib/src/math/psMinimize.c (revision 2738) +++ /trunk/psLib/src/math/psMinimize.c (revision 2739) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.93 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-07 22:28:29 $ + * @version $Revision: 1.94 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-16 23:04:17 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -585,8 +585,7 @@ } + // XXX: Generate code that modifies covar matrix if not-NULL. if (covar != NULL) { PS_IMAGE_CHECK_SIZE(covar, params->n, params->n, NULL); - } else { - covar = psImageAlloc(params->n, params->n, PS_TYPE_F64); } @@ -923,4 +922,5 @@ A = psImageAlloc(polyOrder, polyOrder, PS_TYPE_F64); ALUD = psImageAlloc(polyOrder, polyOrder, PS_TYPE_F64); + B = psVectorAlloc(polyOrder, PS_TYPE_F64); coeffs = psVectorAlloc(polyOrder, PS_TYPE_F64); Index: /trunk/psLib/src/math/psStats.c =================================================================== --- /trunk/psLib/src/math/psStats.c (revision 2738) +++ /trunk/psLib/src/math/psStats.c (revision 2739) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.102 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-12-16 03:01:40 $ + * @version $Revision: 1.103 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-12-16 23:04:18 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -1118,4 +1118,7 @@ XXX: After you fit the polynomial, solve for X analytically. + +XXX: the vectors do not have to be the same length. Must insert the proper +tests to ensure that binNum is within acceptable ranges for both vectors. *****************************************************************************/ float fitQuadraticSearchForYThenReturnX(psVector *xVec, @@ -1126,13 +1129,19 @@ PS_VECTOR_CHECK_NULL(xVec, NAN); PS_VECTOR_CHECK_NULL(yVec, NAN); - PS_VECTOR_CHECK_TYPE(xVec, PS_TYPE_F64, NAN); - PS_VECTOR_CHECK_TYPE(yVec, PS_TYPE_F64, NAN); - PS_VECTOR_CHECK_SIZE_EQUAL(xVec, yVec, NAN); + PS_VECTOR_CHECK_TYPE(xVec, PS_TYPE_F32, NAN); + PS_VECTOR_CHECK_TYPE(yVec, PS_TYPE_F32, NAN); + // PS_VECTOR_CHECK_SIZE_EQUAL(xVec, yVec, NAN); PS_INT_CHECK_RANGE(binNum, 0, (xVec->n - 1), NAN); - - PS_VECTOR_DECLARE_ALLOC_STATIC(x, 3, PS_TYPE_F64); - PS_VECTOR_DECLARE_ALLOC_STATIC(y, 3, PS_TYPE_F64); - PS_VECTOR_DECLARE_ALLOC_STATIC(yErr, 3, PS_TYPE_F64); - PS_POLY_1D_DECLARE_ALLOC_STATIC(myPoly, 2, PS_POLYNOMIAL_ORD); + PS_INT_CHECK_RANGE(binNum, 0, (yVec->n - 1), NAN); + + // PS_VECTOR_DECLARE_ALLOC_STATIC(x, 3, PS_TYPE_F64); + // PS_VECTOR_DECLARE_ALLOC_STATIC(y, 3, PS_TYPE_F64); + // PS_VECTOR_DECLARE_ALLOC_STATIC(yErr, 3, PS_TYPE_F64); + // PS_POLY_1D_DECLARE_ALLOC_STATIC(myPoly, 2, PS_POLYNOMIAL_ORD); + psVector *x = psVectorAlloc(3, PS_TYPE_F64); + psVector *y = psVectorAlloc(3, PS_TYPE_F64); + psVector *yErr = psVectorAlloc(3, PS_TYPE_F64); + psPolynomial1D *myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); + float tmpFloat; @@ -1162,4 +1171,8 @@ false, PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLYNOMIAL_1D_FIT); + psFree(myPoly); + psFree(x); + psFree(y); + psFree(yErr); return(NAN); } @@ -1171,4 +1184,8 @@ false, PS_ERRORTEXT_psStats_STATS_FIT_QUADRATIC_POLY_MEDIAN); + psFree(myPoly); + psFree(x); + psFree(y); + psFree(yErr); return(NAN); } @@ -1192,4 +1209,8 @@ } + psFree(myPoly); + psFree(x); + psFree(y); + psFree(yErr); return(tmpFloat); } @@ -1394,5 +1415,4 @@ myStdev = sqrt((sumSquares - (sumDiffs * sumDiffs / countFloat)) / (countFloat - 1)); - // Using the above (myMean, myStdev) as initial estimates, we fit a // Gaussian to the robustHistogramVector. @@ -1411,4 +1431,5 @@ myParams->data.F32[0] = myMean; myParams->data.F32[1] = myStdev; + psImage *covar = NULL; psMinimizeLMChi2(min, NULL, @@ -1419,4 +1440,5 @@ NULL, (psMinimizeLMChi2Func) psMinimizeLMChi2Gauss1D); + psFree(covar); psFree(min); psFree(myParams); @@ -1429,7 +1451,7 @@ if (fabs((myParams->data.F32[0] - myMean)/myMean) > 0.1) { psLogMsg(__func__, PS_LOG_WARN, - "WARNING: the fitted Gaussian has more than 10% error for the mean.\n"); + "WARNING: the fitted Gaussian has more than 10%% error for the mean.\n"); psLogMsg(__func__, PS_LOG_WARN, - "WARNING: Using the calculated mean (calc, fit) is (%f, %f)\n", + "WARNING: Using the calculated mean instead of Gaussian-fitted mean.(calc, fit) is (%f, %f)\n", myMean, myParams->data.F32[0]); stats->robustMean = myMean; @@ -1446,8 +1468,8 @@ if (fabs((myParams->data.F32[1] - myStdev)/myStdev) > 0.1) { psLogMsg(__func__, PS_LOG_WARN, - "WARNING: the fitted Gaussian has more than 10% error for the mean.\n"); + "WARNING: the fitted Gaussian has more than 10%% error for the stdev.\n"); psLogMsg(__func__, PS_LOG_WARN, - "WARNING: Using the calculated mean (calc, fit) is (%f, %f)\n", - myMean, myParams->data.F32[0]); + "WARNING: Using the calculated stdev instead of Gaussian-fitted stdev.(calc, fit) is (%f, %f)\n", + myStdev, myParams->data.F32[1]); stats->robustStdev = myStdev; } else { @@ -1463,8 +1485,12 @@ // Take a psVector. Fit a polynomial y = f(x) to the 3 data elements // surrounding medianBinNum, then find the x-value corresponding y = sumRobust/2.0. - stats->robustMedian = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds, - cumulativeRobustSums, - medianBinNum, - sumRobust/2.0); + + //XXX: robustHistogram->bounds and cumulativeRobustSums are of different lengths. + //Determine id that is okay. + stats->robustMedian = fitQuadraticSearchForYThenReturnX( + robustHistogram->bounds, + cumulativeRobustSums, + medianBinNum, + sumRobust/2.0); } @@ -1476,12 +1502,18 @@ countFloat = cumulativeRobustSums->data.F32[robustHistogramVector->n - 1]; - stats->robustLQ = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds, - cumulativeRobustSums, - LQBinNum, - countFloat/4.0); - stats->robustUQ = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds, - cumulativeRobustSums, - UQBinNum, - 3.0 * countFloat/4.0); + //XXX: robustHistogram->bounds and cumulativeRobustSums are of different lengths. + //Determine id that is okay. + stats->robustLQ = fitQuadraticSearchForYThenReturnX( + robustHistogram->bounds, + cumulativeRobustSums, + LQBinNum, + countFloat/4.0); + //XXX: robustHistogram->bounds and cumulativeRobustSums are of different lengths. + //Determine id that is okay. + stats->robustUQ = fitQuadraticSearchForYThenReturnX( + robustHistogram->bounds, + cumulativeRobustSums, + UQBinNum, + 3.0 * countFloat/4.0); } // XXX: I think sumNfit == sumN50 here.