IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 28, 2005, 4:13:54 PM (21 years ago)
Author:
gusciora
Message:

I added EAM's new version of LM minimization.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r5113 r5175  
    1010 *  @author EAM, IfA
    1111 *
    12  *  @version $Revision: 1.140 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2005-09-23 23:01:30 $
     12 *  @version $Revision: 1.141 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2005-09-29 02:13:54 $
    1414 *
    1515 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    1616 *
    1717 *  XXX: must follow coding name standards on local functions.
    18  *
     18 *  XXX: put local functions in front.
     19 *
    1920 */
    2021/*****************************************************************************/
    … …  
    3334/*****************************************************************************/
    3435/* DEFINE STATEMENTS                                                         */
    35 /* What are the following macros for?                                        */
    3636/*****************************************************************************/
    37 #define PS_SEG psLib
    38 #define PS_PWD dataManip
    39 #define PS_FILE psMinimize
     37
    4038/*****************************************************************************/
    4139/* TYPE DEFINITIONS                                                          */
    … …  
    4442/*****************************************************************************/
    4543/* GLOBAL VARIABLES                                                          */
     44/* XXX: Do these conform to code standard?         */
    4645/*****************************************************************************/
    4746static psMinimizeChi2PowellFunc Chi2PowellFunc = NULL;
    4847static psVector *myValue;
    4948static psVector *myError;
    50 
    5149/*****************************************************************************/
    5250/* FILE STATIC VARIABLES                                                     */
    … …  
    6462 ******************************************************************************
    6563 *****************************************************************************/
    66 // measure the distance to the minimum assuming a linear model
    67 // XXX: Move this to the pub area.
    68 bool psMinimizeGaussNewtonDelta (psVector *delta,
    69                                  const psVector *params,
    70                                  const psVector *paramMask,
    71                                  const psArray  *x,
    72                                  const psVector *y,
    73                                  const psVector *yErr,
    74                                  psMinimizeLMChi2Func func)
     64// XXX EAM : can we use static copies of LUv, LUm, A?
     65psBool p_psMinLM_GuessABP(
     66    psImage  *Alpha,
     67    psVector *Beta,
     68    psVector *Params,
     69    const psImage  *alpha,
     70    const psVector *beta,
     71    const psVector *params,
     72    const psVector *paramMask,
     73    const psVector *beta_lim,
     74    const psVector *params_min,
     75    const psVector *params_max,
     76    psF64 lambda)
     77{
     78    # define USE_LU_DECOMP 1
     79    # if (USE_LU_DECOMP)
     80        psVector *LUv = NULL;
     81    psImage  *LUm = NULL;
     82    psImage  *A   = NULL;
     83    psF32    det;
     84
     85    // LU decomposition version
     86    psTrace(__func__, 5, "using LUD version\n");
     87
     88    // set new guess values (creates matrix A)
     89    A = psImageCopy(NULL, alpha, PS_TYPE_F64);
     90    for (int j = 0; j < params->n; j++) {
     91        if ((paramMask != NULL) && (paramMask->data.U8[j]))
     92            continue;
     93        A->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
     94    }
     95
     96    // solve A*beta = Beta (Alpha = 1/A)
     97    // these operations do not modify the input values (creates LUm, LUv)
     98    LUm   = psMatrixLUD(NULL, &LUv, A);
     99    Beta  = psMatrixLUSolve(Beta, LUm, beta, LUv);
     100    Alpha = psMatrixInvert(Alpha, A, &det);
     101
     102    # else
     103        // gauss-jordan version
     104        psTrace(__func__, 5, "using Gauss-J version");
     105
     106    // set new guess values (creates matrix A)
     107    Beta = psVectorCopy(Beta, beta, PS_TYPE_F64);
     108    Alpha = psImageCopy(Alpha, alpha, PS_TYPE_F64);
     109    for (int j = 0; j < params->n; j++) {
     110        if ((paramMask != NULL) && (paramMask->data.U8[j]))
     111            continue;
     112        Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
     113    }
     114
     115    psGaussJordan(Alpha, Beta);
     116    psFree(A);
     117    psFree(LUm);
     118    psFree(LUv);
     119    # endif
     120
     121    // apply Beta to get new Params values
     122    for (int j = 0; j < params->n; j++) {
     123        if ((paramMask != NULL) && (paramMask->data.U8[j]))
     124            continue;
     125        // Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j];
     126        // compare Beta to beta limits
     127        if (beta_lim != NULL) {
     128            if (fabs(Beta->data.F64[j]) > fabs(beta_lim->data.F32[j])) {
     129                Beta->data.F64[j] = (Beta->data.F64[j] > 0) ? fabs(beta_lim->data.F32[j]) : -fabs(beta_lim->data.F32[j]);
     130            }
     131        }
     132        Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j];
     133        // compare new params to param limits
     134        if (params_max != NULL) {
     135            Params->data.F32[j] = PS_MIN (Params->data.F32[j], params_max->data.F32[j]);
     136        }
     137        if (params_min != NULL) {
     138            Params->data.F32[j] = PS_MAX (Params->data.F32[j], params_min->data.F32[j]);
     139        }
     140    }
     141    # if (USE_LU_DECOMP)
     142        psFree(A);
     143    psFree(LUm);
     144    psFree(LUv);
     145    # endif
     146
     147    return(true);
     148}
     149
     150
     151bool psMinimizeGaussNewtonDelta(
     152    psVector *delta,
     153    const psVector *params,
     154    const psVector *paramMask,
     155    const psArray  *x,
     156    const psVector *y,
     157    const psVector *yWt,
     158    psMinimizeLMChi2Func func)
    75159{
    76160
    77161    // allocate internal arrays (current vs Guess)
    78     psImage  *alpha  = psImageAlloc  (params->n, params->n, PS_TYPE_F64);
    79     psImage  *Alpha  = psImageAlloc  (params->n, params->n, PS_TYPE_F64);
    80     psVector *beta   = psVectorAlloc (params->n, PS_TYPE_F64);
    81     psVector *Params = psVectorAlloc (params->n, PS_TYPE_F64);
    82     psVector *dy     = psVectorAlloc (y->n, PS_TYPE_F32);
     162    psImage  *alpha  = psImageAlloc (params->n, params->n, PS_TYPE_F64);
     163    psImage  *Alpha  = psImageAlloc (params->n, params->n, PS_TYPE_F64);
     164    psVector *beta   = psVectorAlloc(params->n, PS_TYPE_F64);
     165    psVector *Params = psVectorAlloc(params->n, PS_TYPE_F64);
     166    psVector *dy     = NULL;
    83167
    84168    // the user provides the error or NULL.  we need to convert
    85169    // to appropriate weights
    86     if (yErr != NULL) {
    87         for (int i = 0; i < dy->n; i++) {
    88             dy->data.F32[i] = 1.0 / PS_SQR (yErr->data.F32[i]);
    89         }
     170    if (yWt != NULL) {
     171        dy = (psVector *) yWt;
    90172    } else {
    91         for (int i = 0; i < dy->n; i++) {
    92             dy->data.F32[i] = 1.0;
    93         }
    94     }
    95 
    96     p_psMinLM_SetABX (alpha, beta, params, paramMask, x, y, dy, func);
    97     p_psMinLM_GuessABP (Alpha, delta, Params, alpha, beta, params, paramMask, 0.0);
    98 
    99     psFree (alpha);
    100     psFree (Alpha);
    101     psFree (beta);
    102     psFree (Params);
    103     psFree (dy);
     173        dy = psVectorAlloc(y->n, PS_TYPE_F32);
     174        psVectorInit(dy, 1.0);
     175    }
     176
     177    p_psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func);
     178    p_psMinLM_GuessABP(Alpha, delta, Params, alpha, beta, params, paramMask, NULL, NULL, NULL, 0.0);
     179
     180    psFree(alpha);
     181    psFree(Alpha);
     182    psFree(beta);
     183    psFree(Params);
     184    if (yWt == NULL) {
     185        psFree(dy);
     186    }
    104187    return (true);
    105188}
    106189
    107 
    108190// measure linear model prediction
    109 psF64 p_psMinLM_dLinear (const psVector *Beta, const psVector *beta, psF64 lambda)
     191psF64 p_psMinLM_dLinear(
     192    const psVector *Beta,
     193    const psVector *beta,
     194    psF64 lambda)
    110195{
    111196
    … …  
    117202        dLinear += lambda*PS_SQR(B[i]) + B[i]*b[i];
    118203    }
    119     return (0.5*dLinear);
    120 }
     204    return(0.5*dLinear);
     205}
     206
     207// XXX EAM: this needs to respect the mask on params
     208// alpha, beta, params are already allocated
     209psF64 p_psMinLM_SetABX(
     210    psImage  *alpha,
     211    psVector *beta,
     212    const psVector *params,
     213    const psVector *paramMask,
     214    const psArray  *x,
     215    const psVector *y,
     216    const psVector *dy,
     217    psMinimizeLMChi2Func func)
     218{
     219    PS_ASSERT_IMAGE_NON_NULL(alpha, NAN);
     220    PS_ASSERT_VECTOR_NON_NULL(beta, NAN);
     221    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
     222    PS_ASSERT_PTR_NON_NULL(x, NAN);
     223    PS_ASSERT_VECTOR_NON_NULL(y, NAN);
     224    PS_ASSERT_VECTOR_NON_NULL(dy, NAN);
     225
     226    psF64 chisq;
     227    psF64 delta;
     228    psF64 weight;
     229    psF64 ymodel;
     230    psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32);
     231
     232    // zero alpha and beta for summing below
     233    for (int j = 0; j < params->n; j++) {
     234        for (int k = 0; k < params->n; k++) {
     235            alpha->data.F64[j][k] = 0;
     236        }
     237        beta->data.F64[j] = 0;
     238    }
     239    chisq = 0.0;
     240
     241    // calculate chisq, alpha, beta
     242    for (int i = 0; i < y->n; i++) {
     243        ymodel = func(deriv, params, (psVector *) x->data[i]);
     244
     245        delta = ymodel - y->data.F32[i];
     246        chisq += PS_SQR(delta) * dy->data.F32[i];
     247
     248        for (int j = 0; j < params->n; j++) {
     249            if ((paramMask != NULL) && (paramMask->data.U8[j]))
     250                continue;
     251            weight = deriv->data.F32[j] * dy->data.F32[i];
     252            for (int k = 0; k <= j; k++) {
     253                if ((paramMask != NULL) && (paramMask->data.U8[k]))
     254                    continue;
     255                alpha->data.F64[j][k] += weight * deriv->data.F32[k];
     256            }
     257            beta->data.F64[j] += weight * delta;
     258        }
     259    }
     260
     261    // calculate lower-left half of alpha
     262    for (int j = 1; j < params->n; j++) {
     263        for (int k = 0; k < j; k++) {
     264            alpha->data.F64[k][j] = alpha->data.F64[j][k];
     265        }
     266    }
     267
     268    // fill in pivots if we apply a mask
     269    if (paramMask != NULL) {
     270        for (int j = 0; j < params->n; j++) {
     271            if (paramMask->data.U8[j]) {
     272                alpha->data.F64[j][j] = 1;
     273                beta->data.F64[j] = 1;
     274            }
     275        }
     276    }
     277
     278    psFree(deriv);
     279    return(chisq);
     280}
     281
    121282
    122283/******************************************************************************
    123 psMinimizeLMChi2():  This routine will take an procedure which calculates
    124 an arbitrary function and it's derivative and minimize the chi-squared match
    125 between that function at the specified coords and the specified value at
    126 those coords.
     284psMinimizeLMChi2():  This routine will take an procedure which calculates an
     285arbitrary function and it's derivative and minimize the chi-squared match
     286between that function at the specified coords and the specified value at those
     287coords.
     288 
     289XXX: Put the ASSERTS in.
    127290 
    128291XXX EAM this is my re-implementation of MinLM
    … …  
    135298XXX: Change the whole thing to F64, if input data is F32, convert it.
    136299 *****************************************************************************/
    137 psBool psMinimizeLMChi2(psMinimization *min,
    138                         psImage *covar,
    139                         psVector *params,
    140                         const psVector *paramMask,
    141                         const psArray *x,
    142                         const psVector *y,
    143                         const psVector *yErr,
    144                         psMinimizeLMChi2Func func)
    145 {
     300psBool psMinimizeLMChi2(
     301    psMinimization *min,
     302    psImage *covar,
     303    psVector *params,
     304    const psVector *paramMask,
     305    const psArray *x,
     306    const psVector *y,
     307    const psVector *yWt,
     308    psMinimizeLMChi2Func func)
     309{
     310    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
    146311    PS_ASSERT_PTR_NON_NULL(min, false);
     312    // XXX: If covar not NULL, do asserts...
    147313    PS_ASSERT_VECTOR_NON_NULL(params, false);
    148314    PS_ASSERT_VECTOR_NON_EMPTY(params, false);
    … …  
    162328    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false);
    163329    PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, false);
    164     if (yErr != NULL) {
    165         PS_ASSERT_VECTOR_TYPE(yErr, PS_TYPE_F32, false);
    166         PS_ASSERT_VECTORS_SIZE_EQUAL(y, yErr, false);
     330    if (yWt != NULL) {
     331        PS_ASSERT_VECTOR_TYPE(yWt, PS_TYPE_F32, false);
     332        PS_ASSERT_VECTORS_SIZE_EQUAL(y, yWt, false);
    167333    }
    168334    PS_ASSERT_PTR_NON_NULL(func, false);
    … …  
    173339
    174340    // allocate internal arrays (current vs Guess)
    175     psImage *alpha   = psImageAlloc (params->n, params->n, PS_TYPE_F64);
    176     psImage *Alpha   = psImageAlloc (params->n, params->n, PS_TYPE_F64);
    177     psVector *beta   = psVectorAlloc (params->n, PS_TYPE_F64);
    178     psVector *Beta   = psVectorAlloc (params->n, PS_TYPE_F64);
    179     psVector *Params = psVectorAlloc (params->n, PS_TYPE_F32);
     341    psImage *alpha   = psImageAlloc(params->n, params->n, PS_TYPE_F64);
     342    psImage *Alpha   = psImageAlloc(params->n, params->n, PS_TYPE_F64);
     343    psVector *beta   = psVectorAlloc(params->n, PS_TYPE_F64);
     344    psVector *Beta   = psVectorAlloc(params->n, PS_TYPE_F64);
     345    psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32);
    180346    psVector *dy     = NULL;
    181347    psF64 Chisq = 0.0;
    182348    psF64 lambda = 0.001;
    183 
    184     // XXX EAM: why is this needed here? the value is not used, and the memory
    185     // is allocated above.  However, if I drop it, I get weird answers or
    186     // crashes.
    187     Params = psVectorCopy (Params, params, PS_TYPE_F32);
     349    // XXX: Code this properly.  Don't use mustFree00.
     350    psBool mustFree00 = false;
     351    psVector *beta_lim = NULL;
     352    psVector *param_min = NULL;
     353    psVector *param_max = NULL;
     354
     355    // if we are provided a covar image, we expect to find these three vectors in first three rows
     356    if (covar != NULL) {
     357        mustFree00 = true;
     358        beta_lim  = psVectorAlloc(params->n, PS_TYPE_F32);
     359        param_min = psVectorAlloc(params->n, PS_TYPE_F32);
     360        param_max = psVectorAlloc(params->n, PS_TYPE_F32);
     361        for (int i = 0; i < params->n; i++) {
     362            beta_lim->data.F32[i] = covar->data.F64[0][i];
     363            param_min->data.F32[i] = covar->data.F64[1][i];
     364            param_max->data.F32[i] = covar->data.F64[2][i];
     365        }
     366        psImageRecycle(covar, params->n, params->n, PS_TYPE_F64);
     367    }
     368
     369    // why is this needed here??? the initial guess on params is provided by the user
     370    Params = psVectorCopy(Params, params, PS_TYPE_F32);
    188371
    189372    // the user provides the error or NULL.  we need to convert
    190373    // to appropriate weights
    191     dy = psVectorAlloc (y->n, PS_TYPE_F32);
    192     if (yErr != NULL) {
    193         for (int i = 0; i < dy->n; i++) {
    194             dy->data.F32[i] = 1.0 / PS_SQR (yErr->data.F32[i]);
    195         }
     374    if (yWt != NULL) {
     375        dy = (psVector *) yWt;
    196376    } else {
    197         for (int i = 0; i < dy->n; i++) {
    198             dy->data.F32[i] = 1.0;
    199         }
     377        dy = psVectorAlloc(y->n, PS_TYPE_F32);
     378        psVectorInit(dy, 1.0);
    200379    }
    201380
    202381    // calculate initial alpha and beta, set chisq (min->value)
    203     min->value = p_psMinLM_SetABX (alpha, beta, params, paramMask, x, y, dy, func);
    204     # ifndef PS_NO_TRACE
     382    min->value = p_psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func);
     383    if (isnan(min->value)) {
     384        min->iter = min->maxIter;
     385        return(false);
     386    }
    205387    // dump some useful info if trace is defined
    206     /* XXX: This code is seg faulting.
    207         if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) {
    208             p_psImagePrint  (psTraceGetDestination(), alpha, "alpha guess");
    209             p_psVectorPrint (psTraceGetDestination(), beta, "beta guess");
    210             p_psVectorPrint (psTraceGetDestination(), params, "params guess");
    211         }
    212         if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 4) {
    213             // XXX: Where is this?
    214             // p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess");
    215         }
    216     */
    217     # endif /* PS_NO_TRACE */
    218 
    219     min->lastDelta = min->tol + 1.0;
     388    if (psTraceGetLevel(__func__) >= 6) {
     389        p_psImagePrint(psTraceGetDestination(), alpha, "alpha guess (0)");
     390        p_psVectorPrint(psTraceGetDestination(), beta, "beta guess (0)");
     391        p_psVectorPrint(psTraceGetDestination(), params, "params guess (0)");
     392    }
     393    if (psTraceGetLevel (__func__) >= 6) {
     394        //XXX:  p_psVectorPrintRow(psTraceGetDestination(), Params, "params guess");
     395    }
     396
    220397    // iterate until the tolerance is reached, or give up
    221     while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {
     398    while ((min->iter < min->maxIter) && ((min->lastDelta > min->tol) || !isfinite(min->lastDelta))) {
     399        psTrace(__func__, 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
     400        psTrace(__func__, 5, "Last delta is %f.  Min->tol is %f.\n", min->lastDelta, min->tol);
     401
    222402        // set a new guess for Alpha, Beta, Params
    223         p_psMinLM_GuessABP (Alpha, Beta, Params, alpha, beta, params, paramMask, lambda);
     403        p_psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask,
     404                           beta_lim, param_min, param_max, lambda);
    224405
    225406        // measure linear model prediction
    226         psF64 dLinear = p_psMinLM_dLinear (Beta, beta, lambda);
    227 
    228         # ifndef PS_NO_TRACE
     407        psF64 dLinear = p_psMinLM_dLinear(Beta, beta, lambda);
     408
    229409        // dump some useful info if trace is defined
    230         /* XXX: This code is seg faulting.
    231                 if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) {
    232                     p_psImagePrint  (psTraceGetDestination(), Alpha, "alpha guess");
    233                     p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess");
    234                     p_psVectorPrint (psTraceGetDestination(), Params, "params guess");
    235                 }
    236                 if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 4) {
    237                     // XXX: Where is this?
    238                     // p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess");
    239                 }
    240         */
    241         # endif /* PS_NO_TRACE */
     410        if (psTraceGetLevel(__func__) >= 6) {
     411            p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (1)");
     412            p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (1)");
     413            p_psVectorPrint(psTraceGetDestination(), Params, "params guess (1)");
     414        }
     415        if (psTraceGetLevel(__func__) >= 6) {
     416            //XXX: p_psVectorPrintRow(psTraceGetDestination(), Params, "params guess");
     417        }
    242418
    243419        // calculate Chisq for new guess, update Alpha & Beta
    244         Chisq = p_psMinLM_SetABX (Alpha, Beta, Params, paramMask, x, y, dy, func);
     420        Chisq = p_psMinLM_SetABX(Alpha, Beta, Params, paramMask, x, y, dy, func);
    245421
    246422        // XXX EAM alternate convergence criterion:
    … …  
    250426        psF64 rho = (min->value - Chisq) / dLinear;
    251427
    252         psTrace (".psLib.dataManip.psMinimizeLMChi2", 4, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value, Chisq, min->lastDelta, rho);
    253         # ifndef PS_NO_TRACE
     428        psTrace(__func__, 5, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value,
     429                Chisq, min->lastDelta, rho);
     430
    254431        // dump some useful info if trace is defined
    255         /* XXX: This code is seg faulting.
    256                 if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) {
    257                     p_psImagePrint  (psTraceGetDestination(), Alpha, "alpha guess");
    258                     p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess");
    259                     p_psVectorPrint (psTraceGetDestination(), Params, "params guess");
    260                 }
    261         */
    262         # endif /* PS_NO_TRACE */
     432        if (psTraceGetLevel(__func__) >= 6) {
     433            p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (2)");
     434            p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (2)");
     435            p_psVectorPrint(psTraceGetDestination(), Params, "params guess (2)");
     436        }
    263437
    264438        /* if (Chisq < min->value) {  */
    … …  
    266440            min->lastDelta = (min->value - Chisq) / (dy->n - params->n);
    267441            min->value = Chisq;
    268             alpha  = psImageCopy (alpha, Alpha, PS_TYPE_F64);
    269             beta   = psVectorCopy (beta, Beta, PS_TYPE_F64);
    270             params = psVectorCopy (params, Params, PS_TYPE_F32);
     442            alpha  = psImageCopy(alpha, Alpha, PS_TYPE_F64);
     443            beta   = psVectorCopy(beta, Beta, PS_TYPE_F64);
     444            params = psVectorCopy(params, Params, PS_TYPE_F32);
    271445            lambda *= 0.1;
    272446        } else {
    273447            lambda *= 10.0;
    274448        }
    275         min->iter ++;
    276 
    277         //        printf("CONDITIONS: (%f > %f) && (%d < %d)\n", min->lastDelta, min->tol, min->iter, min->maxIter);
    278     }
    279     psTrace (".psLib.dataManip.psMinimizeLMChi2", 4, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter);
     449        min->iter++;
     450    }
     451    psTrace(__func__, 5, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter);
    280452
    281453    // construct & return the covariance matrix (if requested)
    282454    if (covar != NULL) {
    283         p_psMinLM_GuessABP (covar, Beta, Params, alpha, beta, params, paramMask, 0.0);
     455        p_psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask,
     456                           beta_lim, param_min, param_max, 0.0);
    284457    }
    285458
    286459    // free the internal temporary data
    287     psFree (alpha);
    288     psFree (Alpha);
    289     psFree (beta);
    290     psFree (Beta);
    291     psFree (Params);
    292     psFree (dy);
    293 
     460    psFree(alpha);
     461    psFree(Alpha);
     462    psFree(beta);
     463    psFree(Beta);
     464    psFree(Params);
     465    if (yWt == NULL) {
     466        psFree(dy);
     467    }
     468    if (mustFree00 == true) {
     469        psFree(beta_lim);
     470        psFree(param_min);
     471        psFree(param_max);
     472    }
    294473    if (min->iter == min->maxIter) {
    295         return (false);
    296     }
    297     return (true);
    298 }
    299 
    300 // XXX EAM: this needs to respect the mask on params
    301 // alpha, beta, params are already allocated
    302 psF64 p_psMinLM_SetABX (psImage  *alpha,
    303                         psVector *beta,
    304                         const psVector *params,
    305                         const psVector *paramMask,
    306                         const psArray  *x,
    307                         const psVector *y,
    308                         const psVector *dy,
    309                         psMinimizeLMChi2Func func)
    310 {
    311     PS_ASSERT_IMAGE_NON_NULL(alpha, NAN);
    312     PS_ASSERT_VECTOR_NON_NULL(beta, NAN);
    313     PS_ASSERT_VECTOR_NON_NULL(params, NAN);
    314     PS_ASSERT_PTR_NON_NULL(x, NAN);
    315     PS_ASSERT_VECTOR_NON_NULL(y, NAN);
    316     PS_ASSERT_VECTOR_NON_NULL(dy, NAN);
    317 
    318     psF64 chisq;
    319     psF64 delta;
    320     psF64 weight;
    321     psF64 ymodel;
    322     psVector *deriv = psVectorAlloc (params->n, PS_TYPE_F32);
    323 
    324     // zero alpha and beta for summing below
    325     for (int j = 0; j < params->n; j++) {
    326         for (int k = 0; k < params->n; k++) {
    327             alpha->data.F64[j][k] = 0;
    328         }
    329         beta->data.F64[j] = 0;
    330     }
    331     chisq = 0.0;
    332 
    333     // calculate chisq, alpha, beta
    334     for (int i = 0; i < y->n; i++) {
    335         ymodel = func (deriv, params, (psVector *) x->data[i]);
    336 
    337         delta = ymodel - y->data.F32[i];
    338         chisq += PS_SQR (delta) * dy->data.F32[i];
    339 
    340         for (int j = 0; j < params->n; j++) {
    341             if ((paramMask != NULL) && (paramMask->data.U8[j]))
    342                 continue;
    343             weight = deriv->data.F32[j] * dy->data.F32[i];
    344             for (int k = 0; k <= j; k++) {
    345                 if ((paramMask != NULL) && (paramMask->data.U8[k]))
    346                     continue;
    347                 alpha->data.F64[j][k] += weight * deriv->data.F32[k];
    348             }
    349             beta->data.F64[j] += weight * delta;
    350         }
    351     }
    352 
    353     // calculate lower-left half of alpha
    354     for (int j = 1; j < params->n; j++) {
    355         for (int k = 0; k < j; k++) {
    356             alpha->data.F64[k][j] = alpha->data.F64[j][k];
    357         }
    358     }
    359 
    360     // fill in pivots if we apply a mask
    361     if (paramMask != NULL) {
    362         for (int j = 0; j < params->n; j++) {
    363             if (paramMask->data.U8[j]) {
    364                 alpha->data.F64[j][j] = 1;
    365                 beta->data.F64[j] = 1;
    366             }
    367         }
    368     }
    369 
    370     psFree (deriv);
    371     return (chisq);
    372 }
    373 
    374 // XXX EAM : can we use static copies of LUv, LUm, A?
    375 psBool p_psMinLM_GuessABP (psImage  *Alpha,
    376                            psVector *Beta,
    377                            psVector *Params,
    378                            const psImage  *alpha,
    379                            const psVector *beta,
    380                            const psVector *params,
    381                            const psVector *paramMask,
    382                            psF64 lambda)
    383 {
    384 
    385     # define USE_LU_DECOMP 1
    386     # if (USE_LU_DECOMP)
    387         psVector *LUv = NULL;
    388     psImage  *LUm = NULL;
    389     psImage  *A   = NULL;
    390     psF32    det;
    391 
    392     // LU decomposition version
    393     psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using LUD version\n");
    394 
    395     // set new guess values (creates matrix A)
    396     A = psImageCopy (NULL, alpha, PS_TYPE_F64);
    397     for (int j = 0; j < params->n; j++) {
    398         if ((paramMask != NULL) && (paramMask->data.U8[j]))
    399             continue;
    400         A->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
    401     }
    402 
    403     // solve A*beta = Beta (Alpha = 1/A)
    404     // these operations do not modify the input values (creates LUm, LUv)
    405     LUm   = psMatrixLUD (NULL, &LUv, A);
    406     Beta  = psMatrixLUSolve (Beta, LUm, beta, LUv);
    407     Alpha = psMatrixInvert (Alpha, A, &det);
    408 
    409     # else
    410         // gauss-jordan version
    411         psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using Gauss-J version");
    412 
    413     // set new guess values (creates matrix A)
    414     Beta = psVectorCopy (Beta, beta, PS_TYPE_F64);
    415     Alpha = psImageCopy (Alpha, alpha, PS_TYPE_F64);
    416     for (int j = 0; j < params->n; j++) {
    417         if ((paramMask != NULL) && (paramMask->data.U8[j]))
    418             continue;
    419         Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
    420     }
    421 
    422     psGaussJordan (Alpha, Beta);
    423     # endif
    424 
    425     // apply Beta to get new Params values
    426     for (int j = 0; j < params->n; j++) {
    427         if ((paramMask != NULL) && (paramMask->data.U8[j]))
    428             continue;
    429         Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j];
    430     }
    431 
    432     # if (USE_LU_DECOMP)
    433         psFree (A);
    434     psFree (LUm);
    435     psFree (LUv);
    436     # endif
    437 
    438     return true;
     474        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
     475        return(false);
     476    }
     477    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
     478    return(true);
    439479}
    440480
    441481// XXX EAM : temporary gauss-jordan solver based on gene's
    442482// version based on the Numerical Recipes version
    443 bool psGaussJordan (psImage *a, psVector *b)
     483bool psGaussJordan(
     484    psImage *a,
     485    psVector *b)
    444486{
    445487
    … …  
    455497    vector = b->data.F64;
    456498
    457     indxc = psAlloc (Nx*sizeof(int));
    458     indxr = psAlloc (Nx*sizeof(int));
    459     ipiv  = psAlloc (Nx*sizeof(int));
     499    indxc = psAlloc(Nx*sizeof(int));
     500    indxr = psAlloc(Nx*sizeof(int));
     501    ipiv  = psAlloc(Nx*sizeof(int));
    460502    for (j = 0; j < Nx; j++) {
    461503        ipiv[j] = 0;
    … …  
    476518                    if (ipiv[k] == 0) {
    477519                        if (fabs (matrix[j][k]) >= big) {
    478                             big  = fabs (matrix[j][k]);
     520                            big  = fabs(matrix[j][k]);
    479521                            irow = j;
    480522                            icol = k;
    … …  
    492534        if (irow != icol) {
    493535            for (l = 0; l < Nx; l++) {
    494                 PS_SWAP (matrix[irow][l], matrix[icol][l]);
    495             }
    496             PS_SWAP (vector[irow], vector[icol]);
     536                PS_SWAP(matrix[irow][l], matrix[icol][l]);
     537            }
     538            PS_SWAP(vector[irow], vector[icol]);
    497539        }
    498540        indxr[i] = irow;
    … …  
    524566        if (indxr[l] != indxc[l]) {
    525567            for (k = 0; k < Nx; k++) {
    526                 PS_SWAP (matrix[k][indxr[l]], matrix[k][indxc[l]]);
    527             }
    528         }
    529     }
    530     psFree (ipiv);
    531     psFree (indxr);
    532     psFree (indxc);
    533     return (true);
     568                PS_SWAP(matrix[k][indxr[l]], matrix[k][indxc[l]]);
     569            }
     570        }
     571    }
     572    psFree(ipiv);
     573    psFree(indxr);
     574    psFree(indxc);
     575    return(true);
    534576
    535577fescape:
    536     psFree (ipiv);
    537     psFree (indxr);
    538     psFree (indxc);
    539     return (false);
     578    psFree(ipiv);
     579    psFree(indxr);
     580    psFree(indxc);
     581    return(false);
    540582}
    541583
    … …  
    565607bool psMemCheckMinimization(psPtr ptr)
    566608{
    567     return ( psMemGetDeallocator(ptr) == (psFreeFunc)minimizationFree );
     609    return( psMemGetDeallocator(ptr) == (psFreeFunc)minimizationFree );
    568610}
    569611
Note: See TracChangeset for help on using the changeset viewer.