IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 28, 2006, 4:25:14 PM (19 years ago)
Author:
magnier
Message:

changed psMinimizeLMM limits to a function call (from static vectors); added psMinConstraintMode, psMinimizeLMLimitFunc, converted psMinimizeLMChi2 to F32 only

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/math/psMinimizeLMM.c

    r10192 r10253  
    1010 *  @author EAM, IfA
    1111 *
    12  *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2006-11-26 21:53:57 $
     12 *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2006-11-29 02:25:14 $
    1414 *
    1515 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    4949/*****************************************************************************/
    5050
    51 // XXX EAM : can we use static copies of LUv, LUm, A?
    5251psBool p_psMinLM_GuessABP(
    5352    psImage  *Alpha,
    … …  
    5857    const psVector *params,
    5958    const psVector *paramMask,
    60     const psVector *beta_lim,
    61     const psVector *params_min,
    62     const psVector *params_max,
    63     psF64 lambda)
    64 {
    65     PS_ASSERT_VECTOR_TYPE(Alpha,     PS_TYPE_F64,  false);
    66     PS_ASSERT_VECTOR_TYPE(Beta,      PS_TYPE_F64,  false);
     59    psMinimizeLMLimitFunc checkLimits,
     60    psF32 lambda)
     61{
     62    PS_ASSERT_VECTOR_TYPE(Alpha,     PS_TYPE_F32,  false);
     63    PS_ASSERT_VECTOR_TYPE(Beta,      PS_TYPE_F32,  false);
    6764    PS_ASSERT_VECTOR_TYPE(Params,    PS_TYPE_F32,  false);
    68     PS_ASSERT_VECTOR_TYPE(alpha,     PS_TYPE_F64,  false);
    69     PS_ASSERT_VECTOR_TYPE(beta,      PS_TYPE_F64,  false);
     65    PS_ASSERT_VECTOR_TYPE(alpha,     PS_TYPE_F32,  false);
     66    PS_ASSERT_VECTOR_TYPE(beta,      PS_TYPE_F32,  false);
    7067    PS_ASSERT_VECTOR_TYPE(params,    PS_TYPE_F32,  false);
    7168    if (paramMask) {
    … …  
    7370    }
    7471
     72    // XXX the LU Decomposition code requires F64.  this is now incompatible
     73    // with the rest of this code. drop this segment when the changes are tested.
    7574    # define USE_LU_DECOMP 0
    7675    # if (USE_LU_DECOMP)
    … …  
    115114
    116115    // set new guess values (creates matrix A)
    117     Beta = psVectorCopy(Beta, beta, PS_TYPE_F64);
    118     Alpha = psImageCopy(Alpha, alpha, PS_TYPE_F64);
     116    Beta = psVectorCopy(Beta, beta, PS_TYPE_F32);
     117    Alpha = psImageCopy(Alpha, alpha, PS_TYPE_F32);
    119118    for (int j = 0; j < params->n; j++) {
    120119        if ((paramMask != NULL) && (paramMask->data.U8[j])) {
    121120            continue;
    122121        }
    123         Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
    124     }
    125 
    126     if (false == psMatrixGJSolve(Alpha, Beta)) {
     122        Alpha->data.F32[j][j] = alpha->data.F32[j][j] * (1.0 + lambda);
     123    }
     124
     125    if (false == psMatrixGJSolveF32(Alpha, Beta)) {
    127126        // psError(PS_ERR_UNKNOWN, false, "singular matrix in Guess ABP\n");
    128127        psTrace ("psLib.math", 4, "singular matrix in Guess ABP\n");
    … …  
    137136            continue;
    138137        }
    139         // Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j];
    140         // compare Beta to beta limits
    141         if (beta_lim != NULL) {
    142             if (fabs(Beta->data.F64[j]) > fabs(beta_lim->data.F32[j])) {
    143                 Beta->data.F64[j] = (Beta->data.F64[j] > 0) ? fabs(beta_lim->data.F32[j]) : -fabs(beta_lim->data.F32[j]);
    144             }
    145         }
    146         Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j];
     138        // apply beta limits
     139        if (checkLimits)
     140            checkLimits (PS_MINIMIZE_BETA_LIMIT, j, Params->data.F32, Beta->data.F32);
     141        Params->data.F32[j] = params->data.F32[j] - Beta->data.F32[j];
     142
    147143        // compare new params to param limits
    148         if (params_max != NULL) {
    149             Params->data.F32[j] = PS_MIN (Params->data.F32[j], params_max->data.F32[j]);
    150         }
    151         if (params_min != NULL) {
    152             Params->data.F32[j] = PS_MAX (Params->data.F32[j], params_min->data.F32[j]);
    153         }
    154     }
    155 
    156 
     144        if (checkLimits)
     145            checkLimits (PS_MINIMIZE_PARAM_MIN,  j, Params->data.F32, Beta->data.F32);
     146        if (checkLimits)
     147            checkLimits (PS_MINIMIZE_PARAM_MAX,  j, Params->data.F32, Beta->data.F32);
     148    }
    157149    return(true);
    158150}
    159 
    160151
    161152bool psMinimizeGaussNewtonDelta(
    … …  
    170161    psTrace("psLib.math", 3, "---- begin ----\n");
    171162    // allocate internal arrays (current vs Guess)
    172     psImage  *alpha  = psImageAlloc (params->n, params->n, PS_TYPE_F64);
    173     psImage  *Alpha  = psImageAlloc (params->n, params->n, PS_TYPE_F64);
    174     psVector *beta   = psVectorAlloc(params->n, PS_TYPE_F64);
     163    psImage  *alpha  = psImageAlloc (params->n, params->n, PS_TYPE_F32);
     164    psImage  *Alpha  = psImageAlloc (params->n, params->n, PS_TYPE_F32);
     165    psVector *beta   = psVectorAlloc(params->n, PS_TYPE_F32);
    175166    psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32);
    176167    psVector *dy     = NULL;
    … …  
    186177    }
    187178
    188     psF64 rcF64 = p_psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func);
    189     if (isnan(rcF64)) {
     179    // XXX should we give up if chisq is nan?
     180    psF32 rcF32 = p_psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func);
     181    if (isnan(rcF32)) {
    190182        psTrace ("psLib.math", 5, "p_psMinLM_SetABX() returned a NAN.\n");
    191183        rc = false;
    … …  
    199191    }
    200192
    201     psBool rcBool = p_psMinLM_GuessABP(Alpha, delta, Params, alpha, beta, params, paramMask, NULL, NULL, NULL, 0.0);
     193    psBool rcBool = p_psMinLM_GuessABP(Alpha, delta, Params, alpha, beta, params, paramMask, NULL, 0.0);
    202194    if (rcBool == false) {
    203195        psTrace ("psLib.math", 5, "p_psMinLM_GuessABP() returned FALSE.\n");
    … …  
    223215
    224216// measure linear model prediction
    225 psF64 p_psMinLM_dLinear(
     217psF32 p_psMinLM_dLinear(
    226218    const psVector *Beta,
    227219    const psVector *beta,
    228     psF64 lambda)
     220    psF32 lambda)
    229221{
    230222
    231223    /* get linear model prediction */
    232     psF64 dLinear = 0;
    233     psF64 *B = Beta->data.F64;
    234     psF64 *b = beta->data.F64;
     224    psF32 dLinear = 0;
     225    psF32 *B = Beta->data.F32;
     226    psF32 *b = beta->data.F32;
    235227    for (int i = 0; i < beta->n; i++) {
    236228        dLinear += lambda*PS_SQR(B[i]) + B[i]*b[i];
    … …  
    239231}
    240232
    241 // XXX EAM: this needs to respect the mask on params
    242233// alpha, beta, params are already allocated
    243 psF64 p_psMinLM_SetABX(
     234psF32 p_psMinLM_SetABX(
    244235    psImage  *alpha,
    245236    psVector *beta,
    … …  
    259250    PS_ASSERT_VECTOR_NON_NULL(dy, NAN);
    260251
    261     PS_ASSERT_VECTOR_TYPE(params,    PS_TYPE_F32, false);
     252    PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, false);
    262253    if (paramMask) {
    263254        PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_MASK, false);
    264255    }
    265256
    266     psF64 chisq;
    267     psF64 delta;
    268     psF64 weight;
    269     psF64 ymodel;
     257    psF32 chisq;
     258    psF32 delta;
     259    psF32 weight;
     260    psF32 ymodel;
    270261    psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32);
    271262
    … …  
    273264    for (psS32 j = 0; j < params->n; j++) {
    274265        for (psS32 k = 0; k < params->n; k++) {
    275             alpha->data.F64[j][k] = 0;
    276         }
    277         beta->data.F64[j] = 0;
     266            alpha->data.F32[j][k] = 0;
     267        }
     268        beta->data.F32[j] = 0;
    278269    }
    279270    chisq = 0.0;
    … …  
    301292                    continue;
    302293                }
    303                 alpha->data.F64[j][k] += weight * deriv->data.F32[k];
     294                alpha->data.F32[j][k] += weight * deriv->data.F32[k];
    304295            }
    305             beta->data.F64[j] += weight * delta;
     296            beta->data.F32[j] += weight * delta;
    306297        }
    307298    }
    … …  
    310301    for (psS32 j = 1; j < params->n; j++) {
    311302        for (psS32 k = 0; k < j; k++) {
    312             alpha->data.F64[k][j] = alpha->data.F64[j][k];
     303            alpha->data.F32[k][j] = alpha->data.F32[j][k];
    313304        }
    314305    }
    … …  
    318309        for (psS32 j = 0; j < params->n; j++) {
    319310            if (paramMask->data.U8[j]) {
    320                 alpha->data.F64[j][j] = 1;
    321                 beta->data.F64[j] = 1;
     311                alpha->data.F32[j][j] = 1;
     312                beta->data.F32[j] = 1;
    322313            }
    323314        }
    … …  
    335326coords.
    336327 
    337 XXX: This must work for both F32 and F64.  F32 is currently implemented.
    338      Note: since the LUD routines are only implemented in F64, then we
    339      will have to convert all F32 input vectors to F64 regardless.  So,
    340      the F64 port might be an optimization.
    341  
     328This requires F32 input data; all internal calls use F32.
     329XXX Make an F64 version?
    342330  *****************************************************************************/
    343331psBool psMinimizeLMChi2(
    … …  
    345333    psImage *covar,
    346334    psVector *params,
    347     psMinConstrain *constrain,
     335    psMinConstraint *constraint,
    348336    const psArray *x,
    349337    const psVector *y,
    … …  
    357345    PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, false);
    358346    psVector *paramMask = NULL;
    359     psVector *paramDelta = NULL;
    360     psVector *paramMin = NULL;
    361     psVector *paramMax = NULL;
    362     if (constrain != NULL) {
    363         paramDelta = constrain->paramDelta;
    364         paramMin = constrain->paramMin;
    365         paramMax = constrain->paramMax;
    366         paramMask = constrain->paramMask;
    367         PS_ASSERT_VECTOR_TYPE(paramDelta, PS_TYPE_F32, false);
    368         PS_ASSERT_VECTOR_TYPE(paramMin, PS_TYPE_F32, false);
    369         PS_ASSERT_VECTOR_TYPE(paramMax, PS_TYPE_F32, false);
    370         PS_ASSERT_VECTORS_SIZE_EQUAL(params, paramDelta, false);
    371         PS_ASSERT_VECTORS_SIZE_EQUAL(params, paramMin, false);
    372         PS_ASSERT_VECTORS_SIZE_EQUAL(params, paramMax, false);
     347    if (constraint != NULL) {
     348        paramMask = constraint->paramMask;
    373349        if (paramMask != NULL) {
    374350            PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_U8, false);
    … …  
    392368    PS_ASSERT_PTR_NON_NULL(func, false);
    393369
     370    psMinimizeLMLimitFunc checkLimits = NULL;
     371    if (constraint) {
     372        checkLimits = constraint->checkLimits;
     373    }
     374
    394375    // this function has test and current values for several things
    395376    // the current best value is in lower case
    … …  
    397378
    398379    // allocate internal arrays (current vs Guess)
    399     psImage *alpha   = psImageAlloc(params->n, params->n, PS_TYPE_F64);
    400     psImage *Alpha   = psImageAlloc(params->n, params->n, PS_TYPE_F64);
    401     psVector *beta   = psVectorAlloc(params->n, PS_TYPE_F64);
    402     psVector *Beta   = psVectorAlloc(params->n, PS_TYPE_F64);
     380    psImage *alpha   = psImageAlloc(params->n, params->n, PS_TYPE_F32);
     381    psImage *Alpha   = psImageAlloc(params->n, params->n, PS_TYPE_F32);
     382    psVector *beta   = psVectorAlloc(params->n, PS_TYPE_F32);
     383    psVector *Beta   = psVectorAlloc(params->n, PS_TYPE_F32);
    403384    psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32);
    404385    psVector *dy     = NULL;
    405     psF64 Chisq = 0.0;
    406     psF64 lambda = 0.001;
     386    psF32 Chisq = 0.0;
     387    psF32 lambda = 0.001;
    407388
    408389    // the user provides the error or NULL.  we need to convert
    … …  
    436417
    437418        // set a new guess for Alpha, Beta, Params
    438         if (!p_psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask,
    439                                 paramDelta, paramMin, paramMax, lambda)) {
     419        if (!p_psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda)) {
    440420            min->iter ++;
    441421            lambda *= 10.0;
    … …  
    444424
    445425        // measure linear model prediction
    446         psF64 dLinear = p_psMinLM_dLinear(Beta, beta, lambda);
     426        psF32 dLinear = p_psMinLM_dLinear(Beta, beta, lambda);
    447427
    448428        // dump some useful info if trace is defined
    … …  
    467447        // expected delta from the linear model (dLinear)
    468448        // accept new guess if it is an improvement (rho > 0), or else increase lambda
    469         psF64 rho = (min->value - Chisq) / dLinear;
     449        psF32 rho = (min->value - Chisq) / dLinear;
    470450
    471451        psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value,
    … …  
    482462            min->lastDelta = (min->value - Chisq) / (dy->n - params->n);
    483463            min->value = Chisq;
    484             alpha  = psImageCopy(alpha, Alpha, PS_TYPE_F64);
    485             beta   = psVectorCopy(beta, Beta, PS_TYPE_F64);
     464            alpha  = psImageCopy(alpha, Alpha, PS_TYPE_F32);
     465            beta   = psVectorCopy(beta, Beta, PS_TYPE_F32);
    486466            params = psVectorCopy(params, Params, PS_TYPE_F32);
    487467            lambda *= 0.1;
    … …  
    495475    // construct & return the covariance matrix (if requested)
    496476    if (covar != NULL) {
    497         if (!p_psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask,
    498                                 paramDelta, paramMin, paramMax, 0.0)) {
     477        if (!p_psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0)) {
    499478            psTrace ("psLib.math", 5, "failure to calculate covariance matrix\n");
    500479        }
    … …  
    546525
    547526
    548 static void constrainFree(psMinConstrain *tmp)
    549 {
    550     // There are no dynamically allocated items
    551 }
    552 
    553 psMinConstrain* psMinConstrainAlloc()
    554 {
    555     psMinConstrain *tmp = psAlloc(sizeof(psMinConstrain));
     527static void constraintFree(psMinConstraint *tmp)
     528{
     529    if (tmp == NULL)
     530        return;
     531
     532    psFree (tmp->paramMask);
     533}
     534
     535psMinConstraint* psMinConstraintAlloc()
     536{
     537    psMinConstraint *tmp = psAlloc(sizeof(psMinConstraint));
     538    psMemSetDeallocator(tmp, (psFreeFunc)constraintFree);
    556539    tmp->paramMask = NULL;
    557     tmp->paramMax = NULL;
    558     tmp->paramMin = NULL;
    559     tmp->paramDelta = NULL;
     540    tmp->checkLimits = NULL;
    560541
    561542    return(tmp);
    562543}
    563544
    564 bool psMemCheckConstrain(psPtr tmp)
    565 {
    566     return(psMemGetDeallocator(tmp) == (psFreeFunc) constrainFree);
    567 }
     545bool psMemCheckConstraint(psPtr tmp)
     546{
     547    return(psMemGetDeallocator(tmp) == (psFreeFunc) constraintFree);
     548}
Note: See TracChangeset for help on using the changeset viewer.