IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1185


Ignore:
Timestamp:
Jul 6, 2004, 4:38:32 PM (22 years ago)
Author:
gusciora
Message:

...

Location:
trunk/psLib
Files:
2 added
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/dataManip/Makefile

    r1133 r1185  
    1313           psMatrixVectorArithmetic.o \
    1414           psFFT.o \
    15            psImageIO.o
    16 #           psMinimize.o \
     15           psImageIO.o \
     16           psMinimize.o
    1717
    1818all: $(TARGET_STATIC)

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

    r1176 r1185  
    77
    88#include <gsl/gsl_multifit_nlin.h>
     9#include <gsl/gsl_multimin.h>
    910#include <gsl/gsl_rng.h>
    1011#include <gsl/gsl_randist.h>
    … …  
    3839}
    3940psModelData;
    40 
    4141// The first argument to evalModel() and
    4242// d_evalModel() specifies the data point.
    … …  
    4747
    4848
     49typedef struct
     50{
     51    int paramCount;   // Number of non-masked parameters.
     52    psVector *restrict       initialGuess;
     53    const psVector *restrict  coord;
     54    const psVector *restrict paramMask;
     55    float (*evalModel) (const psVector *, const psVector *);
     56    float (*d_evalModel) (const psVector *, const psVector *, int);
     57}
     58psModelData2;
     59
    4960
    5061/******************************************************************************
    51     This routine must minimize an arbitrary function.
     62p_psMinFunc(*v, *params): This routine calls the user supplied function to
     63be minimized.  The "v" argument of this function corresponds to the
     64parameters of the function to be minimized.
     65 
     66 
     67 *****************************************************************************/
     68double p_psMinFunc(const gsl_vector *v,
     69                   void *params)
     70{
     71    int i;    // Loop index variable.
     72    int j;    // Loop index variable.
     73    float tmpf;    // Temporary floating point variable.
     74    const psVector *restrict coord   = ((psModelData2 *) params)->coord;
     75    const psVector *restrict mask     = ((psModelData2 *) params)->paramMask;
     76    psVector *restrict initialGuess = ((psModelData2 *)params)->initialGuess;
     77    float (*evalModel)(const psVector *, const psVector *) =
     78        ((psModelData2 *) params)->evalModel;
     79    psVector *inputParameterList = NULL;
     80
     81    inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
     82
     83    if (mask != NULL) {
     84        j = 0;
     85        for (i=0;i<mask->n;i++) {
     86            if (mask->data.U8[i] != 0) {
     87                inputParameterList->data.F32[i] = initialGuess->data.F32[i];
     88            } else {
     89                inputParameterList->data.F32[i] = gsl_vector_get(v, j++);
     90            }
     91        }
     92    } else {
     93        for (i=0;i<initialGuess->n;i++) {
     94            inputParameterList->data.F32[i] = gsl_vector_get(v, i);
     95        }
     96    }
     97    //printf("HMMM: (%.1f %.1f %.1f %.1f)\n",
     98    //       inputParameterList->data.F32[0],
     99    //       inputParameterList->data.F32[1],
     100    //       coord->data.F32[0],
     101    //       coord->data.F32[1]);
     102
     103    tmpf = evalModel(inputParameterList, coord);
     104    psFree(inputParameterList);
     105    //printf("Called p_psMinFunc()\n");
     106    return(tmpf);
     107}
     108
     109void p_psMinFuncDeriv(const gsl_vector *v,
     110                      void *params,
     111                      gsl_vector *df)
     112{
     113    int i;    // Loop index variable.
     114    int j;    // Loop index variable.
     115    float tmpf;    // Temporary floating point variable.
     116    const psVector *restrict coord   = ((psModelData2 *) params)->coord;
     117    const psVector *restrict mask     = ((psModelData2 *) params)->paramMask;
     118    psVector *restrict initialGuess = ((psModelData2 *)params)->initialGuess;
     119    float (*d_evalModel)(const psVector *, const psVector *, int) =
     120        ((psModelData2 *) params)->d_evalModel;
     121    psVector *inputParameterList = NULL;
     122
     123    inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
     124
     125    if (mask != NULL) {
     126        j = 0;
     127        for (i=0;i<mask->n;i++) {
     128            if (mask->data.U8[i] != 0) {
     129                inputParameterList->data.F32[i] = initialGuess->data.F32[i];
     130            } else {
     131                inputParameterList->data.F32[i] = gsl_vector_get(v, j++);
     132            }
     133        }
     134    } else {
     135        for (i=0;i<initialGuess->n;i++) {
     136            inputParameterList->data.F32[i] = gsl_vector_get(v, i);
     137        }
     138    }
     139
     140    for (i=0;i<initialGuess->n;i++) {
     141        tmpf = d_evalModel(inputParameterList, coord, i);
     142        gsl_vector_set(df, i, tmpf);
     143    }
     144    psFree(inputParameterList);
     145}
     146
     147/* Compute both f and df together. */
     148void p_psMinFuncFuncDeriv(const gsl_vector *x,
     149                          void *params,
     150                          double *f,
     151                          gsl_vector *df)
     152{
     153    *f = p_psMinFunc(x, params);
     154    p_psMinFuncDeriv(x, params, df);
     155}
     156
     157
     158/******************************************************************************
     159psMinimize(initialGuess, myFunction, myFunctionDeriv, coord, paramMask):
     160 
     161This routine must minimize an arbitrary function; it must determine the set
     162of parameters of that function such that the
     163 
    52164 *****************************************************************************/
    53165psVector *
    54 psMinimize(float (*myFunction)(const psVector *restrict),
    55            psVector *restrict initialGuess,
    56            psVector *restrict paramMask)
    57 {
    58     return(NULL);
    59 }
    60 
    61 
     166psMinimize(psVector *restrict initialGuess,
     167           float (*myFunction)(const psVector *restrict, const psVector *restrict),
     168           float (*myFunctionDeriv)(const psVector *restrict, const psVector *restrict, int),
     169           const psVector *restrict coord,
     170           const psVector *restrict paramMask)
     171{
     172    int status;
     173    int i = 0;
     174    int j = 0;
     175    int iter = 0;
     176    gsl_multimin_function_fdf f;
     177    double *xInit = NULL;
     178    const gsl_multimin_fdfminimizer_type *T;
     179    gsl_multimin_fdfminimizer *s;
     180    psModelData2 inputData;
     181    gsl_vector *x;
     182
     183    inputData.initialGuess = initialGuess;
     184    inputData.coord = coord;
     185    inputData.paramMask = paramMask;
     186    inputData.evalModel = myFunction;
     187    inputData.d_evalModel = myFunctionDeriv;
     188    inputData.paramCount = 0;
     189
     190    printf("psMinimize(): initialGuess->n is %d\n", initialGuess->n);
     191    printf("psMinimize(): coord->n is %d\n", coord->n);
     192    printf("psMinimize(): coord is (%.1f, %.1f)\n", coord->data.F32[0],
     193           coord->data.F32[1]);
     194
     195    if (paramMask != NULL) {
     196        for (i=0;i<paramMask->n;i++) {
     197            if (paramMask->data.U8[i] != 0) {
     198                inputData.paramCount++;
     199            }
     200        }
     201    } else {
     202        inputData.paramCount= initialGuess->n;
     203    }
     204
     205    // The initial guess at the parameters for the function are written into
     206    // the vector inputParameterList.  If the paramMask is not NULL, then those
     207    // parameters are masked out.
     208
     209    xInit = (double *) psAlloc(inputData.paramCount * sizeof(double));
     210    if (paramMask != NULL) {
     211        j = 0;
     212        for (i=0;i<initialGuess->n;i++) {
     213            if (paramMask->data.U8[i] == 0) {
     214                xInit[j++] = initialGuess->data.F32[i];
     215            }
     216        }
     217    } else {
     218        for (i=0;i<initialGuess->n;i++) {
     219            xInit[i] = initialGuess->data.F32[i];
     220        }
     221    }
     222    f.f = &p_psMinFunc;
     223    f.df = &p_psMinFuncDeriv;
     224    f.fdf = &p_psMinFuncFuncDeriv;
     225    f.n = inputData.paramCount;
     226    // GUS: is this correct?
     227    //f.params = xInit;
     228    f.params = &inputData;
     229
     230
     231    printf("inputData.paramCount is %d\n", inputData.paramCount);
     232    x = gsl_vector_alloc(inputData.paramCount);
     233    for (i=0;i<inputData.paramCount;i++) {
     234        gsl_vector_set(x, i, xInit[i]);
     235    }
     236    T = gsl_multimin_fdfminimizer_conjugate_fr;
     237    s = gsl_multimin_fdfminimizer_alloc(T, inputData.paramCount);
     238
     239    gsl_multimin_fdfminimizer_set(s, &f, x, 0.01, 1e-4);
     240
     241    do {
     242        iter++;
     243        status = gsl_multimin_fdfminimizer_iterate(s);
     244
     245        if (status)
     246            break;
     247
     248        status = gsl_multimin_test_gradient(s->gradient, 1e-3);
     249
     250        if (status == GSL_SUCCESS)
     251            printf ("Minimum found at:\n");
     252
     253    } while (status == GSL_CONTINUE && iter < 100);
     254
     255    // In the above steps we had removed the masked elements from the
     256    // the solver.  This next code blocks puts those masked elements
     257    // into the solution.
     258    if (paramMask != NULL) {
     259        j = 0;
     260        for (i=0;i<initialGuess->n;i++) {
     261            if (paramMask->data.U8[i] == 0) {
     262                initialGuess->data.F32[i] = gsl_vector_get(s->x, j++);
     263            } else {
     264                initialGuess->data.F32[i] = initialGuess->data.F32[i];
     265            }
     266        }
     267    } else {
     268        for (i=0;i<initialGuess->n;i++) {
     269            initialGuess->data.F32[i] = gsl_vector_get(s->x, i);
     270        }
     271    }
     272    psFree(xInit);
     273    return(initialGuess);
     274}
     275
     276
     277
     278
     279
     280// The first argument to evalModel() and
     281// d_evalModel() specifies the data point.
     282// It must have the same size as the second
     283// dimension of *domain.
     284// The second argument must have the same size
     285// as *initialGuess and *paramMask.
    62286
    63287/******************************************************************************

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

    r1093 r1185  
    1414/** This routine must minimize a non-linear function */
    1515psVector *
    16 psMinimize(float (*myFunction)(const psVector *restrict), ///< Function to minimize
    17            psVector *restrict initialGuess, ///< Initial guess
    18            psVector *restrict paramMask //!< 1 = fit for parameter, 0 = hold parameter constant
    19           );
     16psMinimize(psVector *restrict initialGuess,
     17           float (*myFunction)(const psVector *restrict, const psVector *restrict),
     18           float (*myFunctionDeriv)(const psVector *restrict, const psVector *restrict, int),
     19           const psVector *restrict coord,
     20           const psVector *restrict paramMask);
    2021
    2122

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

    r1176 r1185  
    77
    88#include <gsl/gsl_multifit_nlin.h>
     9#include <gsl/gsl_multimin.h>
    910#include <gsl/gsl_rng.h>
    1011#include <gsl/gsl_randist.h>
    … …  
    3839}
    3940psModelData;
    40 
    4141// The first argument to evalModel() and
    4242// d_evalModel() specifies the data point.
    … …  
    4747
    4848
     49typedef struct
     50{
     51    int paramCount;   // Number of non-masked parameters.
     52    psVector *restrict       initialGuess;
     53    const psVector *restrict  coord;
     54    const psVector *restrict paramMask;
     55    float (*evalModel) (const psVector *, const psVector *);
     56    float (*d_evalModel) (const psVector *, const psVector *, int);
     57}
     58psModelData2;
     59
    4960
    5061/******************************************************************************
    51     This routine must minimize an arbitrary function.
     62p_psMinFunc(*v, *params): This routine calls the user supplied function to
     63be minimized.  The "v" argument of this function corresponds to the
     64parameters of the function to be minimized.
     65 
     66 
     67 *****************************************************************************/
     68double p_psMinFunc(const gsl_vector *v,
     69                   void *params)
     70{
     71    int i;    // Loop index variable.
     72    int j;    // Loop index variable.
     73    float tmpf;    // Temporary floating point variable.
     74    const psVector *restrict coord   = ((psModelData2 *) params)->coord;
     75    const psVector *restrict mask     = ((psModelData2 *) params)->paramMask;
     76    psVector *restrict initialGuess = ((psModelData2 *)params)->initialGuess;
     77    float (*evalModel)(const psVector *, const psVector *) =
     78        ((psModelData2 *) params)->evalModel;
     79    psVector *inputParameterList = NULL;
     80
     81    inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
     82
     83    if (mask != NULL) {
     84        j = 0;
     85        for (i=0;i<mask->n;i++) {
     86            if (mask->data.U8[i] != 0) {
     87                inputParameterList->data.F32[i] = initialGuess->data.F32[i];
     88            } else {
     89                inputParameterList->data.F32[i] = gsl_vector_get(v, j++);
     90            }
     91        }
     92    } else {
     93        for (i=0;i<initialGuess->n;i++) {
     94            inputParameterList->data.F32[i] = gsl_vector_get(v, i);
     95        }
     96    }
     97    //printf("HMMM: (%.1f %.1f %.1f %.1f)\n",
     98    //       inputParameterList->data.F32[0],
     99    //       inputParameterList->data.F32[1],
     100    //       coord->data.F32[0],
     101    //       coord->data.F32[1]);
     102
     103    tmpf = evalModel(inputParameterList, coord);
     104    psFree(inputParameterList);
     105    //printf("Called p_psMinFunc()\n");
     106    return(tmpf);
     107}
     108
     109void p_psMinFuncDeriv(const gsl_vector *v,
     110                      void *params,
     111                      gsl_vector *df)
     112{
     113    int i;    // Loop index variable.
     114    int j;    // Loop index variable.
     115    float tmpf;    // Temporary floating point variable.
     116    const psVector *restrict coord   = ((psModelData2 *) params)->coord;
     117    const psVector *restrict mask     = ((psModelData2 *) params)->paramMask;
     118    psVector *restrict initialGuess = ((psModelData2 *)params)->initialGuess;
     119    float (*d_evalModel)(const psVector *, const psVector *, int) =
     120        ((psModelData2 *) params)->d_evalModel;
     121    psVector *inputParameterList = NULL;
     122
     123    inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
     124
     125    if (mask != NULL) {
     126        j = 0;
     127        for (i=0;i<mask->n;i++) {
     128            if (mask->data.U8[i] != 0) {
     129                inputParameterList->data.F32[i] = initialGuess->data.F32[i];
     130            } else {
     131                inputParameterList->data.F32[i] = gsl_vector_get(v, j++);
     132            }
     133        }
     134    } else {
     135        for (i=0;i<initialGuess->n;i++) {
     136            inputParameterList->data.F32[i] = gsl_vector_get(v, i);
     137        }
     138    }
     139
     140    for (i=0;i<initialGuess->n;i++) {
     141        tmpf = d_evalModel(inputParameterList, coord, i);
     142        gsl_vector_set(df, i, tmpf);
     143    }
     144    psFree(inputParameterList);
     145}
     146
     147/* Compute both f and df together. */
     148void p_psMinFuncFuncDeriv(const gsl_vector *x,
     149                          void *params,
     150                          double *f,
     151                          gsl_vector *df)
     152{
     153    *f = p_psMinFunc(x, params);
     154    p_psMinFuncDeriv(x, params, df);
     155}
     156
     157
     158/******************************************************************************
     159psMinimize(initialGuess, myFunction, myFunctionDeriv, coord, paramMask):
     160 
     161This routine must minimize an arbitrary function; it must determine the set
     162of parameters of that function such that the
     163 
    52164 *****************************************************************************/
    53165psVector *
    54 psMinimize(float (*myFunction)(const psVector *restrict),
    55            psVector *restrict initialGuess,
    56            psVector *restrict paramMask)
    57 {
    58     return(NULL);
    59 }
    60 
    61 
     166psMinimize(psVector *restrict initialGuess,
     167           float (*myFunction)(const psVector *restrict, const psVector *restrict),
     168           float (*myFunctionDeriv)(const psVector *restrict, const psVector *restrict, int),
     169           const psVector *restrict coord,
     170           const psVector *restrict paramMask)
     171{
     172    int status;
     173    int i = 0;
     174    int j = 0;
     175    int iter = 0;
     176    gsl_multimin_function_fdf f;
     177    double *xInit = NULL;
     178    const gsl_multimin_fdfminimizer_type *T;
     179    gsl_multimin_fdfminimizer *s;
     180    psModelData2 inputData;
     181    gsl_vector *x;
     182
     183    inputData.initialGuess = initialGuess;
     184    inputData.coord = coord;
     185    inputData.paramMask = paramMask;
     186    inputData.evalModel = myFunction;
     187    inputData.d_evalModel = myFunctionDeriv;
     188    inputData.paramCount = 0;
     189
     190    printf("psMinimize(): initialGuess->n is %d\n", initialGuess->n);
     191    printf("psMinimize(): coord->n is %d\n", coord->n);
     192    printf("psMinimize(): coord is (%.1f, %.1f)\n", coord->data.F32[0],
     193           coord->data.F32[1]);
     194
     195    if (paramMask != NULL) {
     196        for (i=0;i<paramMask->n;i++) {
     197            if (paramMask->data.U8[i] != 0) {
     198                inputData.paramCount++;
     199            }
     200        }
     201    } else {
     202        inputData.paramCount= initialGuess->n;
     203    }
     204
     205    // The initial guess at the parameters for the function are written into
     206    // the vector inputParameterList.  If the paramMask is not NULL, then those
     207    // parameters are masked out.
     208
     209    xInit = (double *) psAlloc(inputData.paramCount * sizeof(double));
     210    if (paramMask != NULL) {
     211        j = 0;
     212        for (i=0;i<initialGuess->n;i++) {
     213            if (paramMask->data.U8[i] == 0) {
     214                xInit[j++] = initialGuess->data.F32[i];
     215            }
     216        }
     217    } else {
     218        for (i=0;i<initialGuess->n;i++) {
     219            xInit[i] = initialGuess->data.F32[i];
     220        }
     221    }
     222    f.f = &p_psMinFunc;
     223    f.df = &p_psMinFuncDeriv;
     224    f.fdf = &p_psMinFuncFuncDeriv;
     225    f.n = inputData.paramCount;
     226    // GUS: is this correct?
     227    //f.params = xInit;
     228    f.params = &inputData;
     229
     230
     231    printf("inputData.paramCount is %d\n", inputData.paramCount);
     232    x = gsl_vector_alloc(inputData.paramCount);
     233    for (i=0;i<inputData.paramCount;i++) {
     234        gsl_vector_set(x, i, xInit[i]);
     235    }
     236    T = gsl_multimin_fdfminimizer_conjugate_fr;
     237    s = gsl_multimin_fdfminimizer_alloc(T, inputData.paramCount);
     238
     239    gsl_multimin_fdfminimizer_set(s, &f, x, 0.01, 1e-4);
     240
     241    do {
     242        iter++;
     243        status = gsl_multimin_fdfminimizer_iterate(s);
     244
     245        if (status)
     246            break;
     247
     248        status = gsl_multimin_test_gradient(s->gradient, 1e-3);
     249
     250        if (status == GSL_SUCCESS)
     251            printf ("Minimum found at:\n");
     252
     253    } while (status == GSL_CONTINUE && iter < 100);
     254
     255    // In the above steps we had removed the masked elements from the
     256    // the solver.  This next code blocks puts those masked elements
     257    // into the solution.
     258    if (paramMask != NULL) {
     259        j = 0;
     260        for (i=0;i<initialGuess->n;i++) {
     261            if (paramMask->data.U8[i] == 0) {
     262                initialGuess->data.F32[i] = gsl_vector_get(s->x, j++);
     263            } else {
     264                initialGuess->data.F32[i] = initialGuess->data.F32[i];
     265            }
     266        }
     267    } else {
     268        for (i=0;i<initialGuess->n;i++) {
     269            initialGuess->data.F32[i] = gsl_vector_get(s->x, i);
     270        }
     271    }
     272    psFree(xInit);
     273    return(initialGuess);
     274}
     275
     276
     277
     278
     279
     280// The first argument to evalModel() and
     281// d_evalModel() specifies the data point.
     282// It must have the same size as the second
     283// dimension of *domain.
     284// The second argument must have the same size
     285// as *initialGuess and *paramMask.
    62286
    63287/******************************************************************************

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

    r1093 r1185  
    1414/** This routine must minimize a non-linear function */
    1515psVector *
    16 psMinimize(float (*myFunction)(const psVector *restrict), ///< Function to minimize
    17            psVector *restrict initialGuess, ///< Initial guess
    18            psVector *restrict paramMask //!< 1 = fit for parameter, 0 = hold parameter constant
    19           );
     16psMinimize(psVector *restrict initialGuess,
     17           float (*myFunction)(const psVector *restrict, const psVector *restrict),
     18           float (*myFunctionDeriv)(const psVector *restrict, const psVector *restrict, int),
     19           const psVector *restrict coord,
     20           const psVector *restrict paramMask);
    2021
    2122

  • trunk/psLib/test/dataManip/Makefile

    r1184 r1185  
    33##  Makefile:   test/sysUtils
    44##
    5 ##  $Revision: 1.29 $  $Name: not supported by cvs2svn $
    6 ##  $Date: 2004-07-02 01:08:54 $
     5##  $Revision: 1.30 $  $Name: not supported by cvs2svn $
     6##  $Date: 2004-07-07 02:38:32 $
    77##
    88##  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    4040 tst_psFunc00 \
    4141 tst_psFunc01 \
     42 tst_psMinimize00 \
    4243 tst_psMinimize01 \
    43  tst_psMinimize00 \
     44 tst_psMinimize02 \
     45 tst_psMinimize03 \
    4446 tst_psImageStats00 \
    4547 tst_psImageStats01 \
    … …  
    7779tst_psMinimize00:       tst_psMinimize00.o
    7880tst_psMinimize01:       tst_psMinimize01.o
     81tst_psMinimize02:       tst_psMinimize02.o
     82tst_psMinimize03:       tst_psMinimize03.o
    7983tst_psImageStats00:     tst_psImageStats00.o
    8084tst_psImageStats01:     tst_psImageStats01.o
Note: See TracChangeset for help on using the changeset viewer.