IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 1836


Ignore:
Timestamp:
Sep 20, 2004, 1:16:10 PM (22 years ago)
Author:
gusciora
Message:

Removed a lot of obsolete code required by prior SDRs.

Location:
trunk/psLib/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r1831 r1836  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.40 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-09-18 01:50:45 $
     11 *  @version $Revision: 1.41 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-09-20 23:16:10 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    9494/*****************************************************************************/
    9595
    96 typedef struct
    97 {
    98     size_t n;                   // Number of data points points in domain.
    99     int paramCount;             // Number of non-masked parameters.
    100     psVector* restrict initialGuess;
    101     const psImage* restrict domain;
    102     const psVector* restrict data;
    103     const psVector* restrict errors;
    104     const psVector* restrict paramMask;
    105     psMinimizeFunction evalModel;
    106     psMinimizeFunctionDeriv d_evalModel;
    107 }
    108 psMinChi2Data;
    109 
    110 typedef struct
    111 {
    112     int paramCount;             // Number of non-masked parameters.
    113     psVector* restrict initialGuess;
    114     const psVector* restrict coord;
    115     const psVector* restrict paramMask;
    116     psMinimizeFunction evalModel;
    117     psMinimizeFunctionDeriv d_evalModel;
    118 }
    119 psMinimizeData;
    120 
    12196/*****************************************************************************/
    12297/* GLOBAL VARIABLES                                                          */
    … …  
    137112/* FUNCTION IMPLEMENTATION - LOCAL                                           */
    138113/*****************************************************************************/
    139 
    140 /******************************************************************************
    141 p_psMinFunc(*params, *funcData): We use the GSL procedure
    142 gsl_multimin_fdfminimizer_iterate() to minimize an arbitary function supplied
    143 by the user.  That GSL procedure requires the function to be minimized to be
    144 in a different format than the psLib format.  The purpose of this procedure
    145 is to serve as a GSL-format wrapper for the user-supplied procedure which is
    146 to be minimized.
    147  
    148     params: The parameters of the function to be minimized.  These will be
    149  varied by GSL in order to minimize the function.
    150  
    151     funcData: a private psLib struct which contains the data point to be
    152  minimized, the function and derivative function pointers, an initial
    153  guess at the parameters, an option parameter mask, etc.
    154  *****************************************************************************/
    155 double p_psMinFunc(const gsl_vector * params, void *funcData)
    156 {
    157     int i;                      // Loop index variable.
    158     int j;                      // Loop index variable.
    159     float tmpf;                 // Temporary floating point variable.
    160     const psVector* restrict coord = ((psMinimizeData* ) funcData)->coord;
    161     const psVector* restrict mask = ((psMinimizeData* ) funcData)->paramMask;
    162     psVector* restrict initialGuess = ((psMinimizeData* ) funcData)->initialGuess;
    163     psMinimizeFunction evalModel = ((psMinimizeData* ) funcData)->evalModel;
    164     psVector* inputParameterList = NULL;
    165 
    166     // The GSL routines will call this function with the masked parameters
    167     // removed.  However, the user-supplied function (to be modified) does not
    168     // have those parameters removed.  Here will create a new parameter list
    169     // with the masked parameters added (we expand initialGuess).
    170     inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
    171     if (mask != NULL) {
    172         j = 0;
    173         for (i = 0; i < mask->n; i++) {
    174             if (mask->data.U8[i] != 0) {
    175                 inputParameterList->data.F32[i] = initialGuess->data.F32[i];
    176             } else {
    177                 inputParameterList->data.F32[i] = gsl_vector_get(params, j++);
    178             }
    179         }
    180     } else {
    181         for (i = 0; i < initialGuess->n; i++) {
    182             inputParameterList->data.F32[i] = gsl_vector_get(params, i);
    183         }
    184     }
    185 
    186     // Call the user-supplied function.
    187     tmpf = evalModel(inputParameterList, coord);
    188 
    189     // Free allocated memory and return the value of the function.
    190     psFree(inputParameterList);
    191     return (tmpf);
    192 }
    193 
    194 /******************************************************************************
    195 p_psMinFuncDeriv(*params, *funcData):  a GSL-like wrapper for the
    196 user-supplied procedure which calculates the derviative of the function to be
    197 minimized.
    198  
    199     params: The parameters of the function to be minimized.  These will be
    200  varied by GSL in order to minimize the function.
    201  
    202     funcData: a private psLib struct which contains the data point to be
    203  minimized, the function and derivative function pointers, an initial
    204  guess at the parameters, an option parameter mask, etc.
    205  
    206     df: we calculate the derivative of the function w.r.t. to each parameter
    207  in "params" and return those derivatives in this psVector.
    208  *****************************************************************************/
    209 void p_psMinFuncDeriv(const gsl_vector * params, void *funcData, gsl_vector * df)
    210 {
    211     int i;                      // Loop index variable.
    212     int j;                      // Loop index variable.
    213     float tmpf;                 // Temporary floating point variable.
    214     const psVector* restrict coord = ((psMinimizeData* ) funcData)->coord;
    215     const psVector* restrict mask = ((psMinimizeData* ) funcData)->paramMask;
    216     psVector* restrict initialGuess = ((psMinimizeData* ) funcData)->initialGuess;
    217     psMinimizeFunctionDeriv d_evalModel = ((psMinimizeData* ) funcData)->d_evalModel;
    218     psVector* inputParameterList = NULL;
    219 
    220     // The GSL routines will call this function with the masked parameters
    221     // removed.  However, the user-supplied function (to be modified) does not
    222     // have those parameters removed.  Here will create a new parameter list
    223     // with the masked parameters added (we expand initialGuess).
    224     inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
    225     if (mask != NULL) {
    226         j = 0;
    227         for (i = 0; i < mask->n; i++) {
    228             if (mask->data.U8[i] != 0) {
    229                 inputParameterList->data.F32[i] = initialGuess->data.F32[i];
    230             } else {
    231                 inputParameterList->data.F32[i] = gsl_vector_get(params, j++);
    232             }
    233         }
    234     } else {
    235         for (i = 0; i < initialGuess->n; i++) {
    236             inputParameterList->data.F32[i] = gsl_vector_get(params, i);
    237         }
    238     }
    239 
    240     // Evaluate the derivative w.r.t. each parameter.
    241     // NOTE: we can probably remove the calls for masked parameters.
    242     for (i = 0; i < initialGuess->n; i++) {
    243         tmpf = d_evalModel(inputParameterList, coord, i);
    244         gsl_vector_set(df, i, tmpf);
    245     }
    246 
    247     // Free allocated memory.
    248     psFree(inputParameterList);
    249 }
    250 
    251 /******************************************************************************
    252     Compute both p_psMinFunc and p_psMinFuncDeriv together.
    253  *****************************************************************************/
    254 void p_psMinFuncFuncDeriv(const gsl_vector * params, void *funcData, double *f, gsl_vector * df)
    255 {
    256     *f = p_psMinFunc(params, funcData);
    257     p_psMinFuncDeriv(params, funcData, df);
    258 }
    259 
    260 // The first argument to evalModel() and d_evalModel() specifies the data
    261 // point.  It must have the same size as the second dimension of *domain.
    262 // The second argument must have the same size as *initialGuess and
    263 // *paramMask.
    264 
    265 /******************************************************************************
    266 p_psMinChi2Func(*x, *funcData, *outdata): We use the GSL procedure
    267 gsl_multifit_fdfsolver_iterate() to fit an arbitrary function, supplied by
    268 the user, to a set of data points.  That GSL procedure requires the function
    269 to be fit to be in a different format than the psLib format.  The purpose of
    270 this procedure is to serve as a GSL-format wrapper for the user-supplied
    271 procedure which is to be fit to the data.
    272  
    273     params: These are the parameters which are to be varied by GSL in order
    274   to minimize chi2 over the data set.
    275  
    276     funcData: this data structure contains the input values over which the
    277   function will be evaluated, the expected value of the function at
    278   those points, the amount of error tolerable at those points, a mask
    279   vector which specifies which parameters to the function are to be
    280   constant, and an initial guess at the parameters.
    281  
    282     outData: The function is evaluated at each point, then subtract the
    283   expected value and divide by the error.
    284  *****************************************************************************/
    285 int p_psMinChi2Func(const gsl_vector * params, void *funcData, gsl_vector * outData)
    286 {
    287     int i;                      // Loop index variable.
    288     int j;                      // Loop index variable.
    289     float tmpf;                 // Temporary floating point variable.
    290     const psImage* restrict domain = ((psMinChi2Data* ) funcData)->domain;
    291     const psVector* restrict data = ((psMinChi2Data* ) funcData)->data;
    292     const psVector* restrict errors = ((psMinChi2Data* ) funcData)->errors;
    293     const psVector* restrict mask = ((psMinChi2Data* ) funcData)->paramMask;
    294     psVector* restrict initialGuess = ((psMinChi2Data* ) funcData)->initialGuess;
    295     psMinimizeFunction evalModel = ((psMinChi2Data* ) funcData)->evalModel;
    296     psVector* inputParameterList = NULL;
    297     psVector* tmpVecPtr = NULL;
    298 
    299     tmpVecPtr = psVectorAlloc(domain->numCols, PS_TYPE_F32);
    300 
    301     // The GSL routines will call this function with the masked parameters
    302     // removed.  However, the user-supplied function (to be modified) does not
    303     // have those parameters removed.  Here will create a new parameter list
    304     // with the masked parameters added (we expand initialGuess).
    305 
    306     inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
    307     if (mask != NULL) {
    308         j = 0;
    309         for (i = 0; i < mask->n; i++) {
    310             if (mask->data.U8[i] != 0) {
    311                 inputParameterList->data.F32[i] = initialGuess->data.F32[i];
    312             } else {
    313                 inputParameterList->data.F32[i] = gsl_vector_get(params, j++);
    314             }
    315         }
    316     } else {
    317         for (i = 0; i < initialGuess->n; i++) {
    318             inputParameterList->data.F32[i] = gsl_vector_get(params, i);
    319         }
    320     }
    321 
    322     // Evaluate the function at each data point.
    323     for (i = 0; i < domain->numRows; i++) {
    324         for (j = 0; j < domain->numCols; j++) {
    325             tmpVecPtr->data.F32[j] = domain->data.F32[i][j];
    326         }
    327         tmpf = evalModel(tmpVecPtr, inputParameterList);
    328 
    329         gsl_vector_set(outData, i, (tmpf - data->data.F32[i]) / errors->data.F32[i]);
    330     }
    331 
    332     // Free allocated memory.
    333     psFree(inputParameterList);
    334     psFree(tmpVecPtr);
    335 
    336     return GSL_SUCCESS;
    337 }
    338 
    339 /******************************************************************************
    340 p_psMinChi2FuncDeriv(*x, *funcData, *outdata): a GSL-like wrapper for the
    341 user-supplied procedure which calculates the derviative of the function to be
    342 minimized.
    343     params: These are the parameters which are to be varied by GSL in order
    344   to minimize chi2 over the data set.
    345  
    346     funcData: this data structure contains the input values over which the
    347   function will be evaluated, the expected value of the function at
    348   those points, the amount of error tolerable at those points, a mask
    349   vector which specifies which parameters to the function are to be
    350   constant, and an initial guess at the parameters.
    351  
    352     J: The derivative is evaluated at each point and w.r.t. each parameter
    353  and returned in this data structure.
    354  *****************************************************************************/
    355 int p_psMinChi2FuncDeriv(const gsl_vector * params, void *funcData, gsl_matrix * J)
    356 {
    357     const psImage* restrict domain = ((psMinChi2Data* ) funcData)->domain;
    358     const psVector* restrict errors = ((psMinChi2Data* ) funcData)->errors;
    359     const psVector* restrict mask = ((psMinChi2Data* ) funcData)->paramMask;
    360     psVector* restrict initialGuess = ((psMinChi2Data* ) funcData)->initialGuess;
    361     psVector* inputParameterList = NULL;
    362     psVector* tmpVecPtr = NULL;
    363     psMinimizeFunctionDeriv d_evalModel = ((psMinChi2Data* ) funcData)->d_evalModel;
    364 
    365     size_t i;
    366     int j;
    367     float tmpf;
    368 
    369     tmpVecPtr = psVectorAlloc(domain->numCols, PS_TYPE_F32);
    370 
    371     // The GSL routines will call this functions with the masked parameters
    372     // removed.  However, the user-supplied function (to be modified) does not
    373     // have those parameters removed.  Here will create a new parameter list
    374     // with the masked parameters added (we expand initialGuess).
    375 
    376     inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
    377     if (mask != NULL) {
    378         j = 0;
    379         for (i = 0; i < mask->n; i++) {
    380             if (mask->data.U8[i] != 0) {
    381                 inputParameterList->data.F32[i] = initialGuess->data.F32[i];
    382             } else {
    383                 inputParameterList->data.F32[i] = gsl_vector_get(params, j++);
    384             }
    385         }
    386     } else {
    387         for (i = 0; i < initialGuess->n; i++) {
    388             inputParameterList->data.F32[i] = gsl_vector_get(params, i);
    389         }
    390     }
    391 
    392     // Evaluate the derivtaive at each data point, and w.r.t. each parameter.
    393     for (i = 0; i < domain->numRows; i++) {
    394         for (j = 0; j < tmpVecPtr->n; j++) {
    395             tmpVecPtr->data.F32[j] = domain->data.F32[i][j];
    396         }
    397 
    398         for (j = 0; j < inputParameterList->n; j++) {
    399             tmpf = d_evalModel(tmpVecPtr, inputParameterList, j);
    400             gsl_matrix_set(J, i, j, (tmpf / errors->data.F32[i]));
    401         }
    402     }
    403 
    404     psFree(inputParameterList);
    405     psFree(tmpVecPtr);
    406     return GSL_SUCCESS;
    407 }
    408 
    409 int p_psMinChi2FuncFuncDeriv(const gsl_vector * params, void *funcData, gsl_vector * f, gsl_matrix * J)
    410 {
    411     p_psMinChi2Func(params, funcData, f);
    412     p_psMinChi2FuncDeriv(params, funcData, J);
    413 
    414     return GSL_SUCCESS;
    415 }
    416114
    417115/******************************************************************************
    … …  
    433131
    434132/******************************************************************************
    435 p_psBuildSums1D(x): this routine returns a psVector with "x" elements.  The
     133VectorNormalizeGen(): this routine returns a psVector with "x" elements.  The
    436134values of the vector will be scaled uniformly between -1.0 and 1.0.
    437  *****************************************************************************/
    438 psVector* psBuildImageScalingFactors(int x)
     135*****************************************************************************/
     136psVector* VectorNormalizeGen(int x)
    439137{
    440     int i = 0;                  // loop index variable.
    441     psVector* imageScalingFactors = NULL;
    442 
    443     imageScalingFactors = psVectorAlloc(x, PS_TYPE_F32);
     138    int i = 0;
     139    psVector *tmp = NULL;
     140
     141    tmp = psVectorAlloc(x, PS_TYPE_F32);
    444142
    445143    for (i = 0; i < x; i++) {
    446         imageScalingFactors->data.F32[i] = (((float)2 * i) / ((float)x)) - 1.0;
    447     }
    448 
    449     return (imageScalingFactors);
    450 }
     144        tmp->data.F32[i] = (((float)2 * i) / ((float)x)) - 1.0;
     145    }
     146
     147    return (tmp);
     148}
     149
     150/*****************************************************************************/
     151/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
     152/*****************************************************************************/
    451153
    452154/******************************************************************************
    453 CURRENTLY NOT IN USE.
    454  
    455 p_psPolyOrderCheck(A, N, *indx, *B, polyOrder,*flag) This routine checks if
    456 all polyOrder-th terms in the polyOrder-th order sky background polynomial
    457 defined by the coefficients in the array B[] are consistent with zero.  If
    458 true, then *flag is set to 1.  Otherwise, *flag is set to 0.  The matrix
    459 inversion code in the middle of this procedure draws from Numerical Recipes
    460 in C page 48.
    461 Input:
    462     A       This is the LUD decomposition of the original matrix A.
    463     N       The size of the matrix (plus 1, actually, since offset 1).
    464     indx    misc Numerical Recipes data structure.
    465     B       The coefficients of the sky polynomial.
    466     polyOrder The degree of the sky polynomial.
    467 Output:
    468     *flag   Set this to 1 if we must recalculate the coefficients.
    469  *****************************************************************************/
    470 void p_psPolyOrderCheck(float **A, int N, int *indx, float *B, int polyOrder, int *flag)
    471 {
    472     float **y = NULL;           // This 2-D matrix will hold A^-1
    473     float *col = NULL;          // misc NumerRecipes data structure
    474     float *error = NULL;        // will hold the sqrt() of the
    475 
    476     // diagonal of y[][].
    477     int i = 0;                  // loop-index variable
    478     int j = 0;                  // loop-index variable
    479     int numPolyTerms = 0;       // The number of terms in the
    480 
    481     // polynomial.
    482     int lastTerm = 0;           // The index location of the first
    483 
    484     // n-th order term in array B[].
    485     int firstTerm = 0;          // Index location of last such term.
    486 
    487     // Allocate the necessary data structures for this procedure...
    488     error = (float *)psAlloc((N + 1) * sizeof(float));
    489     col = (float *)psAlloc((N + 1) * sizeof(float));
    490     y = (float **)psAlloc((N + 1) * sizeof(float *));
    491     for (i = 1; i <= N; i++) {
    492         y[i] = (float *)psAlloc((N + 1) * sizeof(float));
    493     }
    494 
    495     // Invert the matrix A and put the result in y[][].  This code is taken
    496     // from Numerical Recipes in C page 48.
    497     for (j = 1; j <= N; j++) {
    498         for (i = 1; i <= N; i++) {
    499             col[i] = 0.0;
    500         }
    501         col[j] = 1.0;
    502         // NOTE: substitue the LUD rotine
    503         // lubksb(A, N, indx, col);
    504         for (i = 1; i <= N; i++) {
    505             y[i][j] = col[i];
    506         }
    507     }
    508 
    509     // Determine where the first n-th order (in this comment, n equals
    510     // polyOrder) polynomial term is stored in the matrix B[], and also were
    511     // the last n-order term is stored.  Then we loop over all the n-order
    512     // terms and check if they are consistent with zero.
    513 
    514     numPolyTerms = (((polyOrder + 1) * (polyOrder + 2)) / 2);
    515     lastTerm = numPolyTerms + 1;
    516     firstTerm = lastTerm - polyOrder;
    517     *flag = 1;
    518     for (i = firstTerm; i <= lastTerm; i++) {
    519         #ifdef DARWIN
    520         error[i] = (float)sqrt(y[i][i]);
    521         #else
    522 
    523         error[i] = sqrtf(y[i][i]);
    524         #endif
    525 
    526         if (!((B[i] <= (2.0f * error[i])) && ((-2.0f * error[i]) <= B[i]))) {
    527             *flag = 0;
    528         }
    529     }
    530 
    531     // Free all memory allocated in this routine.
    532     psFree(error);
    533     psFree(col);
    534     for (j = 1; j <= N; j++) {
    535         psFree(y[j]);
    536     }
    537     psFree(y);
    538 }
    539 
    540 /*****************************************************************************/
    541 /* FUNCTION IMPLEMENTATION - PUBLIC                                          */
    542 /*****************************************************************************/
    543 
    544 /******************************************************************************
    545 psMinimize(initialGuess, myFunction, myFunctionDeriv, coord, paramMask):
    546  
    547 This routine must minimize an arbitrary function; it determines the set of
    548 parameters of that function such that the ...
    549  *****************************************************************************/
    550 psVector* psMinimize(psVector* restrict initialGuess,
    551                      psMinimizeFunction myFunction,
    552                      psMinimizeFunctionDeriv myFunctionDeriv,
    553                      const psVector* restrict coord,
    554                      const psVector* restrict paramMask)
    555 {
    556     int status;
    557     int i = 0;
    558     int j = 0;
    559     int iter = 0;
    560     gsl_multimin_function_fdf f;
    561     const gsl_multimin_fdfminimizer_type *T;
    562     gsl_multimin_fdfminimizer *s;
    563     psMinimizeData inputData;
    564     gsl_vector *x;
    565 
    566     psTrace(".psLib.dataManip.psMinimize", 4,
    567             "---- psMinimize() begin ----\n");
    568 
    569     PS_CHECK_NULL_VECTOR(initialGuess);
    570     PS_CHECK_EMPTY_VECTOR(initialGuess);
    571     PS_CHECK_NULL_VECTOR(coord);
    572     PS_CHECK_EMPTY_VECTOR(coord);
    573     if (paramMask != NULL) {
    574         PS_CHECK_NULL_VECTOR(paramMask);
    575         PS_CHECK_EMPTY_VECTOR(paramMask);
    576         PS_CHECK_VECTOR_SIZE_EQUAL(initialGuess, paramMask);
    577     }
    578 
    579     inputData.initialGuess = initialGuess;
    580     inputData.coord = coord;
    581     inputData.paramMask = paramMask;
    582     inputData.evalModel = myFunction;
    583     inputData.d_evalModel = myFunctionDeriv;
    584     inputData.paramCount = 0;
    585 
    586     // If the user supplied a parameter mask, then count the number of
    587     // non-masked elements.  This will be used later in allocating a vector
    588     // for the parameters.
    589     if (paramMask != NULL) {
    590         for (i = 0; i < paramMask->n; i++) {
    591             if (paramMask->data.U8[i] != 0) {
    592                 inputData.paramCount++;
    593             }
    594         }
    595     } else {
    596         inputData.paramCount = initialGuess->n;
    597     }
    598 
    599     // The initial guess at the parameters for the function are written into
    600     // the vector inputParameterList.  If the paramMask is not NULL, then
    601     // masked parameters are masked out.
    602     x = gsl_vector_alloc(inputData.paramCount);
    603     if (paramMask != NULL) {
    604         j = 0;
    605         for (i = 0; i < initialGuess->n; i++) {
    606             if (paramMask->data.U8[i] == 0) {
    607                 gsl_vector_set(x, j++, initialGuess->data.F32[i]);
    608             }
    609         }
    610     } else {
    611         for (i = 0; i < initialGuess->n; i++) {
    612             gsl_vector_set(x, i, initialGuess->data.F32[i]);
    613         }
    614     }
    615     f.f = &p_psMinFunc;
    616     f.df = &p_psMinFuncDeriv;
    617     f.fdf = &p_psMinFuncFuncDeriv;
    618     f.n = inputData.paramCount;
    619     f.params = &inputData;
    620 
    621     T = gsl_multimin_fdfminimizer_conjugate_fr;
    622     s = gsl_multimin_fdfminimizer_alloc(T, inputData.paramCount);
    623     gsl_multimin_fdfminimizer_set(s, &f, x, 0.01, 1e-4);
    624     do {
    625         iter++;
    626         status = gsl_multimin_fdfminimizer_iterate(s);
    627 
    628         if (status)
    629             break;
    630 
    631         status = gsl_multimin_test_gradient(s->gradient, 1e-3);
    632 
    633         if (status == GSL_SUCCESS)
    634             printf("Minimum found at:\n");
    635 
    636     } while (status == GSL_CONTINUE && iter < MAX_MINIMIZE_ITERATIONS);
    637 
    638     // In the above steps we had removed the masked elements from the
    639     // the solver.  This next code blocks puts those masked elements
    640     // into the solution.
    641     if (paramMask != NULL) {
    642         j = 0;
    643         for (i = 0; i < initialGuess->n; i++) {
    644             if (paramMask->data.U8[i] == 0) {
    645                 initialGuess->data.F32[i] = gsl_vector_get(s->x, j++);
    646             } else {
    647                 initialGuess->data.F32[i] = initialGuess->data.F32[i];
    648             }
    649         }
    650     } else {
    651         for (i = 0; i < initialGuess->n; i++) {
    652             initialGuess->data.F32[i] = gsl_vector_get(s->x, i);
    653         }
    654     }
    655     psTrace(".psLib.dataManip.psMinimize", 4,
    656             "---- psMinimize() end ----\n");
    657     return (initialGuess);
    658 }
    659 
    660 /******************************************************************************
    661     This routine must determine the parameters of an arbitrary function
    662     such that they best fit the supplied data points.
    663  *****************************************************************************/
    664 psVector* psMinimizeChi2(psMinimizeFunction evalModel,
    665                          psMinimizeFunctionDeriv DevalModel,
    666                          const psImage* restrict domain,
    667                          const psVector* restrict data,
    668                          const psVector* restrict errors,
    669                          psVector* restrict initialGuess,
    670                          const psVector* restrict paramMask,
    671                          float *chiSq)
    672 {
    673     int numData = domain->numRows;      // Number of data points
    674     int status;                 // Return status for the GSL solver.
    675     int i = 0;                  // Loop index variable.
    676     int j = 0;                  // Loop index variable.
    677     int iter = 0;               // Iteration counter.
    678     gsl_multifit_function_fdf f;        // GSL structure that contains the
    679 
    680     // functions/derivative to be solved.
    681     double *xInit = NULL;       // The initial guess at the parameters
    682 
    683     // with masked parameters removed.
    684     const gsl_multifit_fdfsolver_type *T;
    685 
    686     // This tells GSL to use the Levenberg-
    687     // Marquardt algorithm for chi2
    688     // minimization.
    689     gsl_multifit_fdfsolver *s;  // GSL data structure.
    690     psMinChi2Data inputData;
    691     float chiSqOld = 0.0;
    692 
    693     psTrace(".psLib.dataManip.psMinimizeChi2", 4,
    694             "---- psMinimizeChi2() begin ----\n");
    695 
    696     PS_CHECK_NULL_IMAGE(domain);
    697     PS_CHECK_EMPTY_IMAGE(domain);
    698     PS_CHECK_NULL_VECTOR(data);
    699     PS_CHECK_EMPTY_VECTOR(data);
    700     PS_CHECK_NULL_VECTOR(errors);
    701     PS_CHECK_EMPTY_VECTOR(errors);
    702     PS_CHECK_NULL_VECTOR(initialGuess);
    703     PS_CHECK_EMPTY_VECTOR(initialGuess);
    704     PS_CHECK_VECTOR_SIZE_EQUAL(data, errors);
    705     if (domain->numRows != data->n) {
    706         psAbort(__func__, "Number of data points and data values not equal.");
    707     }
    708     if (paramMask != NULL) {
    709         PS_CHECK_NULL_VECTOR(paramMask);
    710         PS_CHECK_EMPTY_VECTOR(paramMask);
    711         PS_CHECK_VECTOR_SIZE_EQUAL(initialGuess, paramMask);
    712     }
    713 
    714     inputData.n = numData;
    715     inputData.paramCount = 0;
    716     inputData.initialGuess = initialGuess;
    717     inputData.domain = domain;
    718     inputData.data = data;
    719     inputData.errors = errors;
    720     inputData.paramMask = paramMask;
    721     inputData.evalModel = evalModel;
    722     inputData.d_evalModel = DevalModel;
    723 
    724     // If the user supplied a parameter mask, then count the number of
    725     // non-masked elements.  This will be used later in allocating a vector
    726     // for the parameters.
    727     if (paramMask != NULL) {
    728         for (i = 0; i < paramMask->n; i++) {
    729             if (paramMask->data.U8[i] != 0) {
    730                 inputData.paramCount++;
    731             }
    732         }
    733     } else {
    734         inputData.paramCount = initialGuess->n;
    735     }
    736 
    737     // The initial guess at the parameters for the function are written into
    738     // the vector inputParameterList.  If the paramMask is not NULL, then those
    739     // parameters are masked out.
    740     xInit = (double *)psAlloc(inputData.paramCount * sizeof(double));
    741     if (paramMask != NULL) {
    742         j = 0;
    743         for (i = 0; i < initialGuess->n; i++) {
    744             if (paramMask->data.U8[i] == 0) {
    745                 xInit[j++] = initialGuess->data.F32[i];
    746             }
    747         }
    748     } else {
    749         for (i = 0; i < initialGuess->n; i++) {
    750             xInit[i] = initialGuess->data.F32[i];
    751         }
    752     }
    753 
    754     const gsl_rng_type *type;
    755     gsl_rng *r;
    756 
    757     gsl_rng_env_setup();
    758 
    759     type = gsl_rng_default;
    760     r = gsl_rng_alloc(type);
    761 
    762     // Initialize the main data structure used by the GSL solver.  This will
    763     // contain pointers to the function to be minimized, it's derivative
    764     // function, the number of data points, the number of free parameters,
    765     // and the data structures those functions use.
    766 
    767     f.f = &p_psMinChi2Func;
    768     f.df = &p_psMinChi2FuncDeriv;
    769     f.fdf = &p_psMinChi2FuncFuncDeriv;
    770     f.n = numData;
    771     f.p = inputData.paramCount;
    772     f.params = &inputData;
    773 
    774     gsl_vector_view x = gsl_vector_view_array(xInit, inputData.paramCount);
    775 
    776     T = gsl_multifit_fdfsolver_lmsder;
    777     s = gsl_multifit_fdfsolver_alloc(T, numData, inputData.paramCount);
    778     gsl_multifit_fdfsolver_set(s, &f, &x.vector);
    779     *chiSq = 0.0;
    780     chiSqOld = 0.0;
    781     do {
    782         iter++;
    783         for (i = 0; i < initialGuess->n; i++) {
    784             printf("Iteration %d: parameter %d is %.3f\n", iter, i, gsl_vector_get(s->x, i));
    785         }
    786         // Perform an iteration of the GSL solver.
    787         status = gsl_multifit_fdfsolver_iterate(s);
    788         printf("gsl_multifit_fdfsolver_iterate() status is %s\n", gsl_strerror(status));
    789         for (i = 0; i < initialGuess->n; i++) {
    790             printf("Iteration %d: parameter %d is %.3f\n", iter, i, gsl_vector_get(s->x, i));
    791         }
    792 
    793         // If there was a problem, abort.
    794         if (status) {
    795             psAbort(__func__, "gsl_multifit_fdfsolver_iterate(%s)\n", gsl_strerror(status));
    796         }
    797         // Test if the parameters changed by a small enough amount.
    798         // NOTE: This wasn't working right when the parameters fit exactly.
    799         // Figure out why.
    800         // status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
    801 
    802         // We test for convergence if chiSquared changes by less than 1.0
    803         // as specified in the ADD.
    804         *chiSq = gsl_blas_dnrm2(s->f);
    805         printf("psMinimize.c: chiSq is %.3f\n", *chiSq);
    806         if (fabs(*chiSq - chiSqOld) < 1.0) {
    807             status = GSL_SUCCESS;
    808         } else {
    809             status = GSL_CONTINUE;
    810         }
    811         chiSqOld = *chiSq;
    812 
    813     } while (status == GSL_CONTINUE && iter < MAX_LMM_ITERATIONS);
    814 
    815     // In the above steps we had removed the masked elements from the
    816     // the solver.  This next code blocks puts those masked elements
    817     // into the solution.
    818     if (paramMask != NULL) {
    819         j = 0;
    820         for (i = 0; i < initialGuess->n; i++) {
    821             if (paramMask->data.U8[i] == 0) {
    822                 initialGuess->data.F32[i] = gsl_vector_get(s->x, j++);
    823             } else {
    824                 initialGuess->data.F32[i] = initialGuess->data.F32[i];
    825             }
    826         }
    827     } else {
    828         for (i = 0; i < initialGuess->n; i++) {
    829             initialGuess->data.F32[i] = gsl_vector_get(s->x, i);
    830         }
    831     }
    832 
    833     // Calculate the chi-squared for the derived solution.
    834     *chiSq = gsl_blas_dnrm2(s->f);
    835 
    836     // Free all allocated memory
    837     // NOTE: Free x.
    838     gsl_multifit_fdfsolver_free(s);
    839     psFree(xInit);
    840 
    841     // Bye bye.
    842     psTrace(".psLib.dataManip.psMinimizeChi2", 4,
    843             "---- psMinimizeChi2() end ----\n");
    844     return (initialGuess);
    845 }
    846 
    847 /******************************************************************************
    848 This routine will take an procedure which calculates an arbitrary function
    849 and it's derivative and minimize it.
    850  
    851 XXX: Do this:
    852  After checking that all entries in the paramMask are 1 or 0, when
    853  forming the A matrix from alpha, try this:
    854  
    855      A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
    856  *****************************************************************************/
    857 bool psMinimizeLM(psMinimization *min,
    858                   psImage *covar,
    859                   psVector *params,
    860                   const psVector *paramMask,
    861                   const psArray *coords,
    862                   psMinimizeLMFunc func)
    863 {
    864     psVector *beta = psVectorAlloc(params->n, PS_TYPE_F64);
    865     psVector *perm = psVectorAlloc(params->n, PS_TYPE_F64);
    866     psVector *newParamsF64 = psVectorAlloc(params->n, PS_TYPE_F64);
    867     psVector *newParamsF32 = psVectorAlloc(params->n, PS_TYPE_F32);
    868     psVector *origParams = psVectorAlloc(params->n, PS_TYPE_F32);
    869     psImage *alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32);
    870     psImage *A = psImageAlloc(params->n, params->n, PS_TYPE_F64);
    871     psImage *aOut = psImageAlloc(params->n, params->n, PS_TYPE_F64);
    872     psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32);
    873     psVector *newDeriv = psVectorAlloc(params->n, PS_TYPE_F32);
    874     psVector *NRparams = psVectorAlloc(params->n, PS_TYPE_F32);
    875     int i;
    876     int j;
    877     int k;
    878     float newValue;
    879     float oldValue;
    880     float lamda = 0.00005;
    881 
    882     psTrace(".psLib.dataManip.psMinimizeLM", 4,
    883             "---- psMinimizeLM() begin ----\n");
    884     psTrace(".psLib.dataManip.psMinimizeLM", 6,
    885             "min->maxIter is %d\n", min->maxIter);
    886     psTrace(".psLib.dataManip.psMinimizeLM", 6,
    887             "min->tol is %f\n", min->tol);
    888 
    889     for (i=0;i<params->n;i++) {
    890         origParams->data.F32[i] = params->data.F32[i];
    891     }
    892 
    893     min->lastDelta = HUGE;
    894     min->lastDelta = 12345.0;
    895     min->iter = 0;
    896 
    897     while ((min->lastDelta > min->tol) &&
    898             (min->iter < min->maxIter)) {
    899         psTrace(".psLib.dataManip.psMinimizeLM", 4,
    900                 "------------------------------------------------------\n");
    901         psTrace(".psLib.dataManip.psMinimizeLM", 4,
    902                 "Iteration %d.  Delta is %f\n", min->iter, min->lastDelta);
    903 
    904         min->value = func(deriv, params, coords);
    905 
    906         for (i=0;i<params->n;i++) {
    907             psTrace(".psLib.dataManip.psMinimizeLM", 6,
    908                     "params->data.F32[%d] is %f.\n", i, params->data.F32[i]);
    909         }
    910         for (i=0;i<params->n;i++) {
    911             psTrace(".psLib.dataManip.psMinimizeLM", 6,
    912                     "deriv->data.F32[%d] is %f.\n", i, deriv->data.F32[i]);
    913         }
    914         psTrace(".psLib.dataManip.psMinimizeLM", 6,
    915                 "min->value is (%f)\n", min->value);
    916 
    917         for (i=0;i<params->n;i++) {
    918             if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
    919                 deriv->data.F32[i] = 0.0;
    920             }
    921         }
    922 
    923         // Calculate the BETA vector.
    924         for (i=0;i<params->n;i++) {
    925             beta->data.F64[i] = (float) deriv->data.F32[i];
    926             psTrace(".psLib.dataManip.psMinimizeLM", 6,
    927                     "beta->data.F64[%d] is %f.\n", i, beta->data.F64[i]);
    928         }
    929 
    930         // Calculate the ALPHA matrix.
    931         for (j=0;j<params->n;j++) {
    932             for (k=0;k<params->n;k++) {
    933                 alpha->data.F32[j][k] = deriv->data.F32[k] *
    934                                         deriv->data.F32[j];
    935             }
    936         }
    937 
    938         // Calculate the matrix A.
    939         for (j=0;j<params->n;j++) {
    940             for (k=0;k<params->n;k++) {
    941                 if (j == k) {
    942                     A->data.F64[j][k] =
    943                         (double) ((1.0 + lamda) * alpha->data.F32[j][k]);
    944                 } else {
    945                     A->data.F64[j][k] = (double) alpha->data.F32[j][k];
    946                 }
    947             }
    948         }
    949         for (j=0;j<params->n;j++) {
    950             psTrace(".psLib.dataManip.psMinimizeLM", 6, "Matrix A[][]:\n");
    951             for (k=0;k<params->n;k++) {
    952                 psTrace(".psLib.dataManip.psMinimizeLM", 6, "%f ", A->data.F64[j][k]);
    953             }
    954             psTrace(".psLib.dataManip.psMinimizeLM", 6, "Matrix A[][]:\n");
    955         }
    956 
    957         // Solve A * alpha = Beta
    958         aOut = psMatrixLUD(aOut, perm, A);
    959         newParamsF64 = psMatrixLUSolve(newParamsF64, aOut, beta, perm);
    960 
    961         // Mask any masked parameters.
    962         for (i=0;i<params->n;i++) {
    963             if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
    964                 newParamsF64->data.F64[i] = (double) origParams->data.F32[i];
    965             }
    966             newParamsF32->data.F32[i] = (float) newParamsF64->data.F64[i];
    967 
    968             psTrace(".psLib.dataManip.psMinimizeLM", 6,
    969                     "newParamsF32->data.F32[%d] is %f.\n", i, newParamsF32->data.F32[i]);
    970             NRparams->data.F32[i] = params->data.F32[i] - newParamsF32->data.F32[i];
    971         }
    972 
    973         psTrace(".psLib.dataManip.psMinimizeLM", 6,
    974                 "Calling func() with new parameters:\n");
    975         for (i=0;i<params->n;i++) {
    976             psTrace(".psLib.dataManip.psMinimizeLM", 6,
    977                     "NRparams->data.F32[%d] is %f.\n", i, NRparams->data.F32[i]);
    978         }
    979 
    980         oldValue = min->value;
    981         newValue = func(deriv, NRparams, coords);
    982         psTrace(".psLib.dataManip.psMinimizeLM", 4,
    983                 "old/new values are (%f, %f)\n", oldValue, newValue);
    984         for (i=0;i<params->n;i++) {
    985             if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
    986                 deriv->data.F32[i] = 0.0;
    987             }
    988         }
    989 
    990         if (oldValue > newValue) {
    991             min->lastDelta = oldValue - newValue;
    992             min->value = newValue;
    993 
    994             // No need to check the paramMask here since we already did so
    995             // before the last function evaluation.
    996             for (i=0;i<params->n;i++) {
    997                 //                params->data.F32[i] = (float) newParamsF64->data.F64[i];
    998                 params->data.F32[i] = (float) NRparams->data.F32[i];
    999             }
    1000             min->value = func(deriv, params, coords);
    1001 
    1002             lamda*= 0.1;
    1003         } else {
    1004             lamda*= 10.0;
    1005         }
    1006         psTrace(".psLib.dataManip.psMinimizeLM", 4,
    1007                 "lamda is %f\n", lamda);
    1008         min->iter++;
    1009     }
    1010     psFree(beta);
    1011     psFree(perm);
    1012     psFree(newParamsF64);
    1013     psFree(newParamsF32);
    1014     psFree(origParams);
    1015     psFree(alpha);
    1016     psFree(A);
    1017     psFree(aOut);
    1018     psFree(deriv);
    1019     psFree(newDeriv);
    1020 
    1021     if ((min->iter < min->maxIter) ||
    1022             (min->lastDelta <= min->tol)) {
    1023         return(true);
    1024     }
    1025 
    1026     psTrace(".psLib.dataManip.psMinimizeLM", 4,
    1027             "---- psMinimizeLM() end (false) ----\n");
    1028     return(false);
    1029 }
    1030 
    1031 /******************************************************************************
    1032 This routine will take an procedure which calculates an arbitrary function
    1033 and it's derivative and minimize the chi-squared match between that function
    1034 at the specified coords and the specified value at those coords.
     155psMinimizeLMChi2():  This routine will take an procedure which calculates an
     156arbitrary function and it's derivative and minimize the chi-squared match
     157between that function at the specified coords and the specified value at
     158those coords.
    1035159 
    1036160XXX: Do this:
    … …  
    1265389
    1266390/******************************************************************************
    1267     This routine must fit a polynomial of degree myPoly to the data points
    1268     (x, y) and return the coefficients of that polynomial, as well as the
    1269     error for each data point (yErr).
     391psVectorFitPolynomial1D():  This routine must fit a polynomial of degree
     392myPoly to the data points (x, y) and return the coefficients of that
     393polynomial, as well as the error for each data point (yErr).
    1270394 
    1271395XXX: NOTE: yErr is currently ignored.
    … …  
    1273397psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* myPoly,
    1274398                                        const psVector* restrict x,
    1275                                         const psVector* restrict y, const psVector* restrict yErr)
     399                                        const psVector* restrict y,
     400                                        const psVector* restrict yErr)
    1276401{
    1277402    int polyOrder = myPoly->n;
    … …  
    1381506
    1382507#define STEP_SIZE 0.10
    1383 /******************************************************************************
    1384     This routine takes as input an arbitrary function, and the parameter to
    1385     vary.  This function produces as output a bracket [a, b, c] such that
    1386     f(b) is less than f(a) and f(b).
    1387  
    1388     Algorithm: XXX completely ad hoc: start with the user-supplied starting
    1389     parameter and call that b.  Calculate a/c as a fractional amount
    1390     smaller/larger than b.  Repeat this process until a local minimum is
    1391     found.
    1392  *****************************************************************************/
    1393 psVector *p_psDetermineBracketOld(psVector *params,
    1394                                   int dim,
    1395                                   const psArray *coords,
    1396                                   psMinimizePowellFunc func)
    1397 {
    1398     float a = 0.0;
    1399     float b = 0.0;
    1400     float c = 0.0;
    1401     float fa = 0.0;
    1402     float fb = 0.0;
    1403     float fc = 0.0;
    1404     int iter = 100;
    1405     float aDir = 0.0;
    1406     float cDir = 0.0;
    1407     float new_aDir = 0.0;
    1408     float new_cDir = 0.0;
    1409     psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
    1410     float stepSize = params->data.F32[dim] * STEP_SIZE;
    1411     //    float initialParam = params->data.F32[dim];
    1412 
    1413     if (stepSize == 0.0) {
    1414         stepSize = 1.0;
    1415     }
    1416     a = b = c = params->data.F32[0];
    1417     a-= stepSize;
    1418     c+= stepSize;
    1419 
    1420     params->data.F32[dim] = a;
    1421     fa = func(params, coords);
    1422     params->data.F32[dim] = b;
    1423     fb = func(params, coords);
    1424     params->data.F32[dim] = c;
    1425     fc = func(params, coords);
    1426     if (fa < fb) {
    1427         aDir = -1;
    1428     } else {
    1429         aDir = 1;
    1430     }
    1431 
    1432     if (fc < fb) {
    1433         cDir = -1;
    1434     } else {
    1435         cDir = 1;
    1436     }
    1437 
    1438     while (iter > 0) {
    1439         if ((b < a) && (b < c)) {
    1440             bracket->data.F32[0] = a;
    1441             bracket->data.F32[1] = b;
    1442             bracket->data.F32[2] = c;
    1443             return(bracket);
    1444         }
    1445         stepSize*= (1.0 + stepSize);
    1446         a = a - stepSize;
    1447         c = c + stepSize;
    1448 
    1449         params->data.F32[dim] = a;
    1450         fa = func(params, coords);
    1451         params->data.F32[dim] = c;
    1452         fc = func(params, coords);
    1453 
    1454         //printf("HMMM(%d): (%f %f %f) (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
    1455 
    1456         if (fa < fb) {
    1457             new_aDir = -1;
    1458         } else {
    1459             new_aDir = 1;
    1460         }
    1461 
    1462         if (fc < fb) {
    1463             new_cDir = -1;
    1464         } else {
    1465             new_cDir = 1;
    1466         }
    1467         if ((new_aDir == 1) && (aDir == -1)) {
    1468             bracket->data.F32[0] = a;
    1469             bracket->data.F32[1] = b;
    1470             bracket->data.F32[2] = c;
    1471             return(bracket);
    1472         }
    1473 
    1474         if ((new_cDir == 1) && (cDir == -1)) {
    1475             bracket->data.F32[0] = a;
    1476             bracket->data.F32[1] = b;
    1477             bracket->data.F32[2] = c;
    1478             return(bracket);
    1479         }
    1480         aDir = new_aDir;
    1481         cDir = new_cDir;
    1482         iter--;
    1483     }
    1484     psFree(bracket);
    1485     return(NULL);
    1486 }
    1487 
    1488508/******************************************************************************
    1489509    This routine takes as input an arbitrary function, and the parameter to
    … …  
    1666686}
    1667687
    1668 
    1669 /******************************************************************************
    1670     This routine must minimize a possibly multi-dimensional function
    1671     (several parameters) along a single dimension.
    1672  *****************************************************************************/
    1673 bool psMinimize1DFunc(psMinimization *min,
    1674                       psVector *params,
    1675                       int dim,
    1676                       const psArray *coords,
    1677                       psMinimizePowellFunc func)
    1678 {
    1679     psVector *bracket;
    1680     float a = 0.0;
    1681     float b = 0.0;
    1682     float c = 0.0;
    1683     float n = 0.0;
    1684     float fa = 0.0;
    1685     float fb = 0.0;
    1686     float fc = 0.0;
    1687     float fn = 0.0;
    1688     //    float initialParam = params->data.F32[dim];
    1689 
    1690     bracket = p_psDetermineBracketOld(params, dim, coords, func);
    1691     if (bracket == NULL) {
    1692         psAbort(__func__, "Could not bracket minimum.");
    1693     }
    1694 
    1695     min->iter = 0;
    1696     while (min->iter < min->maxIter) {
    1697         min->iter++;
    1698         //printf("psMinimize1DFunc(): iteration %d\n", min->iter);
    1699         a = bracket->data.F32[0];
    1700         b = bracket->data.F32[1];
    1701         c = bracket->data.F32[2];
    1702 
    1703         params->data.F32[dim] = a;
    1704         fa = func(params, coords);
    1705         params->data.F32[dim] = b;
    1706         fb = func(params, coords);
    1707         params->data.F32[dim] = c;
    1708         fc = func(params, coords);
    1709         //printf("Iteration %d: f(%f %f %f) is (%f %f %f)\n", min->iter, a, b, c, fa, fb, fc);
    1710 
    1711         // We determine which is the biggest segment in [a,b,c] then split
    1712         // that with the point n.
    1713         if ((b-a) > (c-b)) {
    1714             // This is the golden section formula
    1715             params->data.F32[dim] = n = a + (0.69 * (b-a));
    1716             fn = func(params, coords);
    1717             if (fn > fb) {
    1718                 // a = n, b = b, c = c
    1719                 bracket->data.F32[0] = n;
    1720             } else {
    1721                 // a = a, b = n, c = b
    1722                 bracket->data.F32[1] = n;
    1723                 bracket->data.F32[2] = b;
    1724             }
    1725         } else {
    1726             params->data.F32[dim] = n = b + (0.69 * (c-b));
    1727             fn = func(params, coords);
    1728             if (fn > fb) {
    1729                 // a = a, b = b, c = n
    1730                 bracket->data.F32[2] = n;
    1731             } else {
    1732                 // a = b, b = n, c = c
    1733                 bracket->data.F32[0] = b;
    1734                 bracket->data.F32[1] = n;
    1735             }
    1736         }
    1737 
    1738         if ((fabs(a-b) < min->tol) &&
    1739                 (fabs(b-c) < min->tol)) {
    1740             //            psFree(bracket);
    1741             //  XXX: is this line correct?
    1742             params->data.F32[dim] = bracket->data.F32[1];
    1743             min->value = func(params, coords);
    1744             return(true);
    1745         }
    1746     }
    1747 
    1748     //    psFree(bracket);
    1749     return(false);
    1750 }
    1751688
    1752689/******************************************************************************
    … …  
    21681105    return(psMinimizePowell(min, params, paramMask, coords, myPowellChi2Func));
    21691106}
     1107
     1108

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

    r1831 r1836  
    99 *  @author George Gusciora, MHPCC
    1010 *
    11  *  @version $Revision: 1.40 $ $Name: not supported by cvs2svn $
    12  *  @date $Date: 2004-09-18 01:50:45 $
     11 *  @version $Revision: 1.41 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-09-20 23:16:10 $
    1313 *
    1414 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
    … …  
    9494/*****************************************************************************/
    9595
    96 typedef struct
    97 {
    98     size_t n;                   // Number of data points points in domain.
    99     int paramCount;             // Number of non-masked parameters.
    100     psVector* restrict initialGuess;
    101     const psImage* restrict domain;
    102     const psVector* restrict data;
    103     const psVector* restrict errors;
    104     const psVector* restrict paramMask;
    105     psMinimizeFunction evalModel;
    106     psMinimizeFunctionDeriv d_evalModel;
    107 }
    108 psMinChi2Data;
    109 
    110 typedef struct
    111 {
    112     int paramCount;             // Number of non-masked parameters.
    113     psVector* restrict initialGuess;
    114     const psVector* restrict coord;
    115     const psVector* restrict paramMask;
    116     psMinimizeFunction evalModel;
    117     psMinimizeFunctionDeriv d_evalModel;
    118 }
    119 psMinimizeData;
    120 
    12196/*****************************************************************************/
    12297/* GLOBAL VARIABLES                                                          */
    … …  
    137112/* FUNCTION IMPLEMENTATION - LOCAL                                           */
    138113/*****************************************************************************/
    139 
    140 /******************************************************************************
    141 p_psMinFunc(*params, *funcData): We use the GSL procedure
    142 gsl_multimin_fdfminimizer_iterate() to minimize an arbitary function supplied
    143 by the user.  That GSL procedure requires the function to be minimized to be
    144 in a different format than the psLib format.  The purpose of this procedure
    145 is to serve as a GSL-format wrapper for the user-supplied procedure which is
    146 to be minimized.
    147  
    148     params: The parameters of the function to be minimized.  These will be
    149  varied by GSL in order to minimize the function.
    150  
    151     funcData: a private psLib struct which contains the data point to be
    152  minimized, the function and derivative function pointers, an initial
    153  guess at the parameters, an option parameter mask, etc.
    154  *****************************************************************************/
    155 double p_psMinFunc(const gsl_vector * params, void *funcData)
    156 {
    157     int i;                      // Loop index variable.
    158     int j;                      // Loop index variable.
    159     float tmpf;                 // Temporary floating point variable.
    160     const psVector* restrict coord = ((psMinimizeData* ) funcData)->coord;
    161     const psVector* restrict mask = ((psMinimizeData* ) funcData)->paramMask;
    162     psVector* restrict initialGuess = ((psMinimizeData* ) funcData)->initialGuess;
    163     psMinimizeFunction evalModel = ((psMinimizeData* ) funcData)->evalModel;
    164     psVector* inputParameterList = NULL;
    165 
    166     // The GSL routines will call this function with the masked parameters
    167     // removed.  However, the user-supplied function (to be modified) does not
    168     // have those parameters removed.  Here will create a new parameter list
    169     // with the masked parameters added (we expand initialGuess).
    170     inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
    171     if (mask != NULL) {
    172         j = 0;
    173         for (i = 0; i < mask->n; i++) {
    174             if (mask->data.U8[i] != 0) {
    175                 inputParameterList->data.F32[i] = initialGuess->data.F32[i];
    176             } else {
    177                 inputParameterList->data.F32[i] = gsl_vector_get(params, j++);
    178             }
    179         }
    180     } else {
    181         for (i = 0; i < initialGuess->n; i++) {
    182             inputParameterList->data.F32[i] = gsl_vector_get(params, i);
    183         }
    184     }
    185 
    186     // Call the user-supplied function.
    187     tmpf = evalModel(inputParameterList, coord);
    188 
    189     // Free allocated memory and return the value of the function.
    190     psFree(inputParameterList);
    191     return (tmpf);
    192 }
    193 
    194 /******************************************************************************
    195 p_psMinFuncDeriv(*params, *funcData):  a GSL-like wrapper for the
    196 user-supplied procedure which calculates the derviative of the function to be
    197 minimized.
    198  
    199     params: The parameters of the function to be minimized.  These will be
    200  varied by GSL in order to minimize the function.
    201  
    202     funcData: a private psLib struct which contains the data point to be
    203  minimized, the function and derivative function pointers, an initial
    204  guess at the parameters, an option parameter mask, etc.
    205  
    206     df: we calculate the derivative of the function w.r.t. to each parameter
    207  in "params" and return those derivatives in this psVector.
    208  *****************************************************************************/
    209 void p_psMinFuncDeriv(const gsl_vector * params, void *funcData, gsl_vector * df)
    210 {
    211     int i;                      // Loop index variable.
    212     int j;                      // Loop index variable.
    213     float tmpf;                 // Temporary floating point variable.
    214     const psVector* restrict coord = ((psMinimizeData* ) funcData)->coord;
    215     const psVector* restrict mask = ((psMinimizeData* ) funcData)->paramMask;
    216     psVector* restrict initialGuess = ((psMinimizeData* ) funcData)->initialGuess;
    217     psMinimizeFunctionDeriv d_evalModel = ((psMinimizeData* ) funcData)->d_evalModel;
    218     psVector* inputParameterList = NULL;
    219 
    220     // The GSL routines will call this function with the masked parameters
    221     // removed.  However, the user-supplied function (to be modified) does not
    222     // have those parameters removed.  Here will create a new parameter list
    223     // with the masked parameters added (we expand initialGuess).
    224     inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
    225     if (mask != NULL) {
    226         j = 0;
    227         for (i = 0; i < mask->n; i++) {
    228             if (mask->data.U8[i] != 0) {
    229                 inputParameterList->data.F32[i] = initialGuess->data.F32[i];
    230             } else {
    231                 inputParameterList->data.F32[i] = gsl_vector_get(params, j++);
    232             }
    233         }
    234     } else {
    235         for (i = 0; i < initialGuess->n; i++) {
    236             inputParameterList->data.F32[i] = gsl_vector_get(params, i);
    237         }
    238     }
    239 
    240     // Evaluate the derivative w.r.t. each parameter.
    241     // NOTE: we can probably remove the calls for masked parameters.
    242     for (i = 0; i < initialGuess->n; i++) {
    243         tmpf = d_evalModel(inputParameterList, coord, i);
    244         gsl_vector_set(df, i, tmpf);
    245     }
    246 
    247     // Free allocated memory.
    248     psFree(inputParameterList);
    249 }
    250 
    251 /******************************************************************************
    252     Compute both p_psMinFunc and p_psMinFuncDeriv together.
    253  *****************************************************************************/
    254 void p_psMinFuncFuncDeriv(const gsl_vector * params, void *funcData, double *f, gsl_vector * df)
    255 {
    256     *f = p_psMinFunc(params, funcData);
    257     p_psMinFuncDeriv(params, funcData, df);
    258 }
    259 
    260 // The first argument to evalModel() and d_evalModel() specifies the data
    261 // point.  It must have the same size as the second dimension of *domain.
    262 // The second argument must have the same size as *initialGuess and
    263 // *paramMask.
    264 
    265 /******************************************************************************
    266 p_psMinChi2Func(*x, *funcData, *outdata): We use the GSL procedure
    267 gsl_multifit_fdfsolver_iterate() to fit an arbitrary function, supplied by
    268 the user, to a set of data points.  That GSL procedure requires the function
    269 to be fit to be in a different format than the psLib format.  The purpose of
    270 this procedure is to serve as a GSL-format wrapper for the user-supplied
    271 procedure which is to be fit to the data.
    272  
    273     params: These are the parameters which are to be varied by GSL in order
    274   to minimize chi2 over the data set.
    275  
    276     funcData: this data structure contains the input values over which the
    277   function will be evaluated, the expected value of the function at
    278   those points, the amount of error tolerable at those points, a mask
    279   vector which specifies which parameters to the function are to be
    280   constant, and an initial guess at the parameters.
    281  
    282     outData: The function is evaluated at each point, then subtract the
    283   expected value and divide by the error.
    284  *****************************************************************************/
    285 int p_psMinChi2Func(const gsl_vector * params, void *funcData, gsl_vector * outData)
    286 {
    287     int i;                      // Loop index variable.
    288     int j;                      // Loop index variable.
    289     float tmpf;                 // Temporary floating point variable.
    290     const psImage* restrict domain = ((psMinChi2Data* ) funcData)->domain;
    291     const psVector* restrict data = ((psMinChi2Data* ) funcData)->data;
    292     const psVector* restrict errors = ((psMinChi2Data* ) funcData)->errors;
    293     const psVector* restrict mask = ((psMinChi2Data* ) funcData)->paramMask;
    294     psVector* restrict initialGuess = ((psMinChi2Data* ) funcData)->initialGuess;
    295     psMinimizeFunction evalModel = ((psMinChi2Data* ) funcData)->evalModel;
    296     psVector* inputParameterList = NULL;
    297     psVector* tmpVecPtr = NULL;
    298 
    299     tmpVecPtr = psVectorAlloc(domain->numCols, PS_TYPE_F32);
    300 
    301     // The GSL routines will call this function with the masked parameters
    302     // removed.  However, the user-supplied function (to be modified) does not
    303     // have those parameters removed.  Here will create a new parameter list
    304     // with the masked parameters added (we expand initialGuess).
    305 
    306     inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
    307     if (mask != NULL) {
    308         j = 0;
    309         for (i = 0; i < mask->n; i++) {
    310             if (mask->data.U8[i] != 0) {
    311                 inputParameterList->data.F32[i] = initialGuess->data.F32[i];
    312             } else {
    313                 inputParameterList->data.F32[i] = gsl_vector_get(params, j++);
    314             }
    315         }
    316     } else {
    317         for (i = 0; i < initialGuess->n; i++) {
    318             inputParameterList->data.F32[i] = gsl_vector_get(params, i);
    319         }
    320     }
    321 
    322     // Evaluate the function at each data point.
    323     for (i = 0; i < domain->numRows; i++) {
    324         for (j = 0; j < domain->numCols; j++) {
    325             tmpVecPtr->data.F32[j] = domain->data.F32[i][j];
    326         }
    327         tmpf = evalModel(tmpVecPtr, inputParameterList);
    328 
    329         gsl_vector_set(outData, i, (tmpf - data->data.F32[i]) / errors->data.F32[i]);
    330     }
    331 
    332     // Free allocated memory.
    333     psFree(inputParameterList);
    334     psFree(tmpVecPtr);
    335 
    336     return GSL_SUCCESS;
    337 }
    338 
    339 /******************************************************************************
    340 p_psMinChi2FuncDeriv(*x, *funcData, *outdata): a GSL-like wrapper for the
    341 user-supplied procedure which calculates the derviative of the function to be
    342 minimized.
    343     params: These are the parameters which are to be varied by GSL in order
    344   to minimize chi2 over the data set.
    345  
    346     funcData: this data structure contains the input values over which the
    347   function will be evaluated, the expected value of the function at
    348   those points, the amount of error tolerable at those points, a mask
    349   vector which specifies which parameters to the function are to be
    350   constant, and an initial guess at the parameters.
    351  
    352     J: The derivative is evaluated at each point and w.r.t. each parameter
    353  and returned in this data structure.
    354  *****************************************************************************/
    355 int p_psMinChi2FuncDeriv(const gsl_vector * params, void *funcData, gsl_matrix * J)
    356 {
    357     const psImage* restrict domain = ((psMinChi2Data* ) funcData)->domain;
    358     const psVector* restrict errors = ((psMinChi2Data* ) funcData)->errors;
    359     const psVector* restrict mask = ((psMinChi2Data* ) funcData)->paramMask;
    360     psVector* restrict initialGuess = ((psMinChi2Data* ) funcData)->initialGuess;
    361     psVector* inputParameterList = NULL;
    362     psVector* tmpVecPtr = NULL;
    363     psMinimizeFunctionDeriv d_evalModel = ((psMinChi2Data* ) funcData)->d_evalModel;
    364 
    365     size_t i;
    366     int j;
    367     float tmpf;
    368 
    369     tmpVecPtr = psVectorAlloc(domain->numCols, PS_TYPE_F32);
    370 
    371     // The GSL routines will call this functions with the masked parameters
    372     // removed.  However, the user-supplied function (to be modified) does not
    373     // have those parameters removed.  Here will create a new parameter list
    374     // with the masked parameters added (we expand initialGuess).
    375 
    376     inputParameterList = psVectorAlloc(initialGuess->n, PS_TYPE_F32);
    377     if (mask != NULL) {
    378         j = 0;
    379         for (i = 0; i < mask->n; i++) {
    380             if (mask->data.U8[i] != 0) {
    381                 inputParameterList->data.F32[i] = initialGuess->data.F32[i];
    382             } else {
    383                 inputParameterList->data.F32[i] = gsl_vector_get(params, j++);
    384             }
    385         }
    386     } else {
    387         for (i = 0; i < initialGuess->n; i++) {
    388             inputParameterList->data.F32[i] = gsl_vector_get(params, i);
    389         }
    390     }
    391 
    392     // Evaluate the derivtaive at each data point, and w.r.t. each parameter.
    393     for (i = 0; i < domain->numRows; i++) {
    394         for (j = 0; j < tmpVecPtr->n; j++) {
    395             tmpVecPtr->data.F32[j] = domain->data.F32[i][j];
    396         }
    397 
    398         for (j = 0; j < inputParameterList->n; j++) {
    399             tmpf = d_evalModel(tmpVecPtr, inputParameterList, j);
    400             gsl_matrix_set(J, i, j, (tmpf / errors->data.F32[i]));
    401         }
    402     }
    403 
    404     psFree(inputParameterList);
    405     psFree(tmpVecPtr);
    406     return GSL_SUCCESS;
    407 }
    408 
    409 int p_psMinChi2FuncFuncDeriv(const gsl_vector * params, void *funcData, gsl_vector * f, gsl_matrix * J)
    410 {
    411     p_psMinChi2Func(params, funcData, f);
    412     p_psMinChi2FuncDeriv(params, funcData, J);
    413 
    414     return GSL_SUCCESS;
    415 }
    416114
    417115/******************************************************************************
    … …  
    433131
    434132/******************************************************************************
    435 p_psBuildSums1D(x): this routine returns a psVector with "x" elements.  The
     133VectorNormalizeGen(): this routine returns a psVector with "x" elements.  The
    436134values of the vector will be scaled uniformly between -1.0 and 1.0.
    437  *****************************************************************************/
    438 psVector* psBuildImageScalingFactors(int x)
     135*****************************************************************************/
     136psVector* VectorNormalizeGen(int x)
    439137{
    440     int i = 0;                  // loop index variable.
    441     psVector* imageScalingFactors = NULL;
    442 
    443     imageScalingFactors = psVectorAlloc(x, PS_TYPE_F32);
     138    int i = 0;
     139    psVector *tmp = NULL;
     140
     141    tmp = psVectorAlloc(x, PS_TYPE_F32);
    444142
    445143    for (i = 0; i < x; i++) {
    446         imageScalingFactors->data.F32[i] = (((float)2 * i) / ((float)x)) - 1.0;
    447     }
    448 
    449     return (imageScalingFactors);
    450 }
     144        tmp->data.F32[i] = (((float)2 * i) / ((float)x)) - 1.0;
     145    }
     146
     147    return (tmp);
     148}
     149
     150/*****************************************************************************/
     151/* FUNCTION IMPLEMENTATION - PUBLIC                                          */
     152/*****************************************************************************/
    451153
    452154/******************************************************************************
    453 CURRENTLY NOT IN USE.
    454  
    455 p_psPolyOrderCheck(A, N, *indx, *B, polyOrder,*flag) This routine checks if
    456 all polyOrder-th terms in the polyOrder-th order sky background polynomial
    457 defined by the coefficients in the array B[] are consistent with zero.  If
    458 true, then *flag is set to 1.  Otherwise, *flag is set to 0.  The matrix
    459 inversion code in the middle of this procedure draws from Numerical Recipes
    460 in C page 48.
    461 Input:
    462     A       This is the LUD decomposition of the original matrix A.
    463     N       The size of the matrix (plus 1, actually, since offset 1).
    464     indx    misc Numerical Recipes data structure.
    465     B       The coefficients of the sky polynomial.
    466     polyOrder The degree of the sky polynomial.
    467 Output:
    468     *flag   Set this to 1 if we must recalculate the coefficients.
    469  *****************************************************************************/
    470 void p_psPolyOrderCheck(float **A, int N, int *indx, float *B, int polyOrder, int *flag)
    471 {
    472     float **y = NULL;           // This 2-D matrix will hold A^-1
    473     float *col = NULL;          // misc NumerRecipes data structure
    474     float *error = NULL;        // will hold the sqrt() of the
    475 
    476     // diagonal of y[][].
    477     int i = 0;                  // loop-index variable
    478     int j = 0;                  // loop-index variable
    479     int numPolyTerms = 0;       // The number of terms in the
    480 
    481     // polynomial.
    482     int lastTerm = 0;           // The index location of the first
    483 
    484     // n-th order term in array B[].
    485     int firstTerm = 0;          // Index location of last such term.
    486 
    487     // Allocate the necessary data structures for this procedure...
    488     error = (float *)psAlloc((N + 1) * sizeof(float));
    489     col = (float *)psAlloc((N + 1) * sizeof(float));
    490     y = (float **)psAlloc((N + 1) * sizeof(float *));
    491     for (i = 1; i <= N; i++) {
    492         y[i] = (float *)psAlloc((N + 1) * sizeof(float));
    493     }
    494 
    495     // Invert the matrix A and put the result in y[][].  This code is taken
    496     // from Numerical Recipes in C page 48.
    497     for (j = 1; j <= N; j++) {
    498         for (i = 1; i <= N; i++) {
    499             col[i] = 0.0;
    500         }
    501         col[j] = 1.0;
    502         // NOTE: substitue the LUD rotine
    503         // lubksb(A, N, indx, col);
    504         for (i = 1; i <= N; i++) {
    505             y[i][j] = col[i];
    506         }
    507     }
    508 
    509     // Determine where the first n-th order (in this comment, n equals
    510     // polyOrder) polynomial term is stored in the matrix B[], and also were
    511     // the last n-order term is stored.  Then we loop over all the n-order
    512     // terms and check if they are consistent with zero.
    513 
    514     numPolyTerms = (((polyOrder + 1) * (polyOrder + 2)) / 2);
    515     lastTerm = numPolyTerms + 1;
    516     firstTerm = lastTerm - polyOrder;
    517     *flag = 1;
    518     for (i = firstTerm; i <= lastTerm; i++) {
    519         #ifdef DARWIN
    520         error[i] = (float)sqrt(y[i][i]);
    521         #else
    522 
    523         error[i] = sqrtf(y[i][i]);
    524         #endif
    525 
    526         if (!((B[i] <= (2.0f * error[i])) && ((-2.0f * error[i]) <= B[i]))) {
    527             *flag = 0;
    528         }
    529     }
    530 
    531     // Free all memory allocated in this routine.
    532     psFree(error);
    533     psFree(col);
    534     for (j = 1; j <= N; j++) {
    535         psFree(y[j]);
    536     }
    537     psFree(y);
    538 }
    539 
    540 /*****************************************************************************/
    541 /* FUNCTION IMPLEMENTATION - PUBLIC                                          */
    542 /*****************************************************************************/
    543 
    544 /******************************************************************************
    545 psMinimize(initialGuess, myFunction, myFunctionDeriv, coord, paramMask):
    546  
    547 This routine must minimize an arbitrary function; it determines the set of
    548 parameters of that function such that the ...
    549  *****************************************************************************/
    550 psVector* psMinimize(psVector* restrict initialGuess,
    551                      psMinimizeFunction myFunction,
    552                      psMinimizeFunctionDeriv myFunctionDeriv,
    553                      const psVector* restrict coord,
    554                      const psVector* restrict paramMask)
    555 {
    556     int status;
    557     int i = 0;
    558     int j = 0;
    559     int iter = 0;
    560     gsl_multimin_function_fdf f;
    561     const gsl_multimin_fdfminimizer_type *T;
    562     gsl_multimin_fdfminimizer *s;
    563     psMinimizeData inputData;
    564     gsl_vector *x;
    565 
    566     psTrace(".psLib.dataManip.psMinimize", 4,
    567             "---- psMinimize() begin ----\n");
    568 
    569     PS_CHECK_NULL_VECTOR(initialGuess);
    570     PS_CHECK_EMPTY_VECTOR(initialGuess);
    571     PS_CHECK_NULL_VECTOR(coord);
    572     PS_CHECK_EMPTY_VECTOR(coord);
    573     if (paramMask != NULL) {
    574         PS_CHECK_NULL_VECTOR(paramMask);
    575         PS_CHECK_EMPTY_VECTOR(paramMask);
    576         PS_CHECK_VECTOR_SIZE_EQUAL(initialGuess, paramMask);
    577     }
    578 
    579     inputData.initialGuess = initialGuess;
    580     inputData.coord = coord;
    581     inputData.paramMask = paramMask;
    582     inputData.evalModel = myFunction;
    583     inputData.d_evalModel = myFunctionDeriv;
    584     inputData.paramCount = 0;
    585 
    586     // If the user supplied a parameter mask, then count the number of
    587     // non-masked elements.  This will be used later in allocating a vector
    588     // for the parameters.
    589     if (paramMask != NULL) {
    590         for (i = 0; i < paramMask->n; i++) {
    591             if (paramMask->data.U8[i] != 0) {
    592                 inputData.paramCount++;
    593             }
    594         }
    595     } else {
    596         inputData.paramCount = initialGuess->n;
    597     }
    598 
    599     // The initial guess at the parameters for the function are written into
    600     // the vector inputParameterList.  If the paramMask is not NULL, then
    601     // masked parameters are masked out.
    602     x = gsl_vector_alloc(inputData.paramCount);
    603     if (paramMask != NULL) {
    604         j = 0;
    605         for (i = 0; i < initialGuess->n; i++) {
    606             if (paramMask->data.U8[i] == 0) {
    607                 gsl_vector_set(x, j++, initialGuess->data.F32[i]);
    608             }
    609         }
    610     } else {
    611         for (i = 0; i < initialGuess->n; i++) {
    612             gsl_vector_set(x, i, initialGuess->data.F32[i]);
    613         }
    614     }
    615     f.f = &p_psMinFunc;
    616     f.df = &p_psMinFuncDeriv;
    617     f.fdf = &p_psMinFuncFuncDeriv;
    618     f.n = inputData.paramCount;
    619     f.params = &inputData;
    620 
    621     T = gsl_multimin_fdfminimizer_conjugate_fr;
    622     s = gsl_multimin_fdfminimizer_alloc(T, inputData.paramCount);
    623     gsl_multimin_fdfminimizer_set(s, &f, x, 0.01, 1e-4);
    624     do {
    625         iter++;
    626         status = gsl_multimin_fdfminimizer_iterate(s);
    627 
    628         if (status)
    629             break;
    630 
    631         status = gsl_multimin_test_gradient(s->gradient, 1e-3);
    632 
    633         if (status == GSL_SUCCESS)
    634             printf("Minimum found at:\n");
    635 
    636     } while (status == GSL_CONTINUE && iter < MAX_MINIMIZE_ITERATIONS);
    637 
    638     // In the above steps we had removed the masked elements from the
    639     // the solver.  This next code blocks puts those masked elements
    640     // into the solution.
    641     if (paramMask != NULL) {
    642         j = 0;
    643         for (i = 0; i < initialGuess->n; i++) {
    644             if (paramMask->data.U8[i] == 0) {
    645                 initialGuess->data.F32[i] = gsl_vector_get(s->x, j++);
    646             } else {
    647                 initialGuess->data.F32[i] = initialGuess->data.F32[i];
    648             }
    649         }
    650     } else {
    651         for (i = 0; i < initialGuess->n; i++) {
    652             initialGuess->data.F32[i] = gsl_vector_get(s->x, i);
    653         }
    654     }
    655     psTrace(".psLib.dataManip.psMinimize", 4,
    656             "---- psMinimize() end ----\n");
    657     return (initialGuess);
    658 }
    659 
    660 /******************************************************************************
    661     This routine must determine the parameters of an arbitrary function
    662     such that they best fit the supplied data points.
    663  *****************************************************************************/
    664 psVector* psMinimizeChi2(psMinimizeFunction evalModel,
    665                          psMinimizeFunctionDeriv DevalModel,
    666                          const psImage* restrict domain,
    667                          const psVector* restrict data,
    668                          const psVector* restrict errors,
    669                          psVector* restrict initialGuess,
    670                          const psVector* restrict paramMask,
    671                          float *chiSq)
    672 {
    673     int numData = domain->numRows;      // Number of data points
    674     int status;                 // Return status for the GSL solver.
    675     int i = 0;                  // Loop index variable.
    676     int j = 0;                  // Loop index variable.
    677     int iter = 0;               // Iteration counter.
    678     gsl_multifit_function_fdf f;        // GSL structure that contains the
    679 
    680     // functions/derivative to be solved.
    681     double *xInit = NULL;       // The initial guess at the parameters
    682 
    683     // with masked parameters removed.
    684     const gsl_multifit_fdfsolver_type *T;
    685 
    686     // This tells GSL to use the Levenberg-
    687     // Marquardt algorithm for chi2
    688     // minimization.
    689     gsl_multifit_fdfsolver *s;  // GSL data structure.
    690     psMinChi2Data inputData;
    691     float chiSqOld = 0.0;
    692 
    693     psTrace(".psLib.dataManip.psMinimizeChi2", 4,
    694             "---- psMinimizeChi2() begin ----\n");
    695 
    696     PS_CHECK_NULL_IMAGE(domain);
    697     PS_CHECK_EMPTY_IMAGE(domain);
    698     PS_CHECK_NULL_VECTOR(data);
    699     PS_CHECK_EMPTY_VECTOR(data);
    700     PS_CHECK_NULL_VECTOR(errors);
    701     PS_CHECK_EMPTY_VECTOR(errors);
    702     PS_CHECK_NULL_VECTOR(initialGuess);
    703     PS_CHECK_EMPTY_VECTOR(initialGuess);
    704     PS_CHECK_VECTOR_SIZE_EQUAL(data, errors);
    705     if (domain->numRows != data->n) {
    706         psAbort(__func__, "Number of data points and data values not equal.");
    707     }
    708     if (paramMask != NULL) {
    709         PS_CHECK_NULL_VECTOR(paramMask);
    710         PS_CHECK_EMPTY_VECTOR(paramMask);
    711         PS_CHECK_VECTOR_SIZE_EQUAL(initialGuess, paramMask);
    712     }
    713 
    714     inputData.n = numData;
    715     inputData.paramCount = 0;
    716     inputData.initialGuess = initialGuess;
    717     inputData.domain = domain;
    718     inputData.data = data;
    719     inputData.errors = errors;
    720     inputData.paramMask = paramMask;
    721     inputData.evalModel = evalModel;
    722     inputData.d_evalModel = DevalModel;
    723 
    724     // If the user supplied a parameter mask, then count the number of
    725     // non-masked elements.  This will be used later in allocating a vector
    726     // for the parameters.
    727     if (paramMask != NULL) {
    728         for (i = 0; i < paramMask->n; i++) {
    729             if (paramMask->data.U8[i] != 0) {
    730                 inputData.paramCount++;
    731             }
    732         }
    733     } else {
    734         inputData.paramCount = initialGuess->n;
    735     }
    736 
    737     // The initial guess at the parameters for the function are written into
    738     // the vector inputParameterList.  If the paramMask is not NULL, then those
    739     // parameters are masked out.
    740     xInit = (double *)psAlloc(inputData.paramCount * sizeof(double));
    741     if (paramMask != NULL) {
    742         j = 0;
    743         for (i = 0; i < initialGuess->n; i++) {
    744             if (paramMask->data.U8[i] == 0) {
    745                 xInit[j++] = initialGuess->data.F32[i];
    746             }
    747         }
    748     } else {
    749         for (i = 0; i < initialGuess->n; i++) {
    750             xInit[i] = initialGuess->data.F32[i];
    751         }
    752     }
    753 
    754     const gsl_rng_type *type;
    755     gsl_rng *r;
    756 
    757     gsl_rng_env_setup();
    758 
    759     type = gsl_rng_default;
    760     r = gsl_rng_alloc(type);
    761 
    762     // Initialize the main data structure used by the GSL solver.  This will
    763     // contain pointers to the function to be minimized, it's derivative
    764     // function, the number of data points, the number of free parameters,
    765     // and the data structures those functions use.
    766 
    767     f.f = &p_psMinChi2Func;
    768     f.df = &p_psMinChi2FuncDeriv;
    769     f.fdf = &p_psMinChi2FuncFuncDeriv;
    770     f.n = numData;
    771     f.p = inputData.paramCount;
    772     f.params = &inputData;
    773 
    774     gsl_vector_view x = gsl_vector_view_array(xInit, inputData.paramCount);
    775 
    776     T = gsl_multifit_fdfsolver_lmsder;
    777     s = gsl_multifit_fdfsolver_alloc(T, numData, inputData.paramCount);
    778     gsl_multifit_fdfsolver_set(s, &f, &x.vector);
    779     *chiSq = 0.0;
    780     chiSqOld = 0.0;
    781     do {
    782         iter++;
    783         for (i = 0; i < initialGuess->n; i++) {
    784             printf("Iteration %d: parameter %d is %.3f\n", iter, i, gsl_vector_get(s->x, i));
    785         }
    786         // Perform an iteration of the GSL solver.
    787         status = gsl_multifit_fdfsolver_iterate(s);
    788         printf("gsl_multifit_fdfsolver_iterate() status is %s\n", gsl_strerror(status));
    789         for (i = 0; i < initialGuess->n; i++) {
    790             printf("Iteration %d: parameter %d is %.3f\n", iter, i, gsl_vector_get(s->x, i));
    791         }
    792 
    793         // If there was a problem, abort.
    794         if (status) {
    795             psAbort(__func__, "gsl_multifit_fdfsolver_iterate(%s)\n", gsl_strerror(status));
    796         }
    797         // Test if the parameters changed by a small enough amount.
    798         // NOTE: This wasn't working right when the parameters fit exactly.
    799         // Figure out why.
    800         // status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
    801 
    802         // We test for convergence if chiSquared changes by less than 1.0
    803         // as specified in the ADD.
    804         *chiSq = gsl_blas_dnrm2(s->f);
    805         printf("psMinimize.c: chiSq is %.3f\n", *chiSq);
    806         if (fabs(*chiSq - chiSqOld) < 1.0) {
    807             status = GSL_SUCCESS;
    808         } else {
    809             status = GSL_CONTINUE;
    810         }
    811         chiSqOld = *chiSq;
    812 
    813     } while (status == GSL_CONTINUE && iter < MAX_LMM_ITERATIONS);
    814 
    815     // In the above steps we had removed the masked elements from the
    816     // the solver.  This next code blocks puts those masked elements
    817     // into the solution.
    818     if (paramMask != NULL) {
    819         j = 0;
    820         for (i = 0; i < initialGuess->n; i++) {
    821             if (paramMask->data.U8[i] == 0) {
    822                 initialGuess->data.F32[i] = gsl_vector_get(s->x, j++);
    823             } else {
    824                 initialGuess->data.F32[i] = initialGuess->data.F32[i];
    825             }
    826         }
    827     } else {
    828         for (i = 0; i < initialGuess->n; i++) {
    829             initialGuess->data.F32[i] = gsl_vector_get(s->x, i);
    830         }
    831     }
    832 
    833     // Calculate the chi-squared for the derived solution.
    834     *chiSq = gsl_blas_dnrm2(s->f);
    835 
    836     // Free all allocated memory
    837     // NOTE: Free x.
    838     gsl_multifit_fdfsolver_free(s);
    839     psFree(xInit);
    840 
    841     // Bye bye.
    842     psTrace(".psLib.dataManip.psMinimizeChi2", 4,
    843             "---- psMinimizeChi2() end ----\n");
    844     return (initialGuess);
    845 }
    846 
    847 /******************************************************************************
    848 This routine will take an procedure which calculates an arbitrary function
    849 and it's derivative and minimize it.
    850  
    851 XXX: Do this:
    852  After checking that all entries in the paramMask are 1 or 0, when
    853  forming the A matrix from alpha, try this:
    854  
    855      A[i][i] = (1 + lambda*paramask[i]) * alpha[i][i];
    856  *****************************************************************************/
    857 bool psMinimizeLM(psMinimization *min,
    858                   psImage *covar,
    859                   psVector *params,
    860                   const psVector *paramMask,
    861                   const psArray *coords,
    862                   psMinimizeLMFunc func)
    863 {
    864     psVector *beta = psVectorAlloc(params->n, PS_TYPE_F64);
    865     psVector *perm = psVectorAlloc(params->n, PS_TYPE_F64);
    866     psVector *newParamsF64 = psVectorAlloc(params->n, PS_TYPE_F64);
    867     psVector *newParamsF32 = psVectorAlloc(params->n, PS_TYPE_F32);
    868     psVector *origParams = psVectorAlloc(params->n, PS_TYPE_F32);
    869     psImage *alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32);
    870     psImage *A = psImageAlloc(params->n, params->n, PS_TYPE_F64);
    871     psImage *aOut = psImageAlloc(params->n, params->n, PS_TYPE_F64);
    872     psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32);
    873     psVector *newDeriv = psVectorAlloc(params->n, PS_TYPE_F32);
    874     psVector *NRparams = psVectorAlloc(params->n, PS_TYPE_F32);
    875     int i;
    876     int j;
    877     int k;
    878     float newValue;
    879     float oldValue;
    880     float lamda = 0.00005;
    881 
    882     psTrace(".psLib.dataManip.psMinimizeLM", 4,
    883             "---- psMinimizeLM() begin ----\n");
    884     psTrace(".psLib.dataManip.psMinimizeLM", 6,
    885             "min->maxIter is %d\n", min->maxIter);
    886     psTrace(".psLib.dataManip.psMinimizeLM", 6,
    887             "min->tol is %f\n", min->tol);
    888 
    889     for (i=0;i<params->n;i++) {
    890         origParams->data.F32[i] = params->data.F32[i];
    891     }
    892 
    893     min->lastDelta = HUGE;
    894     min->lastDelta = 12345.0;
    895     min->iter = 0;
    896 
    897     while ((min->lastDelta > min->tol) &&
    898             (min->iter < min->maxIter)) {
    899         psTrace(".psLib.dataManip.psMinimizeLM", 4,
    900                 "------------------------------------------------------\n");
    901         psTrace(".psLib.dataManip.psMinimizeLM", 4,
    902                 "Iteration %d.  Delta is %f\n", min->iter, min->lastDelta);
    903 
    904         min->value = func(deriv, params, coords);
    905 
    906         for (i=0;i<params->n;i++) {
    907             psTrace(".psLib.dataManip.psMinimizeLM", 6,
    908                     "params->data.F32[%d] is %f.\n", i, params->data.F32[i]);
    909         }
    910         for (i=0;i<params->n;i++) {
    911             psTrace(".psLib.dataManip.psMinimizeLM", 6,
    912                     "deriv->data.F32[%d] is %f.\n", i, deriv->data.F32[i]);
    913         }
    914         psTrace(".psLib.dataManip.psMinimizeLM", 6,
    915                 "min->value is (%f)\n", min->value);
    916 
    917         for (i=0;i<params->n;i++) {
    918             if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
    919                 deriv->data.F32[i] = 0.0;
    920             }
    921         }
    922 
    923         // Calculate the BETA vector.
    924         for (i=0;i<params->n;i++) {
    925             beta->data.F64[i] = (float) deriv->data.F32[i];
    926             psTrace(".psLib.dataManip.psMinimizeLM", 6,
    927                     "beta->data.F64[%d] is %f.\n", i, beta->data.F64[i]);
    928         }
    929 
    930         // Calculate the ALPHA matrix.
    931         for (j=0;j<params->n;j++) {
    932             for (k=0;k<params->n;k++) {
    933                 alpha->data.F32[j][k] = deriv->data.F32[k] *
    934                                         deriv->data.F32[j];
    935             }
    936         }
    937 
    938         // Calculate the matrix A.
    939         for (j=0;j<params->n;j++) {
    940             for (k=0;k<params->n;k++) {
    941                 if (j == k) {
    942                     A->data.F64[j][k] =
    943                         (double) ((1.0 + lamda) * alpha->data.F32[j][k]);
    944                 } else {
    945                     A->data.F64[j][k] = (double) alpha->data.F32[j][k];
    946                 }
    947             }
    948         }
    949         for (j=0;j<params->n;j++) {
    950             psTrace(".psLib.dataManip.psMinimizeLM", 6, "Matrix A[][]:\n");
    951             for (k=0;k<params->n;k++) {
    952                 psTrace(".psLib.dataManip.psMinimizeLM", 6, "%f ", A->data.F64[j][k]);
    953             }
    954             psTrace(".psLib.dataManip.psMinimizeLM", 6, "Matrix A[][]:\n");
    955         }
    956 
    957         // Solve A * alpha = Beta
    958         aOut = psMatrixLUD(aOut, perm, A);
    959         newParamsF64 = psMatrixLUSolve(newParamsF64, aOut, beta, perm);
    960 
    961         // Mask any masked parameters.
    962         for (i=0;i<params->n;i++) {
    963             if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
    964                 newParamsF64->data.F64[i] = (double) origParams->data.F32[i];
    965             }
    966             newParamsF32->data.F32[i] = (float) newParamsF64->data.F64[i];
    967 
    968             psTrace(".psLib.dataManip.psMinimizeLM", 6,
    969                     "newParamsF32->data.F32[%d] is %f.\n", i, newParamsF32->data.F32[i]);
    970             NRparams->data.F32[i] = params->data.F32[i] - newParamsF32->data.F32[i];
    971         }
    972 
    973         psTrace(".psLib.dataManip.psMinimizeLM", 6,
    974                 "Calling func() with new parameters:\n");
    975         for (i=0;i<params->n;i++) {
    976             psTrace(".psLib.dataManip.psMinimizeLM", 6,
    977                     "NRparams->data.F32[%d] is %f.\n", i, NRparams->data.F32[i]);
    978         }
    979 
    980         oldValue = min->value;
    981         newValue = func(deriv, NRparams, coords);
    982         psTrace(".psLib.dataManip.psMinimizeLM", 4,
    983                 "old/new values are (%f, %f)\n", oldValue, newValue);
    984         for (i=0;i<params->n;i++) {
    985             if ((paramMask != NULL) && (paramMask->data.U8[i] != 0)) {
    986                 deriv->data.F32[i] = 0.0;
    987             }
    988         }
    989 
    990         if (oldValue > newValue) {
    991             min->lastDelta = oldValue - newValue;
    992             min->value = newValue;
    993 
    994             // No need to check the paramMask here since we already did so
    995             // before the last function evaluation.
    996             for (i=0;i<params->n;i++) {
    997                 //                params->data.F32[i] = (float) newParamsF64->data.F64[i];
    998                 params->data.F32[i] = (float) NRparams->data.F32[i];
    999             }
    1000             min->value = func(deriv, params, coords);
    1001 
    1002             lamda*= 0.1;
    1003         } else {
    1004             lamda*= 10.0;
    1005         }
    1006         psTrace(".psLib.dataManip.psMinimizeLM", 4,
    1007                 "lamda is %f\n", lamda);
    1008         min->iter++;
    1009     }
    1010     psFree(beta);
    1011     psFree(perm);
    1012     psFree(newParamsF64);
    1013     psFree(newParamsF32);
    1014     psFree(origParams);
    1015     psFree(alpha);
    1016     psFree(A);
    1017     psFree(aOut);
    1018     psFree(deriv);
    1019     psFree(newDeriv);
    1020 
    1021     if ((min->iter < min->maxIter) ||
    1022             (min->lastDelta <= min->tol)) {
    1023         return(true);
    1024     }
    1025 
    1026     psTrace(".psLib.dataManip.psMinimizeLM", 4,
    1027             "---- psMinimizeLM() end (false) ----\n");
    1028     return(false);
    1029 }
    1030 
    1031 /******************************************************************************
    1032 This routine will take an procedure which calculates an arbitrary function
    1033 and it's derivative and minimize the chi-squared match between that function
    1034 at the specified coords and the specified value at those coords.
     155psMinimizeLMChi2():  This routine will take an procedure which calculates an
     156arbitrary function and it's derivative and minimize the chi-squared match
     157between that function at the specified coords and the specified value at
     158those coords.
    1035159 
    1036160XXX: Do this:
    … …  
    1265389
    1266390/******************************************************************************
    1267     This routine must fit a polynomial of degree myPoly to the data points
    1268     (x, y) and return the coefficients of that polynomial, as well as the
    1269     error for each data point (yErr).
     391psVectorFitPolynomial1D():  This routine must fit a polynomial of degree
     392myPoly to the data points (x, y) and return the coefficients of that
     393polynomial, as well as the error for each data point (yErr).
    1270394 
    1271395XXX: NOTE: yErr is currently ignored.
    … …  
    1273397psPolynomial1D* psVectorFitPolynomial1D(psPolynomial1D* myPoly,
    1274398                                        const psVector* restrict x,
    1275                                         const psVector* restrict y, const psVector* restrict yErr)
     399                                        const psVector* restrict y,
     400                                        const psVector* restrict yErr)
    1276401{
    1277402    int polyOrder = myPoly->n;
    … …  
    1381506
    1382507#define STEP_SIZE 0.10
    1383 /******************************************************************************
    1384     This routine takes as input an arbitrary function, and the parameter to
    1385     vary.  This function produces as output a bracket [a, b, c] such that
    1386     f(b) is less than f(a) and f(b).
    1387  
    1388     Algorithm: XXX completely ad hoc: start with the user-supplied starting
    1389     parameter and call that b.  Calculate a/c as a fractional amount
    1390     smaller/larger than b.  Repeat this process until a local minimum is
    1391     found.
    1392  *****************************************************************************/
    1393 psVector *p_psDetermineBracketOld(psVector *params,
    1394                                   int dim,
    1395                                   const psArray *coords,
    1396                                   psMinimizePowellFunc func)
    1397 {
    1398     float a = 0.0;
    1399     float b = 0.0;
    1400     float c = 0.0;
    1401     float fa = 0.0;
    1402     float fb = 0.0;
    1403     float fc = 0.0;
    1404     int iter = 100;
    1405     float aDir = 0.0;
    1406     float cDir = 0.0;
    1407     float new_aDir = 0.0;
    1408     float new_cDir = 0.0;
    1409     psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
    1410     float stepSize = params->data.F32[dim] * STEP_SIZE;
    1411     //    float initialParam = params->data.F32[dim];
    1412 
    1413     if (stepSize == 0.0) {
    1414         stepSize = 1.0;
    1415     }
    1416     a = b = c = params->data.F32[0];
    1417     a-= stepSize;
    1418     c+= stepSize;
    1419 
    1420     params->data.F32[dim] = a;
    1421     fa = func(params, coords);
    1422     params->data.F32[dim] = b;
    1423     fb = func(params, coords);
    1424     params->data.F32[dim] = c;
    1425     fc = func(params, coords);
    1426     if (fa < fb) {
    1427         aDir = -1;
    1428     } else {
    1429         aDir = 1;
    1430     }
    1431 
    1432     if (fc < fb) {
    1433         cDir = -1;
    1434     } else {
    1435         cDir = 1;
    1436     }
    1437 
    1438     while (iter > 0) {
    1439         if ((b < a) && (b < c)) {
    1440             bracket->data.F32[0] = a;
    1441             bracket->data.F32[1] = b;
    1442             bracket->data.F32[2] = c;
    1443             return(bracket);
    1444         }
    1445         stepSize*= (1.0 + stepSize);
    1446         a = a - stepSize;
    1447         c = c + stepSize;
    1448 
    1449         params->data.F32[dim] = a;
    1450         fa = func(params, coords);
    1451         params->data.F32[dim] = c;
    1452         fc = func(params, coords);
    1453 
    1454         //printf("HMMM(%d): (%f %f %f) (%f %f %f)\n", iter, a, b, c, fa, fb, fc);
    1455 
    1456         if (fa < fb) {
    1457             new_aDir = -1;
    1458         } else {
    1459             new_aDir = 1;
    1460         }
    1461 
    1462         if (fc < fb) {
    1463             new_cDir = -1;
    1464         } else {
    1465             new_cDir = 1;
    1466         }
    1467         if ((new_aDir == 1) && (aDir == -1)) {
    1468             bracket->data.F32[0] = a;
    1469             bracket->data.F32[1] = b;
    1470             bracket->data.F32[2] = c;
    1471             return(bracket);
    1472         }
    1473 
    1474         if ((new_cDir == 1) && (cDir == -1)) {
    1475             bracket->data.F32[0] = a;
    1476             bracket->data.F32[1] = b;
    1477             bracket->data.F32[2] = c;
    1478             return(bracket);
    1479         }
    1480         aDir = new_aDir;
    1481         cDir = new_cDir;
    1482         iter--;
    1483     }
    1484     psFree(bracket);
    1485     return(NULL);
    1486 }
    1487 
    1488508/******************************************************************************
    1489509    This routine takes as input an arbitrary function, and the parameter to
    … …  
    1666686}
    1667687
    1668 
    1669 /******************************************************************************
    1670     This routine must minimize a possibly multi-dimensional function
    1671     (several parameters) along a single dimension.
    1672  *****************************************************************************/
    1673 bool psMinimize1DFunc(psMinimization *min,
    1674                       psVector *params,
    1675                       int dim,
    1676                       const psArray *coords,
    1677                       psMinimizePowellFunc func)
    1678 {
    1679     psVector *bracket;
    1680     float a = 0.0;
    1681     float b = 0.0;
    1682     float c = 0.0;
    1683     float n = 0.0;
    1684     float fa = 0.0;
    1685     float fb = 0.0;
    1686     float fc = 0.0;
    1687     float fn = 0.0;
    1688     //    float initialParam = params->data.F32[dim];
    1689 
    1690     bracket = p_psDetermineBracketOld(params, dim, coords, func);
    1691     if (bracket == NULL) {
    1692         psAbort(__func__, "Could not bracket minimum.");
    1693     }
    1694 
    1695     min->iter = 0;
    1696     while (min->iter < min->maxIter) {
    1697         min->iter++;
    1698         //printf("psMinimize1DFunc(): iteration %d\n", min->iter);
    1699         a = bracket->data.F32[0];
    1700         b = bracket->data.F32[1];
    1701         c = bracket->data.F32[2];
    1702 
    1703         params->data.F32[dim] = a;
    1704         fa = func(params, coords);
    1705         params->data.F32[dim] = b;
    1706         fb = func(params, coords);
    1707         params->data.F32[dim] = c;
    1708         fc = func(params, coords);
    1709         //printf("Iteration %d: f(%f %f %f) is (%f %f %f)\n", min->iter, a, b, c, fa, fb, fc);
    1710 
    1711         // We determine which is the biggest segment in [a,b,c] then split
    1712         // that with the point n.
    1713         if ((b-a) > (c-b)) {
    1714             // This is the golden section formula
    1715             params->data.F32[dim] = n = a + (0.69 * (b-a));
    1716             fn = func(params, coords);
    1717             if (fn > fb) {
    1718                 // a = n, b = b, c = c
    1719                 bracket->data.F32[0] = n;
    1720             } else {
    1721                 // a = a, b = n, c = b
    1722                 bracket->data.F32[1] = n;
    1723                 bracket->data.F32[2] = b;
    1724             }
    1725         } else {
    1726             params->data.F32[dim] = n = b + (0.69 * (c-b));
    1727             fn = func(params, coords);
    1728             if (fn > fb) {
    1729                 // a = a, b = b, c = n
    1730                 bracket->data.F32[2] = n;
    1731             } else {
    1732                 // a = b, b = n, c = c
    1733                 bracket->data.F32[0] = b;
    1734                 bracket->data.F32[1] = n;
    1735             }
    1736         }
    1737 
    1738         if ((fabs(a-b) < min->tol) &&
    1739                 (fabs(b-c) < min->tol)) {
    1740             //            psFree(bracket);
    1741             //  XXX: is this line correct?
    1742             params->data.F32[dim] = bracket->data.F32[1];
    1743             min->value = func(params, coords);
    1744             return(true);
    1745         }
    1746     }
    1747 
    1748     //    psFree(bracket);
    1749     return(false);
    1750 }
    1751688
    1752689/******************************************************************************
    … …  
    21681105    return(psMinimizePowell(min, params, paramMask, coords, myPowellChi2Func));
    21691106}
     1107
     1108
Note: See TracChangeset for help on using the changeset viewer.