IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1079


Ignore:
Timestamp:
Jun 23, 2004, 4:39:34 PM (22 years ago)
Author:
gusciora
Message:

...

Location:
trunk/psLib/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r1022 r1079  
    55#include <float.h>
    66#include <math.h>
     7#include <gsl/gsl_rng.h>
     8#include <gsl/gsl_randist.h>
     9#include <gsl/gsl_vector.h>
     10#include <gsl/gsl_blas.h>
     11#include <gsl/gsl_multifit_nlin.h>
    712
    813#include "psMemory.h"
    … …  
    1924#define MAX_POLYNOMIAL_TERMS  (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2)
    2025int MyInfoLevel = 0;
     26
     27typedef struct
     28{
     29    size_t n;
     30    psVector *initialGuess;
     31    psVector *data;
     32    psVector *errors;
     33    psVector *paramMask;
     34    float (*evalModel) (psVector *, psVector *);
     35    float (*d_evalModel) (psVector *, psVector *);
     36}
     37psmodeldata;
     38
     39
     40
     41
    2142/******************************************************************************
    2243    This routine must minimize an arbitrary function.
    … …  
    3051}
    3152
     53
     54
     55int gsl_function_f(const gsl_vector *x,
     56                   void *params,
     57                   gsl_vector *f)
     58{
     59    size_t n  = ((psmodeldata *)params)->n;
     60    psVector *data   = ((psmodeldata *)params)->data;
     61    psVector *errors = ((psmodeldata *) params)->errors;
     62    psVector *mask   = ((psmodeldata *) params)->paramMask;
     63    psVector *initialguess = ((psmodeldata *)params)->initialguess;
     64    int count = 0;
     65    float (*evalModel) = ((psmodeldata *) params)->evalModel;
     66    size_t i;
     67
     68    for (i=0;i<initialguess->n;i++) {
     69        if (mask->data.U8[i] != 0) {
     70            count++;
     71        }
     72    }
     73
     74    if (mask != NULL) {
     75        j = 0;
     76        inputparameterlist = psVectorAlloc(count, PS_TYPE_F32);
     77        for (i=0;i<initialguess->n;i++) {
     78            if (mask->data.U8[i] == 0) {
     79                inputparameterlist->data.F32[j++] = initialguess->data.F32[i];
     80            }
     81        }
     82    } else {
     83        for (i=0;i<initialguess->n;i++) {
     84            inputparameterlist->data.F32[i] = initialguess->data.F32[i];
     85        }
     86    }
     87    // PAUL: what is this?
     88    evalModel(entiredomain, inputparameterlist);
     89    ???
     90
     91    loop i over domain size n
     92    {
     93        gsl_vector_set (f, i, (modelvalue[i] - data[i])/errors[i]);
     94    }
     95    return GSL_SUCCESS;
     96}
     97
     98int gsl_function_df(const gsl_vector * x,
     99                    void *params,
     100                    gsl_matrix* J)
     101{
     102    size_t n  = ((psmodeldata *)params)->n;
     103    psVector *data   = ((psmodeldata *)params)->data;
     104    psVector *errors = ((psmodeldata *) params)->errors;
     105    psVector *mask   = ((psmodeldata *) params)->paramMask;
     106    psVector *initialguess = ((psmodeldata *)params)->initialguess;
     107    int count = 0;
     108    float (*d_evalModel) = ((psmodeldata *) params)->evalModel;
     109    size_t i;
     110
     111    for (i=0;i<initialguess->n;i++) {
     112        if (mask->data.U8[i] != 0) {
     113            count++;
     114        }
     115    }
     116
     117
     118    if (mask != NULL) {
     119        j = 0;
     120        inputparameterlist = psVectorAlloc(count, PS_TYPE_F32);
     121        for (i=0;i<initialguess->n;i++) {
     122            if (mask->data.U8[i] == 0) {
     123                inputparameterlist->data.F32[j++] = initialguess->data.F32[i];
     124            }
     125        }
     126    } else {
     127        for (i=0;i<initialguess->n;i++) {
     128            inputparameterlist->data.F32[i] = initialguess->data.F32[i];
     129        }
     130    }
     131
     132    d_evalModel(entiredomain, inputparameterlist);
     133    ???
     134    loop i over domain size n
     135    // Jacobian matrix J(i,j) = dfi / dxj, where fi = (Modeli - datai)/errors[i]
     136    {
     137        loop j over the number of coefficients in model (excluding masks)
     138        {
     139            gsl_matrix_set (J, i, j, (d_modelvalue[i][j])/errors[i]);
     140        }
     141    }
     142    return GSL_SUCCESS;
     143}
     144
     145int gsl_function_fdf(const gsl_vector * x,
     146                     void *params,
     147                     gsl_vector * f,
     148                     gsl_matrix * J)
     149{
     150    gsl_function_f (x, params, f);
     151    gsl_function_df (x, params, J);
     152
     153    return GSL_SUCCESS;
     154}
     155
     156
     157
    32158/******************************************************************************
    33159    This routine must minimize an arbitrary function.
    34160 *****************************************************************************/
    35161psVector *
    36 psMinimizeChi2(float (*evalModel)(const psVector *restrict,
    37                                   const psVector *restrict),
     162psMinimizeChi2(float (*evalModel)(const psVector *restrict, const psVector *restrict),
     163               float (*DevalModel)(psVector *, psVector *),
    38164               const psVector *restrict domain,
    39165               const psVector *restrict data,
    … …  
    42168               const psVector *restrict paramMask)
    43169{
    44     return(NULL);
    45 }
     170    size_t n = size(domain);
     171    size_t p = size(initialGuess);
     172    gsl_multifit_fdfsolver_type *T = NULL;
     173    gsl_multifit_fdfsolver *s = NULL;
     174    int status;
     175    size_t i = 0;
     176    size_t iter = 0;
     177    psmodeldata inputdata = { n, initialguess, data, errors, paramMask, evalmodel, d_evalmodel};
     178    gsl_multifit_function_fdf f;
     179    double *x_init = NULL;
     180
     181
     182    //GUS: call psStatsGetnValues, allocate x_init.
     183    /* compress initialguess in light of mask */
     184    j = 0;
     185    loop i over length of initialguess
     186    {
     187        if not masked[i]
     188        {
     189            x_init[j] = initialguess[i];
     190            j++;
     191        } else
     192        { /* reduce p accordingly */
     193            p--;
     194        }
     195}
     196
     197
     198// Creates the vector for x which GSL uses.  Must deallocate.
     199gsl_vector_view x = gsl_vector_view_array (x_init, p);
     200
     201    const gsl_rng_type * type;
     202    gsl_rng *r;
     203    gsl_rng_env_setup();
     204
     205    type = gsl_rng_default;
     206    r = gsl_rng_alloc (type);
     207
     208    f.f = &gsl_function_f;
     209    f.df = &gsl_function_df;
     210    f.fdf = &gsl_function_fdf;
     211    f.n = n;
     212    f.p = p;
     213    f.params = &inputdata;
     214
     215    T = gsl_multifit_fdfsolver_lmsder;
     216    s = gsl_multifit_fdfsolver_alloc (T, n, p);
     217    gsl_multifit_fdfsolver_set (s, &f, &x.vector);
     218
     219    print_state (iter, s);
     220
     221    do {
     222        iter++;
     223        status = gsl_multifit_fdfsolver_iterate (s);
     224
     225        printf ("status = %s\n", gsl_strerror (status));
     226
     227        print_state (iter, s);
     228
     229        // Figure this out
     230        if (status)
     231            break;
     232
     233        // Checks whether the (L2-norm) computed derivative and the difference
     234        // between the real/actual for that test x-vector.  If were close enough
     235        // exit loop.
     236        status = gsl_multifit_test_delta (s->dx, s->x,
     237                                          1e-4, 1e-4);
     238    } while (status == GSL_CONTINUE && iter < 500);
     239    gsl_matrix *covar = gsl_matrix_alloc (p, p);
     240    gsl_multifit_covar (s->J, 0.0, covar);
     241    gsl_matrix_fprintf (stdout, covar, "%g");
     242
     243    /* add in masked values */
     244    j = 0;
     245    loop i over the length of initialguess
     246    {
     247        if mask(i)
     248    {
     249        // do nothing, initialguess shouldn't have changed */
     250    } else
     251    {
     252        initialguess(i)  = gsl_vector_get(s->x, j);  /* output */
     253            uncertainity = sqrt(gsl_matrix_get(covar,j,j));
     254            j++;
     255        }
     256    }
     257    double chi = gsl_blas_dnrm2(s->f);
     258    gsl_multifit_fdfsolver_free (s);
     259    return 0;
     260}
     261
     262// For 3-element vectors
     263int print_state(size_t iter,
     264                gsl_multifit_fdfsolver * s)
     265{
     266    printf ("iter: %3u x = % 15.8f % 15.8f % 15.8f "
     267            "|f(x)| = %g\n",
     268            iter,
     269            gsl_vector_get (s->x, 0),
     270            gsl_vector_get (s->x, 1),
     271            gsl_vector_get (s->x, 2),
     272            gsl_blas_dnrm2 (s->f));
     273}
     274
     275
     276
     277
     278
     279
    46280
    47281/** @brief This procedure calculates various combinations of powers of x and y

  • trunk/psLib/src/image/psImageStats.h

    r1022 r1079  
    1010 *  @author George Gusciora, MHPCC
    1111 *
    12  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2004-06-14 19:40:14 $
     12 *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2004-06-24 02:39:34 $
    1414 *
    1515 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    4848                     );
    4949
     50float psImagePixelInterpolation(psImage *input,
     51                                float x,
     52                                float y);
    5053#endif

  • trunk/psLib/src/imageops/psImageStats.h

    r1022 r1079  
    1010 *  @author George Gusciora, MHPCC
    1111 *
    12  *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
    13  *  @date $Date: 2004-06-14 19:40:14 $
     12 *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
     13 *  @date $Date: 2004-06-24 02:39:34 $
    1414 *
    1515 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    4848                     );
    4949
     50float psImagePixelInterpolation(psImage *input,
     51                                float x,
     52                                float y);
    5053#endif

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

    r1022 r1079  
    55#include <float.h>
    66#include <math.h>
     7#include <gsl/gsl_rng.h>
     8#include <gsl/gsl_randist.h>
     9#include <gsl/gsl_vector.h>
     10#include <gsl/gsl_blas.h>
     11#include <gsl/gsl_multifit_nlin.h>
    712
    813#include "psMemory.h"
    … …  
    1924#define MAX_POLYNOMIAL_TERMS  (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2)
    2025int MyInfoLevel = 0;
     26
     27typedef struct
     28{
     29    size_t n;
     30    psVector *initialGuess;
     31    psVector *data;
     32    psVector *errors;
     33    psVector *paramMask;
     34    float (*evalModel) (psVector *, psVector *);
     35    float (*d_evalModel) (psVector *, psVector *);
     36}
     37psmodeldata;
     38
     39
     40
     41
    2142/******************************************************************************
    2243    This routine must minimize an arbitrary function.
    … …  
    3051}
    3152
     53
     54
     55int gsl_function_f(const gsl_vector *x,
     56                   void *params,
     57                   gsl_vector *f)
     58{
     59    size_t n  = ((psmodeldata *)params)->n;
     60    psVector *data   = ((psmodeldata *)params)->data;
     61    psVector *errors = ((psmodeldata *) params)->errors;
     62    psVector *mask   = ((psmodeldata *) params)->paramMask;
     63    psVector *initialguess = ((psmodeldata *)params)->initialguess;
     64    int count = 0;
     65    float (*evalModel) = ((psmodeldata *) params)->evalModel;
     66    size_t i;
     67
     68    for (i=0;i<initialguess->n;i++) {
     69        if (mask->data.U8[i] != 0) {
     70            count++;
     71        }
     72    }
     73
     74    if (mask != NULL) {
     75        j = 0;
     76        inputparameterlist = psVectorAlloc(count, PS_TYPE_F32);
     77        for (i=0;i<initialguess->n;i++) {
     78            if (mask->data.U8[i] == 0) {
     79                inputparameterlist->data.F32[j++] = initialguess->data.F32[i];
     80            }
     81        }
     82    } else {
     83        for (i=0;i<initialguess->n;i++) {
     84            inputparameterlist->data.F32[i] = initialguess->data.F32[i];
     85        }
     86    }
     87    // PAUL: what is this?
     88    evalModel(entiredomain, inputparameterlist);
     89    ???
     90
     91    loop i over domain size n
     92    {
     93        gsl_vector_set (f, i, (modelvalue[i] - data[i])/errors[i]);
     94    }
     95    return GSL_SUCCESS;
     96}
     97
     98int gsl_function_df(const gsl_vector * x,
     99                    void *params,
     100                    gsl_matrix* J)
     101{
     102    size_t n  = ((psmodeldata *)params)->n;
     103    psVector *data   = ((psmodeldata *)params)->data;
     104    psVector *errors = ((psmodeldata *) params)->errors;
     105    psVector *mask   = ((psmodeldata *) params)->paramMask;
     106    psVector *initialguess = ((psmodeldata *)params)->initialguess;
     107    int count = 0;
     108    float (*d_evalModel) = ((psmodeldata *) params)->evalModel;
     109    size_t i;
     110
     111    for (i=0;i<initialguess->n;i++) {
     112        if (mask->data.U8[i] != 0) {
     113            count++;
     114        }
     115    }
     116
     117
     118    if (mask != NULL) {
     119        j = 0;
     120        inputparameterlist = psVectorAlloc(count, PS_TYPE_F32);
     121        for (i=0;i<initialguess->n;i++) {
     122            if (mask->data.U8[i] == 0) {
     123                inputparameterlist->data.F32[j++] = initialguess->data.F32[i];
     124            }
     125        }
     126    } else {
     127        for (i=0;i<initialguess->n;i++) {
     128            inputparameterlist->data.F32[i] = initialguess->data.F32[i];
     129        }
     130    }
     131
     132    d_evalModel(entiredomain, inputparameterlist);
     133    ???
     134    loop i over domain size n
     135    // Jacobian matrix J(i,j) = dfi / dxj, where fi = (Modeli - datai)/errors[i]
     136    {
     137        loop j over the number of coefficients in model (excluding masks)
     138        {
     139            gsl_matrix_set (J, i, j, (d_modelvalue[i][j])/errors[i]);
     140        }
     141    }
     142    return GSL_SUCCESS;
     143}
     144
     145int gsl_function_fdf(const gsl_vector * x,
     146                     void *params,
     147                     gsl_vector * f,
     148                     gsl_matrix * J)
     149{
     150    gsl_function_f (x, params, f);
     151    gsl_function_df (x, params, J);
     152
     153    return GSL_SUCCESS;
     154}
     155
     156
     157
    32158/******************************************************************************
    33159    This routine must minimize an arbitrary function.
    34160 *****************************************************************************/
    35161psVector *
    36 psMinimizeChi2(float (*evalModel)(const psVector *restrict,
    37                                   const psVector *restrict),
     162psMinimizeChi2(float (*evalModel)(const psVector *restrict, const psVector *restrict),
     163               float (*DevalModel)(psVector *, psVector *),
    38164               const psVector *restrict domain,
    39165               const psVector *restrict data,
    … …  
    42168               const psVector *restrict paramMask)
    43169{
    44     return(NULL);
    45 }
     170    size_t n = size(domain);
     171    size_t p = size(initialGuess);
     172    gsl_multifit_fdfsolver_type *T = NULL;
     173    gsl_multifit_fdfsolver *s = NULL;
     174    int status;
     175    size_t i = 0;
     176    size_t iter = 0;
     177    psmodeldata inputdata = { n, initialguess, data, errors, paramMask, evalmodel, d_evalmodel};
     178    gsl_multifit_function_fdf f;
     179    double *x_init = NULL;
     180
     181
     182    //GUS: call psStatsGetnValues, allocate x_init.
     183    /* compress initialguess in light of mask */
     184    j = 0;
     185    loop i over length of initialguess
     186    {
     187        if not masked[i]
     188        {
     189            x_init[j] = initialguess[i];
     190            j++;
     191        } else
     192        { /* reduce p accordingly */
     193            p--;
     194        }
     195}
     196
     197
     198// Creates the vector for x which GSL uses.  Must deallocate.
     199gsl_vector_view x = gsl_vector_view_array (x_init, p);
     200
     201    const gsl_rng_type * type;
     202    gsl_rng *r;
     203    gsl_rng_env_setup();
     204
     205    type = gsl_rng_default;
     206    r = gsl_rng_alloc (type);
     207
     208    f.f = &gsl_function_f;
     209    f.df = &gsl_function_df;
     210    f.fdf = &gsl_function_fdf;
     211    f.n = n;
     212    f.p = p;
     213    f.params = &inputdata;
     214
     215    T = gsl_multifit_fdfsolver_lmsder;
     216    s = gsl_multifit_fdfsolver_alloc (T, n, p);
     217    gsl_multifit_fdfsolver_set (s, &f, &x.vector);
     218
     219    print_state (iter, s);
     220
     221    do {
     222        iter++;
     223        status = gsl_multifit_fdfsolver_iterate (s);
     224
     225        printf ("status = %s\n", gsl_strerror (status));
     226
     227        print_state (iter, s);
     228
     229        // Figure this out
     230        if (status)
     231            break;
     232
     233        // Checks whether the (L2-norm) computed derivative and the difference
     234        // between the real/actual for that test x-vector.  If were close enough
     235        // exit loop.
     236        status = gsl_multifit_test_delta (s->dx, s->x,
     237                                          1e-4, 1e-4);
     238    } while (status == GSL_CONTINUE && iter < 500);
     239    gsl_matrix *covar = gsl_matrix_alloc (p, p);
     240    gsl_multifit_covar (s->J, 0.0, covar);
     241    gsl_matrix_fprintf (stdout, covar, "%g");
     242
     243    /* add in masked values */
     244    j = 0;
     245    loop i over the length of initialguess
     246    {
     247        if mask(i)
     248    {
     249        // do nothing, initialguess shouldn't have changed */
     250    } else
     251    {
     252        initialguess(i)  = gsl_vector_get(s->x, j);  /* output */
     253            uncertainity = sqrt(gsl_matrix_get(covar,j,j));
     254            j++;
     255        }
     256    }
     257    double chi = gsl_blas_dnrm2(s->f);
     258    gsl_multifit_fdfsolver_free (s);
     259    return 0;
     260}
     261
     262// For 3-element vectors
     263int print_state(size_t iter,
     264                gsl_multifit_fdfsolver * s)
     265{
     266    printf ("iter: %3u x = % 15.8f % 15.8f % 15.8f "
     267            "|f(x)| = %g\n",
     268            iter,
     269            gsl_vector_get (s->x, 0),
     270            gsl_vector_get (s->x, 1),
     271            gsl_vector_get (s->x, 2),
     272            gsl_blas_dnrm2 (s->f));
     273}
     274
     275
     276
     277
     278
     279
    46280
    47281/** @brief This procedure calculates various combinations of powers of x and y
Note: See TracChangeset for help on using the changeset viewer.