Index: trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- trunk/psLib/src/dataManip/psMinimize.c (revision 1192) +++ trunk/psLib/src/dataManip/psMinimize.c (revision 1199) @@ -1,2 +1,20 @@ +/** @file psMinimize.c + * \brief basic minimization functions + * @ingroup Math + * + * This file will contain functions to minimize an arbitrary function at + * a data point, fit an arbitrary function to a set of data points, and + * fit a 1-D polynomial to a set of data points. + * + * @author George Gusciora, MHPCC + * + * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-07-08 20:50:46 $ + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ +/*****************************************************************************/ +/* INCLUDE FILES */ +/*****************************************************************************/ #include #include @@ -23,9 +41,51 @@ #include "psMinimize.h" #include "psMatrix.h" -#include "float.h" -#include - +/*****************************************************************************/ +/* DEFINE STATEMENTS */ +/*****************************************************************************/ #define MAX_LMM_ITERATIONS 100 #define MAX_MINIMIZE_ITERATIONS 100 + +/** Preprocessor macro to generate error on a NULL 1DPolynomial */ +#define PS_CHECK_NULL_1DPOLY(NAME) \ +if (NAME == NULL || NAME->coeff == NULL) { \ + psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME); \ +} + +/** Preprocessor macro to generate error on a NULL vector */ +#define PS_CHECK_NULL_VECTOR(NAME) \ +if (NAME == NULL || NAME->data.V == NULL) { \ + psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME); \ +} + +/** Preprocessor macro to generate error for zero length vector */ +#define PS_CHECK_EMPTY_VECTOR(NAME) \ +if (NAME->n < 1) { \ + psError(__func__,"Invalid operation: %s has zero n value.", #NAME); \ +} + +/** Preprocessor macro to generate error on differing size vectors */ +#define PS_CHECK_VECTOR_SIZE_EQUAL(VEC1, VEC2) \ +if (VEC1->n != VEC2->n) { \ + psError(__func__,"Vector %s has size %d, Vector %s has size %d.", #VEC1, VEC1->n, #VEC2, VEC2->n); \ +} + +/** Preprocessor macro to generate error on a NULL image */ +#define PS_CHECK_NULL_IMAGE(NAME) \ +if (NAME == NULL || NAME->data.V == NULL) { \ + psError(__func__,"Invalid operation: %s or its data is NULL.", #NAME); \ +} + +/** Preprocessor macro to generate error for zero length rows or columns */ +#define PS_CHECK_EMPTY_IMAGE(NAME) \ +if (NAME->numCols < 1 || NAME->numRows < 1) { \ + psError(__func__,"Invalid operation: %s has zero rows or columns (%dx%d).", #NAME, \ + NAME->numCols, NAME->numRows); \ +} + + +/*****************************************************************************/ +/* TYPE DEFINITIONS */ +/*****************************************************************************/ typedef struct { @@ -53,17 +113,34 @@ psMinimizeData; - -/****************************************************************************** -p_psMinFunc(*params, *funcData): We use the GSL-supplied function +/*****************************************************************************/ +/* GLOBAL VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FILE STATIC VARIABLES */ +/*****************************************************************************/ + +// None + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - LOCAL */ +/*****************************************************************************/ + +/****************************************************************************** +p_psMinFunc(*params, *funcData): We use the GSL procedure gsl_multimin_fdfminimizer_iterate() to minimize an arbitary function supplied -by the user. The GSL function requires the user-supplied function to be in -a different format than the psLib format. The purpose of this procedure is -to serve as a GSL-format wrapper for the psLib user-supplied function which -is to be minimized. - *params: The parameters of the function to be minimized. These will be +by the user. That GSL procedure requires the function to be minimized to be +in a different format than the psLib format. The purpose of this procedure +is to serve as a GSL-format wrapper for the user-supplied procedure which is +to be minimized. + + params: The parameters of the function to be minimized. These will be varied by GSL in order to minimize the function. - *funcData: a psLib struct which contains the data point to be minimized, - the function and derivative function pointers, an initial guess at - the parameters, an option parameter mask, etc. + + funcData: a private psLib struct which contains the data point to be + minimized, the function and derivative function pointers, an initial + guess at the parameters, an option parameter mask, etc. *****************************************************************************/ double p_psMinFunc(const gsl_vector *params, @@ -109,16 +186,16 @@ /****************************************************************************** -p_psMinFuncDeriv(*params, *funcData): We use the GSL-supplied function -gsl_multimin_fdfminimizer_iterate() to minimize an arbitary function supplied -by the user. The GSL function requires the user-supplied function to be in -a different format than the psLib format. The purpose of this procedure is -to serve as a GSL-format wrapper for the psLib user-supplied function which -is to be minimized. - *params: The parameters of the function to be minimized. These will be +p_psMinFuncDeriv(*params, *funcData): a GSL-like wrapper for the +user-supplied procedure which calculates the derviative of the function to be +minimized. + + params: The parameters of the function to be minimized. These will be varied by GSL in order to minimize the function. - *funcData: a psLib struct which contains the data point to be minimized, - the function and derivative function pointers, an initial guess at - the parameters, an option parameter mask, etc. - *df: we calculate the derivative of the function w.r.t. to each parameter + + funcData: a private psLib struct which contains the data point to be + minimized, the function and derivative function pointers, an initial + guess at the parameters, an option parameter mask, etc. + + df: we calculate the derivative of the function w.r.t. to each parameter in "params" and return those derivatives in this psVector. *****************************************************************************/ @@ -180,111 +257,4 @@ } - -/****************************************************************************** -psMinimize(initialGuess, myFunction, myFunctionDeriv, coord, paramMask): - -This routine must minimize an arbitrary function; it must determine the set -of parameters of that function such that the - *****************************************************************************/ -psVector * -psMinimize(psVector *restrict initialGuess, - float (*myFunction)(const psVector *restrict, const psVector *restrict), - float (*myFunctionDeriv)(const psVector *restrict, const psVector *restrict, int), - const psVector *restrict coord, - const psVector *restrict paramMask) -{ - int status; - int i = 0; - int j = 0; - int iter = 0; - gsl_multimin_function_fdf f; - const gsl_multimin_fdfminimizer_type *T; - gsl_multimin_fdfminimizer *s; - psMinimizeData inputData; - gsl_vector *x; - - inputData.initialGuess = initialGuess; - inputData.coord = coord; - inputData.paramMask = paramMask; - inputData.evalModel = myFunction; - inputData.d_evalModel = myFunctionDeriv; - inputData.paramCount = 0; - - // If the user supplied a parameter mask, then count the number of - // non-masked elements. This will be used later in allocating a vector - // for the parameters. - if (paramMask != NULL) { - for (i=0;in;i++) { - if (paramMask->data.U8[i] != 0) { - inputData.paramCount++; - } - } - } else { - inputData.paramCount= initialGuess->n; - } - - // The initial guess at the parameters for the function are written into - // the vector inputParameterList. If the paramMask is not NULL, then - // masked parameters are masked out. - x = gsl_vector_alloc(inputData.paramCount); - if (paramMask != NULL) { - j = 0; - for (i=0;in;i++) { - if (paramMask->data.U8[i] == 0) { - gsl_vector_set(x, j++, initialGuess->data.F32[i]); - } - } - } else { - for (i=0;in;i++) { - gsl_vector_set(x, i, initialGuess->data.F32[i]); - } - } - f.f = &p_psMinFunc; - f.df = &p_psMinFuncDeriv; - f.fdf = &p_psMinFuncFuncDeriv; - f.n = inputData.paramCount; - f.params = &inputData; - - T = gsl_multimin_fdfminimizer_conjugate_fr; - s = gsl_multimin_fdfminimizer_alloc(T, inputData.paramCount); - gsl_multimin_fdfminimizer_set(s, &f, x, 0.01, 1e-4); - do { - iter++; - status = gsl_multimin_fdfminimizer_iterate(s); - - if (status) - break; - - status = gsl_multimin_test_gradient(s->gradient, 1e-3); - - if (status == GSL_SUCCESS) - printf ("Minimum found at:\n"); - - } while (status == GSL_CONTINUE && iter < MAX_MINIMIZE_ITERATIONS); - - // In the above steps we had removed the masked elements from the - // the solver. This next code blocks puts those masked elements - // into the solution. - if (paramMask != NULL) { - j = 0; - for (i=0;in;i++) { - if (paramMask->data.U8[i] == 0) { - initialGuess->data.F32[i] = gsl_vector_get(s->x, j++); - } else { - initialGuess->data.F32[i] = initialGuess->data.F32[i]; - } - } - } else { - for (i=0;in;i++) { - initialGuess->data.F32[i] = gsl_vector_get(s->x, i); - } - } - return(initialGuess); -} - - - - - // The first argument to evalModel() and d_evalModel() specifies the data // point. It must have the same size as the second dimension of *domain. @@ -293,19 +263,22 @@ /****************************************************************************** -p_psMinChi2Func(*x, *funcData, *outdata): We use the GSL-supplied function -gsl_multifit_fdfsolver_iterate() to determine the function parameters that -best fit the supllied set of data points. That GSL function requires the -user-supplied function to be in a different format than the psLib format. -The purpose of this procedure is to serve as a GSL-format wrapper for the -psLib user-supplied function which is to be minimized. - x: These are the parameters which are to be varied by GSL in order to - minimized chi2 over the data set. +p_psMinChi2Func(*x, *funcData, *outdata): We use the GSL procedure +gsl_multifit_fdfsolver_iterate() to fit an arbitrary function, supplied by +the user, to a set of data points. That GSL procedure requires the function +to be fit to be in a different format than the psLib format. The purpose of +this procedure is to serve as a GSL-format wrapper for the user-supplied +procedure which is to be fit to the data. + + params: These are the parameters which are to be varied by GSL in order + to minimize chi2 over the data set. + funcData: this data structure contains the input values over which the - function will be evaluated, the expected value of the function at - those points, the amount of error tolerable at those points, a mask - vector which specifies which parameters to the function are to be - constant, and an initial guess at the parameters. + function will be evaluated, the expected value of the function at + those points, the amount of error tolerable at those points, a mask + vector which specifies which parameters to the function are to be + constant, and an initial guess at the parameters. + outData: The function is evaluated at each point, then subtract the - expected value and divide by the error. + expected value and divide by the error. *****************************************************************************/ int p_psMinChi2Func(const gsl_vector *params, @@ -366,4 +339,20 @@ } +/****************************************************************************** +p_psMinChi2FuncDeriv(*x, *funcData, *outdata): a GSL-like wrapper for the +user-supplied procedure which calculates the derviative of the function to be +minimized. + params: These are the parameters which are to be varied by GSL in order + to minimize chi2 over the data set. + + funcData: this data structure contains the input values over which the + function will be evaluated, the expected value of the function at + those points, the amount of error tolerable at those points, a mask + vector which specifies which parameters to the function are to be + constant, and an initial guess at the parameters. + + J: The derivative is evaluated at each point and w.r.t. each parameter + and returned in this data structure. + *****************************************************************************/ int p_psMinChi2FuncDeriv(const gsl_vector *params, void *funcData, @@ -433,4 +422,262 @@ return GSL_SUCCESS; } + + +/****************************************************************************** +p_psBuildSums1D(x, polyOrder, sums): this routine calculates the powers of +input parameter "x" between 0 and input parameter polyOrder. The result is +returned as a psVector sums. + *****************************************************************************/ +void p_psBuildSums1D(double x, + int polyOrder, + psVector *sums) +{ + int i = 0; + double xSum = 0.0; + + xSum = 1.0; + for(i=0;i<=polyOrder;i++) { + sums->data.F64[i] = xSum; + xSum*= x; + } +} + + +/****************************************************************************** +p_psBuildSums1D(x): this routine returns a psVector with "x" elements. The +values of the vector will be scaled uniformly between -1.0 and 1.0. + *****************************************************************************/ +psVector *psBuildImageScalingFactors(int x) + +{ + int i = 0; // loop index variable. + psVector *imageScalingFactors = NULL; + + + imageScalingFactors = psVectorAlloc(x, PS_TYPE_F32); + + for (i=0;idata.F32[i] = (((float) 2*i) / ((float) x)) - 1.0; + } + + return(imageScalingFactors); +} + +/****************************************************************************** +CURRENTLY NOT IN USE. + +p_psPolyOrderCheck(A, N, *indx, *B, polyOrder,*flag) This routine checks if +all polyOrder-th terms in the polyOrder-th order sky background polynomial +defined by the coefficients in the array B[] are consistent with zero. If +true, then *flag is set to 1. Otherwise, *flag is set to 0. The matrix +inversion code in the middle of this procedure draws from Numerical Recipes +in C page 48. +Input: + A This is the LUD decomposition of the original matrix A. + N The size of the matrix (plus 1, actually, since offset 1). + indx misc Numerical Recipes data structure. + B The coefficients of the sky polynomial. + polyOrder The degree of the sky polynomial. +Output: + *flag Set this to 1 if we must recalculate the coefficients. + *****************************************************************************/ +void p_psPolyOrderCheck(float **A, + int N, + int *indx, + float *B, + int polyOrder, + int *flag) +{ + float **y = NULL; // This 2-D matrix will hold A^-1 + float *col = NULL; // misc NumerRecipes data structure + float *error=NULL; // will hold the sqrt() of the + // diagonal of y[][]. + int i=0; // loop-index variable + int j=0; // loop-index variable + int numPolyTerms = 0; // The number of terms in the + // polynomial. + int lastTerm = 0; // The index location of the first + // n-th order term in array B[]. + int firstTerm = 0; // Index location of last such term. + + // Allocate the necessary data structures for this procedure... + error = (float *) psAlloc((N + 1) * sizeof(float)); + col = (float *) psAlloc((N + 1) * sizeof(float)); + y = (float **) psAlloc((N + 1) * sizeof(float *)); + for(i=1;i<=N;i++) { + y[i] = (float *) psAlloc((N + 1) * sizeof(float)); + } + + // Invert the matrix A and put the result in y[][]. This code is taken + // from Numerical Recipes in C page 48. + for(j=1;j<=N;j++) { + for(i=1;i<=N;i++) { + col[i] = 0.0; + } + col[j] = 1.0; + // NOTE: substitue the LUD rotine + // lubksb(A, N, indx, col); + for(i=1;i<=N;i++) { + y[i][j] = col[i]; + } + } + + // Determine where the first n-th order (in this comment, n equals + // polyOrder) polynomial term is stored in the matrix B[], and also were + // the last n-order term is stored. Then we loop over all the n-order + // terms and check if they are consistent with zero. + + numPolyTerms = (((polyOrder+1) * (polyOrder + 2)) / 2); + lastTerm = numPolyTerms + 1; + firstTerm = lastTerm - polyOrder; + *flag = 1; + for (i=firstTerm; i<=lastTerm; i++) { + #ifdef DARWIN + error[i] = (float)sqrt(y[i][i]); + #else + + error[i] = sqrtf(y[i][i]); + #endif + + if (!((B[i] <= (2.0f * error[i])) && + ((-2.0f * error[i]) <= B[i]))) { + *flag = 0; + } + } + + // Free all memory allocated in this routine. + psFree(error); + psFree(col); + for(j=1;j<=N;j++) { + psFree(y[j]); + } + psFree(y); +} + + + + +/*****************************************************************************/ +/* FUNCTION IMPLEMENTATION - PUBLIC */ +/*****************************************************************************/ + + + +/****************************************************************************** +psMinimize(initialGuess, myFunction, myFunctionDeriv, coord, paramMask): + +This routine must minimize an arbitrary function; it determines the set of +parameters of that function such that the. + *****************************************************************************/ +psVector * +psMinimize(psVector *restrict initialGuess, + float (*myFunction)(const psVector *restrict, const psVector *restrict), + float (*myFunctionDeriv)(const psVector *restrict, const psVector *restrict, int), + const psVector *restrict coord, + const psVector *restrict paramMask) +{ + int status; + int i = 0; + int j = 0; + int iter = 0; + gsl_multimin_function_fdf f; + const gsl_multimin_fdfminimizer_type *T; + gsl_multimin_fdfminimizer *s; + psMinimizeData inputData; + gsl_vector *x; + + PS_CHECK_NULL_VECTOR(initialGuess); + PS_CHECK_EMPTY_VECTOR(initialGuess); + PS_CHECK_NULL_VECTOR(coord); + PS_CHECK_EMPTY_VECTOR(coord); + if (paramMask != NULL) { + PS_CHECK_NULL_VECTOR(paramMask); + PS_CHECK_EMPTY_VECTOR(paramMask); + PS_CHECK_VECTOR_SIZE_EQUAL(initialGuess, paramMask); + } + + inputData.initialGuess = initialGuess; + inputData.coord = coord; + inputData.paramMask = paramMask; + inputData.evalModel = myFunction; + inputData.d_evalModel = myFunctionDeriv; + inputData.paramCount = 0; + + // If the user supplied a parameter mask, then count the number of + // non-masked elements. This will be used later in allocating a vector + // for the parameters. + if (paramMask != NULL) { + for (i=0;in;i++) { + if (paramMask->data.U8[i] != 0) { + inputData.paramCount++; + } + } + } else { + inputData.paramCount= initialGuess->n; + } + + // The initial guess at the parameters for the function are written into + // the vector inputParameterList. If the paramMask is not NULL, then + // masked parameters are masked out. + x = gsl_vector_alloc(inputData.paramCount); + if (paramMask != NULL) { + j = 0; + for (i=0;in;i++) { + if (paramMask->data.U8[i] == 0) { + gsl_vector_set(x, j++, initialGuess->data.F32[i]); + } + } + } else { + for (i=0;in;i++) { + gsl_vector_set(x, i, initialGuess->data.F32[i]); + } + } + f.f = &p_psMinFunc; + f.df = &p_psMinFuncDeriv; + f.fdf = &p_psMinFuncFuncDeriv; + f.n = inputData.paramCount; + f.params = &inputData; + + T = gsl_multimin_fdfminimizer_conjugate_fr; + s = gsl_multimin_fdfminimizer_alloc(T, inputData.paramCount); + gsl_multimin_fdfminimizer_set(s, &f, x, 0.01, 1e-4); + do { + iter++; + status = gsl_multimin_fdfminimizer_iterate(s); + + if (status) + break; + + status = gsl_multimin_test_gradient(s->gradient, 1e-3); + + if (status == GSL_SUCCESS) + printf ("Minimum found at:\n"); + + } while (status == GSL_CONTINUE && iter < MAX_MINIMIZE_ITERATIONS); + + // In the above steps we had removed the masked elements from the + // the solver. This next code blocks puts those masked elements + // into the solution. + if (paramMask != NULL) { + j = 0; + for (i=0;in;i++) { + if (paramMask->data.U8[i] == 0) { + initialGuess->data.F32[i] = gsl_vector_get(s->x, j++); + } else { + initialGuess->data.F32[i] = initialGuess->data.F32[i]; + } + } + } else { + for (i=0;in;i++) { + initialGuess->data.F32[i] = gsl_vector_get(s->x, i); + } + } + return(initialGuess); +} + + + + + /****************************************************************************** @@ -454,5 +701,5 @@ gsl_multifit_function_fdf f; // GSL structure that contains the // functions/derivative to be solved. - double *xInit = NULL; // The initial guess at the parameters + double *xInit = NULL; // The initial guess at the parameters // with masked parameters removed. const gsl_multifit_fdfsolver_type *T; @@ -463,4 +710,22 @@ psMinChi2Data inputData; float chiSqOld = 0.0; + + PS_CHECK_NULL_IMAGE(domain); + PS_CHECK_EMPTY_IMAGE(domain); + PS_CHECK_NULL_VECTOR(data); + PS_CHECK_EMPTY_VECTOR(data); + PS_CHECK_NULL_VECTOR(errors); + PS_CHECK_EMPTY_VECTOR(errors); + PS_CHECK_NULL_VECTOR(initialGuess); + PS_CHECK_EMPTY_VECTOR(initialGuess); + PS_CHECK_VECTOR_SIZE_EQUAL(data, errors); + if (domain->numRows != data->n) { + psAbort(__func__,"Number of data points and data values not equal."); + } + if (paramMask != NULL) { + PS_CHECK_NULL_VECTOR(paramMask); + PS_CHECK_EMPTY_VECTOR(paramMask); + PS_CHECK_VECTOR_SIZE_EQUAL(initialGuess, paramMask); + } inputData.n = numData; @@ -589,145 +854,9 @@ /****************************************************************************** -p_psBuildSums1D(x, polyOrder, sums): this routine calculates the powers of -input parameter "x" between 0 and input parameter polyOrder. The result is -returned as a psVector sums. - *****************************************************************************/ -void p_psBuildSums1D(double x, - int polyOrder, - psVector *sums) -{ - int i = 0; - double xSum = 0.0; - - xSum = 1.0; - for(i=0;i<=polyOrder;i++) { - sums->data.F64[i] = xSum; - xSum*= x; - } -} - - -/****************************************************************************** -p_psBuildSums1D(x): this routine returns a psVector with "x" elements. The -values of the vector will be scaled uniformly between -1.0 and 1.0. - *****************************************************************************/ -psVector *psBuildImageScalingFactors(int x) - -{ - int i = 0; // loop index variable. - psVector *imageScalingFactors = NULL; - - - imageScalingFactors = psVectorAlloc(x, PS_TYPE_F32); - - for (i=0;idata.F32[i] = (((float) 2*i) / ((float) x)) - 1.0; - } - - return(imageScalingFactors); -} - -/** @brief This routine checks if all polyOrder-th terms in the polyOrder-th - * order sky background polynomial defined by the coefficients in the array B[] - * are consistent with zero. If true, then *flag is set to 1. Otherwise, - * *flag is set to 0. The matrix inversion code in the middle of this - * procedure draws from Numerical Recipes in C page 48. - * - * Input: - *
    - *
  • A This is the LUD decomposition of the original matrix A. - *
  • N The size of the matrix (plus 1, actually, since offset 1). - *
  • indx misc Numerical Recipes data structure. - *
  • B The coefficients of the sky polynomial. - *
  • polyOrder The degree of the sky polynomial. - *
- * Output: - *
    - *
  • *flag Set this to 1 if we must recalculate the coefficients. - *
- * - * @return error status (PsError) indicating error information, or NULL on - * success. - */ - - - -void polyOrderCheck(float **A, - int N, - int *indx, - float *B, - int polyOrder, - int *flag) -{ - float **y = NULL; // This 2-D matrix will hold A^-1 - float *col = NULL; // misc NumerRecipes data structure - float *error=NULL; // will hold the sqrt() of the - // diagonal of y[][]. - int i=0; // loop-index variable - int j=0; // loop-index variable - int numPolyTerms = 0; // The number of terms in the - // polynomial. - int lastTerm = 0; // The index location of the first - // n-th order term in array B[]. - int firstTerm = 0; // Index location of last such term. - - // Allocate the necessary data structures for this procedure... - error = (float *) psAlloc((N + 1) * sizeof(float)); - col = (float *) psAlloc((N + 1) * sizeof(float)); - y = (float **) psAlloc((N + 1) * sizeof(float *)); - for(i=1;i<=N;i++) { - y[i] = (float *) psAlloc((N + 1) * sizeof(float)); - } - - // Invert the matrix A and put the result in y[][]. This code is taken - // from Numerical Recipes in C page 48. - for(j=1;j<=N;j++) { - for(i=1;i<=N;i++) { - col[i] = 0.0; - } - col[j] = 1.0; - // NOTE: substitue the LUD rotine - // lubksb(A, N, indx, col); - for(i=1;i<=N;i++) { - y[i][j] = col[i]; - } - } - - // Determine where the first n-th order (in this comment, n equals - // polyOrder) polynomial term is stored in the matrix B[], and also were - // the last n-order term is stored. Then we loop over all the n-order - // terms and check if they are consistent with zero. - - numPolyTerms = (((polyOrder+1) * (polyOrder + 2)) / 2); - lastTerm = numPolyTerms + 1; - firstTerm = lastTerm - polyOrder; - *flag = 1; - for (i=firstTerm; i<=lastTerm; i++) { - #ifdef DARWIN - error[i] = (float)sqrt(y[i][i]); - #else - - error[i] = sqrtf(y[i][i]); - #endif - - if (!((B[i] <= (2.0f * error[i])) && - ((-2.0f * error[i]) <= B[i]))) { - *flag = 0; - } - } - - // Free all memory allocated in this routine. - psFree(error); - psFree(col); - for(j=1;j<=N;j++) { - psFree(y[j]); - } - psFree(y); -} - -/****************************************************************************** This routine must fit a polynomial of degree myPoly to the data points (x, y) and return the coefficients of that polynomial, as well as the error for each data poiny (yErr). + +NOTE: yErr is currently ignored. *****************************************************************************/ psPolynomial1D * @@ -749,10 +878,13 @@ psVector *xSums = NULL; - if (x->n != y->n) { - psAbort(__func__, "x and y arguments have different sizes.\n"); - } - if (x->n != yErr->n) { - psAbort(__func__, "y and yErr arguments have different sizes.\n"); - } + PS_CHECK_NULL_1DPOLY(myPoly); + PS_CHECK_NULL_VECTOR(x); + PS_CHECK_EMPTY_VECTOR(x); + PS_CHECK_NULL_VECTOR(y); + PS_CHECK_EMPTY_VECTOR(y); + PS_CHECK_NULL_VECTOR(yErr); + PS_CHECK_EMPTY_VECTOR(yErr); + PS_CHECK_VECTOR_SIZE_EQUAL(x, y); + PS_CHECK_VECTOR_SIZE_EQUAL(y, yErr); A = psImageAlloc(polyOrder, polyOrder, PS_TYPE_F64);