IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1093


Ignore:
Timestamp:
Jun 24, 2004, 8:51:20 PM (22 years ago)
Author:
gusciora
Message:

...

Location:
trunk/psLib/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

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

    r1086 r1093  
    1919#include "psFunctions.h"
    2020#include "psSort.h"
     21#include "psMinimize.h"
    2122#include "float.h"
    2223#include <math.h>
    23 // NOTE: rewrite so there is no maximum order for the polynomials.
    24 #define MAX_POLY_ORDER 10
    25 #define MAX_POLYNOMIAL_TERMS  (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2)
    26 int MyInfoLevel = 0;
    27 
     24
     25#define MAX_LMM_NUM_ITERATIONS 500
    2826typedef struct
    2927{
    3028    size_t n;   // Number of data points points in domain.
    3129    int count;   // Number of non-masked parameters.
    32     psVector *restrict initialGuess;
    33     psImage  *domain;  //
    34     psVector *data;
    35     psVector *errors;
    36     psVector *paramMask;
    37     // The first argument to evalModel() and
    38     // d_evalModel() specifies the data point.
    39     // It must have the same size as the second
    40     // dimension of *domain.
    41     // The second argument must have the same size
    42     // as *initialGuess and *paramMask.
     30    psVector *restrict       initialGuess;
     31    const psImage *restrict  domain;
     32    const psVector *restrict data;
     33    const psVector *restrict  errors;
     34    const psVector *restrict paramMask;
    4335    float (*evalModel) (psVector *, psVector *);
    44     float (*d_evalModel) (psVector *, psVector *);
     36    float (*d_evalModel) (psVector *, psVector *, int);
    4537}
    4638psModelData;
     39// The first argument to evalModel() and
     40// d_evalModel() specifies the data point.
     41// It must have the same size as the second
     42// dimension of *domain.
     43// The second argument must have the same size
     44// as *initialGuess and *paramMask.
    4745
    4846
    … …  
    6563for the user supplied function which is to be minimized.  The GSL
    6664minimization routines have no knowledge of the user-supplied function.  They
    67 only clal this routine, which then calls the user-supplied function.  The
     65only call this routine, which then calls the user-supplied function.  The
    6866arguments are:
    6967    x: These are the parameter which are to be varied by GSL in order to
    7068       minimized chi2 over the data set.
    7169    params:
    72     f
     70    f:
    7371 
    7472 
    … …  
    7977                   gsl_vector *f)
    8078{
    81     psImage  *domain   = ((psModelData *)params)->domain;
    82     psVector *data     = ((psModelData *)params)->data;
    83     psVector *errors   = ((psModelData *) params)->errors;
    84     psVector *mask     = ((psModelData *) params)->paramMask;
    85     psVector *initialGuess = ((psModelData *)params)->initialGuess;
     79    int i;    // Loop index variable.
     80    int j;    // Loop index variable.
     81    float tmpf;    // Temporary floating point variable.
     82    const psImage *restrict  domain   = ((psModelData *)params)->domain;
     83    const psVector *restrict data     = ((psModelData *)params)->data;
     84    const psVector *restrict errors   = ((psModelData *) params)->errors;
     85    const psVector *restrict mask     = ((psModelData *) params)->paramMask;
     86    psVector *restrict initialGuess = ((psModelData *)params)->initialGuess;
     87    float (*evalModel)(psVector *, psVector *) = ((psModelData *) params)->evalModel;
    8688    psVector *inputParameterList = NULL;
    87     //    float (*evalModel)(psVector *, psVector *) = ((psModelData *) params)->evalModel;
    88     size_t i;
    89     int j;
    90     float tmpf;
    91 
    92     // The GSL routines will call this functions with the masked parameters
     89    psVector *tmpVecPtr = NULL;
     90
     91
     92    tmpVecPtr = psVectorAlloc(domain->numCols, PS_TYPE_F32);
     93
     94    // The GSL routines will call this function with the masked parameters
    9395    // removed.  However, the user-supplied function (to be modified) does not
    9496    // have those parameters removed.  Here will create a new parameter list
    … …  
    112114
    113115    for (i=0;i<domain->numRows;i++) {
    114         //GUS        tmpf = evalModel(inputParameterList, domain->data.F32[i]);
     116        for (j=0;j<tmpVecPtr->n;j++) {
     117            tmpVecPtr->data.F32[j] = domain->data.F32[i][j];
     118        }
     119        tmpf = evalModel(inputParameterList, tmpVecPtr);
    115120        gsl_vector_set(f, i, (tmpf - data->data.F32[i])/errors->data.F32[i]);
    116121    }
    117122
     123    psFree(inputParameterList);
     124    psFree(tmpVecPtr);
    118125    return GSL_SUCCESS;
    119126}
    … …  
    123130                    gsl_matrix *J)
    124131{
    125     //    size_t n  = ((psModelData *)params)->n;
    126     psImage  *domain   = ((psModelData *)params)->domain;
    127     //    psVector *data   = ((psModelData *)params)->data;
    128     psVector *errors = ((psModelData *) params)->errors;
    129     psVector *mask   = ((psModelData *) params)->paramMask;
    130     psVector *initialGuess = ((psModelData *)params)->initialGuess;
     132    const psImage *restrict domain   = ((psModelData *)params)->domain;
     133    const psVector *restrict errors   = ((psModelData *) params)->errors;
     134    const psVector *restrict mask     = ((psModelData *) params)->paramMask;
     135    psVector *restrict initialGuess = ((psModelData *)params)->initialGuess;
    131136    psVector *inputParameterList = NULL;
    132     //    int count = 0;
    133     //    float (*d_evalModel) = ((psModelData *) params)->evalModel;
     137    psVector *tmpVecPtr = NULL;
     138    float (*d_evalModel)(psVector *, psVector *, int) = ((psModelData *) params)->d_evalModel;
     139
    134140    size_t i;
    135141    int j;
    136142    float tmpf;
     143
     144    tmpVecPtr = psVectorAlloc(domain->numCols, PS_TYPE_F32);
    137145
    138146    // The GSL routines will call this functions with the masked parameters
    … …  
    157165
    158166    for (i=0;i<domain->numRows;i++) {
     167        for (j=0;j<tmpVecPtr->n;j++) {
     168            tmpVecPtr->data.F32[j] = domain->data.F32[i][j];
     169        }
     170
    159171        for (j=0;j<inputParameterList->n;j++) {
    160             //GUS            tmpf = d_evalModel(inputParameterList, domain->data.F32[i], j);
    161 
     172            tmpf = d_evalModel(inputParameterList, tmpVecPtr, j);
    162173            gsl_matrix_set(J, i, j, (tmpf/errors->data.F32[i]));
    163174        }
    164175    }
    165176
     177    psFree(inputParameterList);
     178    psFree(tmpVecPtr);
    166179    return GSL_SUCCESS;
    167180}
    … …  
    185198psVector *
    186199psMinimizeChi2(float (*evalModel)(const psVector *restrict, const psVector *restrict),
    187                float (*DevalModel)(psVector *, psVector *),
    188                const psVector *restrict domain,
     200               float (*DevalModel)(const psVector *restrict, const psVector *restrict, int),
     201               const psImage *restrict domain,
    189202               const psVector *restrict data,
    190203               const psVector *restrict errors,
    … …  
    193206               float *chiSq)
    194207{
    195     size_t n = domain->n;
     208
     209    int numData = domain->numCols; // Number of data points
     210    gsl_multifit_fdfsolver *s = NULL; // GSL data structure.
     211    int status;    // Return status for the GSL solver.
     212    int i = 0;    // Loop index variable.
     213    int j = 0;    // Loop index variable.
     214    int iter = 0;   // Iteration counter.
     215    psModelData *inputData = NULL; // Contains data that the user-supplied
     216    // function/derivate need.
     217    gsl_multifit_function_fdf f; // GSL structure that contains the
     218    // functions/derivative to be solved.
     219    double *xInit = NULL;  // The initial guess at the parameters
     220    // with masked parameters removed.
    196221    const gsl_multifit_fdfsolver_type *T = gsl_multifit_fdfsolver_lmsder;
    197     gsl_multifit_fdfsolver *s      = NULL;
    198     int status;
    199     size_t i = 0;
    200     int j = 0;
    201     size_t iter = 0;
    202     psModelData *inputData = NULL;
    203     gsl_multifit_function_fdf f;
    204     double *xInit = NULL;
     222    // This tells GSL to use the Levenberg-
     223    // Marquardt algorithm for chi2
     224    // minimization.
    205225
    206226    inputData = (psModelData *) psAlloc(sizeof(psModelData));
    207     inputData->n = n;
     227    inputData->n = numData;
    208228    inputData->count = 0;
    209229    inputData->initialGuess = initialGuess;
    210     //    inputData->domain = domain;
     230    inputData->domain = domain;
    211231    //    inputData->data = data;
    212232    //    inputData->errors = errors;
    … …  
    226246    }
    227247
    228 
    229 
    230248    // The initial guess at the parameters for the function are written into
    231249    // the vector inputParameterList.  If the paramMask is not NULL, then those
    232     // parameters are masked out.  How can this possibly work?  The user-
    233     // supplied function will require a fixed number of parameters.
     250    // parameters are masked out.
    234251
    235252    xInit = (double *) psAlloc(inputData->count * sizeof(double));
    … …  
    257274    r = gsl_rng_alloc(type);
    258275
     276    // Initialize the main data structure used by the GSL solver.  This will
     277    // contain pointers to the function to be minimized, it's derivative
     278    // function, the number of data points, the number of free parameters,
     279    // and the data structures those functions use.
     280
    259281    f.f = &gsl_function_f;
    260282    f.df = &gsl_function_df;
    261283    f.fdf = &gsl_my_function_fdf;
    262     f.n = n;
     284    f.n = numData;
    263285    f.p = inputData->count;
    264286    f.params = &inputData;
    265287
    266     // This tells GSL to use the Levenberg-Marquardt algorithm for chi2
    267     // minimization.
    268     //    T = gsl_multifit_fdfsolver_lmsder;
    269 
    270     // Creates an instance of the GSL solver that we will be iterating on.
    271     s = gsl_multifit_fdfsolver_alloc(T, n, inputData->count);
    272 
    273     // Initialize the GSL minimizer.
     288    // Create an instance of the GSL solver that we will be iterating on.
     289    // It will have numData data points and inputData->count parameters.
     290    s = gsl_multifit_fdfsolver_alloc(T, numData, inputData->count);
     291
     292    // Initialize the GSL minimizer to use function defined by the data
     293    // structure "f" and x.vector as an initial guess for the parameters.
    274294    gsl_multifit_fdfsolver_set(s, &f, &x.vector);
    275295
    276     //    print_state(iter, s);
     296
     297    // Each iteration of the following loop will perform one step in an
     298    // attempt to minimized chi-squared for the function.  The loop exits
     299    // either when the change in parameters is small enough, or when the
     300    // maximum number of iterations is reached.
    277301    do {
    278302        iter++;
    … …  
    280304        status = gsl_multifit_fdfsolver_iterate(s);
    281305
    282         //        printf("status = %s\n", gsl_strerror(status));
    283         //        print_state(iter, s);
    284 
    285306        // If there was a problem, abort.
    286307        if (status) {
    … …  
    288309        }
    289310
    290         // Checks whether the (L2-norm) computed derivative and the difference
    291         // between the real/actual for that test x-vector.  If were close
    292         // enough exit loop.
    293 
    294311        // Test if the parameters changed by a small enough amount.
    295312        status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
    296     } while (status == GSL_CONTINUE && iter < 500);
    297 
    298 
    299 
    300 
     313    } while (status == GSL_CONTINUE && iter < MAX_LMM_NUM_ITERATIONS);
     314
     315
     316    // In the above steps we had removed the masked elements from the
     317    // the solver.  This next code blocks puts those masked elements
     318    // into the solution.
    301319    if (paramMask != NULL) {
    302320        j = 0;
    … …  
    314332    }
    315333
     334    // Calculate the chi-squared for the derived solution.
    316335    *chiSq = gsl_blas_dnrm2(s->f);
     336
     337    // Free all allocated memory
     338    // NOTE: Free x.
    317339    gsl_multifit_fdfsolver_free(s);
    318     // Free all allocated memory
    319340    psFree(xInit);
    320341    psFree(inputData);
    321342
     343    // Bye bye.
    322344    return(initialGuess);
    323345}
    324346
     347
     348
     349
     350
     351
     352
     353
     354
     355
     356
     357
     358
     359
     360
     361// NOTE: rewrite so there is no maximum order for the polynomials.
     362#define MAX_POLY_ORDER 10
     363#define MAX_POLYNOMIAL_TERMS  (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2)
     364int MyInfoLevel = 0;
    325365/** @brief This procedure calculates various combinations of powers of x and y
    326366 *   and stores them in the data structure sums[][].  After it completes:

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

    r974 r1093  
    1212 */
    1313
    14 /** This routine must ninimize a particular non-linear function */
    15 psFloatArray *
    16 psMinimize(float (*myFunction)(const psFloatArray *restrict), ///< Function to minimize
    17            psFloatArray *restrict initialGuess, ///< Initial guess
    18            psIntArray *restrict paramMask //!< 1 = fit for parameter, 0 = hold parameter constant
     14/** This routine must minimize a non-linear function */
     15psVector *
     16psMinimize(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
    1919          );
    2020
    2121
    2222/** Minimize chi^2 for input data */
    23 psFloatArray *
    24 psMinimizeChi2(float (*evalModel)(const psFloatArray *restrict,
    25                                   const psFloatArray *restrict), ///< Model to fit; (domain and params)
    26                const psFloatArray *restrict domain, ///< The domain values for the corresponding measurements
    27                const psFloatArray *restrict data, ///< Data to fit
    28                const psFloatArray *restrict errors, ///< Errors in the data
    29                psFloatArray *restrict initialGuess, ///< Initial guess
    30                const psIntArray *restrict paramMask ///< 1 = fit for parameter, 0 = hold parameter constant
     23psVector *
     24psMinimizeChi2(float (*evalModel)(const psVector *restrict,
     25                                  const psVector *restrict), ///< Model to fit; (domain and params)
     26               float (*DevalModel)(const psVector *restrict,
     27                                   const psVector *restrict,
     28                                   int), ///< Derivative of model to fit; (domain and params)
     29               const psImage *restrict domain, ///< The domain values for the corresponding measurements
     30               const psVector *restrict data, ///< Data to fit
     31               const psVector *restrict errors, ///< Errors in the data
     32               psVector *restrict initialGuess, ///< Initial guess
     33               const psVector *restrict paramMask, ///< 1 = fit for parameter, 0 = hold parameter constant
     34               float *chiSq
    3135              );
    3236
    33 /** Derive a polynomial fit by chi^2 minimisation --- can be done analytically */
     37/** Derive a polynomial fit by chi^2 minimisation (analytically) */
    3438psPolynomial1D *
    3539psGetArrayPolynomial(psPolynomial1D myPoly, ///< Polynomial to fit
    36                      const psFloatArray *restrict x, ///< Ordinates (or NULL to just use the indices)
    37                      const psFloatArray *restrict y, ///< Coordinates
    38                      const psFloatArray *restrict yErr ///< Errors in coordinates, or NULL
     40                     const psVector *restrict x, ///< Ordinates (or NULL to just use the indices)
     41                     const psVector *restrict y, ///< Coordinates
     42                     const psVector *restrict yErr ///< Errors in coordinates, or NULL
    3943                    );
    4044

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

    r1086 r1093  
    1919#include "psFunctions.h"
    2020#include "psSort.h"
     21#include "psMinimize.h"
    2122#include "float.h"
    2223#include <math.h>
    23 // NOTE: rewrite so there is no maximum order for the polynomials.
    24 #define MAX_POLY_ORDER 10
    25 #define MAX_POLYNOMIAL_TERMS  (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2)
    26 int MyInfoLevel = 0;
    27 
     24
     25#define MAX_LMM_NUM_ITERATIONS 500
    2826typedef struct
    2927{
    3028    size_t n;   // Number of data points points in domain.
    3129    int count;   // Number of non-masked parameters.
    32     psVector *restrict initialGuess;
    33     psImage  *domain;  //
    34     psVector *data;
    35     psVector *errors;
    36     psVector *paramMask;
    37     // The first argument to evalModel() and
    38     // d_evalModel() specifies the data point.
    39     // It must have the same size as the second
    40     // dimension of *domain.
    41     // The second argument must have the same size
    42     // as *initialGuess and *paramMask.
     30    psVector *restrict       initialGuess;
     31    const psImage *restrict  domain;
     32    const psVector *restrict data;
     33    const psVector *restrict  errors;
     34    const psVector *restrict paramMask;
    4335    float (*evalModel) (psVector *, psVector *);
    44     float (*d_evalModel) (psVector *, psVector *);
     36    float (*d_evalModel) (psVector *, psVector *, int);
    4537}
    4638psModelData;
     39// The first argument to evalModel() and
     40// d_evalModel() specifies the data point.
     41// It must have the same size as the second
     42// dimension of *domain.
     43// The second argument must have the same size
     44// as *initialGuess and *paramMask.
    4745
    4846
    … …  
    6563for the user supplied function which is to be minimized.  The GSL
    6664minimization routines have no knowledge of the user-supplied function.  They
    67 only clal this routine, which then calls the user-supplied function.  The
     65only call this routine, which then calls the user-supplied function.  The
    6866arguments are:
    6967    x: These are the parameter which are to be varied by GSL in order to
    7068       minimized chi2 over the data set.
    7169    params:
    72     f
     70    f:
    7371 
    7472 
    … …  
    7977                   gsl_vector *f)
    8078{
    81     psImage  *domain   = ((psModelData *)params)->domain;
    82     psVector *data     = ((psModelData *)params)->data;
    83     psVector *errors   = ((psModelData *) params)->errors;
    84     psVector *mask     = ((psModelData *) params)->paramMask;
    85     psVector *initialGuess = ((psModelData *)params)->initialGuess;
     79    int i;    // Loop index variable.
     80    int j;    // Loop index variable.
     81    float tmpf;    // Temporary floating point variable.
     82    const psImage *restrict  domain   = ((psModelData *)params)->domain;
     83    const psVector *restrict data     = ((psModelData *)params)->data;
     84    const psVector *restrict errors   = ((psModelData *) params)->errors;
     85    const psVector *restrict mask     = ((psModelData *) params)->paramMask;
     86    psVector *restrict initialGuess = ((psModelData *)params)->initialGuess;
     87    float (*evalModel)(psVector *, psVector *) = ((psModelData *) params)->evalModel;
    8688    psVector *inputParameterList = NULL;
    87     //    float (*evalModel)(psVector *, psVector *) = ((psModelData *) params)->evalModel;
    88     size_t i;
    89     int j;
    90     float tmpf;
    91 
    92     // The GSL routines will call this functions with the masked parameters
     89    psVector *tmpVecPtr = NULL;
     90
     91
     92    tmpVecPtr = psVectorAlloc(domain->numCols, PS_TYPE_F32);
     93
     94    // The GSL routines will call this function with the masked parameters
    9395    // removed.  However, the user-supplied function (to be modified) does not
    9496    // have those parameters removed.  Here will create a new parameter list
    … …  
    112114
    113115    for (i=0;i<domain->numRows;i++) {
    114         //GUS        tmpf = evalModel(inputParameterList, domain->data.F32[i]);
     116        for (j=0;j<tmpVecPtr->n;j++) {
     117            tmpVecPtr->data.F32[j] = domain->data.F32[i][j];
     118        }
     119        tmpf = evalModel(inputParameterList, tmpVecPtr);
    115120        gsl_vector_set(f, i, (tmpf - data->data.F32[i])/errors->data.F32[i]);
    116121    }
    117122
     123    psFree(inputParameterList);
     124    psFree(tmpVecPtr);
    118125    return GSL_SUCCESS;
    119126}
    … …  
    123130                    gsl_matrix *J)
    124131{
    125     //    size_t n  = ((psModelData *)params)->n;
    126     psImage  *domain   = ((psModelData *)params)->domain;
    127     //    psVector *data   = ((psModelData *)params)->data;
    128     psVector *errors = ((psModelData *) params)->errors;
    129     psVector *mask   = ((psModelData *) params)->paramMask;
    130     psVector *initialGuess = ((psModelData *)params)->initialGuess;
     132    const psImage *restrict domain   = ((psModelData *)params)->domain;
     133    const psVector *restrict errors   = ((psModelData *) params)->errors;
     134    const psVector *restrict mask     = ((psModelData *) params)->paramMask;
     135    psVector *restrict initialGuess = ((psModelData *)params)->initialGuess;
    131136    psVector *inputParameterList = NULL;
    132     //    int count = 0;
    133     //    float (*d_evalModel) = ((psModelData *) params)->evalModel;
     137    psVector *tmpVecPtr = NULL;
     138    float (*d_evalModel)(psVector *, psVector *, int) = ((psModelData *) params)->d_evalModel;
     139
    134140    size_t i;
    135141    int j;
    136142    float tmpf;
     143
     144    tmpVecPtr = psVectorAlloc(domain->numCols, PS_TYPE_F32);
    137145
    138146    // The GSL routines will call this functions with the masked parameters
    … …  
    157165
    158166    for (i=0;i<domain->numRows;i++) {
     167        for (j=0;j<tmpVecPtr->n;j++) {
     168            tmpVecPtr->data.F32[j] = domain->data.F32[i][j];
     169        }
     170
    159171        for (j=0;j<inputParameterList->n;j++) {
    160             //GUS            tmpf = d_evalModel(inputParameterList, domain->data.F32[i], j);
    161 
     172            tmpf = d_evalModel(inputParameterList, tmpVecPtr, j);
    162173            gsl_matrix_set(J, i, j, (tmpf/errors->data.F32[i]));
    163174        }
    164175    }
    165176
     177    psFree(inputParameterList);
     178    psFree(tmpVecPtr);
    166179    return GSL_SUCCESS;
    167180}
    … …  
    185198psVector *
    186199psMinimizeChi2(float (*evalModel)(const psVector *restrict, const psVector *restrict),
    187                float (*DevalModel)(psVector *, psVector *),
    188                const psVector *restrict domain,
     200               float (*DevalModel)(const psVector *restrict, const psVector *restrict, int),
     201               const psImage *restrict domain,
    189202               const psVector *restrict data,
    190203               const psVector *restrict errors,
    … …  
    193206               float *chiSq)
    194207{
    195     size_t n = domain->n;
     208
     209    int numData = domain->numCols; // Number of data points
     210    gsl_multifit_fdfsolver *s = NULL; // GSL data structure.
     211    int status;    // Return status for the GSL solver.
     212    int i = 0;    // Loop index variable.
     213    int j = 0;    // Loop index variable.
     214    int iter = 0;   // Iteration counter.
     215    psModelData *inputData = NULL; // Contains data that the user-supplied
     216    // function/derivate need.
     217    gsl_multifit_function_fdf f; // GSL structure that contains the
     218    // functions/derivative to be solved.
     219    double *xInit = NULL;  // The initial guess at the parameters
     220    // with masked parameters removed.
    196221    const gsl_multifit_fdfsolver_type *T = gsl_multifit_fdfsolver_lmsder;
    197     gsl_multifit_fdfsolver *s      = NULL;
    198     int status;
    199     size_t i = 0;
    200     int j = 0;
    201     size_t iter = 0;
    202     psModelData *inputData = NULL;
    203     gsl_multifit_function_fdf f;
    204     double *xInit = NULL;
     222    // This tells GSL to use the Levenberg-
     223    // Marquardt algorithm for chi2
     224    // minimization.
    205225
    206226    inputData = (psModelData *) psAlloc(sizeof(psModelData));
    207     inputData->n = n;
     227    inputData->n = numData;
    208228    inputData->count = 0;
    209229    inputData->initialGuess = initialGuess;
    210     //    inputData->domain = domain;
     230    inputData->domain = domain;
    211231    //    inputData->data = data;
    212232    //    inputData->errors = errors;
    … …  
    226246    }
    227247
    228 
    229 
    230248    // The initial guess at the parameters for the function are written into
    231249    // the vector inputParameterList.  If the paramMask is not NULL, then those
    232     // parameters are masked out.  How can this possibly work?  The user-
    233     // supplied function will require a fixed number of parameters.
     250    // parameters are masked out.
    234251
    235252    xInit = (double *) psAlloc(inputData->count * sizeof(double));
    … …  
    257274    r = gsl_rng_alloc(type);
    258275
     276    // Initialize the main data structure used by the GSL solver.  This will
     277    // contain pointers to the function to be minimized, it's derivative
     278    // function, the number of data points, the number of free parameters,
     279    // and the data structures those functions use.
     280
    259281    f.f = &gsl_function_f;
    260282    f.df = &gsl_function_df;
    261283    f.fdf = &gsl_my_function_fdf;
    262     f.n = n;
     284    f.n = numData;
    263285    f.p = inputData->count;
    264286    f.params = &inputData;
    265287
    266     // This tells GSL to use the Levenberg-Marquardt algorithm for chi2
    267     // minimization.
    268     //    T = gsl_multifit_fdfsolver_lmsder;
    269 
    270     // Creates an instance of the GSL solver that we will be iterating on.
    271     s = gsl_multifit_fdfsolver_alloc(T, n, inputData->count);
    272 
    273     // Initialize the GSL minimizer.
     288    // Create an instance of the GSL solver that we will be iterating on.
     289    // It will have numData data points and inputData->count parameters.
     290    s = gsl_multifit_fdfsolver_alloc(T, numData, inputData->count);
     291
     292    // Initialize the GSL minimizer to use function defined by the data
     293    // structure "f" and x.vector as an initial guess for the parameters.
    274294    gsl_multifit_fdfsolver_set(s, &f, &x.vector);
    275295
    276     //    print_state(iter, s);
     296
     297    // Each iteration of the following loop will perform one step in an
     298    // attempt to minimized chi-squared for the function.  The loop exits
     299    // either when the change in parameters is small enough, or when the
     300    // maximum number of iterations is reached.
    277301    do {
    278302        iter++;
    … …  
    280304        status = gsl_multifit_fdfsolver_iterate(s);
    281305
    282         //        printf("status = %s\n", gsl_strerror(status));
    283         //        print_state(iter, s);
    284 
    285306        // If there was a problem, abort.
    286307        if (status) {
    … …  
    288309        }
    289310
    290         // Checks whether the (L2-norm) computed derivative and the difference
    291         // between the real/actual for that test x-vector.  If were close
    292         // enough exit loop.
    293 
    294311        // Test if the parameters changed by a small enough amount.
    295312        status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
    296     } while (status == GSL_CONTINUE && iter < 500);
    297 
    298 
    299 
    300 
     313    } while (status == GSL_CONTINUE && iter < MAX_LMM_NUM_ITERATIONS);
     314
     315
     316    // In the above steps we had removed the masked elements from the
     317    // the solver.  This next code blocks puts those masked elements
     318    // into the solution.
    301319    if (paramMask != NULL) {
    302320        j = 0;
    … …  
    314332    }
    315333
     334    // Calculate the chi-squared for the derived solution.
    316335    *chiSq = gsl_blas_dnrm2(s->f);
     336
     337    // Free all allocated memory
     338    // NOTE: Free x.
    317339    gsl_multifit_fdfsolver_free(s);
    318     // Free all allocated memory
    319340    psFree(xInit);
    320341    psFree(inputData);
    321342
     343    // Bye bye.
    322344    return(initialGuess);
    323345}
    324346
     347
     348
     349
     350
     351
     352
     353
     354
     355
     356
     357
     358
     359
     360
     361// NOTE: rewrite so there is no maximum order for the polynomials.
     362#define MAX_POLY_ORDER 10
     363#define MAX_POLYNOMIAL_TERMS  (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2)
     364int MyInfoLevel = 0;
    325365/** @brief This procedure calculates various combinations of powers of x and y
    326366 *   and stores them in the data structure sums[][].  After it completes:

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

    r974 r1093  
    1212 */
    1313
    14 /** This routine must ninimize a particular non-linear function */
    15 psFloatArray *
    16 psMinimize(float (*myFunction)(const psFloatArray *restrict), ///< Function to minimize
    17            psFloatArray *restrict initialGuess, ///< Initial guess
    18            psIntArray *restrict paramMask //!< 1 = fit for parameter, 0 = hold parameter constant
     14/** This routine must minimize a non-linear function */
     15psVector *
     16psMinimize(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
    1919          );
    2020
    2121
    2222/** Minimize chi^2 for input data */
    23 psFloatArray *
    24 psMinimizeChi2(float (*evalModel)(const psFloatArray *restrict,
    25                                   const psFloatArray *restrict), ///< Model to fit; (domain and params)
    26                const psFloatArray *restrict domain, ///< The domain values for the corresponding measurements
    27                const psFloatArray *restrict data, ///< Data to fit
    28                const psFloatArray *restrict errors, ///< Errors in the data
    29                psFloatArray *restrict initialGuess, ///< Initial guess
    30                const psIntArray *restrict paramMask ///< 1 = fit for parameter, 0 = hold parameter constant
     23psVector *
     24psMinimizeChi2(float (*evalModel)(const psVector *restrict,
     25                                  const psVector *restrict), ///< Model to fit; (domain and params)
     26               float (*DevalModel)(const psVector *restrict,
     27                                   const psVector *restrict,
     28                                   int), ///< Derivative of model to fit; (domain and params)
     29               const psImage *restrict domain, ///< The domain values for the corresponding measurements
     30               const psVector *restrict data, ///< Data to fit
     31               const psVector *restrict errors, ///< Errors in the data
     32               psVector *restrict initialGuess, ///< Initial guess
     33               const psVector *restrict paramMask, ///< 1 = fit for parameter, 0 = hold parameter constant
     34               float *chiSq
    3135              );
    3236
    33 /** Derive a polynomial fit by chi^2 minimisation --- can be done analytically */
     37/** Derive a polynomial fit by chi^2 minimisation (analytically) */
    3438psPolynomial1D *
    3539psGetArrayPolynomial(psPolynomial1D myPoly, ///< Polynomial to fit
    36                      const psFloatArray *restrict x, ///< Ordinates (or NULL to just use the indices)
    37                      const psFloatArray *restrict y, ///< Coordinates
    38                      const psFloatArray *restrict yErr ///< Errors in coordinates, or NULL
     40                     const psVector *restrict x, ///< Ordinates (or NULL to just use the indices)
     41                     const psVector *restrict y, ///< Coordinates
     42                     const psVector *restrict yErr ///< Errors in coordinates, or NULL
    3943                    );
    4044
Note: See TracChangeset for help on using the changeset viewer.