Index: /branches/eam_branches/ipp-20130509/psLib/src/math/psMinimizeLMM.c =================================================================== --- /branches/eam_branches/ipp-20130509/psLib/src/math/psMinimizeLMM.c (revision 35622) +++ /branches/eam_branches/ipp-20130509/psLib/src/math/psMinimizeLMM.c (revision 35623) @@ -365,4 +365,29 @@ /****************************************************************************** +psMinimizeLMChi2(): wrapper to call either _Old or _Alt + *****************************************************************************/ +bool psMinimizeLMChi2( + psMinimization *min, + psImage *covar, + psVector *params, + psMinConstraint *constraint, + const psArray *x, + const psVector *y, + const psVector *yWt, + psMinimizeLMChi2Func func) +{ + bool status = psMinimizeLMChi2_Alt( + min, + covar, + params, + constraint, + x, + y, + yWt, + func); + return status; +} + +/****************************************************************************** psMinimizeLMChi2(): This routine will take an procedure which calculates an arbitrary function and it's derivative and minimize the chi-squared match @@ -373,5 +398,5 @@ XXX Make an F64 version? *****************************************************************************/ -bool psMinimizeLMChi2( +bool psMinimizeLMChi2_Old( psMinimization *min, psImage *covar, @@ -510,7 +535,7 @@ psF32 rho = (min->value - Chisq) / dLinear; - psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %f, nDOF: %d\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda, nDOF); - - psTrace("psLib.math.dLinear", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %f\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda); + psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %g, nDOF: %d\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda, nDOF); + + psTrace("psLib.math.dLinear", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %g\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda); // dump some useful info if trace is defined @@ -662,4 +687,5 @@ psF32 Chisq = 0.0; psF32 lambda = 0.001; + psF32 nu = 2.0; psF32 dLinear = 0.0; @@ -706,5 +732,6 @@ min->iter ++; if (min->iter >= min->maxIter) break; - lambda *= 10.0; + // lambda *= 10.0; + lambda *= 2.0; continue; } @@ -732,9 +759,18 @@ // dParDelta = Beta[i]^2 : change (squared) in parameter i float rParSigma = 0.0; - for (int j = 0; j < Beta->n; j++) { - rParSigma += PS_SQR(Beta->data.F32[j]) / Alpha->data.F32[j][j]; + // for (int j = 0; j < Beta->n; j++) { + // rParSigma += PS_SQR(Beta->data.F32[j]) / Alpha->data.F32[j][j]; + // } + for (int j = 0, J = 0; j < Params->n; j++) { + if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[j])) { + continue; + } + // fprintf (stderr, "del: %f - %f vs %f : %f\n", Params->data.F32[j], params->data.F32[j], sqrt(Alpha->data.F32[J][J]), PS_SQR(Params->data.F32[j] - params->data.F32[j]) / Alpha->data.F32[J][J]); + rParSigma += PS_SQR(Params->data.F32[j] - params->data.F32[j]) / Alpha->data.F32[J][J]; + J++; } rParSigma = sqrt(rParSigma); psTrace("psLib.math", 5, "rParSigma: %f, Niter: %d\n", rParSigma, min->iter); + fprintf (stderr, "rParSigma: %f, Niter: %d\n", rParSigma, min->iter); // calculate Chisq for new guess, update Alpha & Beta @@ -743,5 +779,6 @@ min->iter ++; if (min->iter >= min->maxIter) break; - lambda *= 10.0; + // lambda *= 10.0; + lambda *= 2.0; continue; } @@ -753,7 +790,7 @@ psF32 rho = (min->value - Chisq) / dLinear; - psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %f, nDOF: %d\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda, nDOF); - - psTrace("psLib.math.dLinear", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %f\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda); + psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %g, nDOF: %d\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda, nDOF); + + psTrace("psLib.math.dLinear", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %g\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda); // dump some useful info if trace is defined @@ -767,5 +804,8 @@ // rho is positive if the new chisq is smaller; allow for some insignificant change (slight negative rho) + + // XXX the old version of lambda changes: // XXX : Madsen gives suggestion for better use of rho +# if (0) if (rho >= -1e-6) { min->value = Chisq; @@ -777,4 +817,31 @@ lambda *= 10.0; } +# endif + + // rho is positive if the new chisq is smaller + if (rho >= -1e-6) { + min->value = Chisq; + alpha = psImageCopy(alpha, Alpha, PS_TYPE_F32); + beta = psVectorCopy(beta, Beta, PS_TYPE_F32); + params = psVectorCopy(params, Params, PS_TYPE_F32); + } + +# if (0) + // adjust the gain ratio (lambda) based on rho + if (rho < 0.25) { + lambda *= 2.0; + } + if (rho > 0.75) { + lambda *= 0.333; + } +# else + if (rho > 0.0) { + lambda *= PS_MAX(0.33, (1.0 - pow(2.0*rho - 1.0, 3.0))); + nu = 2.0; + } else { + lambda *= nu; + nu *= 2.0; + } +# endif min->iter++; @@ -831,4 +898,5 @@ // if the last improvement was at least as good as maxTol, accept the fit: + // XXX : maybe re-work this to use rParSigma & maxTol if (min->lastDelta <= min->maxTol) { psTrace("psLib.math", 6, "---- end (true) ----\n"); Index: /branches/eam_branches/ipp-20130509/psLib/src/math/psMinimizeLMM.h =================================================================== --- /branches/eam_branches/ipp-20130509/psLib/src/math/psMinimizeLMM.h (revision 35622) +++ /branches/eam_branches/ipp-20130509/psLib/src/math/psMinimizeLMM.h (revision 35623) @@ -134,4 +134,19 @@ ); +/** Minimizes a specified function based on the Levenberg-Marquardt method. + * + * @return bool: True if successful. + */ +bool psMinimizeLMChi2_Old( + psMinimization *min, ///< Minimization specification + psImage *covar, ///< Covariance matrix + psVector *params, ///< "Best Guess" for the parameters that minimize func + psMinConstraint *constraint, ///< Constraints on the parameters + const psArray *x, ///< Measurement ordinates of multiple vectors + const psVector *y, ///< Measurement coordinates + const psVector *yWt, ///< Errors in the measurement coordinates + psMinimizeLMChi2Func func ///< Specified function +); + /** Minimizes a specified function based on the Levenberg-Marquardt method Uses alternative convergence criterion. Index: /branches/eam_branches/ipp-20130509/psLib/test/math/tap_psMinimizeLMM_Trail.c =================================================================== --- /branches/eam_branches/ipp-20130509/psLib/test/math/tap_psMinimizeLMM_Trail.c (revision 35622) +++ /branches/eam_branches/ipp-20130509/psLib/test/math/tap_psMinimizeLMM_Trail.c (revision 35623) @@ -501,5 +501,5 @@ # if (1) - psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5); + // psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5); bool tmpBool = psMinimizeLMChi2_Alt(min, NULL, params, NULL, ordinates, coordinates, errors, (psMinimizeLMChi2Func)fitFunc); # else