IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1819


Ignore:
Timestamp:
Sep 16, 2004, 9:30:43 AM (22 years ago)
Author:
gusciora
Message:

...

Location:
trunk/psLib
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r1798 r1819  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-09-11 03:03:36 $
     11 *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-09-16 19:30:24 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    865865                  psMinimizeLMFunc func)
    866866{
    867     psVector *beta = psVectorAlloc(params->n, PS_TYPE_F32);
    868     psVector *perm = psVectorAlloc(params->n, PS_TYPE_F32);
    869     psVector *newParams = psVectorAlloc(params->n, PS_TYPE_F32);
     867    psVector *beta = psVectorAlloc(params->n, PS_TYPE_F64);
     868    psVector *perm = psVectorAlloc(params->n, PS_TYPE_F64);
     869    psVector *newParamsF64 = psVectorAlloc(params->n, PS_TYPE_F64);
     870    psVector *newParamsF32 = psVectorAlloc(params->n, PS_TYPE_F32);
    870871    psVector *origParams = psVectorAlloc(params->n, PS_TYPE_F32);
    871872    psImage *alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32);
    872     psImage *A = psImageAlloc(params->n, params->n, PS_TYPE_F32);
    873     psImage *aOut = psImageAlloc(params->n, params->n, PS_TYPE_F32);
     873    psImage *A = psImageAlloc(params->n, params->n, PS_TYPE_F64);
     874    psImage *aOut = psImageAlloc(params->n, params->n, PS_TYPE_F64);
    874875    psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32);
    875876    psVector *newDeriv = psVectorAlloc(params->n, PS_TYPE_F32);
     877    psVector *NRparams = psVectorAlloc(params->n, PS_TYPE_F32);
    876878    int i;
    877879    int j;
    … …  
    879881    float newValue;
    880882    float oldValue;
    881     float lamda = 1.0;
     883    float lamda = 0.00005;
     884
     885    psTrace(".psLib.dataManip.psMinimizeLM", 4,
     886            "---- psMinimizeLM() begin ----\n");
     887    psTrace(".psLib.dataManip.psMinimizeLM", 6,
     888            "min->maxIter is %d\n", min->maxIter);
     889    psTrace(".psLib.dataManip.psMinimizeLM", 6,
     890            "min->tol is %f\n", min->tol);
    882891
    883892    for (i=0;i<params->n;i++) {
    … …  
    885894    }
    886895
    887     min->lastDelta = -HUGE;
     896    min->lastDelta = HUGE;
     897    min->lastDelta = 12345.0;
    888898    min->iter = 0;
    889     min->value = func(deriv, params, coords);
    890     for (i=0;i<params->n;i++) {
    891         if ((paramMask != NULL) &&
    892                 (paramMask->data.U8[i] != 0)) {
    893             deriv->data.F32[i] = 0.0;
    894         }
    895     }
    896899
    897900    while ((min->lastDelta > min->tol) &&
    898901            (min->iter < min->maxIter)) {
     902        psTrace(".psLib.dataManip.psMinimizeLM", 4,
     903                "------------------------------------------------------\n");
     904        psTrace(".psLib.dataManip.psMinimizeLM", 4,
     905                "Iteration %d.  Delta is %f\n", min->iter, min->lastDelta);
     906
     907        min->value = func(deriv, params, coords);
    899908
    900909        for (i=0;i<params->n;i++) {
    901             beta->data.F32[i] = deriv->data.F32[i];
    902         }
    903 
     910            psTrace(".psLib.dataManip.psMinimizeLM", 6,
     911                    "params->data.F32[%d] is %f.\n", i, params->data.F32[i]);
     912        }
     913        for (i=0;i<params->n;i++) {
     914            psTrace(".psLib.dataManip.psMinimizeLM", 6,
     915                    "deriv->data.F32[%d] is %f.\n", i, deriv->data.F32[i]);
     916        }
     917        psTrace(".psLib.dataManip.psMinimizeLM", 6,
     918                "min->value is (%f)\n", min->value);
     919
     920        for (i=0;i<params->n;i++) {
     921            if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
     922                deriv->data.F32[i] = 0.0;
     923            }
     924        }
     925
     926        // Calculate the BETA vector.
     927        for (i=0;i<params->n;i++) {
     928            beta->data.F64[i] = (float) deriv->data.F32[i];
     929            psTrace(".psLib.dataManip.psMinimizeLM", 6,
     930                    "beta->data.F64[%d] is %f.\n", i, beta->data.F64[i]);
     931        }
     932
     933        // Calculate the ALPHA matrix.
    904934        for (j=0;j<params->n;j++) {
    905935            for (k=0;k<params->n;k++) {
    906                 alpha->data.F32[j][k] = deriv->data.F32[j] * deriv->data.F32[j];
    907             }
    908         }
     936                alpha->data.F32[j][k] = deriv->data.F32[k] *
     937                                        deriv->data.F32[j];
     938            }
     939        }
     940
     941        // Calculate the matrix A.
    909942        for (j=0;j<params->n;j++) {
    910943            for (k=0;k<params->n;k++) {
    911                 if (j != k) {
    912                     A->data.F32[j][k] = alpha->data.F32[j][k];
     944                if (j == k) {
     945                    A->data.F64[j][k] =
     946                        (double) ((1.0 + lamda) * alpha->data.F32[j][k]);
    913947                } else {
    914                     A->data.F32[j][k] = (1.0 + lamda) * alpha->data.F32[j][k];
     948                    A->data.F64[j][k] = (double) alpha->data.F32[j][k];
    915949                }
    916950            }
    917951        }
     952        for (j=0;j<params->n;j++) {
     953            psTrace(".psLib.dataManip.psMinimizeLM", 6, "Matrix A[][]:\n");
     954            for (k=0;k<params->n;k++) {
     955                psTrace(".psLib.dataManip.psMinimizeLM", 6, "%f ", A->data.F64[j][k]);
     956            }
     957            psTrace(".psLib.dataManip.psMinimizeLM", 6, "Matrix A[][]:\n");
     958        }
     959
     960        // Solve A * alpha = Beta
    918961        aOut = psMatrixLUD(aOut, perm, A);
    919         newParams = psMatrixLUSolve(newParams, aOut, beta, perm);
    920 
     962        newParamsF64 = psMatrixLUSolve(newParamsF64, aOut, beta, perm);
     963
     964        // Mask any masked parameters.
    921965        for (i=0;i<params->n;i++) {
    922             if ((paramMask != NULL) &&
    923                     (paramMask->data.U8[i] != 0)) {
    924                 newParams->data.F32[i] = origParams->data.F32[i];
    925             }
     966            if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
     967                newParamsF64->data.F64[i] = (double) origParams->data.F32[i];
     968            }
     969            newParamsF32->data.F32[i] = (float) newParamsF64->data.F64[i];
     970
     971            psTrace(".psLib.dataManip.psMinimizeLM", 6,
     972                    "newParamsF32->data.F32[%d] is %f.\n", i, newParamsF32->data.F32[i]);
     973            NRparams->data.F32[i] = params->data.F32[i] - newParamsF32->data.F32[i];
     974        }
     975
     976        psTrace(".psLib.dataManip.psMinimizeLM", 6,
     977                "Calling func() with new parameters:\n");
     978        for (i=0;i<params->n;i++) {
     979            psTrace(".psLib.dataManip.psMinimizeLM", 6,
     980                    "NRparams->data.F32[%d] is %f.\n", i, NRparams->data.F32[i]);
    926981        }
    927982
    928983        oldValue = min->value;
    929         newValue = func(newDeriv, newParams, coords);
     984        newValue = func(newDeriv, NRparams, coords);
     985        psTrace(".psLib.dataManip.psMinimizeLM", 4,
     986                "old/new values are (%f, %f)\n", oldValue, newValue);
    930987        for (i=0;i<params->n;i++) {
    931             if ((paramMask != NULL) &&
    932                     (paramMask->data.U8[i] != 0)) {
     988            if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
    933989                deriv->data.F32[i] = 0.0;
    934990            }
    … …  
    942998            // before the last function evaluation.
    943999            for (i=0;i<params->n;i++) {
    944                 params->data.F32[i] = newParams->data.F32[i];
    945                 deriv->data.F32[i] = newDeriv->data.F32[i];
     1000                //                params->data.F32[i] = (float) newParamsF64->data.F64[i];
     1001                params->data.F32[i] = (float) NRparams->data.F32[i];
    9461002            }
    9471003            min->value = func(deriv, params, coords);
    … …  
    9511007            lamda*= 10.0;
    9521008        }
     1009        psTrace(".psLib.dataManip.psMinimizeLM", 4,
     1010                "lamda is %f\n", lamda);
    9531011        min->iter++;
    9541012    }
    9551013    psFree(beta);
    9561014    psFree(perm);
    957     psFree(newParams);
     1015    psFree(newParamsF64);
     1016    psFree(newParamsF32);
    9581017    psFree(origParams);
    9591018    psFree(alpha);
    … …  
    9681027    }
    9691028
     1029    psTrace(".psLib.dataManip.psMinimizeLM", 4,
     1030            "---- psMinimizeLM() end (false) ----\n");
    9701031    return(false);
    9711032}
    … …  
    12991360        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
    13001361                "psDetermineBracket(): iteration %d\n", iter);
    1301         if ((b < a) && (b < c)) {
     1362        if ((fb < fa) && (fb < fc)) {
    13021363            bracket->data.F32[0] = a;
    13031364            bracket->data.F32[1] = b;
    … …  
    17381799
    17391800                psTrace(".psLib.dataManip.psMinimizePowell", 6,
    1740                         "WOW (%f %f %f)\n", u->data.F32[i], Q->data.F32[i], params->data.F32[i]);
     1801                        "u=Q-P (%f = %f - %f)\n", u->data.F32[i],
     1802                        Q->data.F32[i],
     1803                        params->data.F32[i]);
    17411804
    17421805            } else {
    … …  
    17471810        // 5: Move Q only in the direction u, and minimize the function.
    17481811        psTrace(".psLib.dataManip.psMinimizePowell", 6, "HERE 01\n");
    1749         psTrace(".psLib.dataManip.psMinimizePowell", 6,
    1750                 "u[0] is %f\n", u->data.F32[0]);
     1812        for (i=0;i<numDims;i++) {
     1813            psTrace(".psLib.dataManip.psMinimizePowell", 6,
     1814                    "u[i] is %f\n", u->data.F32[i]);
     1815        }
    17511816
    17521817        mul = psLineMin(min, params, u, paramMask, coords, func);

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

    r1798 r1819  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.38 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-09-11 03:03:36 $
     11 *  @version $Revision: 1.39 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-09-16 19:30:24 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    865865                  psMinimizeLMFunc func)
    866866{
    867     psVector *beta = psVectorAlloc(params->n, PS_TYPE_F32);
    868     psVector *perm = psVectorAlloc(params->n, PS_TYPE_F32);
    869     psVector *newParams = psVectorAlloc(params->n, PS_TYPE_F32);
     867    psVector *beta = psVectorAlloc(params->n, PS_TYPE_F64);
     868    psVector *perm = psVectorAlloc(params->n, PS_TYPE_F64);
     869    psVector *newParamsF64 = psVectorAlloc(params->n, PS_TYPE_F64);
     870    psVector *newParamsF32 = psVectorAlloc(params->n, PS_TYPE_F32);
    870871    psVector *origParams = psVectorAlloc(params->n, PS_TYPE_F32);
    871872    psImage *alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32);
    872     psImage *A = psImageAlloc(params->n, params->n, PS_TYPE_F32);
    873     psImage *aOut = psImageAlloc(params->n, params->n, PS_TYPE_F32);
     873    psImage *A = psImageAlloc(params->n, params->n, PS_TYPE_F64);
     874    psImage *aOut = psImageAlloc(params->n, params->n, PS_TYPE_F64);
    874875    psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32);
    875876    psVector *newDeriv = psVectorAlloc(params->n, PS_TYPE_F32);
     877    psVector *NRparams = psVectorAlloc(params->n, PS_TYPE_F32);
    876878    int i;
    877879    int j;
    … …  
    879881    float newValue;
    880882    float oldValue;
    881     float lamda = 1.0;
     883    float lamda = 0.00005;
     884
     885    psTrace(".psLib.dataManip.psMinimizeLM", 4,
     886            "---- psMinimizeLM() begin ----\n");
     887    psTrace(".psLib.dataManip.psMinimizeLM", 6,
     888            "min->maxIter is %d\n", min->maxIter);
     889    psTrace(".psLib.dataManip.psMinimizeLM", 6,
     890            "min->tol is %f\n", min->tol);
    882891
    883892    for (i=0;i<params->n;i++) {
    … …  
    885894    }
    886895
    887     min->lastDelta = -HUGE;
     896    min->lastDelta = HUGE;
     897    min->lastDelta = 12345.0;
    888898    min->iter = 0;
    889     min->value = func(deriv, params, coords);
    890     for (i=0;i<params->n;i++) {
    891         if ((paramMask != NULL) &&
    892                 (paramMask->data.U8[i] != 0)) {
    893             deriv->data.F32[i] = 0.0;
    894         }
    895     }
    896899
    897900    while ((min->lastDelta > min->tol) &&
    898901            (min->iter < min->maxIter)) {
     902        psTrace(".psLib.dataManip.psMinimizeLM", 4,
     903                "------------------------------------------------------\n");
     904        psTrace(".psLib.dataManip.psMinimizeLM", 4,
     905                "Iteration %d.  Delta is %f\n", min->iter, min->lastDelta);
     906
     907        min->value = func(deriv, params, coords);
    899908
    900909        for (i=0;i<params->n;i++) {
    901             beta->data.F32[i] = deriv->data.F32[i];
    902         }
    903 
     910            psTrace(".psLib.dataManip.psMinimizeLM", 6,
     911                    "params->data.F32[%d] is %f.\n", i, params->data.F32[i]);
     912        }
     913        for (i=0;i<params->n;i++) {
     914            psTrace(".psLib.dataManip.psMinimizeLM", 6,
     915                    "deriv->data.F32[%d] is %f.\n", i, deriv->data.F32[i]);
     916        }
     917        psTrace(".psLib.dataManip.psMinimizeLM", 6,
     918                "min->value is (%f)\n", min->value);
     919
     920        for (i=0;i<params->n;i++) {
     921            if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
     922                deriv->data.F32[i] = 0.0;
     923            }
     924        }
     925
     926        // Calculate the BETA vector.
     927        for (i=0;i<params->n;i++) {
     928            beta->data.F64[i] = (float) deriv->data.F32[i];
     929            psTrace(".psLib.dataManip.psMinimizeLM", 6,
     930                    "beta->data.F64[%d] is %f.\n", i, beta->data.F64[i]);
     931        }
     932
     933        // Calculate the ALPHA matrix.
    904934        for (j=0;j<params->n;j++) {
    905935            for (k=0;k<params->n;k++) {
    906                 alpha->data.F32[j][k] = deriv->data.F32[j] * deriv->data.F32[j];
    907             }
    908         }
     936                alpha->data.F32[j][k] = deriv->data.F32[k] *
     937                                        deriv->data.F32[j];
     938            }
     939        }
     940
     941        // Calculate the matrix A.
    909942        for (j=0;j<params->n;j++) {
    910943            for (k=0;k<params->n;k++) {
    911                 if (j != k) {
    912                     A->data.F32[j][k] = alpha->data.F32[j][k];
     944                if (j == k) {
     945                    A->data.F64[j][k] =
     946                        (double) ((1.0 + lamda) * alpha->data.F32[j][k]);
    913947                } else {
    914                     A->data.F32[j][k] = (1.0 + lamda) * alpha->data.F32[j][k];
     948                    A->data.F64[j][k] = (double) alpha->data.F32[j][k];
    915949                }
    916950            }
    917951        }
     952        for (j=0;j<params->n;j++) {
     953            psTrace(".psLib.dataManip.psMinimizeLM", 6, "Matrix A[][]:\n");
     954            for (k=0;k<params->n;k++) {
     955                psTrace(".psLib.dataManip.psMinimizeLM", 6, "%f ", A->data.F64[j][k]);
     956            }
     957            psTrace(".psLib.dataManip.psMinimizeLM", 6, "Matrix A[][]:\n");
     958        }
     959
     960        // Solve A * alpha = Beta
    918961        aOut = psMatrixLUD(aOut, perm, A);
    919         newParams = psMatrixLUSolve(newParams, aOut, beta, perm);
    920 
     962        newParamsF64 = psMatrixLUSolve(newParamsF64, aOut, beta, perm);
     963
     964        // Mask any masked parameters.
    921965        for (i=0;i<params->n;i++) {
    922             if ((paramMask != NULL) &&
    923                     (paramMask->data.U8[i] != 0)) {
    924                 newParams->data.F32[i] = origParams->data.F32[i];
    925             }
     966            if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
     967                newParamsF64->data.F64[i] = (double) origParams->data.F32[i];
     968            }
     969            newParamsF32->data.F32[i] = (float) newParamsF64->data.F64[i];
     970
     971            psTrace(".psLib.dataManip.psMinimizeLM", 6,
     972                    "newParamsF32->data.F32[%d] is %f.\n", i, newParamsF32->data.F32[i]);
     973            NRparams->data.F32[i] = params->data.F32[i] - newParamsF32->data.F32[i];
     974        }
     975
     976        psTrace(".psLib.dataManip.psMinimizeLM", 6,
     977                "Calling func() with new parameters:\n");
     978        for (i=0;i<params->n;i++) {
     979            psTrace(".psLib.dataManip.psMinimizeLM", 6,
     980                    "NRparams->data.F32[%d] is %f.\n", i, NRparams->data.F32[i]);
    926981        }
    927982
    928983        oldValue = min->value;
    929         newValue = func(newDeriv, newParams, coords);
     984        newValue = func(newDeriv, NRparams, coords);
     985        psTrace(".psLib.dataManip.psMinimizeLM", 4,
     986                "old/new values are (%f, %f)\n", oldValue, newValue);
    930987        for (i=0;i<params->n;i++) {
    931             if ((paramMask != NULL) &&
    932                     (paramMask->data.U8[i] != 0)) {
     988            if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
    933989                deriv->data.F32[i] = 0.0;
    934990            }
    … …  
    942998            // before the last function evaluation.
    943999            for (i=0;i<params->n;i++) {
    944                 params->data.F32[i] = newParams->data.F32[i];
    945                 deriv->data.F32[i] = newDeriv->data.F32[i];
     1000                //                params->data.F32[i] = (float) newParamsF64->data.F64[i];
     1001                params->data.F32[i] = (float) NRparams->data.F32[i];
    9461002            }
    9471003            min->value = func(deriv, params, coords);
    … …  
    9511007            lamda*= 10.0;
    9521008        }
     1009        psTrace(".psLib.dataManip.psMinimizeLM", 4,
     1010                "lamda is %f\n", lamda);
    9531011        min->iter++;
    9541012    }
    9551013    psFree(beta);
    9561014    psFree(perm);
    957     psFree(newParams);
     1015    psFree(newParamsF64);
     1016    psFree(newParamsF32);
    9581017    psFree(origParams);
    9591018    psFree(alpha);
    … …  
    9681027    }
    9691028
     1029    psTrace(".psLib.dataManip.psMinimizeLM", 4,
     1030            "---- psMinimizeLM() end (false) ----\n");
    9701031    return(false);
    9711032}
    … …  
    12991360        psTrace(".psLib.dataManip.p_psDetermineBracket", 6,
    13001361                "psDetermineBracket(): iteration %d\n", iter);
    1301         if ((b < a) && (b < c)) {
     1362        if ((fb < fa) && (fb < fc)) {
    13021363            bracket->data.F32[0] = a;
    13031364            bracket->data.F32[1] = b;
    … …  
    17381799
    17391800                psTrace(".psLib.dataManip.psMinimizePowell", 6,
    1740                         "WOW (%f %f %f)\n", u->data.F32[i], Q->data.F32[i], params->data.F32[i]);
     1801                        "u=Q-P (%f = %f - %f)\n", u->data.F32[i],
     1802                        Q->data.F32[i],
     1803                        params->data.F32[i]);
    17411804
    17421805            } else {
    … …  
    17471810        // 5: Move Q only in the direction u, and minimize the function.
    17481811        psTrace(".psLib.dataManip.psMinimizePowell", 6, "HERE 01\n");
    1749         psTrace(".psLib.dataManip.psMinimizePowell", 6,
    1750                 "u[0] is %f\n", u->data.F32[0]);
     1812        for (i=0;i<numDims;i++) {
     1813            psTrace(".psLib.dataManip.psMinimizePowell", 6,
     1814                    "u[i] is %f\n", u->data.F32[i]);
     1815        }
    17511816
    17521817        mul = psLineMin(min, params, u, paramMask, coords, func);

  • trunk/psLib/test/dataManip/Makefile

    r1798 r1819  
    33##  Makefile:   test/sysUtils
    44##
    5 ##  $Revision: 1.38 $  $Name: not supported by cvs2svn $
    6 ##  $Date: 2004-09-11 03:03:18 $
     5##  $Revision: 1.39 $  $Name: not supported by cvs2svn $
     6##  $Date: 2004-09-16 19:30:43 $
    77##
    88##  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    4545tst_psMinimize04 \
    4646tst_psMinimize05 \
     47tst_psMinimize06 \
    4748tst_psStats00 \
    4849tst_psStats01 \
Note: See TracChangeset for help on using the changeset viewer.