IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jul 8, 2004, 10:50:46 AM (22 years ago)
Author:
gusciora
Message:

Added vector/image size tests.

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r1192 r1199  
     1/** @file  psMinimize.c
     2 *  \brief basic minimization functions
     3 *  @ingroup Math
     4 *
     5 *  This file will contain functions to minimize an arbitrary function at
     6 *  a data point, fit an arbitrary function to a set of data points, and
     7 *  fit a 1-D polynomial to a set of data points.
     8 *
     9 *  @author George Gusciora, MHPCC
     10 *
     11 *  @version $Revision: 1.22 $ $Name: not supported by cvs2svn $
     12 *  @date $Date: 2004-07-08 20:50:46 $
     13 *
     14 *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
     15 */
     16/*****************************************************************************/
     17/* INCLUDE FILES            */
     18/*****************************************************************************/
    119#include <stdlib.h>
    220#include <stdio.h>
    … …  
    2341#include "psMinimize.h"
    2442#include "psMatrix.h"
    25 #include "float.h"
    26 #include <math.h>
    27 
     43/*****************************************************************************/
     44/* DEFINE STATEMENTS           */
     45/*****************************************************************************/
    2846#define MAX_LMM_ITERATIONS 100
    2947#define MAX_MINIMIZE_ITERATIONS 100
     48
     49/** Preprocessor macro to generate error on a NULL 1DPolynomial */
     50#define PS_CHECK_NULL_1DPOLY(NAME)                                                          \
     51if (NAME == NULL || NAME->coeff == NULL) {                                                         \
     52    psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME);                          \
     53}
     54
     55/** Preprocessor macro to generate error on a NULL vector */
     56#define PS_CHECK_NULL_VECTOR(NAME)                                                          \
     57if (NAME == NULL || NAME->data.V == NULL) {                                                         \
     58    psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME);                          \
     59}
     60
     61/** Preprocessor macro to generate error for zero length vector */
     62#define PS_CHECK_EMPTY_VECTOR(NAME)                                                          \
     63if (NAME->n < 1) {                                                                                  \
     64    psError(__func__,"Invalid operation: %s has zero n value.", #NAME);                             \
     65}
     66
     67/** Preprocessor macro to generate error on differing size vectors */
     68#define PS_CHECK_VECTOR_SIZE_EQUAL(VEC1, VEC2)                                                          \
     69if (VEC1->n != VEC2->n) {               \
     70    psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \
     71}
     72
     73/** Preprocessor macro to generate error on a NULL image */
     74#define PS_CHECK_NULL_IMAGE(NAME)                                                           \
     75if (NAME == NULL || NAME->data.V == NULL) {                                                         \
     76    psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME);                          \
     77}
     78
     79/** Preprocessor macro to generate error for zero length rows or columns */
     80#define PS_CHECK_EMPTY_IMAGE(NAME)                                                           \
     81if (NAME->numCols < 1 || NAME->numRows < 1) {                                                       \
     82    psError(__func__,"Invalid operation: %s has zero rows or columns (%dx%d).", #NAME,              \
     83            NAME->numCols, NAME->numRows);                                                          \
     84}
     85
     86
     87/*****************************************************************************/
     88/* TYPE DEFINITIONS           */
     89/*****************************************************************************/
    3090typedef struct
    3191{
    … …  
    53113psMinimizeData;
    54114
    55 
    56 /******************************************************************************
    57 p_psMinFunc(*params, *funcData): We use the GSL-supplied function
     115/*****************************************************************************/
     116/* GLOBAL VARIABLES           */
     117/*****************************************************************************/
     118
     119// None
     120
     121/*****************************************************************************/
     122/* FILE STATIC VARIABLES           */
     123/*****************************************************************************/
     124
     125// None
     126
     127/*****************************************************************************/
     128/* FUNCTION IMPLEMENTATION - LOCAL          */
     129/*****************************************************************************/
     130
     131/******************************************************************************
     132p_psMinFunc(*params, *funcData): We use the GSL procedure
    58133gsl_multimin_fdfminimizer_iterate() to minimize an arbitary function supplied
    59 by the user.  The GSL function requires the user-supplied function to be in
    60 a different format than the psLib format.  The purpose of this procedure is
    61 to serve as a GSL-format wrapper for the psLib user-supplied function which
    62 is to be minimized.
    63     *params: The parameters of the function to be minimized.  These will be
     134by the user.  That GSL procedure requires the function to be minimized to be
     135in a different format than the psLib format.  The purpose of this procedure
     136is to serve as a GSL-format wrapper for the user-supplied procedure which is
     137to be minimized.
     138 
     139    params: The parameters of the function to be minimized.  These will be
    64140 varied by GSL in order to minimize the function.
    65     *funcData: a psLib struct which contains the data point to be minimized,
    66  the function and derivative function pointers, an initial guess at
    67  the parameters, an option parameter mask, etc.
     141 
     142    funcData: a private psLib struct which contains the data point to be
     143 minimized, the function and derivative function pointers, an initial
     144 guess at the parameters, an option parameter mask, etc.
    68145 *****************************************************************************/
    69146double p_psMinFunc(const gsl_vector *params,
    … …  
    109186
    110187/******************************************************************************
    111 p_psMinFuncDeriv(*params, *funcData): We use the GSL-supplied function
    112 gsl_multimin_fdfminimizer_iterate() to minimize an arbitary function supplied
    113 by the user.  The GSL function requires the user-supplied function to be in
    114 a different format than the psLib format.  The purpose of this procedure is
    115 to serve as a GSL-format wrapper for the psLib user-supplied function which
    116 is to be minimized.
    117     *params: The parameters of the function to be minimized.  These will be
     188p_psMinFuncDeriv(*params, *funcData):  a GSL-like wrapper for the
     189user-supplied procedure which calculates the derviative of the function to be
     190minimized.
     191 
     192    params: The parameters of the function to be minimized.  These will be
    118193 varied by GSL in order to minimize the function.
    119     *funcData: a psLib struct which contains the data point to be minimized,
    120  the function and derivative function pointers, an initial guess at
    121  the parameters, an option parameter mask, etc.
    122     *df: we calculate the derivative of the function w.r.t. to each parameter
     194 
     195    funcData: a private psLib struct which contains the data point to be
     196 minimized, the function and derivative function pointers, an initial
     197 guess at the parameters, an option parameter mask, etc.
     198 
     199    df: we calculate the derivative of the function w.r.t. to each parameter
    123200 in "params" and return those derivatives in this psVector.
    124201 *****************************************************************************/
    … …  
    180257}
    181258
    182 
    183 /******************************************************************************
    184 psMinimize(initialGuess, myFunction, myFunctionDeriv, coord, paramMask):
    185  
    186 This routine must minimize an arbitrary function; it must determine the set
    187 of parameters of that function such that the
    188  *****************************************************************************/
    189 psVector *
    190 psMinimize(psVector *restrict initialGuess,
    191            float (*myFunction)(const psVector *restrict, const psVector *restrict),
    192            float (*myFunctionDeriv)(const psVector *restrict, const psVector *restrict, int),
    193            const psVector *restrict coord,
    194            const psVector *restrict paramMask)
    195 {
    196     int status;
    197     int i = 0;
    198     int j = 0;
    199     int iter = 0;
    200     gsl_multimin_function_fdf f;
    201     const gsl_multimin_fdfminimizer_type *T;
    202     gsl_multimin_fdfminimizer *s;
    203     psMinimizeData inputData;
    204     gsl_vector *x;
    205 
    206     inputData.initialGuess = initialGuess;
    207     inputData.coord = coord;
    208     inputData.paramMask = paramMask;
    209     inputData.evalModel = myFunction;
    210     inputData.d_evalModel = myFunctionDeriv;
    211     inputData.paramCount = 0;
    212 
    213     // If the user supplied a parameter mask, then count the number of
    214     // non-masked elements.  This will be used later in allocating a vector
    215     // for the parameters.
    216     if (paramMask != NULL) {
    217         for (i=0;i<paramMask->n;i++) {
    218             if (paramMask->data.U8[i] != 0) {
    219                 inputData.paramCount++;
    220             }
    221         }
    222     } else {
    223         inputData.paramCount= initialGuess->n;
    224     }
    225 
    226     // The initial guess at the parameters for the function are written into
    227     // the vector inputParameterList.  If the paramMask is not NULL, then
    228     // masked parameters are masked out.
    229     x = gsl_vector_alloc(inputData.paramCount);
    230     if (paramMask != NULL) {
    231         j = 0;
    232         for (i=0;i<initialGuess->n;i++) {
    233             if (paramMask->data.U8[i] == 0) {
    234                 gsl_vector_set(x, j++, initialGuess->data.F32[i]);
    235             }
    236         }
    237     } else {
    238         for (i=0;i<initialGuess->n;i++) {
    239             gsl_vector_set(x, i, initialGuess->data.F32[i]);
    240         }
    241     }
    242     f.f = &p_psMinFunc;
    243     f.df = &p_psMinFuncDeriv;
    244     f.fdf = &p_psMinFuncFuncDeriv;
    245     f.n = inputData.paramCount;
    246     f.params = &inputData;
    247 
    248     T = gsl_multimin_fdfminimizer_conjugate_fr;
    249     s = gsl_multimin_fdfminimizer_alloc(T, inputData.paramCount);
    250     gsl_multimin_fdfminimizer_set(s, &f, x, 0.01, 1e-4);
    251     do {
    252         iter++;
    253         status = gsl_multimin_fdfminimizer_iterate(s);
    254 
    255         if (status)
    256             break;
    257 
    258         status = gsl_multimin_test_gradient(s->gradient, 1e-3);
    259 
    260         if (status == GSL_SUCCESS)
    261             printf ("Minimum found at:\n");
    262 
    263     } while (status == GSL_CONTINUE && iter < MAX_MINIMIZE_ITERATIONS);
    264 
    265     // In the above steps we had removed the masked elements from the
    266     // the solver.  This next code blocks puts those masked elements
    267     // into the solution.
    268     if (paramMask != NULL) {
    269         j = 0;
    270         for (i=0;i<initialGuess->n;i++) {
    271             if (paramMask->data.U8[i] == 0) {
    272                 initialGuess->data.F32[i] = gsl_vector_get(s->x, j++);
    273             } else {
    274                 initialGuess->data.F32[i] = initialGuess->data.F32[i];
    275             }
    276         }
    277     } else {
    278         for (i=0;i<initialGuess->n;i++) {
    279             initialGuess->data.F32[i] = gsl_vector_get(s->x, i);
    280         }
    281     }
    282     return(initialGuess);
    283 }
    284 
    285 
    286 
    287 
    288 
    289259// The first argument to evalModel() and d_evalModel() specifies the data
    290260// point.  It must have the same size as the second dimension of *domain.
    … …  
    293263
    294264/******************************************************************************
    295 p_psMinChi2Func(*x, *funcData, *outdata): We use the GSL-supplied function
    296 gsl_multifit_fdfsolver_iterate() to determine the function parameters that
    297 best fit the supllied set of data points.  That GSL function requires the
    298 user-supplied function to be in a different format than the psLib format.
    299 The purpose of this procedure is to serve as a GSL-format wrapper for the
    300 psLib user-supplied function which is to be minimized.
    301     x: These are the parameters which are to be varied by GSL in order to
    302  minimized chi2 over the data set.
     265p_psMinChi2Func(*x, *funcData, *outdata): We use the GSL procedure
     266gsl_multifit_fdfsolver_iterate() to fit an arbitrary function, supplied by
     267the user, to a set of data points.  That GSL procedure requires the function
     268to be fit to be in a different format than the psLib format.  The purpose of
     269this procedure is to serve as a GSL-format wrapper for the user-supplied
     270procedure which is to be fit to the data.
     271 
     272    params: These are the parameters which are to be varied by GSL in order
     273  to minimize chi2 over the data set.
     274 
    303275    funcData: this data structure contains the input values over which the
    304  function will be evaluated, the expected value of the function at
    305  those points, the amount of error tolerable at those points, a mask
    306  vector which specifies which parameters to the function are to be
    307  constant, and an initial guess at the parameters.
     276  function will be evaluated, the expected value of the function at
     277  those points, the amount of error tolerable at those points, a mask
     278  vector which specifies which parameters to the function are to be
     279  constant, and an initial guess at the parameters.
     280 
    308281    outData: The function is evaluated at each point, then subtract the
    309  expected value and divide by the error.
     282  expected value and divide by the error.
    310283 *****************************************************************************/
    311284int p_psMinChi2Func(const gsl_vector *params,
    … …  
    366339}
    367340
     341/******************************************************************************
     342p_psMinChi2FuncDeriv(*x, *funcData, *outdata): a GSL-like wrapper for the
     343user-supplied procedure which calculates the derviative of the function to be
     344minimized.
     345    params: These are the parameters which are to be varied by GSL in order
     346  to minimize chi2 over the data set.
     347 
     348    funcData: this data structure contains the input values over which the
     349  function will be evaluated, the expected value of the function at
     350  those points, the amount of error tolerable at those points, a mask
     351  vector which specifies which parameters to the function are to be
     352  constant, and an initial guess at the parameters.
     353 
     354    J: The derivative is evaluated at each point and w.r.t. each parameter
     355 and returned in this data structure.
     356 *****************************************************************************/
    368357int p_psMinChi2FuncDeriv(const gsl_vector *params,
    369358                         void *funcData,
    … …  
    433422    return GSL_SUCCESS;
    434423}
     424
     425
     426/******************************************************************************
     427p_psBuildSums1D(x, polyOrder, sums): this routine calculates the powers of
     428input parameter "x" between 0 and input parameter polyOrder.  The result is
     429returned as a psVector sums.
     430 *****************************************************************************/
     431void p_psBuildSums1D(double x,
     432                     int polyOrder,
     433                     psVector *sums)
     434{
     435    int       i = 0;
     436    double    xSum = 0.0;
     437
     438    xSum = 1.0;
     439    for(i=0;i<=polyOrder;i++) {
     440        sums->data.F64[i] = xSum;
     441        xSum*= x;
     442    }
     443}
     444
     445
     446/******************************************************************************
     447p_psBuildSums1D(x): this routine returns a psVector with "x" elements.  The
     448values of the vector will be scaled uniformly between -1.0 and 1.0.
     449 *****************************************************************************/
     450psVector *psBuildImageScalingFactors(int x)
     451
     452{
     453    int i = 0;                  // loop index variable.
     454    psVector *imageScalingFactors = NULL;
     455
     456
     457    imageScalingFactors = psVectorAlloc(x, PS_TYPE_F32);
     458
     459    for (i=0;i<x;i++) {
     460        imageScalingFactors->data.F32[i] = (((float) 2*i) / ((float) x)) - 1.0;
     461    }
     462
     463    return(imageScalingFactors);
     464}
     465
     466/******************************************************************************
     467CURRENTLY NOT IN USE.
     468 
     469p_psPolyOrderCheck(A, N, *indx, *B, polyOrder,*flag) This routine checks if
     470all polyOrder-th terms in the polyOrder-th order sky background polynomial
     471defined by the coefficients in the array B[] are consistent with zero.  If
     472true, then *flag is set to 1.  Otherwise, *flag is set to 0.  The matrix
     473inversion code in the middle of this procedure draws from Numerical Recipes
     474in C page 48.
     475Input:
     476    A       This is the LUD decomposition of the original matrix A.
     477    N       The size of the matrix (plus 1, actually, since offset 1).
     478    indx    misc Numerical Recipes data structure.
     479    B       The coefficients of the sky polynomial.
     480    polyOrder The degree of the sky polynomial.
     481Output:
     482    *flag   Set this to 1 if we must recalculate the coefficients.
     483 *****************************************************************************/
     484void p_psPolyOrderCheck(float **A,
     485                        int N,
     486                        int *indx,
     487                        float *B,
     488                        int polyOrder,
     489                        int *flag)
     490{
     491    float     **y = NULL;  // This 2-D matrix will hold A^-1
     492    float      *col = NULL;             // misc NumerRecipes data structure
     493    float      *error=NULL;             // will hold the sqrt() of the
     494    // diagonal of y[][].
     495    int         i=0;                    // loop-index variable
     496    int         j=0;                    // loop-index variable
     497    int         numPolyTerms = 0;       // The number of terms in the
     498    // polynomial.
     499    int         lastTerm = 0;           // The index location of the first
     500    // n-th order term in array B[].
     501    int         firstTerm = 0;          // Index location of last such term.
     502
     503    // Allocate the necessary data structures for this procedure...
     504    error = (float *) psAlloc((N + 1) * sizeof(float));
     505    col = (float *) psAlloc((N + 1) * sizeof(float));
     506    y = (float **) psAlloc((N + 1) * sizeof(float *));
     507    for(i=1;i<=N;i++) {
     508        y[i] = (float *) psAlloc((N + 1) * sizeof(float));
     509    }
     510
     511    // Invert the matrix A and put the result in y[][].  This code is taken
     512    // from Numerical Recipes in C page 48.
     513    for(j=1;j<=N;j++) {
     514        for(i=1;i<=N;i++) {
     515            col[i] = 0.0;
     516        }
     517        col[j] = 1.0;
     518        // NOTE: substitue the LUD rotine
     519        //        lubksb(A, N, indx, col);
     520        for(i=1;i<=N;i++) {
     521            y[i][j] = col[i];
     522        }
     523    }
     524
     525    // Determine where the first n-th order (in this comment, n equals
     526    // polyOrder) polynomial term is stored in the matrix B[], and also were
     527    // the last n-order term is stored.  Then we loop over all the n-order
     528    // terms and check if they are consistent with zero.
     529
     530    numPolyTerms = (((polyOrder+1) * (polyOrder + 2)) / 2);
     531    lastTerm = numPolyTerms + 1;
     532    firstTerm = lastTerm - polyOrder;
     533    *flag = 1;
     534    for (i=firstTerm; i<=lastTerm; i++) {
     535        #ifdef DARWIN
     536        error[i] = (float)sqrt(y[i][i]);
     537        #else
     538
     539        error[i] = sqrtf(y[i][i]);
     540        #endif
     541
     542        if (!((B[i]  <= (2.0f * error[i])) &&
     543                ((-2.0f * error[i]) <= B[i]))) {
     544            *flag = 0;
     545        }
     546    }
     547
     548    // Free all memory allocated in this routine.
     549    psFree(error);
     550    psFree(col);
     551    for(j=1;j<=N;j++) {
     552        psFree(y[j]);
     553    }
     554    psFree(y);
     555}
     556
     557
     558
     559
     560/*****************************************************************************/
     561/* FUNCTION IMPLEMENTATION - PUBLIC         */
     562/*****************************************************************************/
     563
     564
     565
     566/******************************************************************************
     567psMinimize(initialGuess, myFunction, myFunctionDeriv, coord, paramMask):
     568 
     569This routine must minimize an arbitrary function; it determines the set of
     570parameters of that function such that the.
     571 *****************************************************************************/
     572psVector *
     573psMinimize(psVector *restrict initialGuess,
     574           float (*myFunction)(const psVector *restrict, const psVector *restrict),
     575           float (*myFunctionDeriv)(const psVector *restrict, const psVector *restrict, int),
     576           const psVector *restrict coord,
     577           const psVector *restrict paramMask)
     578{
     579    int status;
     580    int i = 0;
     581    int j = 0;
     582    int iter = 0;
     583    gsl_multimin_function_fdf f;
     584    const gsl_multimin_fdfminimizer_type *T;
     585    gsl_multimin_fdfminimizer *s;
     586    psMinimizeData inputData;
     587    gsl_vector *x;
     588
     589    PS_CHECK_NULL_VECTOR(initialGuess);
     590    PS_CHECK_EMPTY_VECTOR(initialGuess);
     591    PS_CHECK_NULL_VECTOR(coord);
     592    PS_CHECK_EMPTY_VECTOR(coord);
     593    if (paramMask != NULL) {
     594        PS_CHECK_NULL_VECTOR(paramMask);
     595        PS_CHECK_EMPTY_VECTOR(paramMask);
     596        PS_CHECK_VECTOR_SIZE_EQUAL(initialGuess, paramMask);
     597    }
     598
     599    inputData.initialGuess = initialGuess;
     600    inputData.coord = coord;
     601    inputData.paramMask = paramMask;
     602    inputData.evalModel = myFunction;
     603    inputData.d_evalModel = myFunctionDeriv;
     604    inputData.paramCount = 0;
     605
     606    // If the user supplied a parameter mask, then count the number of
     607    // non-masked elements.  This will be used later in allocating a vector
     608    // for the parameters.
     609    if (paramMask != NULL) {
     610        for (i=0;i<paramMask->n;i++) {
     611            if (paramMask->data.U8[i] != 0) {
     612                inputData.paramCount++;
     613            }
     614        }
     615    } else {
     616        inputData.paramCount= initialGuess->n;
     617    }
     618
     619    // The initial guess at the parameters for the function are written into
     620    // the vector inputParameterList.  If the paramMask is not NULL, then
     621    // masked parameters are masked out.
     622    x = gsl_vector_alloc(inputData.paramCount);
     623    if (paramMask != NULL) {
     624        j = 0;
     625        for (i=0;i<initialGuess->n;i++) {
     626            if (paramMask->data.U8[i] == 0) {
     627                gsl_vector_set(x, j++, initialGuess->data.F32[i]);
     628            }
     629        }
     630    } else {
     631        for (i=0;i<initialGuess->n;i++) {
     632            gsl_vector_set(x, i, initialGuess->data.F32[i]);
     633        }
     634    }
     635    f.f = &p_psMinFunc;
     636    f.df = &p_psMinFuncDeriv;
     637    f.fdf = &p_psMinFuncFuncDeriv;
     638    f.n = inputData.paramCount;
     639    f.params = &inputData;
     640
     641    T = gsl_multimin_fdfminimizer_conjugate_fr;
     642    s = gsl_multimin_fdfminimizer_alloc(T, inputData.paramCount);
     643    gsl_multimin_fdfminimizer_set(s, &f, x, 0.01, 1e-4);
     644    do {
     645        iter++;
     646        status = gsl_multimin_fdfminimizer_iterate(s);
     647
     648        if (status)
     649            break;
     650
     651        status = gsl_multimin_test_gradient(s->gradient, 1e-3);
     652
     653        if (status == GSL_SUCCESS)
     654            printf ("Minimum found at:\n");
     655
     656    } while (status == GSL_CONTINUE && iter < MAX_MINIMIZE_ITERATIONS);
     657
     658    // In the above steps we had removed the masked elements from the
     659    // the solver.  This next code blocks puts those masked elements
     660    // into the solution.
     661    if (paramMask != NULL) {
     662        j = 0;
     663        for (i=0;i<initialGuess->n;i++) {
     664            if (paramMask->data.U8[i] == 0) {
     665                initialGuess->data.F32[i] = gsl_vector_get(s->x, j++);
     666            } else {
     667                initialGuess->data.F32[i] = initialGuess->data.F32[i];
     668            }
     669        }
     670    } else {
     671        for (i=0;i<initialGuess->n;i++) {
     672            initialGuess->data.F32[i] = gsl_vector_get(s->x, i);
     673        }
     674    }
     675    return(initialGuess);
     676}
     677
     678
     679
     680
     681
    435682
    436683/******************************************************************************
    … …  
    454701    gsl_multifit_function_fdf f; // GSL structure that contains the
    455702    // functions/derivative to be solved.
    456     double *xInit = NULL;     // The initial guess at the parameters
     703    double *xInit = NULL;        // The initial guess at the parameters
    457704    // with masked parameters removed.
    458705    const gsl_multifit_fdfsolver_type *T;
    … …  
    463710    psMinChi2Data inputData;
    464711    float chiSqOld = 0.0;
     712
     713    PS_CHECK_NULL_IMAGE(domain);
     714    PS_CHECK_EMPTY_IMAGE(domain);
     715    PS_CHECK_NULL_VECTOR(data);
     716    PS_CHECK_EMPTY_VECTOR(data);
     717    PS_CHECK_NULL_VECTOR(errors);
     718    PS_CHECK_EMPTY_VECTOR(errors);
     719    PS_CHECK_NULL_VECTOR(initialGuess);
     720    PS_CHECK_EMPTY_VECTOR(initialGuess);
     721    PS_CHECK_VECTOR_SIZE_EQUAL(data, errors);
     722    if (domain->numRows != data->n) {
     723        psAbort(__func__,"Number of data points and data values not equal.");
     724    }
     725    if (paramMask != NULL) {
     726        PS_CHECK_NULL_VECTOR(paramMask);
     727        PS_CHECK_EMPTY_VECTOR(paramMask);
     728        PS_CHECK_VECTOR_SIZE_EQUAL(initialGuess, paramMask);
     729    }
    465730
    466731    inputData.n = numData;
    … …  
    589854
    590855/******************************************************************************
    591 p_psBuildSums1D(x, polyOrder, sums): this routine calculates the powers of
    592 input parameter "x" between 0 and input parameter polyOrder.  The result is
    593 returned as a psVector sums.
    594  *****************************************************************************/
    595 void p_psBuildSums1D(double x,
    596                      int polyOrder,
    597                      psVector *sums)
    598 {
    599     int       i = 0;
    600     double    xSum = 0.0;
    601 
    602     xSum = 1.0;
    603     for(i=0;i<=polyOrder;i++) {
    604         sums->data.F64[i] = xSum;
    605         xSum*= x;
    606     }
    607 }
    608 
    609 
    610 /******************************************************************************
    611 p_psBuildSums1D(x): this routine returns a psVector with "x" elements.  The
    612 values of the vector will be scaled uniformly between -1.0 and 1.0.
    613  *****************************************************************************/
    614 psVector *psBuildImageScalingFactors(int x)
    615 
    616 {
    617     int i = 0;                  // loop index variable.
    618     psVector *imageScalingFactors = NULL;
    619 
    620 
    621     imageScalingFactors = psVectorAlloc(x, PS_TYPE_F32);
    622 
    623     for (i=0;i<x;i++) {
    624         imageScalingFactors->data.F32[i] = (((float) 2*i) / ((float) x)) - 1.0;
    625     }
    626 
    627     return(imageScalingFactors);
    628 }
    629 
    630 /** @brief This routine checks if all polyOrder-th terms in the polyOrder-th
    631  * order sky background polynomial defined by the coefficients in the array B[]
    632  * are consistent with zero.  If true, then *flag is set to 1.  Otherwise,
    633  * *flag is set to 0.  The matrix inversion code in the middle of this
    634  * procedure draws from Numerical Recipes in C page 48.
    635  *
    636  *     Input:
    637  *     <ul>
    638  *         <li> A       This is the LUD decomposition of the original matrix A.
    639  *         <li> N       The size of the matrix (plus 1, actually, since offset 1).
    640  *         <li> indx    misc Numerical Recipes data structure.
    641  *         <li> B       The coefficients of the sky polynomial.
    642  *         <li> polyOrder The degree of the sky polynomial.
    643  *     </ul>
    644  *     Output:
    645  *     <ul>
    646  *         <li> *flag   Set this to 1 if we must recalculate the coefficients.
    647  *     </ul>
    648  *
    649  * @return error status (PsError) indicating error information, or NULL on
    650  * success.
    651  */
    652 
    653 
    654 
    655 void polyOrderCheck(float **A,
    656                     int N,
    657                     int *indx,
    658                     float *B,
    659                     int polyOrder,
    660                     int *flag)
    661 {
    662     float     **y = NULL;  // This 2-D matrix will hold A^-1
    663     float      *col = NULL;             // misc NumerRecipes data structure
    664     float      *error=NULL;             // will hold the sqrt() of the
    665     // diagonal of y[][].
    666     int         i=0;                    // loop-index variable
    667     int         j=0;                    // loop-index variable
    668     int         numPolyTerms = 0;       // The number of terms in the
    669     // polynomial.
    670     int         lastTerm = 0;           // The index location of the first
    671     // n-th order term in array B[].
    672     int         firstTerm = 0;          // Index location of last such term.
    673 
    674     // Allocate the necessary data structures for this procedure...
    675     error = (float *) psAlloc((N + 1) * sizeof(float));
    676     col = (float *) psAlloc((N + 1) * sizeof(float));
    677     y = (float **) psAlloc((N + 1) * sizeof(float *));
    678     for(i=1;i<=N;i++) {
    679         y[i] = (float *) psAlloc((N + 1) * sizeof(float));
    680     }
    681 
    682     // Invert the matrix A and put the result in y[][].  This code is taken
    683     // from Numerical Recipes in C page 48.
    684     for(j=1;j<=N;j++) {
    685         for(i=1;i<=N;i++) {
    686             col[i] = 0.0;
    687         }
    688         col[j] = 1.0;
    689         // NOTE: substitue the LUD rotine
    690         //        lubksb(A, N, indx, col);
    691         for(i=1;i<=N;i++) {
    692             y[i][j] = col[i];
    693         }
    694     }
    695 
    696     // Determine where the first n-th order (in this comment, n equals
    697     // polyOrder) polynomial term is stored in the matrix B[], and also were
    698     // the last n-order term is stored.  Then we loop over all the n-order
    699     // terms and check if they are consistent with zero.
    700 
    701     numPolyTerms = (((polyOrder+1) * (polyOrder + 2)) / 2);
    702     lastTerm = numPolyTerms + 1;
    703     firstTerm = lastTerm - polyOrder;
    704     *flag = 1;
    705     for (i=firstTerm; i<=lastTerm; i++) {
    706         #ifdef DARWIN
    707         error[i] = (float)sqrt(y[i][i]);
    708         #else
    709 
    710         error[i] = sqrtf(y[i][i]);
    711         #endif
    712 
    713         if (!((B[i]  <= (2.0f * error[i])) &&
    714                 ((-2.0f * error[i]) <= B[i]))) {
    715             *flag = 0;
    716         }
    717     }
    718 
    719     // Free all memory allocated in this routine.
    720     psFree(error);
    721     psFree(col);
    722     for(j=1;j<=N;j++) {
    723         psFree(y[j]);
    724     }
    725     psFree(y);
    726 }
    727 
    728 /******************************************************************************
    729856    This routine must fit a polynomial of degree myPoly to the data points
    730857    (x, y) and return the coefficients of that polynomial, as well as the
    731858    error for each data poiny (yErr).
     859 
     860NOTE: yErr is currently ignored.
    732861 *****************************************************************************/
    733862psPolynomial1D *
    … …  
    749878    psVector *xSums = NULL;
    750879
    751     if (x->n != y->n) {
    752         psAbort(__func__, "x and y arguments have different sizes.\n");
    753     }
    754     if (x->n != yErr->n) {
    755         psAbort(__func__, "y and yErr arguments have different sizes.\n");
    756     }
     880    PS_CHECK_NULL_1DPOLY(myPoly);
     881    PS_CHECK_NULL_VECTOR(x);
     882    PS_CHECK_EMPTY_VECTOR(x);
     883    PS_CHECK_NULL_VECTOR(y);
     884    PS_CHECK_EMPTY_VECTOR(y);
     885    PS_CHECK_NULL_VECTOR(yErr);
     886    PS_CHECK_EMPTY_VECTOR(yErr);
     887    PS_CHECK_VECTOR_SIZE_EQUAL(x, y);
     888    PS_CHECK_VECTOR_SIZE_EQUAL(y, yErr);
    757889
    758890    A       = psImageAlloc(polyOrder, polyOrder, PS_TYPE_F64);
Note: See TracChangeset for help on using the changeset viewer.