Index: /trunk/psLib/src/astro/psCoord.c =================================================================== --- /trunk/psLib/src/astro/psCoord.c (revision 2196) +++ /trunk/psLib/src/astro/psCoord.c (revision 2197) @@ -10,6 +10,6 @@ * @author George Gusciora, MHPCC * -* @version $Revision: 1.25 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-10-14 01:17:50 $ +* @version $Revision: 1.26 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-10-26 21:24:42 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -58,37 +58,4 @@ /*****************************************************************************/ -static float p_psArg(float x, float y); - -/****************************************************************************** -XXX: Verify this arc tan function. -XXX: macro -XXX: I don't think this arg() is correct. - *****************************************************************************/ - -float p_psArg(float x, - float y) -{ - if (y!= 0.0) { - return(atan(y/x)); - } else { - return(0.5 * M_PI); - } - return (0.0); - /* - if (x > 0) { - return (atan(y / x)); - } else if ((x == 0) && (y == 0)) { - return (0.5 * M_PI); - } else if ((x == 0) && (y == 0)) { - return (-0.5 * M_PI); - } else if ((x == 0) && (y == 0)) { - return (M_PI + atan(y / x)); - } else if ((x == 0) && (y == 0)) { - return (-M_PI + atan(y / x)); - } - psAbort(__func__, "Unacceptable range for (p_psArg(%f, %f).\n", x, y); - */ -} - /*****************************************************************************/ /* FUNCTION IMPLEMENTATION - PUBLIC */ @@ -292,10 +259,10 @@ if (projection->type == PS_PROJ_TAN) { R = sqrt((coord->x * coord->x) + (coord->y * coord->y)); - tmp->d = p_psArg(-coord->y, coord->x); + tmp->d = atan2(-coord->y, coord->x); tmp->r = atan(180.0 / (R * M_PI)); } else if (projection->type == PS_PROJ_SIN) { R = sqrt((coord->x * coord->x) + (coord->y * coord->y)); - tmp->d = p_psArg(-coord->y, coord->x); + tmp->d = atan2(-coord->y, coord->x); tmp->r = acos((R * M_PI) / 180.0); @@ -314,5 +281,5 @@ chu2 *= chu2; chu = sqrt(1.0 - chu1 - chu2); - tmp->d = 2.0 * p_psArg((2.0 * chu * chu) - 1.0, (coord->x * chu * M_PI) / 360.0); + tmp->d = 2.0 * atan2((2.0 * chu * chu) - 1.0, (coord->x * chu * M_PI) / 360.0); tmp->r = asin((coord->y * chu * M_PI) / 180.0); Index: /trunk/psLib/src/astronomy/psCoord.c =================================================================== --- /trunk/psLib/src/astronomy/psCoord.c (revision 2196) +++ /trunk/psLib/src/astronomy/psCoord.c (revision 2197) @@ -10,6 +10,6 @@ * @author George Gusciora, MHPCC * -* @version $Revision: 1.25 $ $Name: not supported by cvs2svn $ -* @date $Date: 2004-10-14 01:17:50 $ +* @version $Revision: 1.26 $ $Name: not supported by cvs2svn $ +* @date $Date: 2004-10-26 21:24:42 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -58,37 +58,4 @@ /*****************************************************************************/ -static float p_psArg(float x, float y); - -/****************************************************************************** -XXX: Verify this arc tan function. -XXX: macro -XXX: I don't think this arg() is correct. - *****************************************************************************/ - -float p_psArg(float x, - float y) -{ - if (y!= 0.0) { - return(atan(y/x)); - } else { - return(0.5 * M_PI); - } - return (0.0); - /* - if (x > 0) { - return (atan(y / x)); - } else if ((x == 0) && (y == 0)) { - return (0.5 * M_PI); - } else if ((x == 0) && (y == 0)) { - return (-0.5 * M_PI); - } else if ((x == 0) && (y == 0)) { - return (M_PI + atan(y / x)); - } else if ((x == 0) && (y == 0)) { - return (-M_PI + atan(y / x)); - } - psAbort(__func__, "Unacceptable range for (p_psArg(%f, %f).\n", x, y); - */ -} - /*****************************************************************************/ /* FUNCTION IMPLEMENTATION - PUBLIC */ @@ -292,10 +259,10 @@ if (projection->type == PS_PROJ_TAN) { R = sqrt((coord->x * coord->x) + (coord->y * coord->y)); - tmp->d = p_psArg(-coord->y, coord->x); + tmp->d = atan2(-coord->y, coord->x); tmp->r = atan(180.0 / (R * M_PI)); } else if (projection->type == PS_PROJ_SIN) { R = sqrt((coord->x * coord->x) + (coord->y * coord->y)); - tmp->d = p_psArg(-coord->y, coord->x); + tmp->d = atan2(-coord->y, coord->x); tmp->r = acos((R * M_PI) / 180.0); @@ -314,5 +281,5 @@ chu2 *= chu2; chu = sqrt(1.0 - chu1 - chu2); - tmp->d = 2.0 * p_psArg((2.0 * chu * chu) - 1.0, (coord->x * chu * M_PI) / 360.0); + tmp->d = 2.0 * atan2((2.0 * chu * chu) - 1.0, (coord->x * chu * M_PI) / 360.0); tmp->r = asin((coord->y * chu * M_PI) / 180.0); Index: /trunk/psLib/src/dataManip/psConstants.h =================================================================== --- /trunk/psLib/src/dataManip/psConstants.h (revision 2196) +++ /trunk/psLib/src/dataManip/psConstants.h (revision 2197) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-14 23:43:12 $ + * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-26 21:24:42 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -107,4 +107,11 @@ } \ +/** Preprocessor macro to generate error on a NULL vector */ +#define PS_CHECK_NULL_VECTOR_RETURN_F0(NAME) \ +if (NAME == NULL || NAME->data.V == NULL) { \ + psError(__func__,"Unallowable operation: %s or its data is NULL.", #NAME); \ + return(0.0); \ +} \ + /** Preprocessor macro to generate error on a NULL poniter */ #define PS_CHECK_NULL_PTR_RETURN_NULL(NAME) \ @@ -113,4 +120,12 @@ return(NULL); \ } + +/** Preprocessor macro to generate error on a NULL poniter */ +#define PS_CHECK_NULL_PTR_RETURN_F0(NAME) \ +if (NAME == NULL) { \ + psError(__func__,"Unallowable operation: %s is NULL.", #NAME); \ + return(0.0); \ +} + /** Preprocessor macro to generate error on a NULL poniter */ #define PS_CHECK_NULL_PTR_RETURN_ZERO(NAME) \ @@ -127,4 +142,11 @@ } \ +/** Preprocessor macro to generate error for zero length vector */ +#define PS_CHECK_EMPTY_VECTOR_RETURN_F0(NAME) \ +if (NAME->n < 1) { \ + psError(__func__,"Unallowable operation: %s has zero n value.", #NAME); \ + return(0.0); \ +} \ + /** Preprocessor macro to generate error on differing size vectors */ #define PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(VEC1, VEC2) \ @@ -132,4 +154,11 @@ psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \ return(NULL); \ +} + +/** Preprocessor macro to generate error on differing size vectors */ +#define PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_F0(VEC1, VEC2) \ +if (VEC1->n != VEC2->n) { \ + psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \ + return(0.0); \ } @@ -199,7 +228,7 @@ #define PS_MAX(A, B) \ -((A > B) ? A : B) \ +(((A) > (B)) ? (A) : (B)) \ #define PS_MIN(A, B) \ -((A < B) ? A : B) \ - +(((A) < (B)) ? (A) : (B)) \ + Index: /trunk/psLib/src/dataManip/psFunctions.c =================================================================== --- /trunk/psLib/src/dataManip/psFunctions.c (revision 2196) +++ /trunk/psLib/src/dataManip/psFunctions.c (revision 2197) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.52 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-14 01:44:48 $ + * @version $Revision: 1.53 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-26 21:24:42 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii Index: /trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- /trunk/psLib/src/dataManip/psMinimize.c (revision 2196) +++ /trunk/psLib/src/dataManip/psMinimize.c (revision 2197) @@ -9,6 +9,6 @@ * @author GLF, MHPCC * - * @version $Revision: 1.58 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-14 01:44:48 $ + * @version $Revision: 1.59 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-26 21:24:42 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -636,4 +636,261 @@ psFree(currValueVec); psFree(newValueVec); + + if ((min->iter < min->maxIter) || + (min->lastDelta <= min->tol)) { + return(true); + } + + psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + "---- psMinimizeLMChi2() end (false) ----\n"); + return(false); +} + +/****************************************************************************** +psMyMinimizeLMChi2(): This routine will take an procedure which calculates +an arbitrary function and it's derivative and minimize the chi-squared match +between that function at the specified coords and the specified value at +those coords. + +XXX: Do this: + After checking that all entries in the paramMask are 1 or 0, when + forming the A matrix from alpha, try this: + + A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i]; + +XXX: This is very different from what is specified in the SDR. Must +coordinate with IfA on new SDR. + +XXX: Do vector/image recycles. + +XXX: probably yErr will be part of the SDR. + +XXX: This must work for both F32 and F64. F32 is currently implemented. + Note: since the LUD routines are only implemented in F64, then we + will have to convert all F32 input vectors to F64 regardless. So, + the F64 port might be. + *****************************************************************************/ +bool psMyMinimizeLMChi2(psMinimization *min, + psImage *covar, + psVector *params, + const psVector *paramMask, + const psArray *x, + const psVector *y, + const psVector *yErr, + psMyMinimizeLMChi2Func func) +{ + PS_CHECK_NULL_PTR_RETURN_NULL(min); + PS_CHECK_NULL_VECTOR_RETURN_NULL(params); + PS_CHECK_EMPTY_VECTOR_RETURN_NULL(params); + PS_CHECK_NULL_PTR_RETURN_NULL(x); + PS_CHECK_NULL_VECTOR_RETURN_NULL(y); + PS_CHECK_EMPTY_VECTOR_RETURN_NULL(y); + PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(x, y); + if (paramMask != NULL) { + PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(params, paramMask); + } + if (yErr != NULL) { + PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(y, yErr); + } + + psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + "---- psMinimizeLMChi2() begin ----\n"); + int numData = y->n; + int numParams = params->n; + int i; + int j; + int k; + int l; + int n; + int p; + psVector *beta = psVectorAlloc(numParams, PS_TYPE_F64); + psVector *perm = psVectorAlloc(numParams, PS_TYPE_F64); + + psVector *paramDeltasF64 = psVectorAlloc(numParams, PS_TYPE_F64); + psVector *origParams = psVectorAlloc(numParams, PS_TYPE_F32); + psVector *newParams = psVectorAlloc(numParams, PS_TYPE_F32); + + psImage *alpha = psImageAlloc(numParams, numParams, PS_TYPE_F32); + psImage *A = psImageAlloc(numParams, numParams, PS_TYPE_F64); + psImage *aOut = psImageAlloc(numParams, numParams, PS_TYPE_F64); + + // psVector **deriv = (psVector **) psAlloc(numData * sizeof(psVector *)); + // for (i=0;imaxIter is %d\n", min->maxIter); + psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + "min->tol is %f\n", min->tol); + + for (p=0;pdata.F32[p] = params->data.F32[p]; + } + + min->lastDelta = HUGE; + min->iter = 0; + + while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) { + psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + "------------------------------------------------------\n"); + psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + "Iteration %d. Delta is %f\n", min->iter, min->lastDelta); + + // + // Calculate the current values and chi-squared of the function. + // + currChi2 = 0.0; + currValueVec = func(deriv, params, x); + for (n=0;ndata.F32[n] * currValueVec->data.F32[n]); + } + + for (p=0;pdata.F32[%d] is %f.\n", p, params->data.F32[p]); + } + psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + "Current chi-squared is (%f)\n", currChi2); + + // + // Mask elements of the derivative for each data point. + // + for (p=0;pdata.U8[p] != 0)) { + for (n=0;ndata.F32[n][p] = 0.0; + } + } + } + + // + // Calculate the BETA vector. + // + for (p=0;pdata.F64[p] = 0.0; + for (n=0;ndata.F64[p])+= + (y->data.F32[n] - currValueVec->data.F32[n]) * + deriv->data.F32[n][p]; + } + // XXX: multiple by -1 here? + (beta->data.F64[p])*= -1.0; + psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + "beta->data.F64[%d] is %f.\n", p, beta->data.F64[p]); + } + psFree(currValueVec); + + // + // Calculate the ALPHA matrix. + // + for (k=0;kdata.F32[k][l] = 0.0; + for (n=0;ndata.F32[k][l]+= deriv->data.F32[n][k] * + deriv->data.F32[n][l]; + } + } + } + + // + // Calculate the matrix A. + // + for (j=0;jdata.F64[j][k] = + (double) ((1.0 + lamda) * alpha->data.F32[j][k]); + } else { + A->data.F64[j][k] = (double) alpha->data.F32[j][k]; + } + } + } + for (j=0;jdata.F64[j][k]); + } + psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n"); + } + + // + // Solve A * alpha = Beta + // + aOut = psMatrixLUD(aOut, perm, A); + paramDeltasF64 = psMatrixLUSolve(paramDeltasF64, aOut, beta, perm); + + // + // Mask any masked parameters. + // + for (i=0;idata.F64[%d] is %f.\n", i, paramDeltasF64->data.F64[i]); + if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) { + newParams->data.F32[i] = origParams->data.F32[i]; + } else { + newParams->data.F32[i] = params->data.F32[i] - + (float) paramDeltasF64->data.F64[i]; + } + } + + psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + "Calling func() with new parameters:\n"); + for (i=0;idata.F32[%d] is %f.\n", i, newParams->data.F32[i]); + } + + + // + // Calculate new function values. + // + newChi2 = 0.0; + newValueVec = func(deriv, newParams, x); + for (n=0;ndata.F32[n] * newValueVec->data.F32[n]); + } + psFree(newValueVec); + + psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + "old/new chi-squareds are (%f, %f)\n", currChi2, newChi2); + + // + // If the new chi-squared is lower, then keep it. + // + if (currChi2 > newChi2) { + min->lastDelta = currChi2 - newChi2; + min->value = newChi2; + + // We already masked params. + for (i=0;idata.F32[i] = (float) newParams->data.F32[i]; + } + lamda*= 0.1; + } else { + lamda*= 10.0; + } + psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + "lamda is %f\n", lamda); + min->iter++; + } + psFree(beta); + psFree(perm); + psFree(paramDeltasF64); + psFree(origParams); + psFree(newParams); + psFree(alpha); + psFree(A); + psFree(aOut); + psFree(deriv); if ((min->iter < min->maxIter) || Index: /trunk/psLib/src/dataManip/psMinimize.h =================================================================== --- /trunk/psLib/src/dataManip/psMinimize.h (revision 2196) +++ /trunk/psLib/src/dataManip/psMinimize.h (revision 2197) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.29 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-06 23:38:49 $ + * @version $Revision: 1.30 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-26 21:24:43 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -84,5 +84,82 @@ +/* + Bug 203: + +typedef float (*psMinimizeLMChi2Func)(psVector *deriv, const psVector *params, +const psVector *x); + +This function takes the current guess for the parameters for which we are +trying to get the best values (params), and a single vector (x) of +conditions from the array of x values fed into the minimiser. + + bool psMinimizeLMChi2(psMinimization *min, + psImage *covar, + psVector *params, + const psVector *paramMask, + const psArray *x, + const psVector *y, + const psVector *yErr, + psMinimizeLMChi2Func func); + +This takes the minimization specs (min), returns the covariance matrix +(covar), takes the best guess of initial parameters (params), the parameter +mask (paramMask), and takes multiple vectors of conditions in an array (x), +the corresponding measured values (y) and errors (yErr) and the function to +fit (func). + +For example, for GRB afterglows, I have flux as a function of time and +frequency. So I stuff into the "psArray *x" all my time and frequency values +(so a whole heap of vectors of size 2), I have the measured values in "y" and +errors in "yErr". Then each of the time-frequency pairs are passed to my model +function with the current parameters, and the model function returns the flux +for that time-frequency pair, and the derivative with respect to each of the +parameters. + +This seems reasonable. The only thing we could change would be to have the +function be defined: + +typedef psVector* (*psMinimizeLMChi2Func)(psImage *deriv, + const psVector *params, + const psArray *x); + +So it would return the model value for each of the measurements at once, and +return for each the derivatives (so that it returns a matrix). + +What do you think? + + + +I'm not sure I understand how LM chi-squared minimization will work with the +following function that you define: + + typedef psVector* (*psMinimizeLMChi2Func)(psImage *deriv, + const psVector *params, + const psArray *x); + +The chi-squared minimization algorithm, as defined in NR, requires that +function to be minimized be evaluated at each data point, and that all +derivatives, with respect to each parameter, be calculated at each data +point. In the above, can I assume that + + x is an array of psVectors, with each vector corresponding to a single + data point. + + The returned value has the same length as x. It contains the value of + the function at each data point in x. + + deriv: an n-by-p matrix where "n" is the number of data points, and "p" + is the number of parameters. The [i][j] element of this matrix + holds the derivative of the function at the i-th data point with + respect to the j-th parameter. + + + +*/ + + // XXX: What if any of these arguments are NULL? + + typedef float (*psMinimizeLMChi2Func)(psVector *deriv, @@ -98,4 +175,19 @@ psMinimizeLMChi2Func func); + +typedef +psVector* (*psMyMinimizeLMChi2Func)(psImage *deriv, + const psVector *params, + const psArray *x); + +bool psMyMinimizeLMChi2(psMinimization *min, + psImage *covar, + psVector *params, + const psVector *paramMask, + const psArray *x, + const psVector *y, + const psVector *yErr, + psMyMinimizeLMChi2Func func); + typedef float (*psMinimizePowellFunc)(const psVector *params, Index: /trunk/psLib/src/dataManip/psStats.c =================================================================== --- /trunk/psLib/src/dataManip/psStats.c (revision 2196) +++ /trunk/psLib/src/dataManip/psStats.c (revision 2197) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.67 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-13 02:40:13 $ + * @version $Revision: 1.68 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-26 21:24:43 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -40,5 +40,5 @@ /*****************************************************************************/ // will use robust statistical methods. -#define GAUSS_WIDTH 20 // The width of the Gaussian or boxcar smoothing. +#define GAUSS_WIDTH 5 // The width of the Gaussian or boxcar smoothing. #define CLIPPED_NUM_ITER_LB 1 #define CLIPPED_NUM_ITER_UB 10 @@ -57,4 +57,5 @@ /* TYPE DEFINITIONS */ /*****************************************************************************/ +psVector* p_psConvertToF32(psVector* in); /*****************************************************************************/ @@ -423,22 +424,4 @@ float rangeMax = 0.0; // Exclude date above this - // Determine if the number of data points exceed a threshold which will - // cause to generate robust stats, as opposed to exact stats. - // XXX: This code is no longer used. We now calculate the median exactly - // regardless of the vector size. - /* - * // Calculate the robust quartiles. stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); - * p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); - * - * // Store the robust quartiles into the sample quartile members. stats->sampleMedian = - * stats2->robustMedian; - * - * // Free temporary data buffers. psFree(stats2); - * - * // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. stats->options = stats->options | - * PS_STAT_ROBUST_FOR_SAMPLE; - * - * return; } */ - // Determine how many data points fit inside this min/max range // and are not masked, IF the maskVector is not NULL> @@ -507,67 +490,66 @@ /****************************************************************************** - This routine smoothes the data in the input robustHistogram with a - Gaussian of width sigma. - - XXX: Must consult with IfA on the proper values for the Gaussian - smoothing. Currently, the Gaussian coefficients are such that only - the center coefficient is non-zero. This is because "e" is being - raised to a very large negative number for all points except the - center point. - - XXX: use a static variable for gaussianCoefs[] and compute them once. - *****************************************************************************/ -psVector* p_psVectorsmoothHistGaussian(psHistogram* robustHistogram, +p_psVectorSmoothHistGaussian(): This routine smoothes the data in the input +robustHistogram with a Gaussian of width sigma. + *****************************************************************************/ +#define GAUSS_WIDTH_OTHER 5.0 +psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram, float sigma) { int i = 0; // Loop index variable int j = 0; // Loop index variable - float denom = 0.0; // Temporary variable - float expo = 0.0; // Temporary variable - float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients - psVector* smooth = psVectorAlloc(robustHistogram->nums->n, PS_TYPE_F32); - - for (i = 0; i < (1 + (2 * GAUSS_WIDTH)); i++) { - if (fabs(sigma) >= FLT_EPSILON) { - // If sigma does not equal zero, then we use Gaussian smoothing. - #ifdef DARWIN - denom = (float)sqrt(2.0 * M_PI * sigma * sigma); - #else - - denom = sqrtf(2.0 * M_PI * sigma * sigma); - #endif - - expo = -(float)((i - GAUSS_WIDTH) * (i - GAUSS_WIDTH)); - expo /= (2.0 * sigma * sigma); - gaussianCoefs[i] = exp(expo / denom); - - //printf("p_psVectorsmoothHistGaussian(): gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]); - // NOTE: Gaussian smoothing just isn't working with low sigma - // values. The problem is that the Gaussian coefficients are - // all zero, except for the middle coefficient, which is exactly - // one. Therefore, I'm using boxcar smoothing. - gaussianCoefs[i] = 1.0 / (1.0 + (2.0 * (float)GAUSS_WIDTH)); - // printf("gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]); - } else { - /* If sigma equals zero (all pixels have the same value) the above code will divide by zero. - * Therefore, we don't need to smooth the data. */ - for (i = 0; i < robustHistogram->nums->n; i++) { - smooth->data.F32[i] = (float)robustHistogram->nums->data.S32[i]; - } - return (smooth); - } - } - - // Perform the actual smoothing. - for (i = 0; i < robustHistogram->nums->n; i++) { + float iMid; + float jMid; + int numBins = robustHistogram->nums->n; + int numBounds = robustHistogram->bounds->n; + psVector* smooth = psVectorAlloc(numBins, PS_TYPE_F32); + int jMin = 0; + int jMax = 0; + float firstBound = robustHistogram->bounds->data.F32[0]; + float lastBound = robustHistogram->bounds->data.F32[numBounds-1]; + psScalar x; + + x.type.type = PS_TYPE_F32; + for (i = 0; i < numBins; i++) { + iMid = (robustHistogram->bounds->data.F32[i] + + robustHistogram->bounds->data.F32[i+1]) / 2.0; + + // YYY: The p_psVectorBinDisect() routine does much of the work of + // the following conditionals, however, it also reports a warning + // message. I don't want the warning message so I reproduce the + // conditionals here. Maybe p_psVectorBinDisect() should not produce + // warnings? + + x.data.F32 = iMid - (GAUSS_WIDTH * sigma); + if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { + jMin = p_psVectorBinDisect(robustHistogram->bounds, &x); + } else if (x.data.F32 <= firstBound) { + jMin = 0; + } else if (x.data.F32 >= lastBound) { + jMin = robustHistogram->bounds->n - 1; + + } + + x.data.F32 = iMid + (GAUSS_WIDTH * sigma); + if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { + jMax = p_psVectorBinDisect(robustHistogram->bounds, &x); + } else if (x.data.F32 <= firstBound) { + jMax = 0; + } else if (x.data.F32 >= lastBound) { + jMax = robustHistogram->bounds->n - 1; + + } + smooth->data.F32[i] = 0.0; - for (j = -GAUSS_WIDTH; j <= +GAUSS_WIDTH; j++) { - if (((j + i) >= 0) && ((j + i) < smooth->n)) { - smooth->data.F32[i] += (gaussianCoefs[j + GAUSS_WIDTH] * - (float)robustHistogram->nums->data.S32[j + i]); - } - } - } - return (smooth); + for (j = jMin ; j <= jMax ; j++) { + jMid = (robustHistogram->bounds->data.F32[j] + + robustHistogram->bounds->data.F32[j+1]) / 2.0; + smooth->data.F32[i] += + ((float) robustHistogram->nums->data.U32[j]) * + psGaussian(jMid, iMid, sigma, true); + } + } + + return(smooth); } @@ -913,5 +895,6 @@ specified data point. *****************************************************************************/ -float p_psGaussian(const psVector* restrict myData, const psVector* restrict myParams) +float p_psGaussian(const psVector* restrict myData, + const psVector* restrict myParams) { float x = myData->data.F32[0]; @@ -930,5 +913,6 @@ calculates the specified partial derivative of the above Gaussian function. *****************************************************************************/ -float p_psGaussianDeriv(const psVector* restrict myData, const psVector* restrict myParams, int whichParam) +float p_psGaussianDeriv(const psVector* restrict myData, + const psVector* restrict myParams, int whichParam) { float x = myData->data.F32[0]; @@ -1036,45 +1020,75 @@ /****************************************************************************** -p_psFitQuadratic(robustHistogram. cumulativeSums, binNum, fitFloat): given -the specified histogram, with the specified pre-calculated cumulativeSums, -this routine fits a quadratic f(x) to the 3 bins surrounding bin "binNum", -and then finds the point x such that f(x) == fitFloat. +fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): This routine +takes psVectors of x/y pairs as input, and fits a quadratic to the 3 points +surrounding element binNum in the vectors (the midpoint between element i +and i+1 is used for x[i]). It then determines for what value x does that +quadratic f(x) = yVal (the input parameter). -XXX: This function is currently not being used. - *****************************************************************************/ -float p_psFitQuadratic(psHistogram* histogram, - psVector* cumulativeSums, - int binNum, - float fitFloat) -{ +// XXX: Use static variables. +*****************************************************************************/ +float fitQuadraticSearchForYThenReturnX(psVector *xVec, + psVector *yVec, + int binNum, + float yVal) +{ + // static psVector* x = NULL; + // static psVector* y = NULL; + // static psVector* yErr = NULL; psVector* x = psVectorAlloc(3, PS_TYPE_F64); psVector* y = psVectorAlloc(3, PS_TYPE_F64); psVector* yErr = psVectorAlloc(3, PS_TYPE_F64); + + /* + if (x == NULL) { + x = psVectorAlloc(3, PS_TYPE_F64); + p_psMemSetPersistent(x, true); + } + if (y == NULL) { + y = psVectorAlloc(3, PS_TYPE_F64); + p_psMemSetPersistent(y, true); + } + if (yErr == NULL) { + yErr = psVectorAlloc(3, PS_TYPE_F64); + p_psMemSetPersistent(yErr, true); + } + */ psPolynomial1D* myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); - - if ((binNum > 0) && (binNum < (histogram->nums->n + 1))) { - x->data.F64[0] = (double)0.5 * - (histogram->bounds->data.F32[binNum - 1] + histogram->bounds->data.F32[binNum]); - x->data.F64[1] = (double)0.5 * - (histogram->bounds->data.F32[binNum] + histogram->bounds->data.F32[binNum + 1]); - x->data.F64[2] = (double)0.5 * - (histogram->bounds->data.F32[binNum + 1] + histogram->bounds->data.F32[binNum + 2]); - - y->data.F64[0] = cumulativeSums->data.F32[binNum - 1]; - y->data.F64[1] = cumulativeSums->data.F32[binNum]; - y->data.F64[2] = cumulativeSums->data.F32[binNum + 1]; - - if (!((y->data.F64[0] <= fitFloat) && (fitFloat <= y->data.F64[2]))) { - psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n"); - } - + float tmpFloat; + + if ((binNum > 0) && (binNum < (yVec->n - 2))) { + x->data.F64[0] = (double) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum])); + x->data.F64[1] = (double) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1])); + x->data.F64[2] = (double) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2])); + y->data.F64[0] = yVec->data.F32[binNum - 1]; + y->data.F64[1] = yVec->data.F32[binNum]; + y->data.F64[2] = yVec->data.F32[binNum + 1]; + + // Ensure that yVal is within the range of the bins we are using. + if (!((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2]))) { + printf("((%f), %f, %f)\n", yVal, y->data.F64[0], y->data.F64[2]); + psError(__func__, "yVal not within y-range\n"); + } yErr->data.F64[0] = 1.0; yErr->data.F64[1] = 1.0; yErr->data.F64[2] = 1.0; + // Determine the coefficients of the polynomial. myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr); - return (p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], fitFloat)); + // Call p_ps1DPolyMedian(), which does a binary search on the + // polynomial, looking for the value x such that f(x) = yVal + tmpFloat = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], yVal); } else { - return (0.5 * (histogram->bounds->data.F32[binNum + 1] + histogram->bounds->data.F32[binNum])); + if (binNum == 0) { + tmpFloat = 0.5 * (xVec->data.F32[binNum] + + xVec->data.F32[binNum + 1]); + } else if (binNum == (xVec->n - 1)) { + // XXX: Is this right? + tmpFloat = xVec->data.F32[binNum]; + } else if (binNum == (xVec->n - 2)) { + // XXX: Is this right? + tmpFloat = 0.5 * (xVec->data.F32[binNum] + + xVec->data.F32[binNum + 1]); + } } @@ -1083,6 +1097,9 @@ psFree(yErr); psFree(myPoly); - return (0.0); -} + return(tmpFloat); +} + +#define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \ +(0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1])) /****************************************************************************** @@ -1111,5 +1128,7 @@ *****************************************************************************/ void p_psVectorRobustStats(const psVector* restrict myVector, - const psVector* restrict maskVector, unsigned int maskVal, psStats* stats) + const psVector* restrict maskVector, + unsigned int maskVal, + psStats* stats) { psHistogram* robustHistogram = NULL; @@ -1118,19 +1137,10 @@ int LQBinNum = -1; // Bin num for lower quartile int UQBinNum = -1; // Bin num for upper quartile + int medianBinNum = -1; int i = 0; // Loop index variable - int maxBinNum = 0; - float maxBinCount = 0.0; + int modeBinNum = 0; + float modeBinCount = 0.0; float dL = 0.0; int numBins = 0; - psStats* tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN); - - // psImage* domain; - // psVector* errors; - // psVector* data; - // psVector* initialGuess; - // psVector* theParams; - // float chiSq=0.0; - // float max = -HUGE; - // float max_pos; float myMean = 0.0; float myStdev = 0.0; @@ -1139,14 +1149,11 @@ float sumSquares = 0.0; float sumDiffs = 0.0; - psVector* x = psVectorAlloc(3, PS_TYPE_F64); - psVector* y = psVectorAlloc(3, PS_TYPE_F64); - psVector* yErr = psVectorAlloc(3, PS_TYPE_F64); - psPolynomial1D* myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); - psVector* cumulativeRobustSumsFullRange = NULL; - psVector* cumulativeRobustSumsDlRange = NULL; + psVector* cumulativeRobustSums = NULL; float sumRobust = 0.0; float sumN50 = 0.0; float sumNfit = 0.0; - float cumulativeMedian = 0.0; + psScalar tmpScalar; + tmpScalar.type.type = PS_TYPE_F32; + psStats* tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN); // Compute the initial bin size of the robust histogram. This is done @@ -1188,4 +1195,5 @@ p_psVectorMax(myVector, maskVector, maskVal, stats); } + // Create the histogram structure. NOTE: we can not specify the bin size // precisely since the argument to psHistogramAlloc() is the number of @@ -1198,7 +1206,6 @@ psVectorHistogram(robustHistogram, myVector, maskVector, maskVal); - // Smooth the histogram. - // XXX: is that the right stdev? - robustHistogramVector = p_psVectorsmoothHistGaussian(robustHistogram, + // Smooth the histogram, Gaussian-style. + robustHistogramVector = p_psVectorSmoothHistGaussian(robustHistogram, tmpStats->clippedStdev / 4.0f); @@ -1208,87 +1215,62 @@ /************************************************************************** - Determine the lower/upper quartiles. + Determine the median/lower/upper quartile bin numbers. + + We define a vector called "cumulativeRobustSums" where the value at + index position i is equal to the sum of bins 0:i. This will be used in + determining the median and lower/upper quartiles. **************************************************************************/ - // We define a vector called "cumulativeRobustSums..." where the value at - // index position i is equal to the sum of bins 0:i. This will be used - // now and later in determining the lower/upper quartiles. - cumulativeRobustSumsFullRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); - cumulativeRobustSumsFullRange->data.F32[0] = robustHistogramVector->data.F32[0]; + cumulativeRobustSums = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); + cumulativeRobustSums->data.F32[0] = robustHistogramVector->data.F32[0]; for (i = 1; i < robustHistogramVector->n; i++) { - cumulativeRobustSumsFullRange->data.F32[i] = - cumulativeRobustSumsFullRange->data.F32[i - 1] + robustHistogramVector->data.F32[i]; - } - sumRobust = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n - 1]; - - // Determine the bin number containing the lower quartile point. - LQBinNum = -1; - for (i = 0; i < cumulativeRobustSumsFullRange->n; i++) { - if (cumulativeRobustSumsFullRange->data.F32[i] >= (sumRobust / 4.0)) { - LQBinNum = i; - break; - } - } - - // Determine the bin number containing the upper quartile point. - UQBinNum = -1; - for (i = cumulativeRobustSumsFullRange->n - 1; i >= 0; i--) { - if (cumulativeRobustSumsFullRange->data.F32[i] <= (3.0 * sumRobust / 4.0)) { - UQBinNum = i; - break; - } - } - - if ((LQBinNum == -1) || (UQBinNum == -1)) { - psAbort(__func__, "Could not determine the robust lower/upper quartiles."); - } - + cumulativeRobustSums->data.F32[i] = cumulativeRobustSums->data.F32[i - 1] + + robustHistogramVector->data.F32[i]; + } + sumRobust = cumulativeRobustSums->data.F32[robustHistogramVector->n - 1]; + + tmpScalar.data.F32 = sumRobust / 4.0; + LQBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar); + tmpScalar.data.F32 = 3.0 * sumRobust / 4.0; + UQBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar); + tmpScalar.data.F32 = sumRobust / 2.0; + medianBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar); + + if ((LQBinNum < 0) || (UQBinNum < 0)) { + psError(__func__, "Could not determine the robust lower/upper quartile bin numbers."); + } + if (medianBinNum < 0) { + psError(__func__, "Could not determine the robust lower/upper quartile bin numbers."); + } /************************************************************************** Determine the mode in the range LQ:UQ. **************************************************************************/ // Determine the bin with the peak value in the range LQ to UQ. - maxBinNum = LQBinNum; - maxBinCount = robustHistogramVector->data.F32[LQBinNum]; - sumN50 = (float)robustHistogram->nums->data.S32[LQBinNum]; + modeBinNum = LQBinNum; + modeBinCount = robustHistogramVector->data.F32[LQBinNum]; + sumN50 = (float)robustHistogram->nums->data.U32[LQBinNum]; for (i = LQBinNum + 1; i <= UQBinNum; i++) { - if (robustHistogramVector->data.F32[i] > maxBinCount) { - maxBinNum = i; - maxBinCount = robustHistogramVector->data.F32[i]; - } - sumN50 += (float)robustHistogram->nums->data.S32[i]; - } - - // XXX: is dL defined as the value at the LQ/UQ, or the bin number? + if (robustHistogramVector->data.F32[i] > modeBinCount) { + modeBinNum = i; + modeBinCount = robustHistogramVector->data.F32[i]; + } + sumN50 += (float)robustHistogram->nums->data.U32[i]; + } + dL = (UQBinNum - LQBinNum) / 4; - - printf("(LQBinNum, UQBinNum, maxBinNum) is (%d, %d, %d)\n", LQBinNum, UQBinNum, maxBinNum); - /************************************************************************** Determine the mean/stdev for the bins in the range mode-dL to mode+dL **************************************************************************/ - cumulativeRobustSumsDlRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); - for (i = 0; i < robustHistogramVector->n; i++) { - cumulativeRobustSumsDlRange->data.F32[i] = 0.0; - } - sumNfit = 0.0; - cumulativeMedian = 0.0; - for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) { - if ((0 <= i) && (i < robustHistogramVector->n)) { - cumulativeRobustSumsDlRange->data.F32[i] = - cumulativeRobustSumsDlRange->data.F32[i - 1] + robustHistogramVector->data.F32[i]; - cumulativeMedian += robustHistogramVector->data.F32[i]; - sumNfit += (float)robustHistogram->nums->data.S32[i]; - } - } - // Calculate the mean of the smoothed robust histogram in the range // mode-dL to mode+dL. We use the midpoint of each bin as the mean for // that bin (this is a non-exact approximation). + sumNfit = 0.0; myMean = 0.0; - for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) { + for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) { if ((0 <= i) && (i < robustHistogramVector->n)) { - myMean += (robustHistogramVector->data.F32[i]) * 0.5 * - (robustHistogram->bounds->data.F32[i + 1] + robustHistogram->bounds->data.F32[i]); + myMean += (robustHistogramVector->data.F32[i]) * PS_BIN_MIDPOINT(robustHistogram, i); countFloat += robustHistogramVector->data.F32[i]; } + + sumNfit += (float)robustHistogram->nums->data.U32[i]; } myMean /= countFloat; @@ -1297,8 +1279,7 @@ // mode-dL to mode+dL. We use the midpoint of each bin as the mean for // that bin. - for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) { + for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) { if ((0 <= i) && (i < robustHistogramVector->n)) { - diff = (0.5 * (robustHistogram->bounds->data.F32[i + 1] + - robustHistogram->bounds->data.F32[i])) - myMean; + diff = PS_BIN_MIDPOINT(robustHistogram, i) - myMean; sumSquares += diff * diff * robustHistogramVector->data.F32[i]; sumDiffs += diff * robustHistogramVector->data.F32[i]; @@ -1315,6 +1296,5 @@ if (stats->options & PS_STAT_ROBUST_MODE) { - stats->robustMode = 0.5 * - (robustHistogram->bounds->data.F32[maxBinNum] + robustHistogram->bounds->data.F32[maxBinNum + 1]); + stats->robustMode = PS_BIN_MIDPOINT(robustHistogram, modeBinNum); } @@ -1322,140 +1302,42 @@ stats->robustStdev = myStdev; } - // To determine the median (and later, the lower/upper quartile), we fit - // a quadratic to the three bins surrounding the bin containing the median. - // The quadratic y=f(x) with x being the midpoint of each bin, and y being - // the cumulative number of data points in all bins up to, and including, - // this bin. We then solve the quadratic for - - // XXX: Create a "p_psQuadraticFitAndSolveForX() function. - + + // To determine the median, we fit a quadratic y=f(x) to the three bins + // surrounding the bin containing the median (x is the midpoint of each + // bin and y is the value of each bin). Then we figure out what value + // of x corresponds to f(x) being the median (half of all points). if (stats->options & PS_STAT_ROBUST_MEDIAN) { - if ((maxBinNum > 0) && (maxBinNum < (robustHistogram->nums->n - 1))) { - x->data.F64[0] = (double)0.5 * - (robustHistogram->bounds->data.F32[maxBinNum - 1] + - robustHistogram->bounds->data.F32[maxBinNum]); - x->data.F64[1] = (double)0.5 * - (robustHistogram->bounds->data.F32[maxBinNum] + - robustHistogram->bounds->data.F32[maxBinNum + 1]); - x->data.F64[2] = (double)0.5 * - (robustHistogram->bounds->data.F32[maxBinNum + 1] + - robustHistogram->bounds->data.F32[maxBinNum + 2]); - - y->data.F64[0] = cumulativeRobustSumsDlRange->data.F32[maxBinNum - 1]; - y->data.F64[1] = cumulativeRobustSumsDlRange->data.F32[maxBinNum]; - y->data.F64[2] = cumulativeRobustSumsDlRange->data.F32[maxBinNum + 1]; - - // Ensure that cumulativeMedian/2 is actually within the range of the bins - // we are using. - cumulativeMedian *= 0.5; - if (!((y->data.F64[0] <= cumulativeMedian) && (cumulativeMedian <= y->data.F64[2]))) { - printf("((%f), %f, %f)\n", cumulativeMedian, y->data.F64[0], y->data.F64[2]); - psAbort(__func__, "p_psVectorRobustStats(1): midpoint not within y-range\n"); - } - // XXX: yErr is not currently used by psVectorFitPolynomial1D(). We - // may have to set this meaningfully later. - yErr->data.F64[0] = 1.0; - yErr->data.F64[1] = 1.0; - yErr->data.F64[2] = 1.0; - - // Determine the coefficients of the polynomial. - myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr); - // Call p_ps1DPolyMedian(), which does a binary search on the - // polynomial, looking for the value x such that - // f(x) = cumulativeMedian. - stats->robustMedian = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], cumulativeMedian); - } else { - // If the mode is the first/last histogram bin, then simply use - // the midpoint of that bin. - stats->robustMedian = 0.5 * (robustHistogram->bounds->data.F32[maxBinNum + 1] + - robustHistogram->bounds->data.F32[maxBinNum]); - } - } - // The lower/upper quartile calculations are very similar to the median - // calculations. We fit a quadratic to the array containing the - // cumulative data points in each bin, then search for x such that - // f(x) equals the lower/upper quartile exactly. - // + // Take a psVector. Fit a polynomial y = f(x) to the 3 data elements + // surrounding medianBinNum, then find the x-value corresponding y = sumRobust/2.0. + stats->robustMedian = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds, + cumulativeRobustSums, + medianBinNum, + sumRobust/2.0); + } + + // To determine the quartiles, we fit a quadratic y=f(x) to the three bins + // surrounding the bin containing LQ/UQ (x is the midpoint of each + // bin and y is the value of each bin). Then we figure out what value + // of x corresponds to f(x) being the LQ/UQ. if (stats->options & PS_STAT_ROBUST_QUARTILE) { - countFloat = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n - 1]; - - if ((LQBinNum > 0) && (LQBinNum < (robustHistogram->nums->n - 1))) { - x->data.F64[0] = (double)0.5 * - (robustHistogram->bounds->data.F32[LQBinNum - 1] + - robustHistogram->bounds->data.F32[LQBinNum]); - x->data.F64[1] = (double)0.5 * - (robustHistogram->bounds->data.F32[LQBinNum] + - robustHistogram->bounds->data.F32[LQBinNum + 1]); - x->data.F64[2] = (double)0.5 * - (robustHistogram->bounds->data.F32[LQBinNum + 1] + - robustHistogram->bounds->data.F32[LQBinNum + 2]); - - y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[LQBinNum - 1]; - y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[LQBinNum]; - y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[LQBinNum + 1]; - - if (!((y->data.F64[0] <= (countFloat / 4.0)) && ((countFloat / 4.0) <= y->data.F64[2]))) { - psAbort(__func__, "p_psVectorRobustStats(2): midpoint not within y-range\n"); - } - - yErr->data.F64[0] = 1.0; - yErr->data.F64[1] = 1.0; - yErr->data.F64[2] = 1.0; - - myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr); - stats->robustLQ = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], countFloat / 4.0); - - } else { - // If the LQ is the first/last histogram bin, then simply use - // the midpoint of that bin. - stats->robustLQ = 0.5 * (robustHistogram->bounds->data.F32[LQBinNum + 1] + - robustHistogram->bounds->data.F32[LQBinNum]); - } - - if ((UQBinNum > 0) && (UQBinNum < (robustHistogram->nums->n - 1))) { - x->data.F64[0] = (double)0.5 * - (robustHistogram->bounds->data.F32[UQBinNum - 1] + - robustHistogram->bounds->data.F32[UQBinNum]); - x->data.F64[1] = (double)0.5 * - (robustHistogram->bounds->data.F32[UQBinNum] + - robustHistogram->bounds->data.F32[UQBinNum + 1]); - x->data.F64[2] = (double)0.5 * - (robustHistogram->bounds->data.F32[UQBinNum + 1] + - robustHistogram->bounds->data.F32[UQBinNum + 2]); - - y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[UQBinNum - 1]; - y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[UQBinNum]; - y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[UQBinNum + 1]; - - if (!((y->data.F64[0] <= (3.0 * countFloat / 4.0)) && - ((3.0 * countFloat / 4.0) <= y->data.F64[2]))) { - psAbort(__func__, "p_psVectorRobustStats(3): midpoint not within y-range\n"); - } - - yErr->data.F64[0] = 1.0; - yErr->data.F64[1] = 1.0; - yErr->data.F64[2] = 1.0; - - myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr); - stats->robustUQ = p_ps1DPolyMedian(myPoly, - x->data.F64[0], x->data.F64[2], 3.0 * countFloat / 4.0); - } else { - // If the UQ is the first/last histogram bin, then simply use - // the midpoint of that bin. - stats->robustUQ = 0.5 * (robustHistogram->bounds->data.F32[UQBinNum + 1] + - robustHistogram->bounds->data.F32[UQBinNum]); - } - } + countFloat = cumulativeRobustSums->data.F32[robustHistogramVector->n - 1]; + + stats->robustLQ = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds, + cumulativeRobustSums, + LQBinNum, + countFloat/4.0); + stats->robustUQ = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds, + cumulativeRobustSums, + UQBinNum, + 3.0 * countFloat/4.0); + } + // XXX: I think sumNfit == sumN50 here. stats->robustNfit = sumNfit; stats->robustN50 = sumN50; - psFree(x); - psFree(y); - psFree(yErr); psFree(tmpStats); psFree(robustHistogram); - psFree(myPoly); - psFree(cumulativeRobustSumsFullRange); - psFree(cumulativeRobustSumsDlRange); + psFree(robustHistogramVector); + psFree(cumulativeRobustSums); } @@ -1468,6 +1350,4 @@ // I am commenting this code out, and replacing it with code which // calculates the mean directly on the robustHistogram. - // - // XXX: Replace this with the LaGrange interpolation stuff. domain = psImageAlloc(1, robustHistogramVector->n, PS_TYPE_F32); @@ -1479,6 +1359,5 @@ max_pos = -1; for (i=0;in;i++) { - domain->data.F32[i][0] = 0.5 * (robustHistogram->bounds->data.F32[i+1] + - robustHistogram->bounds->data.F32[i]); + domain->data.F32[i][0] = PS_BIN_MIDPOINT(robustHistogram, i);bounds->data.F32[modeBinNum + 1]); data->data.F32[i] = (float) robustHistogramVector->data.F32[i]; errors->data.F32[i] = 1.0; @@ -1695,4 +1574,6 @@ psScalar tmpScalar; tmpScalar.type.type = PS_TYPE_F32; + psVector* inF32; + int mustFreeTmp = 1; // NOTE: Verify that this is the correct action. @@ -1706,5 +1587,5 @@ if (out->nums->type.type != PS_TYPE_U32) { - psAbort(__func__, "Only data type PS_TYPE_U32 for output.nums member."); + psAbort(__func__, "Only data type PS_TYPE_U32 for out->nums member."); } // NOTE: Verify that this is the correct action. @@ -1713,15 +1594,17 @@ } - if (in->type.type != PS_TYPE_F32) { - psAbort(__func__, "Only data type PS_TYPE_F32 is currently supported (0x%x).", in->type.type); - } - if (mask != NULL) { if (in->n != mask->n) { - psAbort(__func__, "Vector data and vector mask are of different sizes."); + psError(__func__, "Vector data and vector mask are of different sizes."); } if (mask->type.type != PS_TYPE_U8) { - psAbort(__func__, "Vector mask must be type PS_TYPE_U8"); - } + psError(__func__, "Vector mask must be type PS_TYPE_U8"); + } + } + + inF32 = p_psConvertToF32((psVector *) in); + if (inF32 == NULL) { + inF32 = (psVector *) in; + mustFreeTmp = 0; } @@ -1730,14 +1613,14 @@ numBins = out->nums->n; - for (i = 0; i < in->n; i++) { + for (i = 0; i < inF32->n; i++) { // Check if this pixel is masked, and if so, skip it. if ((mask == NULL) || ((mask != NULL) && (!(mask->data.U8[i] & maskVal)))) { // Check if this pixel is below the minimum value, and if so // count it, then skip it. - if (in->data.F32[i] < out->bounds->data.F32[0]) { + if (inF32->data.F32[i] < out->bounds->data.F32[0]) { out->minNum++; // Check if this pixel is above the maximum value, and if so // count it, then skip it. - } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) { + } else if (inF32->data.F32[i] > out->bounds->data.F32[numBins]) { out->maxNum++; } else { @@ -1746,5 +1629,5 @@ if (out->uniform == true) { binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0]; - binNum = (int)((in->data.F32[i] - out->bounds->data.F32[0]) / binSize); + binNum = (int)((inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize); // NOTE: This next if-statement really shouldn't be necessary. @@ -1755,16 +1638,15 @@ } - (out->nums->data.S32[binNum])++; + (out->nums->data.U32[binNum])++; // If this is a non-uniform histogram, determining the correct // bin number requires a bit more work. } else { - tmpScalar.data.F32 = in->data.F32[i]; - tmp = p_psVectorBinDisect(out->bounds, - &tmpScalar); - if (tmp == -1) { - // XXX: Generate warning message + tmpScalar.data.F32 = inF32->data.F32[i]; + tmp = p_psVectorBinDisect(out->bounds, &tmpScalar); + if (tmp < 0) { + // XXX: Generate warning message } else { - (out->nums->data.S32[tmp])++; + (out->nums->data.U32[tmp])++; } } @@ -1772,19 +1654,22 @@ } } + + if (mustFreeTmp == 1) { + psFree(inF32); + } return (out); } /****************************************************************************** -p_psConvertToF32(myVector, maskVector, maskVal, stats): this is the cheap -way to support a variety of vector data types: we simply convert the input -vector to F32 at the beginning, and write all of our functions in F32. If -the vast majority of all vector stat operations are F32 (or any other single -type), then this is probably the best way to go. Otherwise, when the -algorithms stablize, we will then macro everything and put type support in -the various stat functions. +p_psConvertToF32(in): this is the cheap way to support a variety of vector +data types: we simply convert the input vector to F32 at the beginning, and +write all of our functions in F32. If the vast majority of all vector stat +operations are F32 (or any other single type), then this is probably the +best way to go. Otherwise, when the algorithms stablize, we will then macro +everything and put type support in the various stat functions. XXX: Should the default data type be F64? Since we are buying Athlons... *****************************************************************************/ -psVector* p_psConvertToF32(psStats* stats, psVector* in, psVector* mask, unsigned int maskVal) +psVector* p_psConvertToF32(psVector* in) { int i = 0; @@ -1800,4 +1685,9 @@ for (i = 0; i < in->n; i++) { tmp->data.F32[i] = (float)in->data.U16[i]; + } + } else if (in->type.type == PS_TYPE_U8) { + tmp = psVectorAlloc(in->n, PS_TYPE_F32); + for (i = 0; i < in->n; i++) { + tmp->data.F32[i] = (float)in->data.U8[i]; } } else if (in->type.type == PS_TYPE_F64) { @@ -1839,20 +1729,18 @@ int mustFreeTmp = 1; - // NOTE: Verify that this is the correct action. if (in == NULL) { - return (stats); - } - if (stats == NULL) { - return (NULL); - } - - inF32 = p_psConvertToF32(stats, in, mask, maskVal); + psError(__func__, "in is NULL\n"); + return(stats); + } + PS_CHECK_NULL_PTR_RETURN_NULL(stats); + + inF32 = p_psConvertToF32((psVector *) in); if (inF32 == NULL) { inF32 = in; mustFreeTmp = 0; } - // XXX: Should we abort if (stats->min == stats->max)? if ((stats->options & PS_STAT_USE_RANGE) && (stats->min >= stats->max)) { - psAbort(__func__, "psVectorStats() called with range: %f to %f\n", stats->min, stats->max); + psError(__func__, "psVectorStats() called with range: %f to %f\n", stats->min, stats->max); + return(stats); } Index: /trunk/psLib/src/math/psConstants.h =================================================================== --- /trunk/psLib/src/math/psConstants.h (revision 2196) +++ /trunk/psLib/src/math/psConstants.h (revision 2197) @@ -6,6 +6,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-14 23:43:12 $ + * @version $Revision: 1.13 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-26 21:24:42 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -107,4 +107,11 @@ } \ +/** Preprocessor macro to generate error on a NULL vector */ +#define PS_CHECK_NULL_VECTOR_RETURN_F0(NAME) \ +if (NAME == NULL || NAME->data.V == NULL) { \ + psError(__func__,"Unallowable operation: %s or its data is NULL.", #NAME); \ + return(0.0); \ +} \ + /** Preprocessor macro to generate error on a NULL poniter */ #define PS_CHECK_NULL_PTR_RETURN_NULL(NAME) \ @@ -113,4 +120,12 @@ return(NULL); \ } + +/** Preprocessor macro to generate error on a NULL poniter */ +#define PS_CHECK_NULL_PTR_RETURN_F0(NAME) \ +if (NAME == NULL) { \ + psError(__func__,"Unallowable operation: %s is NULL.", #NAME); \ + return(0.0); \ +} + /** Preprocessor macro to generate error on a NULL poniter */ #define PS_CHECK_NULL_PTR_RETURN_ZERO(NAME) \ @@ -127,4 +142,11 @@ } \ +/** Preprocessor macro to generate error for zero length vector */ +#define PS_CHECK_EMPTY_VECTOR_RETURN_F0(NAME) \ +if (NAME->n < 1) { \ + psError(__func__,"Unallowable operation: %s has zero n value.", #NAME); \ + return(0.0); \ +} \ + /** Preprocessor macro to generate error on differing size vectors */ #define PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(VEC1, VEC2) \ @@ -132,4 +154,11 @@ psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \ return(NULL); \ +} + +/** Preprocessor macro to generate error on differing size vectors */ +#define PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_F0(VEC1, VEC2) \ +if (VEC1->n != VEC2->n) { \ + psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \ + return(0.0); \ } @@ -199,7 +228,7 @@ #define PS_MAX(A, B) \ -((A > B) ? A : B) \ +(((A) > (B)) ? (A) : (B)) \ #define PS_MIN(A, B) \ -((A < B) ? A : B) \ - +(((A) < (B)) ? (A) : (B)) \ + Index: /trunk/psLib/src/math/psMinimize.c =================================================================== --- /trunk/psLib/src/math/psMinimize.c (revision 2196) +++ /trunk/psLib/src/math/psMinimize.c (revision 2197) @@ -9,6 +9,6 @@ * @author GLF, MHPCC * - * @version $Revision: 1.58 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-14 01:44:48 $ + * @version $Revision: 1.59 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-26 21:24:42 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -636,4 +636,261 @@ psFree(currValueVec); psFree(newValueVec); + + if ((min->iter < min->maxIter) || + (min->lastDelta <= min->tol)) { + return(true); + } + + psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + "---- psMinimizeLMChi2() end (false) ----\n"); + return(false); +} + +/****************************************************************************** +psMyMinimizeLMChi2(): This routine will take an procedure which calculates +an arbitrary function and it's derivative and minimize the chi-squared match +between that function at the specified coords and the specified value at +those coords. + +XXX: Do this: + After checking that all entries in the paramMask are 1 or 0, when + forming the A matrix from alpha, try this: + + A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i]; + +XXX: This is very different from what is specified in the SDR. Must +coordinate with IfA on new SDR. + +XXX: Do vector/image recycles. + +XXX: probably yErr will be part of the SDR. + +XXX: This must work for both F32 and F64. F32 is currently implemented. + Note: since the LUD routines are only implemented in F64, then we + will have to convert all F32 input vectors to F64 regardless. So, + the F64 port might be. + *****************************************************************************/ +bool psMyMinimizeLMChi2(psMinimization *min, + psImage *covar, + psVector *params, + const psVector *paramMask, + const psArray *x, + const psVector *y, + const psVector *yErr, + psMyMinimizeLMChi2Func func) +{ + PS_CHECK_NULL_PTR_RETURN_NULL(min); + PS_CHECK_NULL_VECTOR_RETURN_NULL(params); + PS_CHECK_EMPTY_VECTOR_RETURN_NULL(params); + PS_CHECK_NULL_PTR_RETURN_NULL(x); + PS_CHECK_NULL_VECTOR_RETURN_NULL(y); + PS_CHECK_EMPTY_VECTOR_RETURN_NULL(y); + PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(x, y); + if (paramMask != NULL) { + PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(params, paramMask); + } + if (yErr != NULL) { + PS_CHECK_VECTOR_SIZE_EQUAL_RETURN_NULL(y, yErr); + } + + psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + "---- psMinimizeLMChi2() begin ----\n"); + int numData = y->n; + int numParams = params->n; + int i; + int j; + int k; + int l; + int n; + int p; + psVector *beta = psVectorAlloc(numParams, PS_TYPE_F64); + psVector *perm = psVectorAlloc(numParams, PS_TYPE_F64); + + psVector *paramDeltasF64 = psVectorAlloc(numParams, PS_TYPE_F64); + psVector *origParams = psVectorAlloc(numParams, PS_TYPE_F32); + psVector *newParams = psVectorAlloc(numParams, PS_TYPE_F32); + + psImage *alpha = psImageAlloc(numParams, numParams, PS_TYPE_F32); + psImage *A = psImageAlloc(numParams, numParams, PS_TYPE_F64); + psImage *aOut = psImageAlloc(numParams, numParams, PS_TYPE_F64); + + // psVector **deriv = (psVector **) psAlloc(numData * sizeof(psVector *)); + // for (i=0;imaxIter is %d\n", min->maxIter); + psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + "min->tol is %f\n", min->tol); + + for (p=0;pdata.F32[p] = params->data.F32[p]; + } + + min->lastDelta = HUGE; + min->iter = 0; + + while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) { + psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + "------------------------------------------------------\n"); + psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + "Iteration %d. Delta is %f\n", min->iter, min->lastDelta); + + // + // Calculate the current values and chi-squared of the function. + // + currChi2 = 0.0; + currValueVec = func(deriv, params, x); + for (n=0;ndata.F32[n] * currValueVec->data.F32[n]); + } + + for (p=0;pdata.F32[%d] is %f.\n", p, params->data.F32[p]); + } + psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + "Current chi-squared is (%f)\n", currChi2); + + // + // Mask elements of the derivative for each data point. + // + for (p=0;pdata.U8[p] != 0)) { + for (n=0;ndata.F32[n][p] = 0.0; + } + } + } + + // + // Calculate the BETA vector. + // + for (p=0;pdata.F64[p] = 0.0; + for (n=0;ndata.F64[p])+= + (y->data.F32[n] - currValueVec->data.F32[n]) * + deriv->data.F32[n][p]; + } + // XXX: multiple by -1 here? + (beta->data.F64[p])*= -1.0; + psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + "beta->data.F64[%d] is %f.\n", p, beta->data.F64[p]); + } + psFree(currValueVec); + + // + // Calculate the ALPHA matrix. + // + for (k=0;kdata.F32[k][l] = 0.0; + for (n=0;ndata.F32[k][l]+= deriv->data.F32[n][k] * + deriv->data.F32[n][l]; + } + } + } + + // + // Calculate the matrix A. + // + for (j=0;jdata.F64[j][k] = + (double) ((1.0 + lamda) * alpha->data.F32[j][k]); + } else { + A->data.F64[j][k] = (double) alpha->data.F32[j][k]; + } + } + } + for (j=0;jdata.F64[j][k]); + } + psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, "Matrix A[][]:\n"); + } + + // + // Solve A * alpha = Beta + // + aOut = psMatrixLUD(aOut, perm, A); + paramDeltasF64 = psMatrixLUSolve(paramDeltasF64, aOut, beta, perm); + + // + // Mask any masked parameters. + // + for (i=0;idata.F64[%d] is %f.\n", i, paramDeltasF64->data.F64[i]); + if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) { + newParams->data.F32[i] = origParams->data.F32[i]; + } else { + newParams->data.F32[i] = params->data.F32[i] - + (float) paramDeltasF64->data.F64[i]; + } + } + + psTrace(".psLib.dataManip.psMinimizeLMChi2", 6, + "Calling func() with new parameters:\n"); + for (i=0;idata.F32[%d] is %f.\n", i, newParams->data.F32[i]); + } + + + // + // Calculate new function values. + // + newChi2 = 0.0; + newValueVec = func(deriv, newParams, x); + for (n=0;ndata.F32[n] * newValueVec->data.F32[n]); + } + psFree(newValueVec); + + psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + "old/new chi-squareds are (%f, %f)\n", currChi2, newChi2); + + // + // If the new chi-squared is lower, then keep it. + // + if (currChi2 > newChi2) { + min->lastDelta = currChi2 - newChi2; + min->value = newChi2; + + // We already masked params. + for (i=0;idata.F32[i] = (float) newParams->data.F32[i]; + } + lamda*= 0.1; + } else { + lamda*= 10.0; + } + psTrace(".psLib.dataManip.psMinimizeLMChi2", 4, + "lamda is %f\n", lamda); + min->iter++; + } + psFree(beta); + psFree(perm); + psFree(paramDeltasF64); + psFree(origParams); + psFree(newParams); + psFree(alpha); + psFree(A); + psFree(aOut); + psFree(deriv); if ((min->iter < min->maxIter) || Index: /trunk/psLib/src/math/psMinimize.h =================================================================== --- /trunk/psLib/src/math/psMinimize.h (revision 2196) +++ /trunk/psLib/src/math/psMinimize.h (revision 2197) @@ -8,6 +8,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.29 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-06 23:38:49 $ + * @version $Revision: 1.30 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-26 21:24:43 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -84,5 +84,82 @@ +/* + Bug 203: + +typedef float (*psMinimizeLMChi2Func)(psVector *deriv, const psVector *params, +const psVector *x); + +This function takes the current guess for the parameters for which we are +trying to get the best values (params), and a single vector (x) of +conditions from the array of x values fed into the minimiser. + + bool psMinimizeLMChi2(psMinimization *min, + psImage *covar, + psVector *params, + const psVector *paramMask, + const psArray *x, + const psVector *y, + const psVector *yErr, + psMinimizeLMChi2Func func); + +This takes the minimization specs (min), returns the covariance matrix +(covar), takes the best guess of initial parameters (params), the parameter +mask (paramMask), and takes multiple vectors of conditions in an array (x), +the corresponding measured values (y) and errors (yErr) and the function to +fit (func). + +For example, for GRB afterglows, I have flux as a function of time and +frequency. So I stuff into the "psArray *x" all my time and frequency values +(so a whole heap of vectors of size 2), I have the measured values in "y" and +errors in "yErr". Then each of the time-frequency pairs are passed to my model +function with the current parameters, and the model function returns the flux +for that time-frequency pair, and the derivative with respect to each of the +parameters. + +This seems reasonable. The only thing we could change would be to have the +function be defined: + +typedef psVector* (*psMinimizeLMChi2Func)(psImage *deriv, + const psVector *params, + const psArray *x); + +So it would return the model value for each of the measurements at once, and +return for each the derivatives (so that it returns a matrix). + +What do you think? + + + +I'm not sure I understand how LM chi-squared minimization will work with the +following function that you define: + + typedef psVector* (*psMinimizeLMChi2Func)(psImage *deriv, + const psVector *params, + const psArray *x); + +The chi-squared minimization algorithm, as defined in NR, requires that +function to be minimized be evaluated at each data point, and that all +derivatives, with respect to each parameter, be calculated at each data +point. In the above, can I assume that + + x is an array of psVectors, with each vector corresponding to a single + data point. + + The returned value has the same length as x. It contains the value of + the function at each data point in x. + + deriv: an n-by-p matrix where "n" is the number of data points, and "p" + is the number of parameters. The [i][j] element of this matrix + holds the derivative of the function at the i-th data point with + respect to the j-th parameter. + + + +*/ + + // XXX: What if any of these arguments are NULL? + + typedef float (*psMinimizeLMChi2Func)(psVector *deriv, @@ -98,4 +175,19 @@ psMinimizeLMChi2Func func); + +typedef +psVector* (*psMyMinimizeLMChi2Func)(psImage *deriv, + const psVector *params, + const psArray *x); + +bool psMyMinimizeLMChi2(psMinimization *min, + psImage *covar, + psVector *params, + const psVector *paramMask, + const psArray *x, + const psVector *y, + const psVector *yErr, + psMyMinimizeLMChi2Func func); + typedef float (*psMinimizePowellFunc)(const psVector *params, Index: /trunk/psLib/src/math/psPolynomial.c =================================================================== --- /trunk/psLib/src/math/psPolynomial.c (revision 2196) +++ /trunk/psLib/src/math/psPolynomial.c (revision 2197) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.52 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-14 01:44:48 $ + * @version $Revision: 1.53 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-26 21:24:42 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii Index: /trunk/psLib/src/math/psSpline.c =================================================================== --- /trunk/psLib/src/math/psSpline.c (revision 2196) +++ /trunk/psLib/src/math/psSpline.c (revision 2197) @@ -7,6 +7,6 @@ * polynomials. It also contains a Gaussian functions. * - * @version $Revision: 1.52 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-14 01:44:48 $ + * @version $Revision: 1.53 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-26 21:24:42 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii Index: /trunk/psLib/src/math/psStats.c =================================================================== --- /trunk/psLib/src/math/psStats.c (revision 2196) +++ /trunk/psLib/src/math/psStats.c (revision 2197) @@ -9,6 +9,6 @@ * @author GLG, MHPCC * - * @version $Revision: 1.67 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-10-13 02:40:13 $ + * @version $Revision: 1.68 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-10-26 21:24:43 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -40,5 +40,5 @@ /*****************************************************************************/ // will use robust statistical methods. -#define GAUSS_WIDTH 20 // The width of the Gaussian or boxcar smoothing. +#define GAUSS_WIDTH 5 // The width of the Gaussian or boxcar smoothing. #define CLIPPED_NUM_ITER_LB 1 #define CLIPPED_NUM_ITER_UB 10 @@ -57,4 +57,5 @@ /* TYPE DEFINITIONS */ /*****************************************************************************/ +psVector* p_psConvertToF32(psVector* in); /*****************************************************************************/ @@ -423,22 +424,4 @@ float rangeMax = 0.0; // Exclude date above this - // Determine if the number of data points exceed a threshold which will - // cause to generate robust stats, as opposed to exact stats. - // XXX: This code is no longer used. We now calculate the median exactly - // regardless of the vector size. - /* - * // Calculate the robust quartiles. stats2 = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); - * p_psVectorRobustStats(myVector, maskVector, maskVal, stats2); - * - * // Store the robust quartiles into the sample quartile members. stats->sampleMedian = - * stats2->robustMedian; - * - * // Free temporary data buffers. psFree(stats2); - * - * // Set the PS_STAT_ROBUST_FOR_SAMPLE bit in the stats structure. stats->options = stats->options | - * PS_STAT_ROBUST_FOR_SAMPLE; - * - * return; } */ - // Determine how many data points fit inside this min/max range // and are not masked, IF the maskVector is not NULL> @@ -507,67 +490,66 @@ /****************************************************************************** - This routine smoothes the data in the input robustHistogram with a - Gaussian of width sigma. - - XXX: Must consult with IfA on the proper values for the Gaussian - smoothing. Currently, the Gaussian coefficients are such that only - the center coefficient is non-zero. This is because "e" is being - raised to a very large negative number for all points except the - center point. - - XXX: use a static variable for gaussianCoefs[] and compute them once. - *****************************************************************************/ -psVector* p_psVectorsmoothHistGaussian(psHistogram* robustHistogram, +p_psVectorSmoothHistGaussian(): This routine smoothes the data in the input +robustHistogram with a Gaussian of width sigma. + *****************************************************************************/ +#define GAUSS_WIDTH_OTHER 5.0 +psVector* p_psVectorSmoothHistGaussian(psHistogram* robustHistogram, float sigma) { int i = 0; // Loop index variable int j = 0; // Loop index variable - float denom = 0.0; // Temporary variable - float expo = 0.0; // Temporary variable - float gaussianCoefs[1 + (2 * GAUSS_WIDTH)]; // The Gaussian Coefficients - psVector* smooth = psVectorAlloc(robustHistogram->nums->n, PS_TYPE_F32); - - for (i = 0; i < (1 + (2 * GAUSS_WIDTH)); i++) { - if (fabs(sigma) >= FLT_EPSILON) { - // If sigma does not equal zero, then we use Gaussian smoothing. - #ifdef DARWIN - denom = (float)sqrt(2.0 * M_PI * sigma * sigma); - #else - - denom = sqrtf(2.0 * M_PI * sigma * sigma); - #endif - - expo = -(float)((i - GAUSS_WIDTH) * (i - GAUSS_WIDTH)); - expo /= (2.0 * sigma * sigma); - gaussianCoefs[i] = exp(expo / denom); - - //printf("p_psVectorsmoothHistGaussian(): gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]); - // NOTE: Gaussian smoothing just isn't working with low sigma - // values. The problem is that the Gaussian coefficients are - // all zero, except for the middle coefficient, which is exactly - // one. Therefore, I'm using boxcar smoothing. - gaussianCoefs[i] = 1.0 / (1.0 + (2.0 * (float)GAUSS_WIDTH)); - // printf("gaussianCoefs[%d] is %f\n", i, gaussianCoefs[i]); - } else { - /* If sigma equals zero (all pixels have the same value) the above code will divide by zero. - * Therefore, we don't need to smooth the data. */ - for (i = 0; i < robustHistogram->nums->n; i++) { - smooth->data.F32[i] = (float)robustHistogram->nums->data.S32[i]; - } - return (smooth); - } - } - - // Perform the actual smoothing. - for (i = 0; i < robustHistogram->nums->n; i++) { + float iMid; + float jMid; + int numBins = robustHistogram->nums->n; + int numBounds = robustHistogram->bounds->n; + psVector* smooth = psVectorAlloc(numBins, PS_TYPE_F32); + int jMin = 0; + int jMax = 0; + float firstBound = robustHistogram->bounds->data.F32[0]; + float lastBound = robustHistogram->bounds->data.F32[numBounds-1]; + psScalar x; + + x.type.type = PS_TYPE_F32; + for (i = 0; i < numBins; i++) { + iMid = (robustHistogram->bounds->data.F32[i] + + robustHistogram->bounds->data.F32[i+1]) / 2.0; + + // YYY: The p_psVectorBinDisect() routine does much of the work of + // the following conditionals, however, it also reports a warning + // message. I don't want the warning message so I reproduce the + // conditionals here. Maybe p_psVectorBinDisect() should not produce + // warnings? + + x.data.F32 = iMid - (GAUSS_WIDTH * sigma); + if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { + jMin = p_psVectorBinDisect(robustHistogram->bounds, &x); + } else if (x.data.F32 <= firstBound) { + jMin = 0; + } else if (x.data.F32 >= lastBound) { + jMin = robustHistogram->bounds->n - 1; + + } + + x.data.F32 = iMid + (GAUSS_WIDTH * sigma); + if ((x.data.F32 >= firstBound) && (x.data.F32 <= lastBound)) { + jMax = p_psVectorBinDisect(robustHistogram->bounds, &x); + } else if (x.data.F32 <= firstBound) { + jMax = 0; + } else if (x.data.F32 >= lastBound) { + jMax = robustHistogram->bounds->n - 1; + + } + smooth->data.F32[i] = 0.0; - for (j = -GAUSS_WIDTH; j <= +GAUSS_WIDTH; j++) { - if (((j + i) >= 0) && ((j + i) < smooth->n)) { - smooth->data.F32[i] += (gaussianCoefs[j + GAUSS_WIDTH] * - (float)robustHistogram->nums->data.S32[j + i]); - } - } - } - return (smooth); + for (j = jMin ; j <= jMax ; j++) { + jMid = (robustHistogram->bounds->data.F32[j] + + robustHistogram->bounds->data.F32[j+1]) / 2.0; + smooth->data.F32[i] += + ((float) robustHistogram->nums->data.U32[j]) * + psGaussian(jMid, iMid, sigma, true); + } + } + + return(smooth); } @@ -913,5 +895,6 @@ specified data point. *****************************************************************************/ -float p_psGaussian(const psVector* restrict myData, const psVector* restrict myParams) +float p_psGaussian(const psVector* restrict myData, + const psVector* restrict myParams) { float x = myData->data.F32[0]; @@ -930,5 +913,6 @@ calculates the specified partial derivative of the above Gaussian function. *****************************************************************************/ -float p_psGaussianDeriv(const psVector* restrict myData, const psVector* restrict myParams, int whichParam) +float p_psGaussianDeriv(const psVector* restrict myData, + const psVector* restrict myParams, int whichParam) { float x = myData->data.F32[0]; @@ -1036,45 +1020,75 @@ /****************************************************************************** -p_psFitQuadratic(robustHistogram. cumulativeSums, binNum, fitFloat): given -the specified histogram, with the specified pre-calculated cumulativeSums, -this routine fits a quadratic f(x) to the 3 bins surrounding bin "binNum", -and then finds the point x such that f(x) == fitFloat. +fitQuadraticSearchForYThenReturnX(*xVec, *yVec, binNum, yVal): This routine +takes psVectors of x/y pairs as input, and fits a quadratic to the 3 points +surrounding element binNum in the vectors (the midpoint between element i +and i+1 is used for x[i]). It then determines for what value x does that +quadratic f(x) = yVal (the input parameter). -XXX: This function is currently not being used. - *****************************************************************************/ -float p_psFitQuadratic(psHistogram* histogram, - psVector* cumulativeSums, - int binNum, - float fitFloat) -{ +// XXX: Use static variables. +*****************************************************************************/ +float fitQuadraticSearchForYThenReturnX(psVector *xVec, + psVector *yVec, + int binNum, + float yVal) +{ + // static psVector* x = NULL; + // static psVector* y = NULL; + // static psVector* yErr = NULL; psVector* x = psVectorAlloc(3, PS_TYPE_F64); psVector* y = psVectorAlloc(3, PS_TYPE_F64); psVector* yErr = psVectorAlloc(3, PS_TYPE_F64); + + /* + if (x == NULL) { + x = psVectorAlloc(3, PS_TYPE_F64); + p_psMemSetPersistent(x, true); + } + if (y == NULL) { + y = psVectorAlloc(3, PS_TYPE_F64); + p_psMemSetPersistent(y, true); + } + if (yErr == NULL) { + yErr = psVectorAlloc(3, PS_TYPE_F64); + p_psMemSetPersistent(yErr, true); + } + */ psPolynomial1D* myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); - - if ((binNum > 0) && (binNum < (histogram->nums->n + 1))) { - x->data.F64[0] = (double)0.5 * - (histogram->bounds->data.F32[binNum - 1] + histogram->bounds->data.F32[binNum]); - x->data.F64[1] = (double)0.5 * - (histogram->bounds->data.F32[binNum] + histogram->bounds->data.F32[binNum + 1]); - x->data.F64[2] = (double)0.5 * - (histogram->bounds->data.F32[binNum + 1] + histogram->bounds->data.F32[binNum + 2]); - - y->data.F64[0] = cumulativeSums->data.F32[binNum - 1]; - y->data.F64[1] = cumulativeSums->data.F32[binNum]; - y->data.F64[2] = cumulativeSums->data.F32[binNum + 1]; - - if (!((y->data.F64[0] <= fitFloat) && (fitFloat <= y->data.F64[2]))) { - psAbort(__func__, "p_psVectorRobustStats(0): midpoint not within y-range\n"); - } - + float tmpFloat; + + if ((binNum > 0) && (binNum < (yVec->n - 2))) { + x->data.F64[0] = (double) (0.5 * (xVec->data.F32[binNum - 1] + xVec->data.F32[binNum])); + x->data.F64[1] = (double) (0.5 * (xVec->data.F32[binNum] + xVec->data.F32[binNum+1])); + x->data.F64[2] = (double) (0.5 * (xVec->data.F32[binNum+1] + xVec->data.F32[binNum+2])); + y->data.F64[0] = yVec->data.F32[binNum - 1]; + y->data.F64[1] = yVec->data.F32[binNum]; + y->data.F64[2] = yVec->data.F32[binNum + 1]; + + // Ensure that yVal is within the range of the bins we are using. + if (!((y->data.F64[0] <= yVal) && (yVal <= y->data.F64[2]))) { + printf("((%f), %f, %f)\n", yVal, y->data.F64[0], y->data.F64[2]); + psError(__func__, "yVal not within y-range\n"); + } yErr->data.F64[0] = 1.0; yErr->data.F64[1] = 1.0; yErr->data.F64[2] = 1.0; + // Determine the coefficients of the polynomial. myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr); - return (p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], fitFloat)); + // Call p_ps1DPolyMedian(), which does a binary search on the + // polynomial, looking for the value x such that f(x) = yVal + tmpFloat = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], yVal); } else { - return (0.5 * (histogram->bounds->data.F32[binNum + 1] + histogram->bounds->data.F32[binNum])); + if (binNum == 0) { + tmpFloat = 0.5 * (xVec->data.F32[binNum] + + xVec->data.F32[binNum + 1]); + } else if (binNum == (xVec->n - 1)) { + // XXX: Is this right? + tmpFloat = xVec->data.F32[binNum]; + } else if (binNum == (xVec->n - 2)) { + // XXX: Is this right? + tmpFloat = 0.5 * (xVec->data.F32[binNum] + + xVec->data.F32[binNum + 1]); + } } @@ -1083,6 +1097,9 @@ psFree(yErr); psFree(myPoly); - return (0.0); -} + return(tmpFloat); +} + +#define PS_BIN_MIDPOINT(HISTOGRAM, BIN_NUM) \ +(0.5 * (HISTOGRAM->bounds->data.F32[(BIN_NUM)] + HISTOGRAM->bounds->data.F32[(BIN_NUM)+1])) /****************************************************************************** @@ -1111,5 +1128,7 @@ *****************************************************************************/ void p_psVectorRobustStats(const psVector* restrict myVector, - const psVector* restrict maskVector, unsigned int maskVal, psStats* stats) + const psVector* restrict maskVector, + unsigned int maskVal, + psStats* stats) { psHistogram* robustHistogram = NULL; @@ -1118,19 +1137,10 @@ int LQBinNum = -1; // Bin num for lower quartile int UQBinNum = -1; // Bin num for upper quartile + int medianBinNum = -1; int i = 0; // Loop index variable - int maxBinNum = 0; - float maxBinCount = 0.0; + int modeBinNum = 0; + float modeBinCount = 0.0; float dL = 0.0; int numBins = 0; - psStats* tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN); - - // psImage* domain; - // psVector* errors; - // psVector* data; - // psVector* initialGuess; - // psVector* theParams; - // float chiSq=0.0; - // float max = -HUGE; - // float max_pos; float myMean = 0.0; float myStdev = 0.0; @@ -1139,14 +1149,11 @@ float sumSquares = 0.0; float sumDiffs = 0.0; - psVector* x = psVectorAlloc(3, PS_TYPE_F64); - psVector* y = psVectorAlloc(3, PS_TYPE_F64); - psVector* yErr = psVectorAlloc(3, PS_TYPE_F64); - psPolynomial1D* myPoly = psPolynomial1DAlloc(2, PS_POLYNOMIAL_ORD); - psVector* cumulativeRobustSumsFullRange = NULL; - psVector* cumulativeRobustSumsDlRange = NULL; + psVector* cumulativeRobustSums = NULL; float sumRobust = 0.0; float sumN50 = 0.0; float sumNfit = 0.0; - float cumulativeMedian = 0.0; + psScalar tmpScalar; + tmpScalar.type.type = PS_TYPE_F32; + psStats* tmpStats = psStatsAlloc(PS_STAT_CLIPPED_STDEV | PS_STAT_CLIPPED_MEAN); // Compute the initial bin size of the robust histogram. This is done @@ -1188,4 +1195,5 @@ p_psVectorMax(myVector, maskVector, maskVal, stats); } + // Create the histogram structure. NOTE: we can not specify the bin size // precisely since the argument to psHistogramAlloc() is the number of @@ -1198,7 +1206,6 @@ psVectorHistogram(robustHistogram, myVector, maskVector, maskVal); - // Smooth the histogram. - // XXX: is that the right stdev? - robustHistogramVector = p_psVectorsmoothHistGaussian(robustHistogram, + // Smooth the histogram, Gaussian-style. + robustHistogramVector = p_psVectorSmoothHistGaussian(robustHistogram, tmpStats->clippedStdev / 4.0f); @@ -1208,87 +1215,62 @@ /************************************************************************** - Determine the lower/upper quartiles. + Determine the median/lower/upper quartile bin numbers. + + We define a vector called "cumulativeRobustSums" where the value at + index position i is equal to the sum of bins 0:i. This will be used in + determining the median and lower/upper quartiles. **************************************************************************/ - // We define a vector called "cumulativeRobustSums..." where the value at - // index position i is equal to the sum of bins 0:i. This will be used - // now and later in determining the lower/upper quartiles. - cumulativeRobustSumsFullRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); - cumulativeRobustSumsFullRange->data.F32[0] = robustHistogramVector->data.F32[0]; + cumulativeRobustSums = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); + cumulativeRobustSums->data.F32[0] = robustHistogramVector->data.F32[0]; for (i = 1; i < robustHistogramVector->n; i++) { - cumulativeRobustSumsFullRange->data.F32[i] = - cumulativeRobustSumsFullRange->data.F32[i - 1] + robustHistogramVector->data.F32[i]; - } - sumRobust = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n - 1]; - - // Determine the bin number containing the lower quartile point. - LQBinNum = -1; - for (i = 0; i < cumulativeRobustSumsFullRange->n; i++) { - if (cumulativeRobustSumsFullRange->data.F32[i] >= (sumRobust / 4.0)) { - LQBinNum = i; - break; - } - } - - // Determine the bin number containing the upper quartile point. - UQBinNum = -1; - for (i = cumulativeRobustSumsFullRange->n - 1; i >= 0; i--) { - if (cumulativeRobustSumsFullRange->data.F32[i] <= (3.0 * sumRobust / 4.0)) { - UQBinNum = i; - break; - } - } - - if ((LQBinNum == -1) || (UQBinNum == -1)) { - psAbort(__func__, "Could not determine the robust lower/upper quartiles."); - } - + cumulativeRobustSums->data.F32[i] = cumulativeRobustSums->data.F32[i - 1] + + robustHistogramVector->data.F32[i]; + } + sumRobust = cumulativeRobustSums->data.F32[robustHistogramVector->n - 1]; + + tmpScalar.data.F32 = sumRobust / 4.0; + LQBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar); + tmpScalar.data.F32 = 3.0 * sumRobust / 4.0; + UQBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar); + tmpScalar.data.F32 = sumRobust / 2.0; + medianBinNum = p_psVectorBinDisect(cumulativeRobustSums, &tmpScalar); + + if ((LQBinNum < 0) || (UQBinNum < 0)) { + psError(__func__, "Could not determine the robust lower/upper quartile bin numbers."); + } + if (medianBinNum < 0) { + psError(__func__, "Could not determine the robust lower/upper quartile bin numbers."); + } /************************************************************************** Determine the mode in the range LQ:UQ. **************************************************************************/ // Determine the bin with the peak value in the range LQ to UQ. - maxBinNum = LQBinNum; - maxBinCount = robustHistogramVector->data.F32[LQBinNum]; - sumN50 = (float)robustHistogram->nums->data.S32[LQBinNum]; + modeBinNum = LQBinNum; + modeBinCount = robustHistogramVector->data.F32[LQBinNum]; + sumN50 = (float)robustHistogram->nums->data.U32[LQBinNum]; for (i = LQBinNum + 1; i <= UQBinNum; i++) { - if (robustHistogramVector->data.F32[i] > maxBinCount) { - maxBinNum = i; - maxBinCount = robustHistogramVector->data.F32[i]; - } - sumN50 += (float)robustHistogram->nums->data.S32[i]; - } - - // XXX: is dL defined as the value at the LQ/UQ, or the bin number? + if (robustHistogramVector->data.F32[i] > modeBinCount) { + modeBinNum = i; + modeBinCount = robustHistogramVector->data.F32[i]; + } + sumN50 += (float)robustHistogram->nums->data.U32[i]; + } + dL = (UQBinNum - LQBinNum) / 4; - - printf("(LQBinNum, UQBinNum, maxBinNum) is (%d, %d, %d)\n", LQBinNum, UQBinNum, maxBinNum); - /************************************************************************** Determine the mean/stdev for the bins in the range mode-dL to mode+dL **************************************************************************/ - cumulativeRobustSumsDlRange = psVectorAlloc(robustHistogramVector->n, PS_TYPE_F32); - for (i = 0; i < robustHistogramVector->n; i++) { - cumulativeRobustSumsDlRange->data.F32[i] = 0.0; - } - sumNfit = 0.0; - cumulativeMedian = 0.0; - for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) { - if ((0 <= i) && (i < robustHistogramVector->n)) { - cumulativeRobustSumsDlRange->data.F32[i] = - cumulativeRobustSumsDlRange->data.F32[i - 1] + robustHistogramVector->data.F32[i]; - cumulativeMedian += robustHistogramVector->data.F32[i]; - sumNfit += (float)robustHistogram->nums->data.S32[i]; - } - } - // Calculate the mean of the smoothed robust histogram in the range // mode-dL to mode+dL. We use the midpoint of each bin as the mean for // that bin (this is a non-exact approximation). + sumNfit = 0.0; myMean = 0.0; - for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) { + for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) { if ((0 <= i) && (i < robustHistogramVector->n)) { - myMean += (robustHistogramVector->data.F32[i]) * 0.5 * - (robustHistogram->bounds->data.F32[i + 1] + robustHistogram->bounds->data.F32[i]); + myMean += (robustHistogramVector->data.F32[i]) * PS_BIN_MIDPOINT(robustHistogram, i); countFloat += robustHistogramVector->data.F32[i]; } + + sumNfit += (float)robustHistogram->nums->data.U32[i]; } myMean /= countFloat; @@ -1297,8 +1279,7 @@ // mode-dL to mode+dL. We use the midpoint of each bin as the mean for // that bin. - for (i = maxBinNum - dL; i <= maxBinNum + dL; i++) { + for (i = modeBinNum - dL; i <= modeBinNum + dL; i++) { if ((0 <= i) && (i < robustHistogramVector->n)) { - diff = (0.5 * (robustHistogram->bounds->data.F32[i + 1] + - robustHistogram->bounds->data.F32[i])) - myMean; + diff = PS_BIN_MIDPOINT(robustHistogram, i) - myMean; sumSquares += diff * diff * robustHistogramVector->data.F32[i]; sumDiffs += diff * robustHistogramVector->data.F32[i]; @@ -1315,6 +1296,5 @@ if (stats->options & PS_STAT_ROBUST_MODE) { - stats->robustMode = 0.5 * - (robustHistogram->bounds->data.F32[maxBinNum] + robustHistogram->bounds->data.F32[maxBinNum + 1]); + stats->robustMode = PS_BIN_MIDPOINT(robustHistogram, modeBinNum); } @@ -1322,140 +1302,42 @@ stats->robustStdev = myStdev; } - // To determine the median (and later, the lower/upper quartile), we fit - // a quadratic to the three bins surrounding the bin containing the median. - // The quadratic y=f(x) with x being the midpoint of each bin, and y being - // the cumulative number of data points in all bins up to, and including, - // this bin. We then solve the quadratic for - - // XXX: Create a "p_psQuadraticFitAndSolveForX() function. - + + // To determine the median, we fit a quadratic y=f(x) to the three bins + // surrounding the bin containing the median (x is the midpoint of each + // bin and y is the value of each bin). Then we figure out what value + // of x corresponds to f(x) being the median (half of all points). if (stats->options & PS_STAT_ROBUST_MEDIAN) { - if ((maxBinNum > 0) && (maxBinNum < (robustHistogram->nums->n - 1))) { - x->data.F64[0] = (double)0.5 * - (robustHistogram->bounds->data.F32[maxBinNum - 1] + - robustHistogram->bounds->data.F32[maxBinNum]); - x->data.F64[1] = (double)0.5 * - (robustHistogram->bounds->data.F32[maxBinNum] + - robustHistogram->bounds->data.F32[maxBinNum + 1]); - x->data.F64[2] = (double)0.5 * - (robustHistogram->bounds->data.F32[maxBinNum + 1] + - robustHistogram->bounds->data.F32[maxBinNum + 2]); - - y->data.F64[0] = cumulativeRobustSumsDlRange->data.F32[maxBinNum - 1]; - y->data.F64[1] = cumulativeRobustSumsDlRange->data.F32[maxBinNum]; - y->data.F64[2] = cumulativeRobustSumsDlRange->data.F32[maxBinNum + 1]; - - // Ensure that cumulativeMedian/2 is actually within the range of the bins - // we are using. - cumulativeMedian *= 0.5; - if (!((y->data.F64[0] <= cumulativeMedian) && (cumulativeMedian <= y->data.F64[2]))) { - printf("((%f), %f, %f)\n", cumulativeMedian, y->data.F64[0], y->data.F64[2]); - psAbort(__func__, "p_psVectorRobustStats(1): midpoint not within y-range\n"); - } - // XXX: yErr is not currently used by psVectorFitPolynomial1D(). We - // may have to set this meaningfully later. - yErr->data.F64[0] = 1.0; - yErr->data.F64[1] = 1.0; - yErr->data.F64[2] = 1.0; - - // Determine the coefficients of the polynomial. - myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr); - // Call p_ps1DPolyMedian(), which does a binary search on the - // polynomial, looking for the value x such that - // f(x) = cumulativeMedian. - stats->robustMedian = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], cumulativeMedian); - } else { - // If the mode is the first/last histogram bin, then simply use - // the midpoint of that bin. - stats->robustMedian = 0.5 * (robustHistogram->bounds->data.F32[maxBinNum + 1] + - robustHistogram->bounds->data.F32[maxBinNum]); - } - } - // The lower/upper quartile calculations are very similar to the median - // calculations. We fit a quadratic to the array containing the - // cumulative data points in each bin, then search for x such that - // f(x) equals the lower/upper quartile exactly. - // + // Take a psVector. Fit a polynomial y = f(x) to the 3 data elements + // surrounding medianBinNum, then find the x-value corresponding y = sumRobust/2.0. + stats->robustMedian = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds, + cumulativeRobustSums, + medianBinNum, + sumRobust/2.0); + } + + // To determine the quartiles, we fit a quadratic y=f(x) to the three bins + // surrounding the bin containing LQ/UQ (x is the midpoint of each + // bin and y is the value of each bin). Then we figure out what value + // of x corresponds to f(x) being the LQ/UQ. if (stats->options & PS_STAT_ROBUST_QUARTILE) { - countFloat = cumulativeRobustSumsFullRange->data.F32[robustHistogramVector->n - 1]; - - if ((LQBinNum > 0) && (LQBinNum < (robustHistogram->nums->n - 1))) { - x->data.F64[0] = (double)0.5 * - (robustHistogram->bounds->data.F32[LQBinNum - 1] + - robustHistogram->bounds->data.F32[LQBinNum]); - x->data.F64[1] = (double)0.5 * - (robustHistogram->bounds->data.F32[LQBinNum] + - robustHistogram->bounds->data.F32[LQBinNum + 1]); - x->data.F64[2] = (double)0.5 * - (robustHistogram->bounds->data.F32[LQBinNum + 1] + - robustHistogram->bounds->data.F32[LQBinNum + 2]); - - y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[LQBinNum - 1]; - y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[LQBinNum]; - y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[LQBinNum + 1]; - - if (!((y->data.F64[0] <= (countFloat / 4.0)) && ((countFloat / 4.0) <= y->data.F64[2]))) { - psAbort(__func__, "p_psVectorRobustStats(2): midpoint not within y-range\n"); - } - - yErr->data.F64[0] = 1.0; - yErr->data.F64[1] = 1.0; - yErr->data.F64[2] = 1.0; - - myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr); - stats->robustLQ = p_ps1DPolyMedian(myPoly, x->data.F64[0], x->data.F64[2], countFloat / 4.0); - - } else { - // If the LQ is the first/last histogram bin, then simply use - // the midpoint of that bin. - stats->robustLQ = 0.5 * (robustHistogram->bounds->data.F32[LQBinNum + 1] + - robustHistogram->bounds->data.F32[LQBinNum]); - } - - if ((UQBinNum > 0) && (UQBinNum < (robustHistogram->nums->n - 1))) { - x->data.F64[0] = (double)0.5 * - (robustHistogram->bounds->data.F32[UQBinNum - 1] + - robustHistogram->bounds->data.F32[UQBinNum]); - x->data.F64[1] = (double)0.5 * - (robustHistogram->bounds->data.F32[UQBinNum] + - robustHistogram->bounds->data.F32[UQBinNum + 1]); - x->data.F64[2] = (double)0.5 * - (robustHistogram->bounds->data.F32[UQBinNum + 1] + - robustHistogram->bounds->data.F32[UQBinNum + 2]); - - y->data.F64[0] = cumulativeRobustSumsFullRange->data.F32[UQBinNum - 1]; - y->data.F64[1] = cumulativeRobustSumsFullRange->data.F32[UQBinNum]; - y->data.F64[2] = cumulativeRobustSumsFullRange->data.F32[UQBinNum + 1]; - - if (!((y->data.F64[0] <= (3.0 * countFloat / 4.0)) && - ((3.0 * countFloat / 4.0) <= y->data.F64[2]))) { - psAbort(__func__, "p_psVectorRobustStats(3): midpoint not within y-range\n"); - } - - yErr->data.F64[0] = 1.0; - yErr->data.F64[1] = 1.0; - yErr->data.F64[2] = 1.0; - - myPoly = psVectorFitPolynomial1D(myPoly, x, y, yErr); - stats->robustUQ = p_ps1DPolyMedian(myPoly, - x->data.F64[0], x->data.F64[2], 3.0 * countFloat / 4.0); - } else { - // If the UQ is the first/last histogram bin, then simply use - // the midpoint of that bin. - stats->robustUQ = 0.5 * (robustHistogram->bounds->data.F32[UQBinNum + 1] + - robustHistogram->bounds->data.F32[UQBinNum]); - } - } + countFloat = cumulativeRobustSums->data.F32[robustHistogramVector->n - 1]; + + stats->robustLQ = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds, + cumulativeRobustSums, + LQBinNum, + countFloat/4.0); + stats->robustUQ = fitQuadraticSearchForYThenReturnX(robustHistogram->bounds, + cumulativeRobustSums, + UQBinNum, + 3.0 * countFloat/4.0); + } + // XXX: I think sumNfit == sumN50 here. stats->robustNfit = sumNfit; stats->robustN50 = sumN50; - psFree(x); - psFree(y); - psFree(yErr); psFree(tmpStats); psFree(robustHistogram); - psFree(myPoly); - psFree(cumulativeRobustSumsFullRange); - psFree(cumulativeRobustSumsDlRange); + psFree(robustHistogramVector); + psFree(cumulativeRobustSums); } @@ -1468,6 +1350,4 @@ // I am commenting this code out, and replacing it with code which // calculates the mean directly on the robustHistogram. - // - // XXX: Replace this with the LaGrange interpolation stuff. domain = psImageAlloc(1, robustHistogramVector->n, PS_TYPE_F32); @@ -1479,6 +1359,5 @@ max_pos = -1; for (i=0;in;i++) { - domain->data.F32[i][0] = 0.5 * (robustHistogram->bounds->data.F32[i+1] + - robustHistogram->bounds->data.F32[i]); + domain->data.F32[i][0] = PS_BIN_MIDPOINT(robustHistogram, i);bounds->data.F32[modeBinNum + 1]); data->data.F32[i] = (float) robustHistogramVector->data.F32[i]; errors->data.F32[i] = 1.0; @@ -1695,4 +1574,6 @@ psScalar tmpScalar; tmpScalar.type.type = PS_TYPE_F32; + psVector* inF32; + int mustFreeTmp = 1; // NOTE: Verify that this is the correct action. @@ -1706,5 +1587,5 @@ if (out->nums->type.type != PS_TYPE_U32) { - psAbort(__func__, "Only data type PS_TYPE_U32 for output.nums member."); + psAbort(__func__, "Only data type PS_TYPE_U32 for out->nums member."); } // NOTE: Verify that this is the correct action. @@ -1713,15 +1594,17 @@ } - if (in->type.type != PS_TYPE_F32) { - psAbort(__func__, "Only data type PS_TYPE_F32 is currently supported (0x%x).", in->type.type); - } - if (mask != NULL) { if (in->n != mask->n) { - psAbort(__func__, "Vector data and vector mask are of different sizes."); + psError(__func__, "Vector data and vector mask are of different sizes."); } if (mask->type.type != PS_TYPE_U8) { - psAbort(__func__, "Vector mask must be type PS_TYPE_U8"); - } + psError(__func__, "Vector mask must be type PS_TYPE_U8"); + } + } + + inF32 = p_psConvertToF32((psVector *) in); + if (inF32 == NULL) { + inF32 = (psVector *) in; + mustFreeTmp = 0; } @@ -1730,14 +1613,14 @@ numBins = out->nums->n; - for (i = 0; i < in->n; i++) { + for (i = 0; i < inF32->n; i++) { // Check if this pixel is masked, and if so, skip it. if ((mask == NULL) || ((mask != NULL) && (!(mask->data.U8[i] & maskVal)))) { // Check if this pixel is below the minimum value, and if so // count it, then skip it. - if (in->data.F32[i] < out->bounds->data.F32[0]) { + if (inF32->data.F32[i] < out->bounds->data.F32[0]) { out->minNum++; // Check if this pixel is above the maximum value, and if so // count it, then skip it. - } else if (in->data.F32[i] > out->bounds->data.F32[numBins]) { + } else if (inF32->data.F32[i] > out->bounds->data.F32[numBins]) { out->maxNum++; } else { @@ -1746,5 +1629,5 @@ if (out->uniform == true) { binSize = out->bounds->data.F32[1] - out->bounds->data.F32[0]; - binNum = (int)((in->data.F32[i] - out->bounds->data.F32[0]) / binSize); + binNum = (int)((inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize); // NOTE: This next if-statement really shouldn't be necessary. @@ -1755,16 +1638,15 @@ } - (out->nums->data.S32[binNum])++; + (out->nums->data.U32[binNum])++; // If this is a non-uniform histogram, determining the correct // bin number requires a bit more work. } else { - tmpScalar.data.F32 = in->data.F32[i]; - tmp = p_psVectorBinDisect(out->bounds, - &tmpScalar); - if (tmp == -1) { - // XXX: Generate warning message + tmpScalar.data.F32 = inF32->data.F32[i]; + tmp = p_psVectorBinDisect(out->bounds, &tmpScalar); + if (tmp < 0) { + // XXX: Generate warning message } else { - (out->nums->data.S32[tmp])++; + (out->nums->data.U32[tmp])++; } } @@ -1772,19 +1654,22 @@ } } + + if (mustFreeTmp == 1) { + psFree(inF32); + } return (out); } /****************************************************************************** -p_psConvertToF32(myVector, maskVector, maskVal, stats): this is the cheap -way to support a variety of vector data types: we simply convert the input -vector to F32 at the beginning, and write all of our functions in F32. If -the vast majority of all vector stat operations are F32 (or any other single -type), then this is probably the best way to go. Otherwise, when the -algorithms stablize, we will then macro everything and put type support in -the various stat functions. +p_psConvertToF32(in): this is the cheap way to support a variety of vector +data types: we simply convert the input vector to F32 at the beginning, and +write all of our functions in F32. If the vast majority of all vector stat +operations are F32 (or any other single type), then this is probably the +best way to go. Otherwise, when the algorithms stablize, we will then macro +everything and put type support in the various stat functions. XXX: Should the default data type be F64? Since we are buying Athlons... *****************************************************************************/ -psVector* p_psConvertToF32(psStats* stats, psVector* in, psVector* mask, unsigned int maskVal) +psVector* p_psConvertToF32(psVector* in) { int i = 0; @@ -1800,4 +1685,9 @@ for (i = 0; i < in->n; i++) { tmp->data.F32[i] = (float)in->data.U16[i]; + } + } else if (in->type.type == PS_TYPE_U8) { + tmp = psVectorAlloc(in->n, PS_TYPE_F32); + for (i = 0; i < in->n; i++) { + tmp->data.F32[i] = (float)in->data.U8[i]; } } else if (in->type.type == PS_TYPE_F64) { @@ -1839,20 +1729,18 @@ int mustFreeTmp = 1; - // NOTE: Verify that this is the correct action. if (in == NULL) { - return (stats); - } - if (stats == NULL) { - return (NULL); - } - - inF32 = p_psConvertToF32(stats, in, mask, maskVal); + psError(__func__, "in is NULL\n"); + return(stats); + } + PS_CHECK_NULL_PTR_RETURN_NULL(stats); + + inF32 = p_psConvertToF32((psVector *) in); if (inF32 == NULL) { inF32 = in; mustFreeTmp = 0; } - // XXX: Should we abort if (stats->min == stats->max)? if ((stats->options & PS_STAT_USE_RANGE) && (stats->min >= stats->max)) { - psAbort(__func__, "psVectorStats() called with range: %f to %f\n", stats->min, stats->max); + psError(__func__, "psVectorStats() called with range: %f to %f\n", stats->min, stats->max); + return(stats); } Index: /trunk/psLib/test/dataManip/tst_psMinimize06.c =================================================================== --- /trunk/psLib/test/dataManip/tst_psMinimize06.c (revision 2196) +++ /trunk/psLib/test/dataManip/tst_psMinimize06.c (revision 2197) @@ -74,6 +74,96 @@ } - -int main() +psVector *myFunc2(psImage *myDeriv, + psVector *myParams, + psArray *myCoords) +{ + psVector *sum = psVectorAlloc(myCoords->n, PS_TYPE_F32); + int i; + int j; + + + if (myDeriv == NULL) { + myDeriv = psImageAlloc(myParams->n, myCoords->n, PS_TYPE_F32); + psError(__func__, "myDeriv is NULL.\n"); + } + + for (i=0;idata.F32[i] = MIN_VALUE; + for (j=0;jdata.F32[i]+= (myParams->data.F32[j] - expectedParm[j]) * + (myParams->data.F32[j] - expectedParm[j]); + + myDeriv->data.F32[i][j] = (2.0 * myParams->data.F32[j]) - + (2.0 * expectedParm[j]); + } + } + + return(sum); +} + +int t01() +{ + int currentId = psMemGetId(); + int memLeaks = 0; + int i = 0; + psArray *myCoords; + psVector *myParams; + psImage *myCovar; + psMinimization *min; + psVector *y; + + psTraceSetLevel(".psLib", 0); + /************************************************************************** + *************************************************************************/ + min = psMinimizationAlloc(NUM_ITERATIONS, ERR_TOL); + myCovar = psImageAlloc(NUM_PARAMS, NUM_PARAMS, PS_TYPE_F32); + myParams = psVectorAlloc(NUM_PARAMS, PS_TYPE_F32); + myCoords = psArrayAlloc(N); + y = psVectorAlloc(N, PS_TYPE_F32); + + for (i=0;idata[i] = (psPTR *) psVectorAlloc(2, PS_TYPE_F32); + ((psVector *) (myCoords->data[i]))->data.F32[0] = (float) (i+10); + ((psVector *) (myCoords->data[i]))->data.F32[1] = (float) (i+3); + y->data.F32[i] = (float) i; + } + for (i=0;idata.F32[i] = (float) i; + myParams->data.F32[i] = expectedParm[i] * 1.3; + myParams->data.F32[i] = (float) (5 + i); + myParams->data.F32[i] = 0.0; + } + + psMyMinimizeLMChi2(min, + myCovar, + myParams, + NULL, + myCoords, + y, + NULL, + (psMyMinimizeLMChi2Func) myFunc2); + + printf("\nThe chi-squared is %f\n", min->value); + for (i=0;idata.F32[i], expectedParm[i]); + } + + psFree(min); + psFree(myCovar); + psFree(myParams); + psFree(myCoords); + psFree(y); + + psMemCheckCorruption(1); + memLeaks = psMemCheckLeaks(currentId,NULL,NULL); + if (0 != memLeaks) { + psAbort(__func__,"Memory Leaks! (%d leaks)", memLeaks); + } + return (!testStatus); +} + +int t02() { int currentId = psMemGetId(); @@ -85,8 +175,8 @@ psImage *myCovar; psMinimization *min; - psVector *x = psVectorAlloc(N, PS_TYPE_F32); psVector *y = psVectorAlloc(N, PS_TYPE_F32); psTraceSetLevel(".psLib", 0); + t02(); /************************************************************************** *************************************************************************/ @@ -101,5 +191,4 @@ ((psVector *) (myCoords->data[i]))->data.F32[0] = (float) (i+10); ((psVector *) (myCoords->data[i]))->data.F32[1] = (float) (i+3); - x->data.F32[i] = (float) i; y->data.F32[i] = (float) i; } @@ -131,5 +220,4 @@ psFree(myParamMask); psFree(min); - psFree(x); psFree(y); psFree(myCovar); @@ -142,2 +230,8 @@ return (!testStatus); } + +int main() +{ + t01(); + // t02(); +} Index: /trunk/psLib/test/dataManip/tst_psStats07.c =================================================================== --- /trunk/psLib/test/dataManip/tst_psStats07.c (revision 2196) +++ /trunk/psLib/test/dataManip/tst_psStats07.c (revision 2197) @@ -10,10 +10,10 @@ #include -#define N 200 +#define N 90 #define MEAN 32.0 #define STDEV 2.0 #define ERROR_TOLERANCE 0.10 -int main() +int t00() { psStats * myStats = NULL; @@ -32,18 +32,11 @@ float realMedianNoMask = MEAN; float realModeNoMask = MEAN; - float realStdevNoMask = STDEV * 0.33; + float realStdevNoMask = STDEV * 0.20; float realLQNoMask = MEAN - ( 0.6 * STDEV ); float realUQNoMask = MEAN + ( 0.6 * STDEV ); int realN50NoMask = N / 4; int realNfitNoMask = N / 4; - float realMeanWithMask = MEAN; - float realMedianWithMask = MEAN; - float realModeWithMask = MEAN; - float realStdevWithMask = STDEV; - float realLQWithMask = MEAN; - float realUQWithMask = MEAN; - int realN50WithMask = N / 4; - int realNfitWithMask = N / 4; - + + psTraceSetLevel(".psLib.dataManip.psStats", 10); /*************************************************************************/ @@ -232,7 +225,78 @@ "PS_STAT_ROBUST_STATS: robust Nfit: no vector mask", testStatus ); - - - return ( 0 ); + return(testStatus); + + /*************************************************************************/ + /* Deallocate data structures */ + /*************************************************************************/ + printPositiveTestHeader( stdout, + "psStats functions", + "psStats(): deallocating memory" ); + + psFree( myStats ); + psFree( myVector ); + psFree( maskVector ); + + psMemCheckCorruption( 1 ); + memLeaks = psMemCheckLeaks(currentId,NULL,stderr); + if ( 0 != memLeaks ) { + psAbort( __func__, "Memory Leaks! (%d leaks)", memLeaks ); + } + + printFooter( stdout, + "psVector functions", + "psStats(): deallocating memory", + testStatus ); + + return ( !globalTestStatus ); +} + + +int t01() +{ + psStats * myStats = NULL; + int testStatus = true; + int globalTestStatus = true; + int i = 0; + psVector *myVector = NULL; + psVector *maskVector = NULL; + // NOTE: These values were calculated by running the function on the data. + // A: They must be changed if we adjust the number of data points. + // B: We don't really know that they are correct. + int count = 0; + int currentId = psMemGetId(); + int memLeaks = 0; + float realMeanWithMask = MEAN; + float realMedianWithMask = MEAN; + float realModeWithMask = MEAN; + float realStdevWithMask = STDEV * 0.20; + float realLQWithMask = MEAN; + float realUQWithMask = MEAN; + int realN50WithMask = N / 4; + int realNfitWithMask = N / 4; + + psTraceSetLevel(".psLib.dataManip.psStats", 10); + + /*************************************************************************/ + /* Allocate and initialize data structures */ + /*************************************************************************/ + myStats = psStatsAlloc( PS_STAT_ROBUST_MEAN | + PS_STAT_ROBUST_MEDIAN | + PS_STAT_ROBUST_MODE | + PS_STAT_ROBUST_STDEV | + PS_STAT_ROBUST_QUARTILE ); + + maskVector = psVectorAlloc( N, PS_TYPE_U8 ); + maskVector->n = N; + myVector = psGaussianDev( MEAN, STDEV, N ); + // Set the mask vector and calculate the expected maximum. + for ( i = 0;i < N;i++ ) { + if ( i < ( N / 2 ) ) { + maskVector->data.U8[ i ] = 0; + count++; + } else { + maskVector->data.U8[ i ] = 1; + } + } /*************************************************************************/ @@ -422,2 +486,8 @@ return ( !globalTestStatus ); } + +int main() +{ + t00(); + t01(); +}