Index: /trunk/psLib/src/math/psMinimizeLMM.c =================================================================== --- /trunk/psLib/src/math/psMinimizeLMM.c (revision 10252) +++ /trunk/psLib/src/math/psMinimizeLMM.c (revision 10253) @@ -10,6 +10,6 @@ * @author EAM, IfA * - * @version $Revision: 1.26 $ $Name: not supported by cvs2svn $ - * @date $Date: 2006-11-26 21:53:57 $ + * @version $Revision: 1.27 $ $Name: not supported by cvs2svn $ + * @date $Date: 2006-11-29 02:25:14 $ * * Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii @@ -49,5 +49,4 @@ /*****************************************************************************/ -// XXX EAM : can we use static copies of LUv, LUm, A? psBool p_psMinLM_GuessABP( psImage *Alpha, @@ -58,14 +57,12 @@ const psVector *params, const psVector *paramMask, - const psVector *beta_lim, - const psVector *params_min, - const psVector *params_max, - psF64 lambda) -{ - PS_ASSERT_VECTOR_TYPE(Alpha, PS_TYPE_F64, false); - PS_ASSERT_VECTOR_TYPE(Beta, PS_TYPE_F64, false); + psMinimizeLMLimitFunc checkLimits, + psF32 lambda) +{ + PS_ASSERT_VECTOR_TYPE(Alpha, PS_TYPE_F32, false); + PS_ASSERT_VECTOR_TYPE(Beta, PS_TYPE_F32, false); PS_ASSERT_VECTOR_TYPE(Params, PS_TYPE_F32, false); - PS_ASSERT_VECTOR_TYPE(alpha, PS_TYPE_F64, false); - PS_ASSERT_VECTOR_TYPE(beta, PS_TYPE_F64, false); + PS_ASSERT_VECTOR_TYPE(alpha, PS_TYPE_F32, false); + PS_ASSERT_VECTOR_TYPE(beta, PS_TYPE_F32, false); PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, false); if (paramMask) { @@ -73,4 +70,6 @@ } + // XXX the LU Decomposition code requires F64. this is now incompatible + // with the rest of this code. drop this segment when the changes are tested. # define USE_LU_DECOMP 0 # if (USE_LU_DECOMP) @@ -115,14 +114,14 @@ // set new guess values (creates matrix A) - Beta = psVectorCopy(Beta, beta, PS_TYPE_F64); - Alpha = psImageCopy(Alpha, alpha, PS_TYPE_F64); + Beta = psVectorCopy(Beta, beta, PS_TYPE_F32); + Alpha = psImageCopy(Alpha, alpha, PS_TYPE_F32); for (int j = 0; j < params->n; j++) { if ((paramMask != NULL) && (paramMask->data.U8[j])) { continue; } - Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda); - } - - if (false == psMatrixGJSolve(Alpha, Beta)) { + Alpha->data.F32[j][j] = alpha->data.F32[j][j] * (1.0 + lambda); + } + + if (false == psMatrixGJSolveF32(Alpha, Beta)) { // psError(PS_ERR_UNKNOWN, false, "singular matrix in Guess ABP\n"); psTrace ("psLib.math", 4, "singular matrix in Guess ABP\n"); @@ -137,25 +136,17 @@ continue; } - // Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j]; - // compare Beta to beta limits - if (beta_lim != NULL) { - if (fabs(Beta->data.F64[j]) > fabs(beta_lim->data.F32[j])) { - Beta->data.F64[j] = (Beta->data.F64[j] > 0) ? fabs(beta_lim->data.F32[j]) : -fabs(beta_lim->data.F32[j]); - } - } - Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j]; + // apply beta limits + if (checkLimits) + checkLimits (PS_MINIMIZE_BETA_LIMIT, j, Params->data.F32, Beta->data.F32); + Params->data.F32[j] = params->data.F32[j] - Beta->data.F32[j]; + // compare new params to param limits - if (params_max != NULL) { - Params->data.F32[j] = PS_MIN (Params->data.F32[j], params_max->data.F32[j]); - } - if (params_min != NULL) { - Params->data.F32[j] = PS_MAX (Params->data.F32[j], params_min->data.F32[j]); - } - } - - + if (checkLimits) + checkLimits (PS_MINIMIZE_PARAM_MIN, j, Params->data.F32, Beta->data.F32); + if (checkLimits) + checkLimits (PS_MINIMIZE_PARAM_MAX, j, Params->data.F32, Beta->data.F32); + } return(true); } - bool psMinimizeGaussNewtonDelta( @@ -170,7 +161,7 @@ psTrace("psLib.math", 3, "---- begin ----\n"); // allocate internal arrays (current vs Guess) - psImage *alpha = psImageAlloc (params->n, params->n, PS_TYPE_F64); - psImage *Alpha = psImageAlloc (params->n, params->n, PS_TYPE_F64); - psVector *beta = psVectorAlloc(params->n, PS_TYPE_F64); + psImage *alpha = psImageAlloc (params->n, params->n, PS_TYPE_F32); + psImage *Alpha = psImageAlloc (params->n, params->n, PS_TYPE_F32); + psVector *beta = psVectorAlloc(params->n, PS_TYPE_F32); psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32); psVector *dy = NULL; @@ -186,6 +177,7 @@ } - psF64 rcF64 = p_psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func); - if (isnan(rcF64)) { + // XXX should we give up if chisq is nan? + psF32 rcF32 = p_psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func); + if (isnan(rcF32)) { psTrace ("psLib.math", 5, "p_psMinLM_SetABX() returned a NAN.\n"); rc = false; @@ -199,5 +191,5 @@ } - psBool rcBool = p_psMinLM_GuessABP(Alpha, delta, Params, alpha, beta, params, paramMask, NULL, NULL, NULL, 0.0); + psBool rcBool = p_psMinLM_GuessABP(Alpha, delta, Params, alpha, beta, params, paramMask, NULL, 0.0); if (rcBool == false) { psTrace ("psLib.math", 5, "p_psMinLM_GuessABP() returned FALSE.\n"); @@ -223,14 +215,14 @@ // measure linear model prediction -psF64 p_psMinLM_dLinear( +psF32 p_psMinLM_dLinear( const psVector *Beta, const psVector *beta, - psF64 lambda) + psF32 lambda) { /* get linear model prediction */ - psF64 dLinear = 0; - psF64 *B = Beta->data.F64; - psF64 *b = beta->data.F64; + psF32 dLinear = 0; + psF32 *B = Beta->data.F32; + psF32 *b = beta->data.F32; for (int i = 0; i < beta->n; i++) { dLinear += lambda*PS_SQR(B[i]) + B[i]*b[i]; @@ -239,7 +231,6 @@ } -// XXX EAM: this needs to respect the mask on params // alpha, beta, params are already allocated -psF64 p_psMinLM_SetABX( +psF32 p_psMinLM_SetABX( psImage *alpha, psVector *beta, @@ -259,13 +250,13 @@ PS_ASSERT_VECTOR_NON_NULL(dy, NAN); - PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, false); + PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, false); if (paramMask) { PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_MASK, false); } - psF64 chisq; - psF64 delta; - psF64 weight; - psF64 ymodel; + psF32 chisq; + psF32 delta; + psF32 weight; + psF32 ymodel; psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32); @@ -273,7 +264,7 @@ for (psS32 j = 0; j < params->n; j++) { for (psS32 k = 0; k < params->n; k++) { - alpha->data.F64[j][k] = 0; - } - beta->data.F64[j] = 0; + alpha->data.F32[j][k] = 0; + } + beta->data.F32[j] = 0; } chisq = 0.0; @@ -301,7 +292,7 @@ continue; } - alpha->data.F64[j][k] += weight * deriv->data.F32[k]; + alpha->data.F32[j][k] += weight * deriv->data.F32[k]; } - beta->data.F64[j] += weight * delta; + beta->data.F32[j] += weight * delta; } } @@ -310,5 +301,5 @@ for (psS32 j = 1; j < params->n; j++) { for (psS32 k = 0; k < j; k++) { - alpha->data.F64[k][j] = alpha->data.F64[j][k]; + alpha->data.F32[k][j] = alpha->data.F32[j][k]; } } @@ -318,6 +309,6 @@ for (psS32 j = 0; j < params->n; j++) { if (paramMask->data.U8[j]) { - alpha->data.F64[j][j] = 1; - beta->data.F64[j] = 1; + alpha->data.F32[j][j] = 1; + beta->data.F32[j] = 1; } } @@ -335,9 +326,6 @@ coords. -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 an optimization. - +This requires F32 input data; all internal calls use F32. +XXX Make an F64 version? *****************************************************************************/ psBool psMinimizeLMChi2( @@ -345,5 +333,5 @@ psImage *covar, psVector *params, - psMinConstrain *constrain, + psMinConstraint *constraint, const psArray *x, const psVector *y, @@ -357,18 +345,6 @@ PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, false); psVector *paramMask = NULL; - psVector *paramDelta = NULL; - psVector *paramMin = NULL; - psVector *paramMax = NULL; - if (constrain != NULL) { - paramDelta = constrain->paramDelta; - paramMin = constrain->paramMin; - paramMax = constrain->paramMax; - paramMask = constrain->paramMask; - PS_ASSERT_VECTOR_TYPE(paramDelta, PS_TYPE_F32, false); - PS_ASSERT_VECTOR_TYPE(paramMin, PS_TYPE_F32, false); - PS_ASSERT_VECTOR_TYPE(paramMax, PS_TYPE_F32, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(params, paramDelta, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(params, paramMin, false); - PS_ASSERT_VECTORS_SIZE_EQUAL(params, paramMax, false); + if (constraint != NULL) { + paramMask = constraint->paramMask; if (paramMask != NULL) { PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_U8, false); @@ -392,4 +368,9 @@ PS_ASSERT_PTR_NON_NULL(func, false); + psMinimizeLMLimitFunc checkLimits = NULL; + if (constraint) { + checkLimits = constraint->checkLimits; + } + // this function has test and current values for several things // the current best value is in lower case @@ -397,12 +378,12 @@ // allocate internal arrays (current vs Guess) - psImage *alpha = psImageAlloc(params->n, params->n, PS_TYPE_F64); - psImage *Alpha = psImageAlloc(params->n, params->n, PS_TYPE_F64); - psVector *beta = psVectorAlloc(params->n, PS_TYPE_F64); - psVector *Beta = psVectorAlloc(params->n, PS_TYPE_F64); + psImage *alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32); + psImage *Alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32); + psVector *beta = psVectorAlloc(params->n, PS_TYPE_F32); + psVector *Beta = psVectorAlloc(params->n, PS_TYPE_F32); psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32); psVector *dy = NULL; - psF64 Chisq = 0.0; - psF64 lambda = 0.001; + psF32 Chisq = 0.0; + psF32 lambda = 0.001; // the user provides the error or NULL. we need to convert @@ -436,6 +417,5 @@ // set a new guess for Alpha, Beta, Params - if (!p_psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, - paramDelta, paramMin, paramMax, lambda)) { + if (!p_psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda)) { min->iter ++; lambda *= 10.0; @@ -444,5 +424,5 @@ // measure linear model prediction - psF64 dLinear = p_psMinLM_dLinear(Beta, beta, lambda); + psF32 dLinear = p_psMinLM_dLinear(Beta, beta, lambda); // dump some useful info if trace is defined @@ -467,5 +447,5 @@ // expected delta from the linear model (dLinear) // accept new guess if it is an improvement (rho > 0), or else increase lambda - psF64 rho = (min->value - Chisq) / dLinear; + psF32 rho = (min->value - Chisq) / dLinear; psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value, @@ -482,6 +462,6 @@ min->lastDelta = (min->value - Chisq) / (dy->n - params->n); min->value = Chisq; - alpha = psImageCopy(alpha, Alpha, PS_TYPE_F64); - beta = psVectorCopy(beta, Beta, PS_TYPE_F64); + alpha = psImageCopy(alpha, Alpha, PS_TYPE_F32); + beta = psVectorCopy(beta, Beta, PS_TYPE_F32); params = psVectorCopy(params, Params, PS_TYPE_F32); lambda *= 0.1; @@ -495,6 +475,5 @@ // construct & return the covariance matrix (if requested) if (covar != NULL) { - if (!p_psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask, - paramDelta, paramMin, paramMax, 0.0)) { + if (!p_psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0)) { psTrace ("psLib.math", 5, "failure to calculate covariance matrix\n"); } @@ -546,22 +525,24 @@ -static void constrainFree(psMinConstrain *tmp) -{ - // There are no dynamically allocated items -} - -psMinConstrain* psMinConstrainAlloc() -{ - psMinConstrain *tmp = psAlloc(sizeof(psMinConstrain)); +static void constraintFree(psMinConstraint *tmp) +{ + if (tmp == NULL) + return; + + psFree (tmp->paramMask); +} + +psMinConstraint* psMinConstraintAlloc() +{ + psMinConstraint *tmp = psAlloc(sizeof(psMinConstraint)); + psMemSetDeallocator(tmp, (psFreeFunc)constraintFree); tmp->paramMask = NULL; - tmp->paramMax = NULL; - tmp->paramMin = NULL; - tmp->paramDelta = NULL; + tmp->checkLimits = NULL; return(tmp); } -bool psMemCheckConstrain(psPtr tmp) -{ - return(psMemGetDeallocator(tmp) == (psFreeFunc) constrainFree); -} +bool psMemCheckConstraint(psPtr tmp) +{ + return(psMemGetDeallocator(tmp) == (psFreeFunc) constraintFree); +}