IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Aug 8, 2005, 11:42:07 AM (21 years ago)
Author:
gusciora
Message:

I incorporated EAMs modifications. There might still be bugs with this code.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r4580 r4730  
    99 *  @author GLG, MHPCC
    1010 *
    11  *  @version $Revision: 1.127 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2005-07-19 02:55:54 $
     11 *  @version $Revision: 1.128 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2005-08-08 21:42:07 $
    1313 *
    1414 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    5454/* FUNCTION IMPLEMENTATION - LOCAL                                           */
    5555/*****************************************************************************/
     56
     57// measure the distance to the minimum assuming a linear model
     58bool psMinimizeGaussNewtonDelta (psVector *delta,
     59                                 const psVector *params,
     60                                 const psVector *paramMask,
     61                                 const psArray  *x,
     62                                 const psVector *y,
     63                                 const psVector *yErr,
     64                                 psMinimizeLMChi2Func func)
     65{
     66
     67    // allocate internal arrays (current vs Guess)
     68    psImage  *alpha  = psImageAlloc  (params->n, params->n, PS_TYPE_F64);
     69    psImage  *Alpha  = psImageAlloc  (params->n, params->n, PS_TYPE_F64);
     70    psVector *beta   = psVectorAlloc (params->n, PS_TYPE_F64);
     71    psVector *Params = psVectorAlloc (params->n, PS_TYPE_F64);
     72    psVector *dy     = psVectorAlloc (y->n, PS_TYPE_F32);
     73
     74    // the user provides the error or NULL.  we need to convert
     75    // to appropriate weights
     76    if (yErr != NULL) {
     77        for (int i = 0; i < dy->n; i++) {
     78            dy->data.F32[i] = 1.0 / PS_SQR (yErr->data.F32[i]);
     79        }
     80    } else {
     81        for (int i = 0; i < dy->n; i++) {
     82            dy->data.F32[i] = 1.0;
     83        }
     84    }
     85
     86    p_psMinLM_SetABX (alpha, beta, params, paramMask, x, y, dy, func);
     87    p_psMinLM_GuessABP (Alpha, delta, Params, alpha, beta, params, paramMask, 0.0);
     88
     89    psFree (alpha);
     90    psFree (Alpha);
     91    psFree (beta);
     92    psFree (Params);
     93    psFree (dy);
     94    return (true);
     95}
     96
     97
     98// measure linear model prediction
     99psF64 p_psMinLM_dLinear (const psVector *Beta, const psVector *beta, psF64 lambda)
     100{
     101
     102    /* get linear model prediction */
     103    psF64 dLinear = 0;
     104    psF64 *B = Beta->data.F64;
     105    psF64 *b = beta->data.F64;
     106    for (int i = 0; i < beta->n; i++) {
     107        dLinear += lambda*PS_SQR(B[i]) + B[i]*b[i];
     108    }
     109    return (0.5*dLinear);
     110}
    56111
    57112/******************************************************************************
    … …  
    642697
    643698
     699
     700
    644701/******************************************************************************
    645702XXX: We assume unnormalized gaussians.
    … …  
    769826}
    770827
     828//XXX: What's this for?
    771829psF64 p_psImageGetElementF64(psImage *a, int i, int j);
    772 
     830/******************************************************************************
     831psMinimizeLMChi2():  This routine will take an procedure which calculates
     832an arbitrary function and it's derivative and minimize the chi-squared match
     833between that function at the specified coords and the specified value at
     834those coords.
     835 
    773836// XXX EAM this is my re-implementation of MinLM
    774 bool psMinimizeLMChi2(psMinimization *min,
    775                       psImage *covar,
    776                       psVector *params,
    777                       const psVector *paramMask,
    778                       const psArray *x,
    779                       const psVector *y,
    780                       const psVector *yErr,
    781                       psMinimizeLMChi2Func func)
     837 
     838XXX: This must work for both F32 and F64.  F32 is currently implemented.
     839     Note: since the LUD routines are only implemented in F64, then we
     840     will have to convert all F32 input vectors to F64 regardless.  So,
     841     the F64 port might be.
     842 *****************************************************************************/
     843psBool psMinimizeLMChi2(psMinimization *min,
     844                        psImage *covar,
     845                        psVector *params,
     846                        const psVector *paramMask,
     847                        const psArray *x,
     848                        const psVector *y,
     849                        const psVector *yErr,
     850                        psMinimizeLMChi2Func func)
    782851{
    783852    PS_ASSERT_PTR_NON_NULL(min, NULL);
    … …  
    799868    psVector *beta   = psVectorAlloc (params->n, PS_TYPE_F64);
    800869    psVector *Beta   = psVectorAlloc (params->n, PS_TYPE_F64);
    801     psVector *Params = psVectorAlloc (params->n, PS_TYPE_F64);
     870    psVector *Params = psVectorAlloc (params->n, PS_TYPE_F32);
    802871    psVector *dy     = NULL;
    803     psF64 chisq = 0.0;
    804872    psF64 Chisq = 0.0;
    805873    psF64 lambda = 0.001;
    806874
    807     // the initial guess on params is provided by the user
     875    // XXX EAM: why is this needed here? the value is not used, and the memory
     876    // is allocated above.  However, if I drop it, I get weird answers or
     877    // crashes.
    808878    Params = psVectorCopy (Params, params, PS_TYPE_F32);
    809879
    … …  
    822892
    823893    // calculate initial alpha and beta, set chisq (min->value)
    824     min->value = p_psMinLM_SetABX (alpha, beta, params, x, y, dy, func);
     894    min->value = p_psMinLM_SetABX (alpha, beta, params, paramMask, x, y, dy, func);
    825895    # ifndef PS_NO_TRACE
    826896    // dump some useful info if trace is defined
    827     if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") > 4) {
     897    if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) {
    828898        p_psImagePrint  (psTraceGetDestination(), alpha, "alpha guess");
    829899        p_psVectorPrint (psTraceGetDestination(), beta, "beta guess");
    830900        p_psVectorPrint (psTraceGetDestination(), params, "params guess");
    831901    }
     902    if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") == 4) {
     903        // XXX: Where is this?
     904        // p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess");
     905    }
    832906    # endif /* PS_NO_TRACE */
    833907
    834 
     908    min->lastDelta = min->tol + 1.0;
    835909    // iterate until the tolerance is reached, or give up
    836910    while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {
    837 
    838911        // set a new guess for Alpha, Beta, Params
    839         p_psMinLM_GuessABP (Alpha, Beta, Params, alpha, beta, params, lambda);
     912        p_psMinLM_GuessABP (Alpha, Beta, Params, alpha, beta, params, paramMask, lambda);
     913
     914        // measure linear model prediction
     915        psF64 dLinear = p_psMinLM_dLinear (Beta, beta, lambda);
    840916
    841917        # ifndef PS_NO_TRACE
    842918        // dump some useful info if trace is defined
    843         if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") > 4) {
     919        if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) {
    844920            p_psImagePrint  (psTraceGetDestination(), Alpha, "alpha guess");
    845921            p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess");
    846922            p_psVectorPrint (psTraceGetDestination(), Params, "params guess");
    847923        }
     924        if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") == 4) {
     925            // XXX: Where is this?
     926            //            p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess");
     927        }
    848928        # endif /* PS_NO_TRACE */
    849929
    850930        // calculate Chisq for new guess, update Alpha & Beta
    851         Chisq = p_psMinLM_SetABX (Alpha, Beta, Params, x, y, dy, func);
    852         psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, Chisq %f, delta: %f\n", chisq, Chisq, min->lastDelta);
    853 
     931        Chisq = p_psMinLM_SetABX (Alpha, Beta, Params, paramMask, x, y, dy, func);
     932
     933        // XXX EAM alternate convergence criterion:
     934        // compare the delta (min->value - Chisq) with the
     935        // expected delta from the linear model (dLinear)
    854936        // accept new guess (if improvement), or increase lambda
    855         if (Chisq < min->value) {
     937        psF64 rho = (min->value - Chisq) / dLinear;
     938
     939        psTrace (".psLib.dataManip.psMinimizeLMChi2", 4, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value, Chisq, min->lastDelta, rho);
     940        # ifndef PS_NO_TRACE
     941        // dump some useful info if trace is defined
     942        if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) {
     943            p_psImagePrint  (psTraceGetDestination(), Alpha, "alpha guess");
     944            p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess");
     945            p_psVectorPrint (psTraceGetDestination(), Params, "params guess");
     946        }
     947        # endif /* PS_NO_TRACE */
     948
     949        /* if (Chisq < min->value) {  */
     950        if (rho > 0.0) {
    856951            min->lastDelta = (min->value - Chisq) / (dy->n - params->n);
    857952            min->value = Chisq;
    … …  
    864959        }
    865960        min->iter ++;
    866     }
    867     psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, Chisq %f, delta: %f\n", chisq, Chisq, min->lastDelta);
     961
     962        //        printf("CONDITIONS: (%f > %f) && (%d < %d)\n", min->lastDelta, min->tol, min->iter, min->maxIter);
     963    }
     964    psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter);
     965
     966    // construct & return the covariance matrix (if requested)
     967    if (covar != NULL) {
     968        p_psMinLM_GuessABP (covar, Beta, Params, alpha, beta, params, paramMask, 0.0);
     969    }
    868970
    869971    // free the internal temporary data
    … …  
    874976    psFree (Params);
    875977    psFree (dy);
     978
     979    if (min->iter == min->maxIter) {
     980        return (false);
     981    }
    876982    return (true);
    877983}
    … …  
    882988psF64 p_psMinLM_SetABX (psImage  *alpha,
    883989                        psVector *beta,
    884                         psVector *params,
     990                        const psVector *params,
     991                        const psVector *paramMask,
    885992                        const psArray  *x,
    886993                        const psVector *y,
    … …  
    8981005    for (int j = 0; j < params->n; j++) {
    8991006        for (int k = 0; k < params->n; k++) {
    900             alpha->data.F64[j][k] = 0.0;
    901         }
    902         beta->data.F64[j] = 0.0;
     1007            alpha->data.F64[j][k] = 0;
     1008        }
     1009        beta->data.F64[j] = 0;
    9031010    }
    9041011    chisq = 0.0;
    … …  
    9121019
    9131020        for (int j = 0; j < params->n; j++) {
     1021            if ((paramMask != NULL) && (paramMask->data.U8[j]))
     1022                continue;
    9141023            weight = deriv->data.F32[j] * dy->data.F32[i];
    9151024            for (int k = 0; k <= j; k++) {
     1025                if ((paramMask != NULL) && (paramMask->data.U8[k]))
     1026                    continue;
    9161027                alpha->data.F64[j][k] += weight * deriv->data.F32[k];
    9171028            }
    … …  
    9261037        }
    9271038    }
     1039
     1040    // fill in pivots if we apply a mask
     1041    if (paramMask != NULL) {
     1042        for (int j = 0; j < params->n; j++) {
     1043            if (paramMask->data.U8[j]) {
     1044                alpha->data.F64[j][j] = 1;
     1045                beta->data.F64[j] = 1;
     1046            }
     1047        }
     1048    }
     1049
    9281050    psFree (deriv);
    9291051    return (chisq);
    … …  
    9341056                           psVector *Beta,
    9351057                           psVector *Params,
    936                            psImage  *alpha,
    937                            psVector *beta,
    938                            psVector *params,
     1058                           const psImage  *alpha,
     1059                           const psVector *beta,
     1060                           const psVector *params,
     1061                           const psVector *paramMask,
    9391062                           psF64 lambda)
    9401063{
    … …  
    9481071
    9491072    // LU decomposition version
    950     psTrace (".pslib.dataManip.psMinLM_GuessABP", 3, "using LUD version");
     1073    psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using LUD version\n");
    9511074
    9521075    // set new guess values (creates matrix A)
    9531076    A = psImageCopy (NULL, alpha, PS_TYPE_F64);
    9541077    for (int j = 0; j < params->n; j++) {
     1078        if ((paramMask != NULL) && (paramMask->data.U8[j]))
     1079            continue;
    9551080        A->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
    9561081    }
    … …  
    9641089    # else
    9651090        // gauss-jordan version
    966         psTrace (".pslib.dataManip.psMinLM_GuessABP", 3, "using Gauss-J version");
     1091        psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using Gauss-J version");
    9671092
    9681093    // set new guess values (creates matrix A)
    … …  
    9701095    Alpha = psImageCopy (Alpha, alpha, PS_TYPE_F64);
    9711096    for (int j = 0; j < params->n; j++) {
     1097        if ((paramMask != NULL) && (paramMask->data.U8[j]))
     1098            continue;
    9721099        Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
    9731100    }
    … …  
    9761103    # endif
    9771104
    978     // apply beta to get new params values
     1105    // apply Beta to get new Params values
    9791106    for (int j = 0; j < params->n; j++) {
     1107        if ((paramMask != NULL) && (paramMask->data.U8[j]))
     1108            continue;
    9801109        Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j];
    9811110    }
    … …  
    10951224    psFree (indxc);
    10961225    return (false);
    1097 }
    1098 
    1099 /******************************************************************************
    1100 psMinimizeLMChi2():  This routine will take an procedure which calculates
    1101 an arbitrary function and it's derivative and minimize the chi-squared match
    1102 between that function at the specified coords and the specified value at
    1103 those coords.
    1104  
    1105 XXX: Do this:
    1106  After checking that all entries in the paramMask are 1 or 0, when
    1107  forming the A matrix from alpha, try this:
    1108  
    1109      A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
    1110  
    1111 XXX: This is very different from what is specified in the SDR.  Must
    1112 coordinate with IfA on new SDR.
    1113  
    1114 XXX: Do vector/image recycles.
    1115  
    1116 XXX: probably yErr will be part of the SDR.
    1117  
    1118 XXX: This must work for both F32 and F64.  F32 is currently implemented.
    1119      Note: since the LUD routines are only implemented in F64, then we
    1120      will have to convert all F32 input vectors to F64 regardless.  So,
    1121      the F64 port might be.
    1122  
    1123 XXX: Must update the covar matrix.
    1124  *****************************************************************************/
    1125 psBool psMinimizeLMChi2Old(psMinimization *min,
    1126                            psImage *covar,
    1127                            psVector *params,
    1128                            const psVector *paramMask,
    1129                            const psArray *x,
    1130                            const psVector *y,
    1131                            const psVector *yErr,
    1132                            psMinimizeLMChi2Func func)
    1133 {
    1134     PS_ASSERT_PTR_NON_NULL(min, NULL);
    1135     PS_ASSERT_VECTOR_NON_NULL(params, NULL);
    1136     PS_ASSERT_VECTOR_NON_EMPTY(params, NULL);
    1137     PS_ASSERT_PTR_NON_NULL(x, NULL);
    1138     PS_ASSERT_VECTOR_NON_NULL(y, NULL);
    1139     PS_ASSERT_VECTOR_NON_EMPTY(y, NULL);
    1140     PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, NULL);
    1141     PS_ASSERT_PTR_NON_NULL(func, NULL);
    1142 
    1143     if (paramMask != NULL) {
    1144         PS_ASSERT_VECTORS_SIZE_EQUAL(params, paramMask, NULL);
    1145     }
    1146     if (yErr != NULL) {
    1147         PS_ASSERT_VECTORS_SIZE_EQUAL(y, yErr, NULL);
    1148     }
    1149     if (covar != NULL) {
    1150         PS_ASSERT_IMAGE_SIZE(covar, params->n, params->n, NULL);
    1151     }
    1152 
    1153     psTrace(".psLib.dataManip.psMinimize", 4,
    1154             "---- psMinimizeLMChi2() begin ----\n");
    1155     psS32 numData = y->n;
    1156     psS32 numParams = params->n;
    1157     psS32 i;
    1158     psS32 j;
    1159     psS32 k;
    1160     psS32 l;
    1161     psS32 n;
    1162     psS32 p;
    1163     psVector *beta = psVectorAlloc(numParams, PS_TYPE_F64);
    1164     psVector *perm = NULL;
    1165 
    1166     psVector *paramDeltasF64 = psVectorAlloc(numParams, PS_TYPE_F64);
    1167     psVector *origParams = psVectorAlloc(numParams, PS_TYPE_F32);
    1168     psVector *newParams = psVectorAlloc(numParams, PS_TYPE_F32);
    1169 
    1170     psImage *alpha = psImageAlloc(numParams, numParams, PS_TYPE_F32);
    1171     psImage *A = psImageAlloc(numParams, numParams, PS_TYPE_F64);
    1172     psImage *aOut = psImageAlloc(numParams, numParams, PS_TYPE_F64);
    1173     psImage *deriv = psImageAlloc(numParams, numData, PS_TYPE_F32);
    1174     psVector *currValueVec = NULL;
    1175     psVector *newValueVec = NULL;
    1176     psF32 currChi2 = 0.0;
    1177     psF32 newChi2 = 0.0;
    1178     psF32 lamda = 0.00005;  // XXX EAM : this starting value is VERY small (lamda is mis-spelt)
    1179     lamda = 0.05;  // XXX EAM : this starting value is quite large (lamda is mis-spelt)
    1180 
    1181     psTrace(".psLib.dataManip.psMinimize", 6,
    1182             "min->maxIter is %d\n", min->maxIter);
    1183     psTrace(".psLib.dataManip.psMinimize", 6,
    1184             "min->tol is %f\n", min->tol);
    1185 
    1186     for (p=0;p<numParams;p++) {
    1187         origParams->data.F32[p] = params->data.F32[p];
    1188     }
    1189 
    1190     min->lastDelta = PS_MAX_F32;
    1191     min->iter = 0;
    1192 
    1193     while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {
    1194         psTrace(".psLib.dataManip.psMinimize", 4,
    1195                 "------------------------------------------------------\n");
    1196         psTrace(".psLib.dataManip.psMinimize", 4,
    1197                 "Iteration %d.  Delta is %f\n", min->iter, min->lastDelta);
    1198 
    1199         //
    1200         // Calculate the current values and chi-squared of the function.
    1201         //
    1202         currChi2 = 0.0;
    1203         // currValueVec = func(deriv, params, x);
    1204 
    1205         // XXX EAM: use BinaryOp ?
    1206         // t1 = BinaryOp (NULL, currValueVec, "-", y);
    1207         // t1 = BinaryOp (t1, t1, "*", t1);
    1208 
    1209         // XXX EAM: this ignores yErr
    1210         for (n=0;n<numData;n++) {
    1211             currChi2+= (currValueVec->data.F32[n] - y->data.F32[n]) *
    1212                        (currValueVec->data.F32[n] - y->data.F32[n]);
    1213             psTrace(".psLib.dataManip.psMinimize", 6,
    1214                     "data[%d], chi2 calculation+= (%f - %f)^2\n", n,
    1215                     currValueVec->data.F32[n], y->data.F32[n]);
    1216         }
    1217 
    1218         // XXX EAM: this is just for tracing
    1219         for (p=0;p<numParams;p++) {
    1220             psTrace(".psLib.dataManip.psMinimize", 6,
    1221                     "params->data.F32[%d] is %f.\n", p, params->data.F32[p]);
    1222         }
    1223         psTrace(".psLib.dataManip.psMinimize", 6,
    1224                 "Current chi-squared is (%f)\n", currChi2);
    1225 
    1226         //
    1227         // Mask elements of the derivative for each data point.
    1228         // XXX EAM : is this necessary?  probably not...
    1229         for (p=0;p<numParams;p++) {
    1230             if ((paramMask != NULL) && (paramMask->data.U8[p] != 0)) {
    1231                 for (n=0;n<numData;n++) {
    1232                     deriv->data.F32[n][p] = 0.0;
    1233                 }
    1234             }
    1235         }
    1236 
    1237         //
    1238         // Calculate the BETA vector.
    1239         // XXX EAM: I think this is wrong
    1240         for (p=0;p<numParams;p++) {
    1241             if ((paramMask != NULL) && (paramMask->data.U8[p] != 0)) {
    1242                 continue;
    1243             }
    1244             beta->data.F64[p] = 0.0;
    1245             for (n=0;n<numData;n++) {
    1246                 (beta->data.F64[p])+=
    1247                     (y->data.F32[n] - currValueVec->data.F32[n]) *
    1248                     deriv->data.F32[n][p];
    1249             }
    1250             // XXX: multiply by -1 here?
    1251             (beta->data.F64[p])*= -1.0;
    1252             psTrace(".psLib.dataManip.psMinimize", 6,
    1253                     "beta->data.F64[%d] is %f.\n", p, beta->data.F64[p]);
    1254         }
    1255         psFree(currValueVec);
    1256 
    1257         //
    1258         // Calculate the ALPHA matrix.
    1259         // XXX EAM: also wrong? (missing yErr)
    1260         for (k=0;k<numParams;k++) {
    1261             for (l=0;l<numParams;l++) {
    1262                 alpha->data.F32[k][l] = 0.0;
    1263                 for (n=0;n<numData;n++) {
    1264                     alpha->data.F32[k][l]+= deriv->data.F32[n][k] *
    1265                                             deriv->data.F32[n][l];
    1266                 }
    1267             }
    1268         }
    1269 
    1270         //
    1271         // Calculate the matrix A.
    1272         //
    1273         for (j=0;j<numParams;j++) {
    1274             for (k=0;k<numParams;k++) {
    1275                 if (j == k) {
    1276                     A->data.F64[j][k] =
    1277                         (psF64) ((1.0 + lamda) * alpha->data.F32[j][k]);
    1278                 } else {
    1279                     A->data.F64[j][k] = (psF64) alpha->data.F32[j][k];
    1280                 }
    1281             }
    1282         }
    1283         for (j=0;j<numParams;j++) {
    1284             psTrace(".psLib.dataManip.psMinimize", 6, "Matrix A[][]:\n");
    1285             for (k=0;k<numParams;k++) {
    1286                 psTrace(".psLib.dataManip.psMinimize", 6, "%f ", A->data.F64[j][k]);
    1287             }
    1288             psTrace(".psLib.dataManip.psMinimize", 6, "Matrix A[][]:\n");
    1289         }
    1290 
    1291         //
    1292         // Solve A * alpha = Beta
    1293         //
    1294         // XXX: How do we know if these functions were successful?
    1295         //
    1296         aOut = psMatrixLUD(aOut, &perm, A);
    1297         paramDeltasF64 = psMatrixLUSolve(paramDeltasF64, aOut, beta, perm);
    1298 
    1299         //
    1300         // Mask any masked parameters.
    1301         //
    1302         for (i=0;i<numParams;i++) {
    1303             psTrace(".psLib.dataManip.psMinimize", 6,
    1304                     "paramDeltasF64->data.F64[%d] is %f.\n", i, paramDeltasF64->data.F64[i]);
    1305             if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
    1306                 newParams->data.F32[i] = origParams->data.F32[i];
    1307             } else {
    1308                 newParams->data.F32[i] = params->data.F32[i] -
    1309                                          (psF32) paramDeltasF64->data.F64[i];
    1310             }
    1311         }
    1312 
    1313         psTrace(".psLib.dataManip.psMinimize", 6,
    1314                 "Calling func() with new parameters:\n");
    1315         for (i=0;i<numParams;i++) {
    1316             psTrace(".psLib.dataManip.psMinimize", 6,
    1317                     "newParams->data.F32[%d] is %f.\n", i, newParams->data.F32[i]);
    1318         }
    1319 
    1320 
    1321         //
    1322         // Calculate new function values.
    1323         //
    1324         newChi2 = 0.0;
    1325         // newValueVec = func(deriv, newParams, x);
    1326         for (n=0;n<numData;n++) {
    1327             newChi2+= (newValueVec->data.F32[n] - y->data.F32[n]) *
    1328                       (newValueVec->data.F32[n] - y->data.F32[n]);
    1329 
    1330         }
    1331         psFree(newValueVec);
    1332 
    1333         psTrace(".psLib.dataManip.psMinimize", 4,
    1334                 "old/new chi-squareds are (%f, %f)\n", currChi2, newChi2);
    1335 
    1336         //
    1337         // If the new chi-squared is lower, then keep it.
    1338         //
    1339         if (currChi2 > newChi2) {
    1340             min->lastDelta = (currChi2 - newChi2)/currChi2;
    1341             min->value = newChi2;
    1342 
    1343             // We already masked params.
    1344             for (i=0;i<numParams;i++) {
    1345                 params->data.F32[i] = (psF32) newParams->data.F32[i];
    1346             }
    1347             lamda*= 0.1;
    1348             psTrace(".psLib.dataManip.psMinimize", 4, "*** Reducing lamda by factor of 10\n");
    1349         } else {
    1350             lamda*= 10.0;
    1351             psTrace(".psLib.dataManip.psMinimize", 4, "*** Increasing lamda by factor of 10\n");
    1352         }
    1353         psTrace(".psLib.dataManip.psMinimize", 4,
    1354                 "lamda is %f\n", lamda);
    1355         min->iter++;
    1356     }
    1357     psFree(beta);
    1358     psFree(perm);
    1359     psFree(paramDeltasF64);
    1360     psFree(origParams);
    1361     psFree(newParams);
    1362     psFree(alpha);
    1363     psFree(A);
    1364     psFree(aOut);
    1365     psFree(deriv);
    1366 
    1367     if ((min->iter < min->maxIter) ||
    1368             (min->lastDelta <= min->tol)) {
    1369         return(true);
    1370     }
    1371 
    1372     psTrace(".psLib.dataManip.psMinimize", 4,
    1373             "---- psMinimizeLMChi2() end (false) ----\n");
    1374     return(false);
    13751226}
    13761227
Note: See TracChangeset for help on using the changeset viewer.