Index: trunk/psLib/src/astro/psCoord.c =================================================================== --- trunk/psLib/src/astro/psCoord.c (revision 10831) +++ trunk/psLib/src/astro/psCoord.c (revision 10848) @@ -10,6 +10,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.128 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-12-22 21:19:47 $ +* @version $Revision: 1.129 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-12-29 04:38:42 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -905,6 +905,7 @@ } - trans->x = psVectorFitPolynomial2D(trans->x, NULL, 0, xOut, NULL, xIn, yIn); - trans->y = psVectorFitPolynomial2D(trans->y, NULL, 0, yOut, NULL, xIn, yIn); + bool result = true; + result &= psVectorFitPolynomial2D(trans->x, NULL, 0, xOut, NULL, xIn, yIn); + result &= psVectorFitPolynomial2D(trans->y, NULL, 0, yOut, NULL, xIn, yIn); psFree(xIn); psFree(yIn); @@ -912,5 +913,5 @@ psFree(yOut); - if ((trans->x == NULL) || (trans->y == NULL)) { + if (!result) { psError( PS_ERR_UNKNOWN, true, "psVectorFitPolynomial2D() returned NULL: could not fit a 2-D polynomial to the data.\n"); return(false); @@ -997,6 +998,7 @@ } - myPT->x = psVectorFitPolynomial2D(myPT->x, NULL, 0, xOut, NULL, xIn, yIn); - myPT->y = psVectorFitPolynomial2D(myPT->y, NULL, 0, yOut, NULL, xIn, yIn); + bool result = true; + result &= psVectorFitPolynomial2D(myPT->x, NULL, 0, xOut, NULL, xIn, yIn); + result &= psVectorFitPolynomial2D(myPT->y, NULL, 0, yOut, NULL, xIn, yIn); psFree(inCoord); @@ -1007,5 +1009,5 @@ psFree(yOut); - if ((myPT->x == NULL) || (myPT->y == NULL)) { + if (!result) { psError( PS_ERR_UNKNOWN, true, "psVectorFitPolynomial2D() returned NULL: could not fit a 2-D polynomial to the data.\n"); psFree(out); Index: trunk/psLib/src/math/psMinimizePolyFit.c =================================================================== --- trunk/psLib/src/math/psMinimizePolyFit.c (revision 10831) +++ trunk/psLib/src/math/psMinimizePolyFit.c (revision 10848) @@ -10,6 +10,6 @@ * @author EAM, IfA * - * @version $Revision: 1.26 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-12-17 09:43:48 $ + * @version $Revision: 1.27 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-12-29 04:38:42 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -55,4 +55,8 @@ }\ }\ + +// free a local temporary F64 vector (TEMP) which is a copy of a non-F64 vector (ORIG) +# define PS_FREE_TEMP_F64_VECTOR(ORIG, TEMP) \ +if ((ORIG != NULL) && (ORIG->type.type != PS_TYPE_F64)) { psFree(TEMP); } /*****************************************************************************/ @@ -281,5 +285,4 @@ *****************************************************************************/ - /****************************************************************************** ****************************************************************************** @@ -288,12 +291,4 @@ *****************************************************************************/ -static psPolynomial1D* VectorFitPolynomial1DOrd( - psPolynomial1D* myPoly, - const psVector *mask, - psMaskType maskValue, - const psVector *f, - const psVector *fErr, - const psVector *x); - /****************************************************************************** vectorFitPolynomial1DCheb(): This routine will fit a Chebyshev @@ -301,5 +296,5 @@ coefficients of that polynomial. *****************************************************************************/ -static psPolynomial1D *vectorFitPolynomial1DCheb( +static bool vectorFitPolynomial1DCheb( psPolynomial1D* myPoly, const psVector *mask, @@ -404,6 +399,5 @@ if (!psMatrixGJSolve(A, B)) { psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations. Returning NULL.\n"); - psFree(myPoly); - myPoly = NULL; + goto escape_GJ; } else { // the first nTerm entries in B correspond directly to the desired @@ -427,12 +421,10 @@ if (ALUD == NULL) { psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix. Returning NULL.\n"); - psFree(myPoly); - myPoly = NULL; + goto escape_LUD; } else { coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); if (coeffs == NULL) { psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix. Returning NULL.\n"); - psFree(myPoly); - myPoly = NULL; + goto escape_LUD; } else { for (psS32 k = 0; k < numTerms; k++) { @@ -441,15 +433,22 @@ } } - - psFree(ALUD); psFree(coeffs); psFree(outPerm); } - psFree(A); psFree(B); - - return(myPoly); + return true; + +escape_LUD: + // XXX drop the LUD version! + psFree(A); + psFree(B); + return false; + +escape_GJ: + psFree(A); + psFree(B); + return false; } @@ -459,6 +458,9 @@ x and fErr vectors may be NULL. All non-NULL vectors must be of type PS_TYPE_F64. - *****************************************************************************/ -static psPolynomial1D* VectorFitPolynomial1DOrd( + +XXX EAM : since this is a private function, can we drop the ASSERTS? +XXX EAM : can we drop the LUD version? it does not calculate coeffErr values!! +*****************************************************************************/ +static bool VectorFitPolynomial1DOrd( psPolynomial1D* myPoly, const psVector *mask, @@ -469,18 +471,18 @@ { psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); - PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL); + PS_ASSERT_POLY_NON_NULL(myPoly, false); + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, false); if (mask) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); } if (x) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, false); } if (fErr) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, false); } @@ -510,9 +512,8 @@ if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) { psError(PS_ERR_UNKNOWN, false, "Could initialize data structures A, B. Returning NULL.\n"); - psFree(myPoly); psFree(A); psFree(B); psTrace("psLib.math", 4, "---- %s() End ----\n", __func__); - return(NULL); + return false; } @@ -596,6 +597,5 @@ if (!psMatrixGJSolve(A, B)) { psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations. Returning NULL.\n"); - psFree(myPoly); - myPoly = NULL; + goto escape_GJ; } else { // the first nTerm entries in B correspond directly to the desired @@ -617,12 +617,10 @@ if (ALUD == NULL) { psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix. Returning NULL.\n"); - psFree(myPoly); - myPoly = NULL; + goto escape_LUD; } else { coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); if (coeffs == NULL) { psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix. Returning NULL.\n"); - psFree(myPoly); - myPoly = NULL; + goto escape_LUD; } else { for (psS32 k = 0; k < nTerm; k++) { @@ -632,16 +630,23 @@ } } - psFree(ALUD); psFree(coeffs); psFree(outPerm); } - - psFree(A); psFree(B); psTrace("psLib.math", 4, "---- %s() End ----\n", __func__); - return (myPoly); + return true; + +escape_LUD: + psFree(A); + psFree(B); + return false; + +escape_GJ: + psFree(A); + psFree(B); + return false; } @@ -653,5 +658,5 @@ conversion only. *****************************************************************************/ -psPolynomial1D *psVectorFitPolynomial1D( +bool psVectorFitPolynomial1D( psPolynomial1D *poly, const psVector *mask, @@ -661,41 +666,40 @@ const psVector *x) { - // Internal pointers for possibly NULL or mis-typed vectors. + PS_ASSERT_POLY_NON_NULL(poly, false); + PS_ASSERT_INT_NONNEGATIVE(poly->nX, false); + + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_NON_EMPTY(f, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false); + if (mask != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); + } + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, false); + } + if (x != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, false); + } + + // Convert input vectors to F64 if necessary. + psVector *f64 = (f->type.type == PS_TYPE_F64) ? (psVector *) f : psVectorCopy (NULL, f, PS_TYPE_F64); psVector *x64 = NULL; - psVector *f64 = NULL; + if (x != NULL) { + x64 = (x->type.type == PS_TYPE_F64) ? (psVector *) x : psVectorCopy (NULL, x, PS_TYPE_F64); + } psVector *fErr64 = NULL; - - PS_ASSERT_POLY_NON_NULL(poly, NULL); - //PS_ASSERT_INT_NONNEGATIVE(poly->nX, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_NON_EMPTY(f, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); - if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - if (x != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(x, NULL); - } if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL); - } - - f64 = (f->type.type == PS_TYPE_F64) ? (psVector *) f : psVectorCopy(NULL, f, PS_TYPE_F64); - - if (x != NULL) { - x64 = (x->type.type == PS_TYPE_F64) ? (psVector *) x : psVectorCopy(NULL, x, PS_TYPE_F64); - } - - if (fErr != NULL) { - fErr64 = (fErr->type.type == PS_TYPE_F64) ? (psVector *) fErr : psVectorCopy(NULL, fErr, PS_TYPE_F64); - } + fErr64 = (fErr->type.type == PS_TYPE_F64) ? (psVector *) fErr : psVectorCopy (NULL, fErr, PS_TYPE_F64); + } + + bool result = true; switch (poly->type) { case PS_POLYNOMIAL_ORD: - poly = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, x64); - if (poly == NULL) { + result = VectorFitPolynomial1DOrd(poly, mask, maskValue, f64, fErr64, x64); + if (!result) { psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial. Returning NULL.\n"); } @@ -703,5 +707,5 @@ case PS_POLYNOMIAL_CHEB: if (mask != NULL) { - // psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n"); + psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n"); } if (fErr != NULL) { @@ -713,5 +717,8 @@ } - poly = vectorFitPolynomial1DCheb(poly, mask, maskValue, f64, fErr64, x64); + result = vectorFitPolynomial1DCheb(poly, mask, maskValue, f64, fErr64, x64); + if (!result) { + psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial. Returning NULL.\n"); + } if (x == NULL) { @@ -721,26 +728,18 @@ default: psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type (%d). Returning NULL.\n", poly->type); - poly = NULL; + result = false; break; } // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if ((x != NULL) && (x->type.type != PS_TYPE_F64)) { - psFree(x64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } - - return(poly); + PS_FREE_TEMP_F64_VECTOR (f, f64); + PS_FREE_TEMP_F64_VECTOR (x, x64); + PS_FREE_TEMP_F64_VECTOR (fErr, fErr64); + + return result; } // This function accepts F32 and F64 input vectors. -psPolynomial1D *psVectorClipFitPolynomial1D( +bool psVectorClipFitPolynomial1D( psPolynomial1D *poly, psStats *stats, @@ -752,14 +751,15 @@ { psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_PTR_NON_NULL(stats, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); - PS_ASSERT_VECTOR_NON_NULL(mask, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(mask, f, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); + PS_ASSERT_POLY_NON_NULL(poly, false); + PS_ASSERT_PTR_NON_NULL(stats, false); + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false); + PS_ASSERT_VECTOR_NON_NULL(mask, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(mask, f, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); + if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, f, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(fErr, f, false); + PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false); } @@ -767,6 +767,6 @@ psVector *x = NULL; if (xIn != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(xIn, f, NULL); - PS_ASSERT_VECTOR_TYPE(xIn, f->type.type, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(xIn, f, false); + PS_ASSERT_VECTOR_TYPE(xIn, f->type.type, false); x = (psVector *) xIn; } else { @@ -781,6 +781,21 @@ } else { psError(PS_ERR_UNKNOWN, true, "Error, bad poly type.\n"); - return(NULL); - } + return false; + } + } + + // the user supplies one of various stats option pairs, + // determine the desired mean and stdev STATS options: + // XXX enforce consistency? + // XXX psStatsGetValue() probably has inverted precedence + psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2); + psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2); + if (!meanOption) { + psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected"); + return false; + } + if (!stdevOption) { + psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected"); + return false; } @@ -798,14 +813,6 @@ minClipSigma = fabs(stats->clipSigma); } - - psVector *fit = NULL; psVector *resid = psVectorAlloc(f->n, PS_TYPE_F64); - // eventual expansion: user supplies one of various stats option pairs, - // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to - // evaluate the clipping sigma - // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used - stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter); psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma); @@ -822,18 +829,17 @@ } } - poly = psVectorFitPolynomial1D(poly, mask, maskValue, f, fErr, x); - if (poly == NULL) { - psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial. Returning NULL.\n"); + if (!psVectorFitPolynomial1D(poly, mask, maskValue, f, fErr, x)) { + psError(PS_ERR_UNKNOWN, false, "Could not fit polynomial. Returning false.\n"); if (xIn == NULL) { psFree(x); } - return(NULL); - } - - fit = psPolynomial1DEvalVector(poly, x); + return false; + } + + psVector *fit = psPolynomial1DEvalVector(poly, x); if (fit == NULL) { - psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector(). Returning NULL.\n"); + psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector(). Returning false.\n"); psFree(resid); - return(NULL); + return false; } for (psS32 i = 0 ; i < f->n ; i++) { @@ -856,23 +862,17 @@ if (!psVectorStats(stats, resid, NULL, mask, maskValue)) { - psError(PS_ERR_UNKNOWN, false, "Could not compute statistics on the resid vector. Returning NULL.\n"); - psFree(resid) - psFree(fit) - return(NULL); - } - # if (USE_ROBUST_STATS_FOR_CLIPPING) - psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian); - psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev); - psF32 minClipValue = -minClipSigma*stats->robustStdev; - psF32 maxClipValue = +maxClipSigma*stats->robustStdev; - psF32 clipMedian = stats->robustMedian; - # else - - psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian); - psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev); - psF32 minClipValue = -minClipSigma*stats->robustStdev; - psF32 maxClipValue = +maxClipSigma*stats->robustStdev; - psF32 clipMedian = stats->sampleMedian; - # endif + psError(PS_ERR_UNKNOWN, false, "Could not compute statistics on the resid vector. Returning false.\n"); + psFree(resid); + psFree(fit); + return false; + } + + double meanValue = psStatsGetValue (stats, meanOption); + double stdevValue = psStatsGetValue (stats, stdevOption); + + psTrace("psLib.math", 5, "Mean is %f\n", meanValue); + psTrace("psLib.math", 5, "Stdev is %f\n", stdevValue); + psF32 minClipValue = -minClipSigma*stdevValue; + psF32 maxClipValue = +maxClipSigma*stdevValue; // set mask if pts are not valid @@ -884,5 +884,5 @@ } - if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue)) { + if ((resid->data.F64[i] - meanValue > maxClipValue) || (resid->data.F64[i] - meanValue < minClipValue)) { if (f->type.type == PS_TYPE_F64) { psTrace("psLib.math", 6, "Masking element %d (%f). resid->data.F64[%d] is %f\n", @@ -906,4 +906,5 @@ // psTrace("psLib.math", 6, "keeping %d of %ld pts for fit\n", Nkeep, x->n); + stats->clippedNvalues = Nkeep; psFree(fit); } @@ -917,5 +918,5 @@ psTrace("psLib.math", 3, "---- %s() end ----\n", __func__); - return (poly); + return true; } @@ -933,5 +934,5 @@ *****************************************************************************/ -static psPolynomial2D* VectorFitPolynomial2DOrd( +static bool VectorFitPolynomial2DOrd( psPolynomial2D* myPoly, const psVector* mask, @@ -943,22 +944,22 @@ { psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); - PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL); + PS_ASSERT_POLY_NON_NULL(myPoly, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, false); + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, false); if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, false); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, false); + } + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false); if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); } @@ -974,9 +975,8 @@ if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) { psError(PS_ERR_UNKNOWN, false, "Could initialize data structures A, B. Returning NULL.\n"); - psFree(myPoly); psFree(A); psFree(B); psTrace("psLib.math", 4, "---- %s() End ----\n", __func__); - return(NULL); + return false; } @@ -1040,21 +1040,21 @@ if (!psMatrixGJSolve(A, B)) { psError(PS_ERR_UNKNOWN, false, "Could not solve linear equations. Returning NULL.\n"); - psFree(myPoly); - myPoly = NULL; - } else { - // select the appropriate solution entries - for (int i = 0; i < nTerm; i++) { - int l = i / nYterm; // x index - int m = i % nYterm; // y index - myPoly->coeff[l][m] = B->data.F64[i]; - myPoly->coeffErr[l][m] = sqrt(A->data.F64[i][i]); - } - } - + psFree(A); + psFree(B); + return false; + } + + // select the appropriate solution entries + for (int i = 0; i < nTerm; i++) { + int l = i / nYterm; // x index + int m = i % nYterm; // y index + myPoly->coeff[l][m] = B->data.F64[i]; + myPoly->coeffErr[l][m] = sqrt(A->data.F64[i][i]); + } psFree(A); psFree(B); psTrace("psLib.math", 4, "---- %s() end ----\n", __func__); - return (myPoly); + return true; } @@ -1065,5 +1065,5 @@ vector conversion only. *****************************************************************************/ -psPolynomial2D *psVectorFitPolynomial2D( +bool psVectorFitPolynomial2D( psPolynomial2D *poly, const psVector *mask, @@ -1074,132 +1074,64 @@ const psVector *y) { - // Internal pointers for possibly NULL or mis-typed vectors. - psVector *x64 = NULL; - psVector *y64 = NULL; - psVector *f64 = NULL; + PS_ASSERT_POLY_NON_NULL(poly, false); + PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false); + + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false); + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false); + if (mask != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); + } + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, false); + } + + // Convert input vectors to F64 if necessary. + psVector *f64 = (f->type.type == PS_TYPE_F64) ? (psVector *) f : psVectorCopy(NULL, f, PS_TYPE_F64); + psVector *x64 = (x->type.type == PS_TYPE_F64) ? (psVector *) x : psVectorCopy(NULL, x, PS_TYPE_F64); + psVector *y64 = (y->type.type == PS_TYPE_F64) ? (psVector *) y : psVectorCopy(NULL, y, PS_TYPE_F64); + psVector *fErr64 = NULL; - - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); - - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); - - if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL); - } - - // - // Convert input vectors to F64 if necessary. - // - if (f->type.type != PS_TYPE_F64) { - f64 = psVectorCopy(NULL, f, PS_TYPE_F64); - } else { - f64 = (psVector *) f; - } - - if (x->type.type != PS_TYPE_F64) { - x64 = psVectorCopy(NULL, x, PS_TYPE_F64); - } else { - x64 = (psVector *) x; - } - - if (y->type.type != PS_TYPE_F64) { - y64 = psVectorCopy(NULL, y, PS_TYPE_F64); - } else { - y64 = (psVector *) y; - } - - // - // fErr - // - if (fErr != NULL) { - if (fErr->type.type != PS_TYPE_F64) { - fErr64 = psVectorCopy(NULL, fErr, PS_TYPE_F64); - } else { - fErr64 = (psVector *) fErr; - } - } - - if (poly->type == PS_POLYNOMIAL_ORD) { - poly = VectorFitPolynomial2DOrd(poly, mask, maskValue, f64, fErr64, x64, y64); - if (poly == NULL) { + fErr64 = (fErr->type.type == PS_TYPE_F64) ? (psVector *) fErr : psVectorCopy(NULL, fErr, PS_TYPE_F64); + } + + bool result = true; + + switch (poly->type) { + case PS_POLYNOMIAL_ORD: + result = VectorFitPolynomial2DOrd(poly, mask, maskValue, f64, fErr64, x64, y64); + if (!result) { psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial. Returning NULL.\n"); - // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if (x->type.type != PS_TYPE_F64) { - psFree(x64); - } - - if (y->type.type != PS_TYPE_F64) { - psFree(y64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } - return(NULL); - } - } else if (poly->type == PS_POLYNOMIAL_CHEB) { + } + break; + case PS_POLYNOMIAL_CHEB: if (mask != NULL) { psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n"); } - psLogMsg(__func__, PS_LOG_WARN, "WARNING: 2-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); - psFree(poly); - poly = NULL; - } else { - // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if (x->type.type != PS_TYPE_F64) { - psFree(x64); - } - - if (y->type.type != PS_TYPE_F64) { - psFree(y64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } + psError(PS_ERR_UNKNOWN, true, "2-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); + result = false; + break; + default: psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type. Returning NULL.\n"); - } - + result = false; + break; + } // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if (x->type.type != PS_TYPE_F64) { - psFree(x64); - } - - if (y->type.type != PS_TYPE_F64) { - psFree(y64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } - - return(poly); + PS_FREE_TEMP_F64_VECTOR (f, f64); + PS_FREE_TEMP_F64_VECTOR (x, x64); + PS_FREE_TEMP_F64_VECTOR (y, y64); + PS_FREE_TEMP_F64_VECTOR (fErr, fErr64); + + return result; } -psPolynomial2D *psVectorClipFitPolynomial2D( +bool psVectorClipFitPolynomial2D( psPolynomial2D *poly, psStats *stats, @@ -1212,25 +1144,40 @@ { psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); - PS_ASSERT_PTR_NON_NULL(stats, NULL); - PS_ASSERT_VECTOR_NON_NULL(mask, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); - if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_TYPE(x, f->type.type, NULL); - - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_TYPE(y, f->type.type, NULL); + PS_ASSERT_POLY_NON_NULL(poly, false); + PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false); + PS_ASSERT_PTR_NON_NULL(stats, false); + PS_ASSERT_VECTOR_NON_NULL(mask, false); + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false); + + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_TYPE(x, f->type.type, false); + + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false); + PS_ASSERT_VECTOR_TYPE(y, f->type.type, false); + + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false); + PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false); + } + + // the user supplies one of various stats option pairs, + // determine the desired mean and stdev STATS options: + // XXX enforce consistency? + // XXX psStatsGetValue() probably has inverted precedence + psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2); + psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2); + if (!meanOption) { + psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected"); + return false; + } + if (!stdevOption) { + psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected"); + return false; } @@ -1250,10 +1197,4 @@ psVector *resid = psVectorAlloc(f->n, PS_TYPE_F64); - // eventual expansion: user supplies one of various stats option pairs, - // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to - // evaluate the clipping sigma - // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used - stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter); psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma); @@ -1270,9 +1211,8 @@ } - poly = psVectorFitPolynomial2D(poly, mask, maskValue, f, fErr, x, y); - if (poly == NULL) { - psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the data. Returning NULL.\n"); + if (!psVectorFitPolynomial2D(poly, mask, maskValue, f, fErr, x, y)) { + psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the data. Returning false.\n"); psFree(resid) - return(NULL); + return false; } @@ -1281,5 +1221,5 @@ psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector(). Returning NULL.\n"); psFree(resid) - return(NULL); + return false; } @@ -1307,24 +1247,14 @@ psFree(resid) psFree(fit) - return(NULL); - } - # if (USE_ROBUST_STATS_FOR_CLIPPING) - psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian); - psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev); - psTrace("psLib.math", 5, "Sample Median is %f\n", stats->sampleMedian); - psTrace("psLib.math", 5, "Sample Stdev is %f\n", stats->sampleStdev); - psF32 minClipValue = -minClipSigma*stats->robustStdev; - psF32 maxClipValue = +maxClipSigma*stats->robustStdev; - psF32 clipMedian = stats->robustMedian; - # else - - psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian); - psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev); - psTrace("psLib.math", 5, "Robust Median is %f\n", stats->robustMedian); - psTrace("psLib.math", 5, "Robust Stdev is %f\n", stats->robustStdev); - psF32 minClipValue = -minClipSigma*stats->robustStdev; - psF32 maxClipValue = +maxClipSigma*stats->robustStdev; - psF32 clipMedian = stats->sampleMedian; - # endif + return false; + } + + double meanValue = psStatsGetValue (stats, meanOption); + double stdevValue = psStatsGetValue (stats, stdevOption); + + psTrace("psLib.math", 5, "Mean is %f\n", meanValue); + psTrace("psLib.math", 5, "Stdev is %f\n", stdevValue); + psF32 minClipValue = -minClipSigma*stdevValue; + psF32 maxClipValue = +maxClipSigma*stdevValue; // set mask if pts are not valid @@ -1336,5 +1266,5 @@ } - if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue)) { + if ((resid->data.F64[i] - meanValue > maxClipValue) || (resid->data.F64[i] - meanValue < minClipValue)) { if (fit->type.type == PS_TYPE_F64) { psTrace("psLib.math", 6, "Masking element %d (%f). resid->data.F64[%d] is %f\n", @@ -1352,6 +1282,6 @@ Nkeep++; } - psTrace("psLib.math", 4, "keeping %d of %ld pts for fit\n", Nkeep, x->n); + stats->clippedNvalues = Nkeep; psFree(fit); } @@ -1359,11 +1289,6 @@ psFree(resid); - if (poly == NULL) { - psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data. Returning NULL.\n"); - return(NULL); - } - psTrace("psLib.math", 3, "---- %s() end ----\n", __func__); - return(poly); + return true; } @@ -1381,5 +1306,5 @@ *****************************************************************************/ -static psPolynomial3D* VectorFitPolynomial3DOrd( +static bool VectorFitPolynomial3DOrd( psPolynomial3D* myPoly, const psVector* mask, @@ -1392,29 +1317,28 @@ { psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); - PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, NULL); - - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL); + PS_ASSERT_POLY_NON_NULL(myPoly, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, false); + + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, false); if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, false); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, false); + } + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false); + PS_ASSERT_VECTOR_NON_NULL(z, false); + PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, false); if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); + } int nXterm = 1 + myPoly->nX; // Number of x terms @@ -1429,9 +1353,8 @@ if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) { psError(PS_ERR_UNKNOWN, false, "Could initialize data structures A, B. Returning NULL.\n"); - psFree(myPoly); psFree(A); psFree(B); psTrace("psLib.math", 4, "---- %s() End ----\n", __func__); - return(NULL); + return false; } @@ -1511,8 +1434,6 @@ // The matrices were overflowing, so I switched to LUD. if (!psMatrixGJSolve(A, B)) { - psFree(A); - psFree(B); psError(PS_ERR_UNKNOWN, false, "Failed to perform GaussJordan elimination.\n"); - return(NULL); + goto escape_GJ; } // select the appropriate solution entries @@ -1534,12 +1455,10 @@ if (ALUD == NULL) { psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix. Returning NULL.\n"); - psFree(myPoly); - myPoly = NULL; + goto escape_LUD; } else { coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); if (coeffs == NULL) { psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix. Returning NULL.\n"); - psFree(myPoly); - myPoly = NULL; + goto escape_LUD; } else { // select the appropriate solution entries @@ -1565,5 +1484,16 @@ psTrace("psLib.math", 4, "---- %s() end ----\n", __func__); - return (myPoly); + return true; + +escape_LUD: + psFree(A); + psFree(B); + return false; + +escape_GJ: + + psFree(A); + psFree(B); + return false; } @@ -1574,5 +1504,5 @@ vector conversion only. *****************************************************************************/ -psPolynomial3D *psVectorFitPolynomial3D( +bool psVectorFitPolynomial3D( psPolynomial3D *poly, const psVector *mask, @@ -1584,147 +1514,68 @@ const psVector *z) { - // Internal pointers for possibly NULL or mis-typed vectors. - psVector *x64 = NULL; - psVector *y64 = NULL; - psVector *z64 = NULL; - psVector *f64 = NULL; + PS_ASSERT_POLY_NON_NULL(poly, false); + PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false); + + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false); + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false); + PS_ASSERT_VECTOR_NON_NULL(z, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, false); + if (mask != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); + } + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, false); + } + + // Convert input vectors to F64 if necessary. + psVector *f64 = (f->type.type == PS_TYPE_F64) ? (psVector *) f : psVectorCopy(NULL, f, PS_TYPE_F64); + psVector *x64 = (x->type.type == PS_TYPE_F64) ? (psVector *) x : psVectorCopy(NULL, x, PS_TYPE_F64); + psVector *y64 = (y->type.type == PS_TYPE_F64) ? (psVector *) y : psVectorCopy(NULL, y, PS_TYPE_F64); + psVector *z64 = (z->type.type == PS_TYPE_F64) ? (psVector *) z : psVectorCopy(NULL, z, PS_TYPE_F64); + psVector *fErr64 = NULL; - - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); - - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); - if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - // PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, NULL); - } - - // - // Convert input vectors to F64 if necessary. - // - if (f->type.type != PS_TYPE_F64) { - f64 = psVectorCopy(NULL, f, PS_TYPE_F64); - } else { - f64 = (psVector *) f; - } - if (x->type.type != PS_TYPE_F64) { - x64 = psVectorCopy(NULL, x, PS_TYPE_F64); - } else { - x64 = (psVector *) x; - } - if (y->type.type != PS_TYPE_F64) { - y64 = psVectorCopy(NULL, y, PS_TYPE_F64); - } else { - y64 = (psVector *) y; - } - - if (z->type.type != PS_TYPE_F64) { - z64 = psVectorCopy(NULL, z, PS_TYPE_F64); - } else { - z64 = (psVector *) z; - } - - if (fErr != NULL) { - if (fErr->type.type != PS_TYPE_F64) { - fErr64 = psVectorCopy(NULL, fErr, PS_TYPE_F64); - } else { - fErr64 = (psVector *) fErr; - } - } - - if (poly->type == PS_POLYNOMIAL_ORD) { - poly = VectorFitPolynomial3DOrd(poly, mask, maskValue, f64, fErr64, x64, y64, z64); - if (poly == NULL) { + fErr64 = (fErr->type.type == PS_TYPE_F64) ? (psVector *) fErr : psVectorCopy(NULL, fErr, PS_TYPE_F64); + } + + bool result = true; + + switch (poly->type) { + case PS_POLYNOMIAL_ORD: + result = VectorFitPolynomial3DOrd(poly, mask, maskValue, f64, fErr64, x64, y64, z64); + if (!result) { psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial. Returning NULL.\n"); - // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if (x->type.type != PS_TYPE_F64) { - psFree(x64); - } - - if (y->type.type != PS_TYPE_F64) { - psFree(y64); - } - - if (z->type.type != PS_TYPE_F64) { - psFree(z64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } - return(NULL); - } - } else if (poly->type == PS_POLYNOMIAL_CHEB) { + } + break; + case PS_POLYNOMIAL_CHEB: if (mask != NULL) { psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n"); } - psLogMsg(__func__, PS_LOG_WARN, "WARNING: 3-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); - psFree(poly); - poly = NULL; - } else { - // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if (x->type.type != PS_TYPE_F64) { - psFree(x64); - } - - if (y->type.type != PS_TYPE_F64) { - psFree(y64); - } - - if (z->type.type != PS_TYPE_F64) { - psFree(z64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } + psError(PS_ERR_UNKNOWN, true, "3-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); + result = false; + break; + default: psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type. Returning NULL.\n"); - } - + result = false; + break; + } // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if (x->type.type != PS_TYPE_F64) { - psFree(x64); - } - - if (y->type.type != PS_TYPE_F64) { - psFree(y64); - } - - if (z->type.type != PS_TYPE_F64) { - psFree(z64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } - - return(poly); + PS_FREE_TEMP_F64_VECTOR (f, f64); + PS_FREE_TEMP_F64_VECTOR (x, x64); + PS_FREE_TEMP_F64_VECTOR (y, y64); + PS_FREE_TEMP_F64_VECTOR (z, z64); + PS_FREE_TEMP_F64_VECTOR (fErr, fErr64); + + return result; } -psPolynomial3D *psVectorClipFitPolynomial3D( +bool psVectorClipFitPolynomial3D( psPolynomial3D *poly, psStats *stats, @@ -1738,29 +1589,44 @@ { psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); - PS_ASSERT_PTR_NON_NULL(stats, NULL); - PS_ASSERT_VECTOR_NON_NULL(mask, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); - if (mask != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_TYPE(x, f->type.type, NULL); - - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_TYPE(y, f->type.type, NULL); - - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); - PS_ASSERT_VECTOR_TYPE(z, f->type.type, NULL); + PS_ASSERT_POLY_NON_NULL(poly, false); + PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false); + PS_ASSERT_PTR_NON_NULL(stats, false); + PS_ASSERT_VECTOR_NON_NULL(mask, false); + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false); + + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_TYPE(x, f->type.type, false); + + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false); + PS_ASSERT_VECTOR_TYPE(y, f->type.type, false); + + PS_ASSERT_VECTOR_NON_NULL(z, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, false); + PS_ASSERT_VECTOR_TYPE(z, f->type.type, false); + + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false); + PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false); + } + + // the user supplies one of various stats option pairs, + // determine the desired mean and stdev STATS options: + // XXX enforce consistency? + // XXX psStatsGetValue() probably has inverted precedence + psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2); + psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2); + if (!meanOption) { + psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected"); + return false; + } + if (!stdevOption) { + psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected"); + return false; } @@ -1778,13 +1644,6 @@ minClipSigma = fabs(stats->clipSigma); } - psVector *fit = NULL; psVector *resid = psVectorAlloc(f->n, PS_TYPE_F64); - // eventual expansion: user supplies one of various stats option pairs, - // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to - // evaluate the clipping sigma - // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used - stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter); psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma); @@ -1801,16 +1660,14 @@ } - poly = psVectorFitPolynomial3D(poly, mask, maskValue, f, fErr, x, y, z); - if (poly == NULL) { + if (!psVectorFitPolynomial3D(poly, mask, maskValue, f, fErr, x, y, z)) { psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the data. Returning NULL.\n"); psFree(resid) - psFree(fit) - return(NULL); - } - fit = psPolynomial3DEvalVector(poly, x, y, z); + return false; + } + psVector *fit = psPolynomial3DEvalVector(poly, x, y, z); if (fit == NULL) { psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial3DEvalVector(). Returning NULL.\n"); psFree(resid) - return(NULL); + return false; } for (psS32 i = 0 ; i < f->n ; i++) { @@ -1837,21 +1694,14 @@ psFree(resid) psFree(fit) - return(NULL); - } - - # if (USE_ROBUST_STATS_FOR_CLIPPING) - psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian); - psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev); - psF32 minClipValue = -minClipSigma*stats->robustStdev; - psF32 maxClipValue = +maxClipSigma*stats->robustStdev; - psF32 clipMedian = stats->robustMedian; - # else - - psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian); - psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev); - psF32 minClipValue = -minClipSigma*stats->robustStdev; - psF32 maxClipValue = +maxClipSigma*stats->robustStdev; - psF32 clipMedian = stats->sampleMedian; - # endif + return false; + } + + double meanValue = psStatsGetValue (stats, meanOption); + double stdevValue = psStatsGetValue (stats, stdevOption); + + psTrace("psLib.math", 5, "Mean is %f\n", meanValue); + psTrace("psLib.math", 5, "Stdev is %f\n", stdevValue); + psF32 minClipValue = -minClipSigma*stdevValue; + psF32 maxClipValue = +maxClipSigma*stdevValue; // set mask if pts are not valid @@ -1863,5 +1713,5 @@ } - if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue)) { + if ((resid->data.F64[i] - meanValue > maxClipValue) || (resid->data.F64[i] - meanValue < minClipValue)) { if (f->type.type == PS_TYPE_F64) { psTrace("psLib.math", 6, "Masking element %d (%f). resid->data.F64[%d] is %f\n", @@ -1879,6 +1729,6 @@ Nkeep++; } - psTrace("psLib.math", 6, "keeping %d of %ld pts for fit\n", Nkeep, x->n); + stats->clippedNvalues = Nkeep; psFree(fit); } @@ -1886,11 +1736,6 @@ psFree(resid); - if (poly == NULL) { - psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data. Returning NULL.\n"); - return(NULL); - } - psTrace("psLib.math", 3, "---- %s() end ----\n", __func__); - return(poly); + return true; } @@ -1906,5 +1751,5 @@ *****************************************************************************/ -static psPolynomial4D* VectorFitPolynomial4DOrd( +static bool VectorFitPolynomial4DOrd( psPolynomial4D* myPoly, const psVector* mask, @@ -1918,30 +1763,30 @@ { psTrace("psLib.math", 4, "---- %s() begin ----\n", __func__); - PS_ASSERT_POLY_NON_NULL(myPoly, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, NULL); - PS_ASSERT_INT_NONNEGATIVE(myPoly->nT, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, NULL); + PS_ASSERT_POLY_NON_NULL(myPoly, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nX, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nY, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nZ, false); + PS_ASSERT_INT_NONNEGATIVE(myPoly->nT, false); + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE(f, PS_TYPE_F64, false); if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); - PS_ASSERT_VECTOR_NON_NULL(t, NULL); - PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F64, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(y, fErr, false); + PS_ASSERT_VECTOR_TYPE(fErr, PS_TYPE_F64, false); + } + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false); + PS_ASSERT_VECTOR_NON_NULL(z, false); + PS_ASSERT_VECTOR_TYPE(z, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, false); + PS_ASSERT_VECTOR_NON_NULL(t, false); + PS_ASSERT_VECTOR_TYPE(t, PS_TYPE_F64, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, false); if (mask) { - PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(y, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); } @@ -1959,9 +1804,8 @@ if (!psImageInit(A, 0.0) || !psVectorInit(B, 0.0)) { psError(PS_ERR_UNKNOWN, false, "Could initialize data structures A, B. Returning NULL.\n"); - psFree(myPoly); psFree(A); psFree(B); psTrace("psLib.math", 4, "---- %s() End ----\n", __func__); - return(NULL); + return false; } @@ -2051,8 +1895,6 @@ // The GaussJordan version was overflowing, so I'm using LUD. if (!psMatrixGJSolve(A, B)) { - psFree(A); - psFree(B); psError(PS_ERR_UNKNOWN, false, "Failed to perform GaussJordan elimination.\n"); - return(NULL); + goto escape_GJ; } @@ -2076,12 +1918,10 @@ if (ALUD == NULL) { psError(PS_ERR_UNKNOWN, false, "Could not do LUD decomposition on matrix. Returning NULL.\n"); - psFree(myPoly); - myPoly = NULL; + goto escape_LUD; } else { coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); if (coeffs == NULL) { psError(PS_ERR_UNKNOWN, false, "Could not solve LUD matrix. Returning NULL.\n"); - psFree(myPoly); - myPoly = NULL; + goto escape_LUD; } else { // select the appropriate solution entries @@ -2100,16 +1940,23 @@ } } - psFree(ALUD); psFree(coeffs); psFree(outPerm); - - } - + } psFree(A); psFree(B); psTrace("psLib.math", 4, "---- %s() end ----\n", __func__); - return (myPoly); + return true; + +escape_LUD: + psFree(A); + psFree(B); + return false; + +escape_GJ: + psFree(A); + psFree(B); + return false; } @@ -2120,5 +1967,5 @@ via vector conversion only. *****************************************************************************/ -psPolynomial4D *psVectorFitPolynomial4D( +bool psVectorFitPolynomial4D( psPolynomial4D *poly, const psVector *mask, @@ -2131,170 +1978,73 @@ const psVector *t) { - // Internal pointers for possibly NULL or mis-typed vectors. - psVector *x64 = NULL; - psVector *y64 = NULL; - psVector *z64 = NULL; - psVector *t64 = NULL; - psVector *f64 = NULL; + PS_ASSERT_POLY_NON_NULL(poly, false); + PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false); + + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false); + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false); + PS_ASSERT_VECTOR_NON_NULL(z, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, false); + PS_ASSERT_VECTOR_NON_NULL(t, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, false); + if (mask) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); + } + if (fErr != NULL) { + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, false); + } + + // Convert input vectors to F64 if necessary. + psVector *f64 = (f->type.type == PS_TYPE_F64) ? (psVector *) f : psVectorCopy(NULL, f, PS_TYPE_F64); + psVector *x64 = (x->type.type == PS_TYPE_F64) ? (psVector *) x : psVectorCopy(NULL, x, PS_TYPE_F64); + psVector *y64 = (y->type.type == PS_TYPE_F64) ? (psVector *) y : psVectorCopy(NULL, y, PS_TYPE_F64); + psVector *z64 = (z->type.type == PS_TYPE_F64) ? (psVector *) z : psVectorCopy(NULL, z, PS_TYPE_F64); + psVector *t64 = (t->type.type == PS_TYPE_F64) ? (psVector *) t : psVectorCopy(NULL, t, PS_TYPE_F64); + psVector *fErr64 = NULL; - - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); - - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); - if (mask) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); - PS_ASSERT_VECTOR_NON_NULL(t, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL); if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(fErr, NULL); - } - - - // - // Convert input vector to F64 if necessary. - // - if (f->type.type != PS_TYPE_F64) { - f64 = psVectorCopy(NULL, f, PS_TYPE_F64); - } else { - f64 = (psVector *) f; - } - if (x->type.type != PS_TYPE_F64) { - x64 = psVectorCopy(NULL, x, PS_TYPE_F64); - } else { - x64 = (psVector *) x; - } - if (y->type.type != PS_TYPE_F64) { - y64 = psVectorCopy(NULL, y, PS_TYPE_F64); - } else { - y64 = (psVector *) y; - } - if (z->type.type != PS_TYPE_F64) { - z64 = psVectorCopy(NULL, z, PS_TYPE_F64); - } else { - z64 = (psVector *) z; - } - if (t->type.type != PS_TYPE_F64) { - t64 = psVectorCopy(NULL, t, PS_TYPE_F64); - } else { - t64 = (psVector *) t; - } - // - // fErr - // - if (fErr != NULL) { - if (fErr->type.type != PS_TYPE_F64) { - fErr64 = psVectorCopy(NULL, fErr, PS_TYPE_F64); - } else { - fErr64 = (psVector *) fErr; - } - } - - if (poly->type == PS_POLYNOMIAL_ORD) { - poly = VectorFitPolynomial4DOrd(poly, mask, maskValue, f64, fErr64, x64, y64, z64, t64); - if (poly == NULL) { + fErr64 = (fErr->type.type == PS_TYPE_F64) ? (psVector *) fErr : psVectorCopy(NULL, fErr, PS_TYPE_F64); + } + + bool result = true; + + switch (poly->type) { + case PS_POLYNOMIAL_ORD: + result = VectorFitPolynomial4DOrd(poly, mask, maskValue, f64, fErr64, x64, y64, z64, t64); + if (!result) { psError(PS_ERR_UNKNOWN, true, "Could not fit polynomial. Returning NULL.\n"); - // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if (x->type.type != PS_TYPE_F64) { - psFree(x64); - } - - if (y->type.type != PS_TYPE_F64) { - psFree(y64); - } - - if (z->type.type != PS_TYPE_F64) { - psFree(z64); - } - - if (t->type.type != PS_TYPE_F64) { - psFree(t64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } - return(NULL); - } - } else if (poly->type == PS_POLYNOMIAL_CHEB) { + } + break; + case PS_POLYNOMIAL_CHEB: if (mask != NULL) { psLogMsg(__func__, PS_LOG_WARN, "WARNING: ignoring mask and maskValue with Chebyshev polynomials.\n"); } - psLogMsg(__func__, PS_LOG_WARN, "WARNING: 4-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); - psFree(poly); - poly = NULL; - } else { - // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if (x->type.type != PS_TYPE_F64) { - psFree(x64); - } - - if (y->type.type != PS_TYPE_F64) { - psFree(y64); - } - - if (z->type.type != PS_TYPE_F64) { - psFree(z64); - } - - if (t->type.type != PS_TYPE_F64) { - psFree(t64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } + psError(PS_ERR_UNKNOWN, true, "4-D Chebyshev polynomial vector fitting has not been implemented. Returning NULL.\n"); + result = false; + break; + default: psError(PS_ERR_UNKNOWN, true, "Incorrect polynomial type. Returning NULL.\n"); - } - + result = false; + break; + } // Free psVectors that were created for NULL arguments. - if (f->type.type != PS_TYPE_F64) { - psFree(f64); - } - - if (x->type.type != PS_TYPE_F64) { - psFree(x64); - } - - if (y->type.type != PS_TYPE_F64) { - psFree(y64); - } - - if (z->type.type != PS_TYPE_F64) { - psFree(z64); - } - - if (t->type.type != PS_TYPE_F64) { - psFree(t64); - } - - if ((fErr != NULL) && (fErr->type.type != PS_TYPE_F64)) { - psFree(fErr64); - } - - return(poly); + PS_FREE_TEMP_F64_VECTOR (f, f64); + PS_FREE_TEMP_F64_VECTOR (x, x64); + PS_FREE_TEMP_F64_VECTOR (y, y64); + PS_FREE_TEMP_F64_VECTOR (z, z64); + PS_FREE_TEMP_F64_VECTOR (t, t64); + PS_FREE_TEMP_F64_VECTOR (fErr, fErr64); + + return result; } -psPolynomial4D *psVectorClipFitPolynomial4D( +bool psVectorClipFitPolynomial4D( psPolynomial4D *poly, psStats *stats, @@ -2309,33 +2059,48 @@ { psTrace("psLib.math", 3, "---- %s() begin ----\n", __func__); - PS_ASSERT_POLY_NON_NULL(poly, NULL); - PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, NULL); - PS_ASSERT_PTR_NON_NULL(stats, NULL); - PS_ASSERT_VECTOR_NON_NULL(mask, NULL); - PS_ASSERT_VECTOR_NON_NULL(f, NULL); - PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, NULL); - if (mask) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, NULL); - PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, NULL); - } - PS_ASSERT_VECTOR_NON_NULL(x, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, NULL); - PS_ASSERT_VECTOR_TYPE(x, f->type.type, NULL); - - PS_ASSERT_VECTOR_NON_NULL(y, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, NULL); - PS_ASSERT_VECTOR_TYPE(y, f->type.type, NULL); - - PS_ASSERT_VECTOR_NON_NULL(z, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, NULL); - PS_ASSERT_VECTOR_TYPE(z, f->type.type, NULL); - - PS_ASSERT_VECTOR_NON_NULL(t, NULL); - PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, NULL); - PS_ASSERT_VECTOR_TYPE(t, f->type.type, NULL); + PS_ASSERT_POLY_NON_NULL(poly, false); + PS_ASSERT_POLY_TYPE(poly, PS_POLYNOMIAL_ORD, false); + PS_ASSERT_PTR_NON_NULL(stats, false); + PS_ASSERT_VECTOR_NON_NULL(mask, false); + PS_ASSERT_VECTOR_NON_NULL(f, false); + PS_ASSERT_VECTOR_TYPE_F32_OR_F64(f, false); + + PS_ASSERT_VECTOR_NON_NULL(x, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, x, false); + PS_ASSERT_VECTOR_TYPE(x, f->type.type, false); + + PS_ASSERT_VECTOR_NON_NULL(y, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, y, false); + PS_ASSERT_VECTOR_TYPE(y, f->type.type, false); + + PS_ASSERT_VECTOR_NON_NULL(z, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, z, false); + PS_ASSERT_VECTOR_TYPE(z, f->type.type, false); + + PS_ASSERT_VECTOR_NON_NULL(t, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, t, false); + PS_ASSERT_VECTOR_TYPE(t, f->type.type, false); + + PS_ASSERT_VECTORS_SIZE_EQUAL(f, mask, false); + PS_ASSERT_VECTOR_TYPE(mask, PS_TYPE_U8, false); if (fErr != NULL) { - PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, NULL); - PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(f, fErr, false); + PS_ASSERT_VECTOR_TYPE(fErr, f->type.type, false); + } + + // the user supplies one of various stats option pairs, + // determine the desired mean and stdev STATS options: + // XXX enforce consistency? + // XXX psStatsGetValue() probably has inverted precedence + psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2); + psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2); + if (!meanOption) { + psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected"); + return false; + } + if (!stdevOption) { + psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected"); + return false; } @@ -2353,13 +2118,6 @@ minClipSigma = fabs(stats->clipSigma); } - psVector *fit = NULL; psVector *resid = psVectorAlloc(f->n, PS_TYPE_F64); - // eventual expansion: user supplies one of various stats option pairs, - // eg (SAMPLE_MEAN | SAMPLE_STDEV) and the correct pair is used to - // evaluate the clipping sigma - // for now, for the SAMPLE_MEDIAN and SAMPLE_STDEV to be used - stats->options |= (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); - stats->options |= (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); psTrace("psLib.math", 4, "stats->clipIter is %d\n", stats->clipIter); psTrace("psLib.math", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma); @@ -2376,17 +2134,15 @@ } - poly = psVectorFitPolynomial4D (poly, mask, maskValue, f, fErr, x, y, z, t); - if (poly == NULL) { + if (!psVectorFitPolynomial4D (poly, mask, maskValue, f, fErr, x, y, z, t)) { psError(PS_ERR_UNKNOWN, false, "Could not fit a polynomial to the data. Returning NULL.\n"); psFree(resid) - psFree(fit) - return(NULL); - } - - fit = psPolynomial4DEvalVector (poly, x, y, z, t); + return false; + } + + psVector *fit = psPolynomial4DEvalVector (poly, x, y, z, t); if (fit == NULL) { psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial4DEvalVector(). Returning NULL.\n"); psFree(resid) - return(NULL); + return false; } for (psS32 i = 0 ; i < f->n ; i++) { @@ -2413,20 +2169,14 @@ psFree(resid) psFree(fit) - return(NULL); - } - # if (USE_ROBUST_STATS_FOR_CLIPPING) - psTrace("psLib.math", 5, "Median is %f\n", stats->robustMedian); - psTrace("psLib.math", 5, "Stdev is %f\n", stats->robustStdev); - psF32 minClipValue = -minClipSigma*stats->robustStdev; - psF32 maxClipValue = +maxClipSigma*stats->robustStdev; - psF32 clipMedian = stats->robustMedian; - # else - - psTrace("psLib.math", 5, "Median is %f\n", stats->sampleMedian); - psTrace("psLib.math", 5, "Stdev is %f\n", stats->sampleStdev); - psF32 minClipValue = -minClipSigma*stats->robustStdev; - psF32 maxClipValue = +maxClipSigma*stats->robustStdev; - psF32 clipMedian = stats->sampleMedian; - # endif + return false; + } + + double meanValue = psStatsGetValue (stats, meanOption); + double stdevValue = psStatsGetValue (stats, stdevOption); + + psTrace("psLib.math", 5, "Mean is %f\n", meanValue); + psTrace("psLib.math", 5, "Stdev is %f\n", stdevValue); + psF32 minClipValue = -minClipSigma*stdevValue; + psF32 maxClipValue = +maxClipSigma*stdevValue; // set mask if pts are not valid @@ -2438,5 +2188,5 @@ } - if ((resid->data.F64[i] - clipMedian > maxClipValue) || (resid->data.F64[i] - clipMedian < minClipValue)) { + if ((resid->data.F64[i] - meanValue > maxClipValue) || (resid->data.F64[i] - meanValue < minClipValue)) { if (f->type.type == PS_TYPE_F64) { psTrace("psLib.math", 6, "Masking element %d (%f). resid->data.F64[%d] is %f\n", @@ -2452,9 +2202,8 @@ continue; } - Nkeep++; } - psTrace("psLib.math", 6, "keeping %d of %ld pts for fit\n", Nkeep, x->n); + stats->clippedNvalues = Nkeep; psFree (fit); } @@ -2462,10 +2211,5 @@ psFree (resid); - if (poly == NULL) { - psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data. Returning NULL.\n"); - return(NULL); - } - psTrace("psLib.math", 3, "---- %s() end ----\n", __func__); - return(poly); + return true; } Index: trunk/psLib/src/math/psMinimizePolyFit.h =================================================================== --- trunk/psLib/src/math/psMinimizePolyFit.h (revision 10831) +++ trunk/psLib/src/math/psMinimizePolyFit.h (revision 10848) @@ -10,6 +10,6 @@ * XXX: Must Doxygenate. * - * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-01-23 20:44:29 $ + * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-12-29 04:38:42 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -54,5 +54,5 @@ */ -psPolynomial1D *psVectorFitPolynomial1D( +bool psVectorFitPolynomial1D( psPolynomial1D *poly, const psVector *mask, @@ -63,5 +63,5 @@ ); -psPolynomial2D *psVectorFitPolynomial2D( +bool psVectorFitPolynomial2D( psPolynomial2D *poly, const psVector *mask, @@ -73,5 +73,5 @@ ); -psPolynomial3D *psVectorFitPolynomial3D( +bool psVectorFitPolynomial3D( psPolynomial3D *poly, const psVector *mask, @@ -84,5 +84,5 @@ ); -psPolynomial4D *psVectorFitPolynomial4D( +bool psVectorFitPolynomial4D( psPolynomial4D *poly, const psVector *mask, @@ -97,5 +97,5 @@ -psPolynomial1D *psVectorClipFitPolynomial1D( +bool psVectorClipFitPolynomial1D( psPolynomial1D *poly, psStats *stats, @@ -107,5 +107,5 @@ ); -psPolynomial2D *psVectorClipFitPolynomial2D( +bool psVectorClipFitPolynomial2D( psPolynomial2D *poly, psStats *stats, @@ -118,5 +118,5 @@ ); -psPolynomial3D *psVectorClipFitPolynomial3D( +bool psVectorClipFitPolynomial3D( psPolynomial3D *poly, psStats *stats, @@ -130,5 +130,5 @@ ); -psPolynomial4D *psVectorClipFitPolynomial4D( +bool psVectorClipFitPolynomial4D( psPolynomial4D *poly, psStats *stats, Index: trunk/psLib/src/math/psPolynomial.c =================================================================== --- trunk/psLib/src/math/psPolynomial.c (revision 10831) +++ trunk/psLib/src/math/psPolynomial.c (revision 10848) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * -* @version $Revision: 1.155 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-12-10 17:28:46 $ +* @version $Revision: 1.156 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-12-29 04:38:42 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -581,5 +581,6 @@ unsigned int nX) { - //PS_ASSERT_INT_NONNEGATIVE(n, NULL); + PS_ASSERT_POLY_VALID_TYPE(type, NULL); + psU32 nOrder = nX; psPolynomial1D *newPoly = (psPolynomial1D* ) psAlloc(sizeof(psPolynomial1D)); @@ -605,6 +606,5 @@ unsigned int nY) { - //PS_ASSERT_INT_NONNEGATIVE(nX, NULL); - //PS_ASSERT_INT_NONNEGATIVE(nY, NULL); + PS_ASSERT_POLY_VALID_TYPE(type, NULL); unsigned int x = 0; @@ -711,4 +711,6 @@ unsigned int nZ) { + PS_ASSERT_POLY_VALID_TYPE(type, NULL); + //PS_ASSERT_INT_NONNEGATIVE(nX, NULL); //PS_ASSERT_INT_NONNEGATIVE(nY, NULL); @@ -761,4 +763,6 @@ unsigned int nT) { + PS_ASSERT_POLY_VALID_TYPE(type, NULL); + //PS_ASSERT_INT_NONNEGATIVE(nX, NULL); //PS_ASSERT_INT_NONNEGATIVE(nY, NULL); Index: trunk/psLib/src/math/psPolynomial.h =================================================================== --- trunk/psLib/src/math/psPolynomial.h (revision 10831) +++ trunk/psLib/src/math/psPolynomial.h (revision 10848) @@ -11,6 +11,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.65 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-12-09 14:19:56 $ + * @version $Revision: 1.66 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-12-29 04:38:42 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -329,4 +329,12 @@ } \ +#define PS_ASSERT_POLY_VALID_TYPE(TYPE, RVAL) \ +if ((TYPE != PS_POLYNOMIAL_ORD) && \ + (TYPE != PS_POLYNOMIAL_CHEB)) { \ + psError(PS_ERR_BAD_PARAMETER_TYPE, true, \ + "Unallowable operation: invalid type %d for polynomial", TYPE); \ + return(RVAL); \ +} \ + #define PS_POLY_PRINT_1D(NAME) \ printf("Poly %s: (nX) is (%d)\n", #NAME, NAME->nX);\ Index: trunk/psLib/src/math/psPolynomialUtils.c =================================================================== --- trunk/psLib/src/math/psPolynomialUtils.c (revision 10831) +++ trunk/psLib/src/math/psPolynomialUtils.c (revision 10848) @@ -13,5 +13,5 @@ #include "psPolynomialUtils.h" -psPolynomial4D *psVectorChiClipFitPolynomial4D( +bool psVectorChiClipFitPolynomial4D( psPolynomial4D *poly, psStats *stats, @@ -75,11 +75,24 @@ int Nkeep = 0; - poly = psVectorFitPolynomial4D (poly, mask, maskValue, f, fErr, x, y, z, t); + if (!psVectorFitPolynomial4D (poly, mask, maskValue, f, fErr, x, y, z, t)) { + psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data. Returning NULL.\n"); + psFree (resid); + return false; + } + fit = psPolynomial4DEvalVector (poly, x, y, z, t); + if (fit == NULL) { + psError(PS_ERR_UNKNOWN, false, "Could not call psPolynomial4DEvalVector(). Returning NULL.\n"); + psFree(resid) + return false; + } + resid = (psVector *) psBinaryOp (resid, (void *) f, "-", (void *) fit); if (!psVectorStats (stats, resid, NULL, mask, maskValue)) { psError(PS_ERR_UNKNOWN, false, "failed to measure vector stats"); - return NULL; + psFree (fit); + psFree (resid); + return false; } psTrace (__func__, 5, "resid stats: %f +/- %f\n", stats->sampleMedian, stats->sampleStdev); @@ -108,7 +121,5 @@ } - psTrace (__func__, 4, "keeping %d of %ld pts for fit\n", - Nkeep, x->n); - + psTrace (__func__, 4, "keeping %d of %ld pts for fit\n", Nkeep, x->n); stats->clippedNvalues = Nkeep; psFree (fit); @@ -117,9 +128,5 @@ psFree (resid); - if (poly == NULL) { - psError(PS_ERR_UNKNOWN, true, "Could not fit a polynomial to the data. Returning NULL.\n"); - return(NULL); - } - return(poly); + return true; } Index: trunk/psLib/src/math/psPolynomialUtils.h =================================================================== --- trunk/psLib/src/math/psPolynomialUtils.h (revision 10831) +++ trunk/psLib/src/math/psPolynomialUtils.h (revision 10848) @@ -14,5 +14,5 @@ // perform vector clip-fit based on significance of deviations -psPolynomial4D *psVectorChiClipFitPolynomial4D( +bool psVectorChiClipFitPolynomial4D( psPolynomial4D *poly, // Polynomial to fit psStats *stats, // Statistics to use in clipping Index: trunk/psLib/src/math/psStats.c =================================================================== --- trunk/psLib/src/math/psStats.c (revision 10831) +++ trunk/psLib/src/math/psStats.c (revision 10848) @@ -13,6 +13,6 @@ * use ->min and ->max (PS_STAT_USE_RANGE) * - * @version $Revision: 1.196 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-12-16 05:33:48 $ + * @version $Revision: 1.197 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-12-29 04:38:42 $ * * Copyright 2006 IfA, University of Hawaii @@ -1856,6 +1856,5 @@ // Determine the coefficients of the polynomial. psPolynomial1D *myPoly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD, 2); - myPoly = psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x); - if (myPoly == NULL) { + if (!psVectorFitPolynomial1D(myPoly, NULL, 0, y, NULL, x)) { psError(PS_ERR_UNEXPECTED_NULL, false, _("Failed to fit a 1-dimensional polynomial to the three specified data points. Returning NAN.")); Index: trunk/psLib/test/math/.cvsignore =================================================================== --- trunk/psLib/test/math/.cvsignore (revision 10831) +++ trunk/psLib/test/math/.cvsignore (revision 10848) @@ -4,52 +4,43 @@ Makefile Makefile.in -tst_psFunc01 -tst_psHist00 -tst_psHist01 -tst_psHist02 -tst_psHist03 -tst_psMatrix01 -tst_psMatrix02 -tst_psMatrix03 -tst_psMatrix04 -tst_psMatrix05 -tst_psMatrix06 -tst_psMatrix07 -tst_psMatrixVectorArithmetic01 -tst_psMatrixVectorArithmetic02 -tst_psMatrixVectorArithmetic03 -tst_psMatrixVectorArithmetic04 -tst_psRandom -tst_psStats00 -tst_psStats01 -tst_psStats02 -tst_psStats03 -tst_psStats05 -tst_psStats06 -tst_psStats07 -tst_psStats08 -tst_psStats09 seed_msglog1.txt seed_msglog2.txt -tst_psSpline1D -tst_psMathUtils -tst_psMinimizeLMM -tst_psMinimizePowell -tst_psPolyFit1D -tst_psPolyFit2D -tst_psPolyFit3D -tst_psPolyFit4D -tst_psPolynomial -tst_psPolynomialEval1D -tst_psPolynomialEval2D -tst_psPolynomialEval3D -tst_psPolynomialEval4D *.bb *.bbg *.da gmon.out + +tap_psHist00 +tap_psHist01 +tap_psHist02 +tap_psHist03 +tap_psMD5 +tap_psMatrix01 +tap_psMatrix02 +tap_psMatrix03 +tap_psMatrix04 +tap_psMatrix05 +tap_psMatrix06 +tap_psMatrix07 +tap_psPolyFit1D +tap_psPolyFit2D +tap_psPolyFit3D +tap_psPolyFit4D +tap_psPolynomial +tap_psPolynomialEval1D +tap_psPolynomialEval2D +tap_psPolynomialEval3D +tap_psPolynomialEval4D +tap_psPolynomialUtils_Derivatives tap_psSparse +tap_psStats00 +tap_psStats01 +tap_psStats02 +tap_psStats03 +tap_psStats05 +tap_psStats06 +tap_psStats07 +tap_psStats08 +tap_psStats09 tap_psStatsTiming -tap_psMD5 tap_psStats_Sample_01 -tap_psPolynomialUtils_Derivatives Index: trunk/psLib/test/math/Makefile.am =================================================================== --- trunk/psLib/test/math/Makefile.am (revision 10831) +++ trunk/psLib/test/math/Makefile.am (revision 10848) @@ -13,9 +13,38 @@ TEST_PROGS = \ + tap_psHist00 \ + tap_psHist01 \ + tap_psHist02 \ + tap_psHist03 \ tap_psMD5 \ + tap_psMatrix01 \ + tap_psMatrix02 \ + tap_psMatrix03 \ + tap_psMatrix04 \ + tap_psMatrix05 \ + tap_psMatrix06 \ + tap_psMatrix07 \ + tap_psPolyFit1D \ + tap_psPolyFit2D \ + tap_psPolyFit3D \ + tap_psPolyFit4D \ + tap_psPolynomial \ + tap_psPolynomialEval1D \ + tap_psPolynomialEval2D \ + tap_psPolynomialEval3D \ + tap_psPolynomialEval4D \ + tap_psPolynomialUtils_Derivatives \ tap_psSparse \ + tap_psStats00 \ + tap_psStats01 \ + tap_psStats02 \ + tap_psStats03 \ + tap_psStats05 \ + tap_psStats06 \ + tap_psStats07 \ + tap_psStats08 \ + tap_psStats09 \ tap_psStatsTiming \ - tap_psStats_Sample_01 \ - tap_psPolynomialUtils_Derivatives + tap_psStats_Sample_01 if BUILD_TESTS @@ -29,2 +58,3 @@ test: check + $(top_srcdir)/test/test.pl Index: trunk/psLib/test/math/tap_psPolyFit1D.c =================================================================== --- trunk/psLib/test/math/tap_psPolyFit1D.c (revision 10831) +++ trunk/psLib/test/math/tap_psPolyFit1D.c (revision 10848) @@ -272,5 +272,5 @@ } - psPolynomial1D *rc = NULL; + bool rc = false; if (flags & TS00_CLIP_FIT) { rc = psVectorClipFitPolynomial1D(myPoly, stats, mask, MASK_VALUE, f, fErr, x); @@ -279,5 +279,5 @@ } - if (rc == NULL) { + if (!rc) { if (expectedRC == true) { diag("TEST ERROR: the 1D polynomial fitting function returned NULL.\n"); Index: trunk/psLib/test/math/tap_psPolyFit2D.c =================================================================== --- trunk/psLib/test/math/tap_psPolyFit2D.c (revision 10831) +++ trunk/psLib/test/math/tap_psPolyFit2D.c (revision 10848) @@ -279,5 +279,5 @@ } - psPolynomial2D *rc = NULL; + bool rc = false; if (flags & TS00_CLIP_FIT) { rc = psVectorClipFitPolynomial2D(myPoly, stats, mask, MASK_VALUE, f, fErr, x, y); @@ -286,5 +286,5 @@ } - if (rc == NULL) { + if (!rc) { if (expectedRC == true) { diag("TEST ERROR: the 2D polynomial fitting function returned NULL.\n"); Index: trunk/psLib/test/math/tap_psPolyFit3D.c =================================================================== --- trunk/psLib/test/math/tap_psPolyFit3D.c (revision 10831) +++ trunk/psLib/test/math/tap_psPolyFit3D.c (revision 10848) @@ -324,5 +324,5 @@ } - psPolynomial3D *rc = NULL; + bool rc = false; if (flags & TS00_CLIP_FIT) { rc = psVectorClipFitPolynomial3D(myPoly, stats, mask, MASK_VALUE, f, fErr, x, y, z); @@ -331,5 +331,5 @@ } - if (rc == NULL) { + if (!rc) { if (expectedRC == true) { diag("TEST ERROR: the 3D polynomial fitting function returned NULL.\n"); Index: trunk/psLib/test/math/tap_psPolyFit4D.c =================================================================== --- trunk/psLib/test/math/tap_psPolyFit4D.c (revision 10831) +++ trunk/psLib/test/math/tap_psPolyFit4D.c (revision 10848) @@ -377,5 +377,5 @@ } - psPolynomial4D *rc = NULL; + bool rc = false; if (flags & TS00_CLIP_FIT) { rc = psVectorClipFitPolynomial4D(myPoly, stats, mask, MASK_VALUE, f, fErr, x, y, z, t); @@ -384,5 +384,5 @@ } - if (rc == NULL) { + if (!rc) { if (expectedRC == true) { diag("TEST ERROR: the 4D polynomial fitting function returned NULL.\n"); Index: trunk/psLib/test/math/tap_psPolynomialEval2D.c =================================================================== --- trunk/psLib/test/math/tap_psPolynomialEval2D.c (revision 10831) +++ trunk/psLib/test/math/tap_psPolynomialEval2D.c (revision 10848) @@ -4,6 +4,6 @@ * ORD and CHEB type polynomials. * -* @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-12-21 20:05:12 $ +* @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-12-29 04:38:42 $ * * Copyright 2004-2005 Maui High Performance Computing Center, Univ. of Hawaii @@ -148,6 +148,6 @@ psMemId id = psMemGetId(); psPolynomial2D* polyOrd = psPolynomial2DAlloc(99, TERMS-1, TERMS-1); - ok(polyOrd != NULL, "Ordinary polynomial allocation successful"); - skip_start(polyOrd == NULL, 2, "Skipping tests because psPolynomial2DAlloc() failed"); + ok(polyOrd == NULL, "Ordinary polynomial allocation successful"); + skip_start(polyOrd == NULL, 1, "Skipping tests because psPolynomial2DAlloc() failed"); // Attempt to evaluation invalid polynomial type Index: trunk/psLib/test/math/tap_psPolynomialEval3D.c =================================================================== --- trunk/psLib/test/math/tap_psPolynomialEval3D.c (revision 10831) +++ trunk/psLib/test/math/tap_psPolynomialEval3D.c (revision 10848) @@ -4,6 +4,6 @@ * ORD and CHEB type polynomials. * -* @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-12-21 20:05:12 $ +* @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-12-29 04:38:42 $ * * Copyright 2004-2005 Maui High Performance Computing Center, Univ. of Hawaii @@ -207,6 +207,6 @@ psMemId id = psMemGetId(); psPolynomial3D* polyOrd = psPolynomial3DAlloc(99, TERMS-1, TERMS-1, TERMS-1); - ok(polyOrd != NULL, "Ordinary polynomial allocation successful"); - skip_start(polyOrd == NULL, 2, "Skipping tests because psPolynomial3DAlloc() failed"); + ok(polyOrd == NULL, "Ordinary polynomial allocation successful"); + skip_start(polyOrd == NULL, 1, "Skipping tests because psPolynomial3DAlloc() failed"); // Attempt to evaluation invalid polynomial type Index: trunk/psLib/test/math/tap_psPolynomialEval4D.c =================================================================== --- trunk/psLib/test/math/tap_psPolynomialEval4D.c (revision 10831) +++ trunk/psLib/test/math/tap_psPolynomialEval4D.c (revision 10848) @@ -4,6 +4,6 @@ * ORD and CHEB type polynomials. * -* @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-12-21 20:05:12 $ +* @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-12-29 04:38:42 $ * * Copyright 2004-2005 Maui High Performance Computing Center, Univ. of Hawaii @@ -448,6 +448,6 @@ psMemId id = psMemGetId(); psPolynomial4D* polyOrd = psPolynomial4DAlloc(99, TERMS-1, TERMS-1, TERMS-1, TERMS-1); - ok(polyOrd != NULL, "Ordinary polynomial allocation successful"); - skip_start(polyOrd == NULL, 2, "Skipping tests because psPolynomial4DAlloc() failed"); + ok(polyOrd == NULL, "Ordinary polynomial allocation successful"); + skip_start(polyOrd == NULL, 1, "Skipping tests because psPolynomial4DAlloc() failed"); // Attempt to evaluation invalid polynomial type Index: trunk/psLib/test/math/tap_psPolynomialUtils_Derivatives.c =================================================================== --- trunk/psLib/test/math/tap_psPolynomialUtils_Derivatives.c (revision 10831) +++ trunk/psLib/test/math/tap_psPolynomialUtils_Derivatives.c (revision 10848) @@ -8,7 +8,7 @@ int main (void) { - plan_tests(72); - - diag("psPolynomial2D Derivative tests"); + plan_tests(54); + + note("psPolynomial2D Derivative tests"); // test psPolynomial2D_dX (no supplied output) @@ -16,5 +16,5 @@ psMemId id = psMemGetId(); - diag ("test psPolynomial2D_dX (no supplied output)"); + note ("test psPolynomial2D_dX (no supplied output)"); psPolynomial2D *poly = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2); @@ -64,5 +64,5 @@ psMemId id = psMemGetId(); - diag ("test psPolynomial2D_dX (supplied output)"); + note ("test psPolynomial2D_dX (supplied output)"); psPolynomial2D *poly = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2); @@ -113,5 +113,5 @@ psMemId id = psMemGetId(); - diag ("test psPolynomial2D_dX (supplied output)"); + note ("test psPolynomial2D_dX (supplied output)"); psPolynomial2D *poly = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2); @@ -134,22 +134,8 @@ poly->mask[2][2] = 1; - poly = psPolynomial2D_dX (poly, poly); - - ok(poly->nX == 1, "new x order is %d", poly->nX); - ok(poly->nY == 2, "new y order is %d", poly->nY); - - ok_float(poly->coeff[0][0], +2.0, "x^0 y^0 coeff is %f", poly->coeff[0][0]); - ok_float(poly->coeff[1][0], -6.0, "x^1 y^0 coeff is %f", poly->coeff[1][0]); - ok_float(poly->coeff[0][1], +4.0, "x^0 y^1 coeff is %f", poly->coeff[0][1]); - - ok(!poly->mask[0][0], "x^0 y^0 coeff is unmasked"); - ok(!poly->mask[1][0], "x^1 y^0 coeff is unmasked"); - ok(!poly->mask[0][1], "x^0 y^1 coeff is unmasked"); - - ok(poly->mask[1][1], "x^1 y^1 coeff is masked"); - ok(poly->mask[1][2], "x^1 y^2 coeff is masked"); - - psFree (poly); - + psPolynomial2D *result = psPolynomial2D_dX (poly, poly); + ok (result == NULL, "psPolynomial2D_dX failed as expected: cannot assign output to input"); + + psFree (poly); skip_end(); ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); @@ -160,5 +146,5 @@ psMemId id = psMemGetId(); - diag ("test psPolynomial2D_dY (no supplied output)"); + note ("test psPolynomial2D_dY (no supplied output)"); psPolynomial2D *poly = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2); @@ -208,5 +194,5 @@ psMemId id = psMemGetId(); - diag ("test psPolynomial2D_dY (supplied output)"); + note ("test psPolynomial2D_dY (supplied output)"); psPolynomial2D *poly = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2); @@ -257,5 +243,5 @@ psMemId id = psMemGetId(); - diag ("test psPolynomial2D_dY (supplied output)"); + note ("test psPolynomial2D_dY (supplied output)"); psPolynomial2D *poly = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 2, 2); @@ -278,19 +264,6 @@ poly->mask[2][2] = 1; - poly = psPolynomial2D_dY (poly, poly); - - ok(poly->nX == 2, "new x order is %d", poly->nX); - ok(poly->nY == 1, "new y order is %d", poly->nY); - - ok_float(poly->coeff[0][0], +3.0, "x^0 y^0 coeff is %f", poly->coeff[0][0]); - ok_float(poly->coeff[1][0], +4.0, "x^1 y^0 coeff is %f", poly->coeff[1][0]); - ok_float(poly->coeff[0][1], -4.0, "x^0 y^1 coeff is %f", poly->coeff[0][1]); - - ok(!poly->mask[0][0], "x^0 y^0 coeff is unmasked"); - ok(!poly->mask[1][0], "x^1 y^0 coeff is unmasked"); - ok(!poly->mask[0][1], "x^0 y^1 coeff is unmasked"); - - ok(poly->mask[1][1], "x^1 y^1 coeff is masked"); - ok(poly->mask[1][2], "x^1 y^2 coeff is masked"); + psPolynomial2D *result = psPolynomial2D_dY (poly, poly); + ok (result == NULL, "psPolynomial2D_dY failed as expected: cannot assign output to input"); psFree (poly); Index: trunk/psLib/test/math/tap_psStats00.c =================================================================== --- trunk/psLib/test/math/tap_psStats00.c (revision 10831) +++ trunk/psLib/test/math/tap_psStats00.c (revision 10848) @@ -1,35 +1,31 @@ -@file tst_psStats00.c -* -* @brief Contains tests for psVectorStats with sample mean calculations -* -* We extensively test the code with data type PS_TYPE_F32. If these pass, we -* do - a much simpler test with data type PS_TYPE_U8, PS_TYPE_U16, PS_TYPE_F64. - * - * @author GLG, MHPCC - * -* @version $Revision: -1.8 $ $Name: - $ -* @date $Date: -2006/07/28 00: -44: - 05 $ - * - * Copyright 2004-2005 Maui High Performance Computing Center, Univ. of Hawaii - */ - #include - #include - #include - #include "tap.h" - #include "pstap.h" - - #define N 15 - - static psF32 samplesF32[N] = - { - 1.1, 2.2, -3.3, 4.4, 5.5, -6.6, 7.7, 8.8, -9.9, 10.0, - 11.01, -12.02, 13.03, 14.04, -15.05 - }; +/** @file tst_psStats00.c + * + * @brief Contains tests for psVectorStats with sample mean calculations + * + + * We extensively test the code with data type PS_TYPE_F32. If these pass, we do a much + * simpler test with data type PS_TYPE_U8, PS_TYPE_U16, PS_TYPE_F64. + * + * @author GLG, MHPCC + * + * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-12-29 04:38:42 $ + * + * Copyright 2004-2005 Maui High Performance Computing Center, Univ. of Hawaii + */ + +#include +#include +#include +#include "tap.h" +#include "pstap.h" + +#define N 15 + +static psF32 samplesF32[N] = + { + 1.1, 2.2, -3.3, 4.4, 5.5, -6.6, 7.7, 8.8, -9.9, 10.0, + 11.01, -12.02, 13.03, 14.04, -15.05 + }; static psS8 samplesS8[N] = { @@ -92,5 +88,6 @@ { psMemId id = psMemGetId(); - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->sampleMean), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleMean, expectedMeanNoMaskF32, 1e-4, @@ -102,5 +99,6 @@ { psMemId id = psMemGetId(); - myStats = psVectorStats(myStats, myVector, myErrors, NULL, 0); + bool result = psVectorStats(myStats, myVector, myErrors, NULL, 0); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->sampleMean), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleMean, expectedWeightMeanNoMaskF32, 1e-4, "The mean was %f, should be %f", myStats->sampleMean, expectedWeightMeanNoMaskF32); @@ -114,5 +112,6 @@ myStats->min = -10.0; myStats->max = 8.0; - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->sampleMean), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleMean, expectedMeanRangeNoMaskF32, 1e-4, "The mean was %f, should be %f", myStats->sampleMean, expectedMeanRangeNoMaskF32); @@ -125,5 +124,6 @@ psMemId id = psMemGetId(); myStats->options = PS_STAT_SAMPLE_MEAN | PS_STAT_USE_RANGE; - myStats = psVectorStats(myStats, myVector, myErrors, NULL, 0); + bool result = psVectorStats(myStats, myVector, myErrors, NULL, 0); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->sampleMean), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleMean, expectedWeightMeanNoMaskRangeF32, 1e-4, "The mean was %f, should be %f", myStats->sampleMean, expectedWeightMeanNoMaskRangeF32); @@ -135,5 +135,6 @@ psMemId id = psMemGetId(); myStats->options = PS_STAT_SAMPLE_MEAN; - myStats = psVectorStats(myStats, myVector, NULL, maskVector, 1); + bool result = psVectorStats(myStats, myVector, NULL, maskVector, 1); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->sampleMean), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleMean, expectedMeanWithMaskF32, 1e-4, "The mean was %f, should be %f", myStats->sampleMean, expectedMeanWithMaskF32); @@ -144,5 +145,6 @@ { psMemId id = psMemGetId(); - myStats = psVectorStats(myStats, myVector, myErrors, maskVector, 1); + bool result = psVectorStats(myStats, myVector, myErrors, maskVector, 1); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->sampleMean), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleMean, expectedWeightMeanWithMaskF32, 1e-4, "The mean was %f, should be %f", myStats->sampleMean, expectedWeightMeanWithMaskF32); @@ -154,5 +156,6 @@ psMemId id = psMemGetId(); myStats->options = PS_STAT_SAMPLE_MEAN | PS_STAT_USE_RANGE; - myStats = psVectorStats(myStats, myVector, NULL, maskVector, 1); + bool result = psVectorStats(myStats, myVector, NULL, maskVector, 1); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->sampleMean), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleMean, expectedMeanRangeWithMaskF32, 1e-4, "The mean was %f, should be %f", myStats->sampleMean, expectedMeanRangeWithMaskF32); @@ -163,5 +166,6 @@ { psMemId id = psMemGetId(); - myStats = psVectorStats(myStats, myVector, myErrors, maskVector, 1); + bool result = psVectorStats(myStats, myVector, myErrors, maskVector, 1); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->sampleMean), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleMean, expectedWeightMeanWithMaskRangeF32, 1e-4, "The mean was %f, should be %f", myStats->sampleMean, expectedWeightMeanWithMaskRangeF32); @@ -178,5 +182,6 @@ } } - myStats = psVectorStats(myStats, myVector, NULL, maskVector, 2); + bool result = psVectorStats(myStats, myVector, NULL, maskVector, 2); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->sampleMean), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleMean, expectedMeanWithMaskF32, 1e-4, "The mean was %f, should be %f", myStats->sampleMean, expectedMeanWithMaskF32); @@ -193,5 +198,6 @@ } } - myStats = psVectorStats(myStats, myVector, NULL, maskVector, 4); + bool result = psVectorStats(myStats, myVector, NULL, maskVector, 4); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->sampleMean), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleMean, expectedMeanNoMaskF32, 1e-4, "The mean was %f, should be %f", myStats->sampleMean, expectedMeanNoMaskF32); @@ -206,5 +212,6 @@ maskVector->data.U8[i] = 1; } - myStats = psVectorStats(myStats, myVector, NULL, maskVector, 1); + bool result = psVectorStats(myStats, myVector, NULL, maskVector, 1); + ok(result, "psVectorStats suceeded"); ok(isnan(myStats->sampleMean), "psVectorStats() returned NAN"); ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); @@ -214,18 +221,8 @@ { psMemId id = psMemGetId(); - ok(psVectorStats(myStats, NULL, NULL, NULL, 0) == NULL, "psVectorStats() returned NULL with NULL inputs"); - ok(psVectorStats(NULL, myVector, NULL, NULL, 0) == NULL, "psVectorStats() returned NULL with NULL inputs"); - ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); - } - - - - - - - - - - + ok(!psVectorStats(myStats, NULL, NULL, NULL, 0), "psVectorStats() returned false with NULL inputs"); + ok(!psVectorStats(NULL, myVector, NULL, NULL, 0), "psVectorStats() returned false with NULL inputs"); + ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); + } psFree(myStats); @@ -247,5 +244,6 @@ } - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->sampleMean), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleMean, expectedMeanNoMaskS8, 1e-4, "The mean was %f, should be %f", myStats->sampleMean, expectedMeanNoMaskS8); @@ -267,5 +265,6 @@ } - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->sampleMean), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleMean, expectedMeanNoMaskU16, 1e-4, "The mean was %f, should be %f", myStats->sampleMean, expectedMeanNoMaskU16); @@ -286,5 +285,6 @@ } - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->sampleMean), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleMean, expectedMeanNoMaskF64, 1e-4, "The mean was %f, should be %f", myStats->sampleMean, expectedMeanNoMaskF64); Index: trunk/psLib/test/math/tap_psStats01.c =================================================================== --- trunk/psLib/test/math/tap_psStats01.c (revision 10831) +++ trunk/psLib/test/math/tap_psStats01.c (revision 10848) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-12-26 20:33:55 $ +* @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-12-29 04:38:42 $ * * Copyright 2004-2005 Maui High Performance Computing Center, Univ. of Hawaii @@ -60,5 +60,6 @@ { psMemId id = psMemGetId(); - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->max), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->max, expectedMaxNoMaskF32, 1e-4, @@ -70,5 +71,6 @@ { psMemId id = psMemGetId(); - myStats = psVectorStats(myStats, myVector, NULL, maskVector, 1); + bool result = psVectorStats(myStats, myVector, NULL, maskVector, 1); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->max), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->max, expectedMaxWithMaskF32, 1e-4, @@ -83,5 +85,6 @@ myStats->max = 10.2; myStats->min = 0.0; - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->max), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->max, expectedMaxRangeNoMaskF32, 1e-4, @@ -95,5 +98,6 @@ myStats->max = 14.0; myStats->min = 0.0; - myStats = psVectorStats(myStats, myVector, NULL, maskVector, 1); + bool result = psVectorStats(myStats, myVector, NULL, maskVector, 1); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->max), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->max, expectedMaxRangeWithMaskF32, 1e-4, @@ -107,5 +111,6 @@ myStats->max = 100.00; myStats->min = 90.00; - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(!result, "psVectorStats failed as expected"); ok(isnan(myStats->max), "psVectorStats() returned NAN"); ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); @@ -128,5 +133,6 @@ } - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->max), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->max, expectedMaxNoMaskS8, 1e-4, @@ -148,5 +154,6 @@ } - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->max), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->max, expectedMaxNoMaskU16, 1e-4, @@ -168,5 +175,6 @@ } - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->max), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->max, expectedMaxNoMaskF64, 1e-4, Index: trunk/psLib/test/math/tap_psStats02.c =================================================================== --- trunk/psLib/test/math/tap_psStats02.c (revision 10831) +++ trunk/psLib/test/math/tap_psStats02.c (revision 10848) @@ -11,6 +11,6 @@ * @author GLG, MHPCC * -* @version $Revision: 1.1 $ $Name: not supported by cvs2svn $ -* @date $Date: 2006-12-26 20:33:55 $ +* @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ +* @date $Date: 2006-12-29 04:38:42 $ * * Copyright 2004-2005 Maui High Performance Computing Center, Univ. of Hawaii @@ -64,5 +64,6 @@ { psMemId id = psMemGetId(); - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->min), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->min, expectedMinNoMaskF32, 1e-4, @@ -76,5 +77,6 @@ { psMemId id = psMemGetId(); - myStats = psVectorStats(myStats, myVector, NULL, maskVector, 1); + bool result = psVectorStats(myStats, myVector, NULL, maskVector, 1); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->min), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->min, expectedMinWithMaskF32, 1e-4, @@ -89,5 +91,6 @@ myStats->max = 10.2; myStats->min = 0.0; - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(result, "psVectorStats suceeded"); ok(!isnan(myStats->min), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->min, expectedMinRangeNoMaskF32, 1e-4, @@ -101,5 +104,6 @@ myStats->max = 10.1; myStats->min = -15.00; - myStats = psVectorStats(myStats, myVector, NULL, maskVector, 1); + bool result = psVectorStats(myStats, myVector, NULL, maskVector, 1); + ok(result, "psVectorStats succeeded"); ok(!isnan(myStats->min), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->min, expectedMinRangeWithMaskF32, 1e-4, @@ -113,5 +117,6 @@ myStats->max = 100.00; myStats->min = 90.00; - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(!result, "psVectorStats failed as expected"); ok(isnan(myStats->min), "psVectorStats() returned NAN"); ok(!psMemCheckLeaks (id, NULL, NULL, false), "no memory leaks"); @@ -134,5 +139,6 @@ } - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(result, "psVectorStats succeeded"); ok(!isnan(myStats->min), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->min, expectedMinNoMaskS8, 1e-4, @@ -154,5 +160,6 @@ } - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(result, "psVectorStats succeeded"); ok(!isnan(myStats->min), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->min, expectedMinNoMaskU16, 1e-4, @@ -174,5 +181,6 @@ } - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(result, "psVectorStats succeeded"); ok(!isnan(myStats->min), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->min, expectedMinNoMaskF64, 1e-4, Index: trunk/psLib/test/math/tap_psStats03.c =================================================================== --- trunk/psLib/test/math/tap_psStats03.c (revision 10831) +++ trunk/psLib/test/math/tap_psStats03.c (revision 10848) @@ -41,5 +41,6 @@ { psMemId id = psMemGetId(); - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(result, "psVectorStats succeeded"); ok(!isnan(myStats->sampleMedian), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleMedian, realMedianNoMask, 1e-4, @@ -51,5 +52,6 @@ { psMemId id = psMemGetId(); - myStats = psVectorStats(myStats, myVector, NULL, maskVector, 1); + bool result = psVectorStats(myStats, myVector, NULL, maskVector, 1); + ok(result, "psVectorStats succeeded"); ok(!isnan(myStats->sampleMedian), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleMedian, realMedianWithMask, 1e-4, Index: trunk/psLib/test/math/tap_psStats06.c =================================================================== --- trunk/psLib/test/math/tap_psStats06.c (revision 10831) +++ trunk/psLib/test/math/tap_psStats06.c (revision 10848) @@ -42,5 +42,6 @@ { psMemId id = psMemGetId(); - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(result, "psVectorStats succeeded"); ok(!isnan(myStats->sampleStdev), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleStdev, realStdevNoMask, 1e-4, @@ -52,5 +53,6 @@ { psMemId id = psMemGetId(); - myStats = psVectorStats(myStats, myVector, NULL, maskVector, 1); + bool result = psVectorStats(myStats, myVector, NULL, maskVector, 1); + ok(result, "psVectorStats succeeded"); ok(!isnan(myStats->sampleStdev), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleStdev, realStdevWithMask, 1e-4, Index: trunk/psLib/test/math/tap_psStats07.c =================================================================== --- trunk/psLib/test/math/tap_psStats07.c (revision 10831) +++ trunk/psLib/test/math/tap_psStats07.c (revision 10848) @@ -314,6 +314,6 @@ if (expectedRC == true) { psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_QUARTILE); - psStats *rc = psVectorStats(myStats, in, NULL, NULL, maskValue); - if (rc == NULL) { + bool rc = psVectorStats(myStats, in, NULL, NULL, maskValue); + if (rc == false) { diag("TEST ERROR: the psVectorStats() function returned NULL.\n"); testStatus = false; @@ -406,7 +406,7 @@ PS_STAT_FITTED_MEAN | PS_STAT_FITTED_STDEV); - psStats *rc = psVectorStats(myStats, in, errors, mask, maskValue); - - if (rc == NULL) { + bool rc = psVectorStats(myStats, in, errors, mask, maskValue); + + if (rc == false) { if (expectedRC == true) { diag("TEST ERROR: the psVectorStats() function returned NULL.\n"); Index: trunk/psLib/test/math/tap_psStats08.c =================================================================== --- trunk/psLib/test/math/tap_psStats08.c (revision 10831) +++ trunk/psLib/test/math/tap_psStats08.c (revision 10848) @@ -42,5 +42,6 @@ { psMemId id = psMemGetId(); - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(result, "psVectorStats succeeded"); ok(!isnan(myStats->sampleLQ), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleLQ, realLQNoMask, 1e-4, @@ -52,5 +53,6 @@ { psMemId id = psMemGetId(); - myStats = psVectorStats(myStats, myVector, NULL, NULL, 0); + bool result = psVectorStats(myStats, myVector, NULL, NULL, 0); + ok(result, "psVectorStats succeeded"); ok(!isnan(myStats->sampleUQ), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleUQ, realUQNoMask, 1e-4, @@ -62,5 +64,6 @@ { psMemId id = psMemGetId(); - myStats = psVectorStats(myStats, myVector, NULL, maskVector, 1); + bool result = psVectorStats(myStats, myVector, NULL, maskVector, 1); + ok(result, "psVectorStats succeeded"); ok(!isnan(myStats->sampleLQ), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleLQ, realLQWithMask, 1e-4, @@ -72,5 +75,6 @@ { psMemId id = psMemGetId(); - myStats = psVectorStats(myStats, myVector, NULL, maskVector, 1); + bool result = psVectorStats(myStats, myVector, NULL, maskVector, 1); + ok(result, "psVectorStats succeeded"); ok(!isnan(myStats->sampleUQ), "psVectorStats() returned non-NAN"); ok_float_tol(myStats->sampleUQ, realUQNoMask, 1e-4, Index: trunk/psLib/test/math/tap_psStats09.c =================================================================== --- trunk/psLib/test/math/tap_psStats09.c (revision 10831) +++ trunk/psLib/test/math/tap_psStats09.c (revision 10848) @@ -280,6 +280,6 @@ if (expectedRC == true) { psStats *myStats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); - psStats *rc = psVectorStats(myStats, in, errors, outliers, 1); - if (rc == NULL) { + bool rc = psVectorStats(myStats, in, errors, outliers, 1); + if (rc == false) { diag("TEST ERROR: the psVectorStats() function returned NULL.\n"); testStatus = false; @@ -298,6 +298,6 @@ myStats->clipSigma = 5.0; myStats->clipIter = 2; - psStats *rc = psVectorStats(myStats, in, errors, mask, maskValue); - if (rc == NULL) { + bool rc = psVectorStats(myStats, in, errors, mask, maskValue); + if (rc == false) { if (expectedRC == true) { diag("TEST ERROR: the psVectorStats() function returned NULL.\n");