IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Sep 27, 2007, 2:35:20 PM (19 years ago)
Author:
eugene
Message:

cleanups; Alpha, Beta now only include unmasked parameters

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r14320 r15046  
    1010 *  @author EAM, IfA
    1111 *
    12  *  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2007-07-20 00:22:24 $
     12 *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2007-09-28 00:35:17 $
    1414 *
    1515 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
    … …  
    5454/*****************************************************************************/
    5555
     56// Alpha & Beta only represent unmasked values
    5657bool psMinLM_GuessABP(
    5758    psImage  *Alpha,
    … …  
    7677    }
    7778
    78     // XXX the LU Decomposition code requires F64.  this is now incompatible
    79     // with the rest of this code. drop this segment when the changes are tested.
    80     # define USE_LU_DECOMP 0
    81     # if (USE_LU_DECOMP)
    82         psVector *LUv = NULL;
    83     psImage  *LUm = NULL;
    84     psImage  *A   = NULL;
    85     psF32    det;
    86     // LU decomposition version
    87     psTrace("psLib.math", 5, "using LUD version\n");
    88 
    89     // set new guess values (creates matrix A)
    90     A = psImageCopy(NULL, alpha, PS_TYPE_F64);
    91     for (int j = 0; j < params->n; j++) {
    92         if ((paramMask != NULL) && (paramMask->data.U8[j])) {
    93             continue;
    94         }
    95         A->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
    96     }
    97 
    98     // solve A*beta = Beta (Alpha = 1/A)
    99     // these operations do not modify the input values (creates LUm, LUv)
    100     LUm   = psMatrixLUD(NULL, &LUv, A);
    101     if (LUm == NULL) {
    102         psError(PS_ERR_UNKNOWN, false, "psMatrixLUD() returned NULL\n");
    103     }
    104     Beta  = psMatrixLUSolve(Beta, LUm, beta, LUv);
    105     if (Beta == NULL) {
    106         psError(PS_ERR_UNKNOWN, false, "psMatrixLUSolve() returned NULL\n");
    107     }
    108     Alpha = psMatrixInvert(Alpha, A, &det);
    109     psFree(LUm);
    110     psFree(LUv);
    111     psFree(A);
    112 
    113     if (Alpha == NULL) {
    114         psError(PS_ERR_UNKNOWN, false, "psMatrixInvert() returned NULL\n");
    115         return(false);
    116     }
    117     # else
    118         // gauss-jordan version
    119         psTrace("psLib.math", 5, "using Gauss-J version");
    120 
    121     // set new guess values (creates matrix A)
     79    assert (alpha->numCols == beta->n);
     80    assert (alpha->numCols == alpha->numRows);
     81
     82    // set new guess values, applying (1+lambda) scaling to pivots
    12283    Beta = psVectorCopy(Beta, beta, PS_TYPE_F32);
    12384    Alpha = psImageCopy(Alpha, alpha, PS_TYPE_F32);
    124     for (int j = 0; j < params->n; j++) {
    125         if ((paramMask != NULL) && (paramMask->data.U8[j])) {
    126             continue;
    127         }
     85    for (int j = 0; j < Alpha->numCols; j++) {
    12886        Alpha->data.F32[j][j] = alpha->data.F32[j][j] * (1.0 + lambda);
    12987    }
    13088
     89    // error and clear above if kept?
    13190    if (false == psMatrixGJSolveF32(Alpha, Beta)) {
    132         // psError(PS_ERR_UNKNOWN, false, "singular matrix in Guess ABP\n");
    13391        psTrace ("psLib.math", 4, "singular matrix in Guess ABP\n");
    13492        return(false);
    13593    }
    136     # endif
    13794
    13895    // measure linear model prediction
    13996    // (we must do this before truncating Beta below)
    14097    if (dLinear) {
    141         *dLinear = psMinLM_dLinear(Beta, beta, paramMask, lambda);
     98        *dLinear = psMinLM_dLinear(Beta, beta, lambda);
     99    }
     100
     101    // full-length Beta for checkLimits functions
     102    psVector *tmpBeta = psVectorAlloc(params->n, PS_TYPE_F32);
     103    psVectorInit (tmpBeta, 0.0);
     104
     105    // set tmpBeta values which are not masked
     106    for (int j = 0, n = 0; j < params->n; j++) {
     107        if (paramMask && (paramMask->data.U8[j])) continue;
     108        tmpBeta->data.F32[j] = Beta->data.F32[n];
     109        n++;
    142110    }
    143111
    144112    // apply Beta to get new Params values
    145113    for (int j = 0; j < params->n; j++) {
    146         if ((paramMask != NULL) && (paramMask->data.U8[j])) {
     114        if (paramMask && (paramMask->data.U8[j])) {
    147115            Params->data.F32[j] = params->data.F32[j];
    148116            continue;
    149117        }
    150118        // apply beta limits
    151         if (checkLimits)
    152             checkLimits (PS_MINIMIZE_BETA_LIMIT, j, Params->data.F32, Beta->data.F32);
    153         Params->data.F32[j] = params->data.F32[j] - Beta->data.F32[j];
     119        if (checkLimits) {
     120            checkLimits (PS_MINIMIZE_BETA_LIMIT, j, Params->data.F32, tmpBeta->data.F32);
     121        }
     122
     123        Params->data.F32[j] = params->data.F32[j] - tmpBeta->data.F32[j];
    154124
    155125        // compare new params to param limits
    156         if (checkLimits)
    157             checkLimits (PS_MINIMIZE_PARAM_MIN,  j, Params->data.F32, Beta->data.F32);
    158         if (checkLimits)
    159             checkLimits (PS_MINIMIZE_PARAM_MAX,  j, Params->data.F32, Beta->data.F32);
    160     }
     126        if (checkLimits) {
     127            checkLimits (PS_MINIMIZE_PARAM_MIN,  j, Params->data.F32, tmpBeta->data.F32);
     128            checkLimits (PS_MINIMIZE_PARAM_MAX,  j, Params->data.F32, tmpBeta->data.F32);
     129        }
     130    }
     131
     132    // apply tmpBeta after limits have been checked
     133    for (int j = 0, n = 0; j < params->n; j++) {
     134        if (paramMask && (paramMask->data.U8[j])) continue;
     135        Beta->data.F32[n] = tmpBeta->data.F32[j];
     136        n++;
     137    }
     138
     139    psFree (tmpBeta);
    161140    return(true);
    162141}
    … …  
    172151{
    173152    psTrace("psLib.math", 3, "---- begin ----\n");
     153
    174154    // allocate internal arrays (current vs Guess)
    175     psImage  *alpha  = psImageAlloc (params->n, params->n, PS_TYPE_F32);
    176     psImage  *Alpha  = psImageAlloc (params->n, params->n, PS_TYPE_F32);
    177     psVector *beta   = psVectorAlloc(params->n, PS_TYPE_F32);
     155    psImage *Alpha = NULL;
     156    psVector *Beta = NULL;
     157
     158    psVectorInit (delta, 0.0);
     159
     160    // Alpha & Beta only contain elements to represent the unmasked parameters
     161    // if none are available, return false
     162    if (!psMinLM_AllocAB (&Alpha, &Beta, params, paramMask)) {
     163        return false;
     164    }
     165
     166    psImage *alpha   = psImageAlloc(Alpha->numCols, Alpha->numRows, PS_TYPE_F32);
    178167    psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32);
     168
    179169    psVector *dy     = NULL;
    180     bool rc = true;
     170    bool retValue = true;
    181171
    182172    // the user provides the error or NULL.  we need to convert
    … …  
    190180
    191181    // XXX should we give up if chisq is nan?
    192     psF32 rcF32 = psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func);
    193     if (isnan(rcF32)) {
    194         psTrace ("psLib.math", 5, "psMinLM_SetABX() returned a NAN.\n");
    195         rc = false;
    196     }
     182    psF32 chisq = psMinLM_SetABX(alpha, Beta, params, paramMask, x, y, dy, func);
     183    if (isnan(chisq)) {
     184        psTrace ("psLib.math", 5, "psMinLM_SetABX() returned a NAN chisq.\n");
     185        retValue = false;
     186    }
     187
    197188    psTrace("psLib.math", 5, "psMinLM_SetABX() was succesful\n");
    198189    // dump some useful info if trace is defined
    199190    if (psTraceGetLevel("psLib.math") >= 6) {
    200191        p_psImagePrint(psTraceGetDestination(), alpha, "alpha guess (0)");
    201         p_psVectorPrint(psTraceGetDestination(), beta, "beta guess (0)");
     192        p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (0)");
    202193        p_psVectorPrint(psTraceGetDestination(), params, "params guess (0)");
    203194    }
    204195
    205     bool rcBool = psMinLM_GuessABP(Alpha, delta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL);
    206     if (rcBool == false) {
     196    bool status = psMinLM_GuessABP(Alpha, delta, Params, alpha, Beta, params, paramMask, NULL, 0.0, NULL);
     197    if (!status) {
    207198        psTrace ("psLib.math", 5, "psMinLM_GuessABP() returned FALSE.\n");
    208         rc = false;
     199        retValue = false;
    209200    }
    210201    psTrace("psLib.math", 5, "psMinLM_GuessABP() was succesful\n");
    … …  
    217208    psFree(alpha);
    218209    psFree(Alpha);
    219     psFree(beta);
     210    psFree(Beta);
    220211    psFree(Params);
    221212    if (yWt == NULL) {
    … …  
    223214    }
    224215    psTrace("psLib.math", 3, "---- end ----\n");
    225     return(rc);
     216    return(retValue);
    226217}
    227218
    … …  
    230221    const psVector *Beta,
    231222    const psVector *beta,
    232     const psVector *paramMask,
    233223    psF32 lambda)
    234224{
    … …  
    240230
    241231    float dh = 0.0, sh = 0.0, Sh = 0.0;
     232
     233    // beta only counts unmasked parameters
    242234    for (int i = 0; i < beta->n; i++) {
    243         if ((paramMask != NULL) && (paramMask->data.U8[i])) continue;
    244235        dh = lambda*B[i] + b[i];
    245236        sh = 0.5*B[i]*dh;
    246237        Sh += sh;
    247         // fprintf (stderr, "%f, %f  * %f %f %f * :  %f : %f + %f\n", B[i], b[i], dh, sh, Sh, lambda*PS_SQR(B[i]) + B[i]*b[i], lambda*PS_SQR(B[i]), B[i]*b[i]);
    248238        dLinear += lambda*PS_SQR(B[i]) + B[i]*b[i];
    249239    }
    … …  
    262252    psMinimizeLMChi2Func func)
    263253{
    264     // XXX: Check vector sizes.
    265254    PS_ASSERT_IMAGE_NON_NULL(alpha, NAN);
    266255    PS_ASSERT_VECTOR_NON_NULL(beta, NAN);
    … …  
    281270    psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32);
    282271
    283     // zero alpha and beta for summing below
    284     for (psS32 j = 0; j < params->n; j++) {
    285         for (psS32 k = 0; k < params->n; k++) {
    286             alpha->data.F32[j][k] = 0;
    287         }
    288         beta->data.F32[j] = 0;
    289     }
     272    // zero alpha, beta, and chisq for summing below
     273    psImageInit (alpha, 0.0);
     274    psVectorInit (beta, 0.0);
    290275    chisq = 0.0;
    291276
    292     // calculate chisq, alpha, beta
     277    // calculate chisq, alpha, beta. alpha & beta only represent unmasked parameters; skip
     278    // masked ones
    293279    for (psS32 i = 0; i < y->n; i++) {
    294280        ymodel = func(deriv, params, (psVector *) x->data[i]);
    … …  
    296282        delta = ymodel - y->data.F32[i];
    297283        chisq += PS_SQR(delta) * dy->data.F32[i];
    298         if (isnan(dy->data.F32[i]))
    299             psAbort("nan in weights");
    300         if (isnan(delta))
    301             psAbort("nan in delta");
    302         if (isnan(chisq))
    303             psAbort("nan in chisq");
    304 
    305         for (psS32 j = 0; j < params->n; j++) {
    306             if ((paramMask != NULL) && (paramMask->data.U8[j])) {
    307                 continue;
     284        assert (!isnan(dy->data.F32[i]));
     285        assert (!isnan(delta));
     286        assert (!isnan(chisq));
     287
     288        // we track alpha,beta and params,deriv separately
     289        for (int j = 0, J = 0; j < params->n; j++) {
     290            if (paramMask && (paramMask->data.U8[j])) continue;
     291
     292            weight = deriv->data.F32[j] * dy->data.F32[i];
     293
     294            for (int k = 0, K = 0; k <= j; k++) {
     295                if (paramMask && (paramMask->data.U8[k])) continue;
     296                alpha->data.F32[J][K] += weight * deriv->data.F32[k];
     297                K++;
    308298            }
    309             weight = deriv->data.F32[j] * dy->data.F32[i];
    310             for (psS32 k = 0; k <= j; k++) {
    311                 if ((paramMask != NULL) && (paramMask->data.U8[k])) {
    312                     continue;
    313                 }
    314                 alpha->data.F32[j][k] += weight * deriv->data.F32[k];
    315             }
    316             beta->data.F32[j] += weight * delta;
     299            beta->data.F32[J] += weight * delta;
     300            J++;
    317301        }
    318302    }
    319303
    320304    // calculate lower-left half of alpha
    321     for (psS32 j = 1; j < params->n; j++) {
    322         for (psS32 k = 0; k < j; k++) {
     305    for (int j = 1; j < alpha->numCols; j++) {
     306        for (int k = 0; k < j; k++) {
    323307            alpha->data.F32[k][j] = alpha->data.F32[j][k];
    324         }
    325     }
    326 
    327     // fill in pivots if we apply a mask
    328     if (paramMask != NULL) {
    329         for (psS32 j = 0; j < params->n; j++) {
    330             if (paramMask->data.U8[j]) {
    331                 alpha->data.F32[j][j] = 1;
    332                 beta->data.F32[j] = 1;
    333             }
    334308        }
    335309    }
    … …  
    398372
    399373    // allocate internal arrays (current vs Guess)
    400     psImage *alpha   = psImageAlloc(params->n, params->n, PS_TYPE_F32);
    401     psImage *Alpha   = psImageAlloc(params->n, params->n, PS_TYPE_F32);
    402     psVector *beta   = psVectorAlloc(params->n, PS_TYPE_F32);
    403     psVector *Beta   = psVectorAlloc(params->n, PS_TYPE_F32);
     374    psImage *Alpha = NULL;
     375    psVector *Beta = NULL;
     376
     377    // Alpha & Beta only contain elements to represent the unmasked parameters
     378    if (!psMinLM_AllocAB (&Alpha, &Beta, params, paramMask)) {
     379        psAbort ("programming error: no unmasked parameters to be fit\n");
     380    }
     381   
     382    psImage *alpha   = psImageAlloc(Alpha->numCols, Alpha->numRows, PS_TYPE_F32);
     383    psVector *beta   = psVectorAlloc(Beta->n, PS_TYPE_F32);
    404384    psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32);
     385
    405386    psVector *dy     = NULL;
    406387    psF32 Chisq = 0.0;
    … …  
    494475    // construct & return the covariance matrix (if requested)
    495476    if (covar != NULL) {
    496         if (!psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) {
     477        if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) {
    497478            psTrace ("psLib.math", 5, "failure to calculate covariance matrix\n");
    498479        }
     480        // set covar values which are not masked
     481        psImageInit (covar, 0.0);
     482        for (int j = 0, J = 0; j < params->n; j++) {
     483            if (paramMask && (paramMask->data.U8[j])) {
     484                covar->data.F32[j][j] = 1.0;
     485                continue;
     486            }
     487            for (int k = 0, K = 0; k < params->n; k++) {
     488                if (paramMask && (paramMask->data.U8[k])) continue;
     489                covar->data.F32[j][k] = Alpha->data.F32[J][K];
     490                K++;
     491            }
     492            J++;
     493        }
    499494    }
    500495
    … …  
    516511}
    517512
     513bool psMinLM_AllocAB (psImage **Alpha, psVector **Beta, const psVector *params, const psVector *paramMask) {
     514
     515    assert (Alpha);
     516    assert (Beta);
     517    assert (params);
     518
     519    int nParams = params->n;
     520
     521    // count unmasked parameters
     522    if (paramMask) {
     523        nParams = 0;
     524        for (int i = 0; i < paramMask->n; i++) {
     525            if (paramMask->data.U8[i]) continue;
     526            nParams ++;
     527        }
     528    }
     529
     530    if (nParams == 0) {
     531        return false;
     532    }
     533
     534    *Alpha = psImageAlloc(nParams, nParams, PS_TYPE_F32);
     535    *Beta  = psVectorAlloc(nParams, PS_TYPE_F32);
     536    return true;
     537}
     538
    518539static void minimizationFree(psMinimization *min)
    519540{
Note: See TracChangeset for help on using the changeset viewer.