Index: /trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- /trunk/psLib/src/dataManip/psMinimize.c (revision 1092) +++ /trunk/psLib/src/dataManip/psMinimize.c (revision 1093) @@ -19,30 +19,28 @@ #include "psFunctions.h" #include "psSort.h" +#include "psMinimize.h" #include "float.h" #include -// NOTE: rewrite so there is no maximum order for the polynomials. -#define MAX_POLY_ORDER 10 -#define MAX_POLYNOMIAL_TERMS (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2) -int MyInfoLevel = 0; - + +#define MAX_LMM_NUM_ITERATIONS 500 typedef struct { size_t n; // Number of data points points in domain. int count; // Number of non-masked parameters. - psVector *restrict initialGuess; - psImage *domain; // - psVector *data; - psVector *errors; - psVector *paramMask; - // The first argument to evalModel() and - // d_evalModel() specifies the data point. - // It must have the same size as the second - // dimension of *domain. - // The second argument must have the same size - // as *initialGuess and *paramMask. + psVector *restrict initialGuess; + const psImage *restrict domain; + const psVector *restrict data; + const psVector *restrict errors; + const psVector *restrict paramMask; float (*evalModel) (psVector *, psVector *); - float (*d_evalModel) (psVector *, psVector *); + float (*d_evalModel) (psVector *, psVector *, int); } psModelData; +// The first argument to evalModel() and +// d_evalModel() specifies the data point. +// It must have the same size as the second +// dimension of *domain. +// The second argument must have the same size +// as *initialGuess and *paramMask. @@ -65,10 +63,10 @@ for the user supplied function which is to be minimized. The GSL minimization routines have no knowledge of the user-supplied function. They -only clal this routine, which then calls the user-supplied function. The +only call this routine, which then calls the user-supplied function. The arguments are: x: These are the parameter which are to be varied by GSL in order to minimized chi2 over the data set. params: - f + f: @@ -79,16 +77,20 @@ gsl_vector *f) { - psImage *domain = ((psModelData *)params)->domain; - psVector *data = ((psModelData *)params)->data; - psVector *errors = ((psModelData *) params)->errors; - psVector *mask = ((psModelData *) params)->paramMask; - psVector *initialGuess = ((psModelData *)params)->initialGuess; + int i; // Loop index variable. + int j; // Loop index variable. + float tmpf; // Temporary floating point variable. + const psImage *restrict domain = ((psModelData *)params)->domain; + const psVector *restrict data = ((psModelData *)params)->data; + const psVector *restrict errors = ((psModelData *) params)->errors; + const psVector *restrict mask = ((psModelData *) params)->paramMask; + psVector *restrict initialGuess = ((psModelData *)params)->initialGuess; + float (*evalModel)(psVector *, psVector *) = ((psModelData *) params)->evalModel; psVector *inputParameterList = NULL; - // float (*evalModel)(psVector *, psVector *) = ((psModelData *) params)->evalModel; - size_t i; - int j; - float tmpf; - - // The GSL routines will call this functions with the masked parameters + psVector *tmpVecPtr = NULL; + + + tmpVecPtr = psVectorAlloc(domain->numCols, PS_TYPE_F32); + + // The GSL routines will call this function with the masked parameters // removed. However, the user-supplied function (to be modified) does not // have those parameters removed. Here will create a new parameter list @@ -112,8 +114,13 @@ for (i=0;inumRows;i++) { - //GUS tmpf = evalModel(inputParameterList, domain->data.F32[i]); + for (j=0;jn;j++) { + tmpVecPtr->data.F32[j] = domain->data.F32[i][j]; + } + tmpf = evalModel(inputParameterList, tmpVecPtr); gsl_vector_set(f, i, (tmpf - data->data.F32[i])/errors->data.F32[i]); } + psFree(inputParameterList); + psFree(tmpVecPtr); return GSL_SUCCESS; } @@ -123,16 +130,17 @@ gsl_matrix *J) { - // size_t n = ((psModelData *)params)->n; - psImage *domain = ((psModelData *)params)->domain; - // psVector *data = ((psModelData *)params)->data; - psVector *errors = ((psModelData *) params)->errors; - psVector *mask = ((psModelData *) params)->paramMask; - psVector *initialGuess = ((psModelData *)params)->initialGuess; + const psImage *restrict domain = ((psModelData *)params)->domain; + const psVector *restrict errors = ((psModelData *) params)->errors; + const psVector *restrict mask = ((psModelData *) params)->paramMask; + psVector *restrict initialGuess = ((psModelData *)params)->initialGuess; psVector *inputParameterList = NULL; - // int count = 0; - // float (*d_evalModel) = ((psModelData *) params)->evalModel; + psVector *tmpVecPtr = NULL; + float (*d_evalModel)(psVector *, psVector *, int) = ((psModelData *) params)->d_evalModel; + size_t i; int j; float tmpf; + + tmpVecPtr = psVectorAlloc(domain->numCols, PS_TYPE_F32); // The GSL routines will call this functions with the masked parameters @@ -157,11 +165,16 @@ for (i=0;inumRows;i++) { + for (j=0;jn;j++) { + tmpVecPtr->data.F32[j] = domain->data.F32[i][j]; + } + for (j=0;jn;j++) { - //GUS tmpf = d_evalModel(inputParameterList, domain->data.F32[i], j); - + tmpf = d_evalModel(inputParameterList, tmpVecPtr, j); gsl_matrix_set(J, i, j, (tmpf/errors->data.F32[i])); } } + psFree(inputParameterList); + psFree(tmpVecPtr); return GSL_SUCCESS; } @@ -185,6 +198,6 @@ psVector * psMinimizeChi2(float (*evalModel)(const psVector *restrict, const psVector *restrict), - float (*DevalModel)(psVector *, psVector *), - const psVector *restrict domain, + float (*DevalModel)(const psVector *restrict, const psVector *restrict, int), + const psImage *restrict domain, const psVector *restrict data, const psVector *restrict errors, @@ -193,20 +206,27 @@ float *chiSq) { - size_t n = domain->n; + + int numData = domain->numCols; // Number of data points + gsl_multifit_fdfsolver *s = NULL; // GSL data structure. + int status; // Return status for the GSL solver. + int i = 0; // Loop index variable. + int j = 0; // Loop index variable. + int iter = 0; // Iteration counter. + psModelData *inputData = NULL; // Contains data that the user-supplied + // function/derivate need. + gsl_multifit_function_fdf f; // GSL structure that contains the + // functions/derivative to be solved. + double *xInit = NULL; // The initial guess at the parameters + // with masked parameters removed. const gsl_multifit_fdfsolver_type *T = gsl_multifit_fdfsolver_lmsder; - gsl_multifit_fdfsolver *s = NULL; - int status; - size_t i = 0; - int j = 0; - size_t iter = 0; - psModelData *inputData = NULL; - gsl_multifit_function_fdf f; - double *xInit = NULL; + // This tells GSL to use the Levenberg- + // Marquardt algorithm for chi2 + // minimization. inputData = (psModelData *) psAlloc(sizeof(psModelData)); - inputData->n = n; + inputData->n = numData; inputData->count = 0; inputData->initialGuess = initialGuess; - // inputData->domain = domain; + inputData->domain = domain; // inputData->data = data; // inputData->errors = errors; @@ -226,10 +246,7 @@ } - - // The initial guess at the parameters for the function are written into // the vector inputParameterList. If the paramMask is not NULL, then those - // parameters are masked out. How can this possibly work? The user- - // supplied function will require a fixed number of parameters. + // parameters are masked out. xInit = (double *) psAlloc(inputData->count * sizeof(double)); @@ -257,22 +274,29 @@ r = gsl_rng_alloc(type); + // Initialize the main data structure used by the GSL solver. This will + // contain pointers to the function to be minimized, it's derivative + // function, the number of data points, the number of free parameters, + // and the data structures those functions use. + f.f = &gsl_function_f; f.df = &gsl_function_df; f.fdf = &gsl_my_function_fdf; - f.n = n; + f.n = numData; f.p = inputData->count; f.params = &inputData; - // This tells GSL to use the Levenberg-Marquardt algorithm for chi2 - // minimization. - // T = gsl_multifit_fdfsolver_lmsder; - - // Creates an instance of the GSL solver that we will be iterating on. - s = gsl_multifit_fdfsolver_alloc(T, n, inputData->count); - - // Initialize the GSL minimizer. + // Create an instance of the GSL solver that we will be iterating on. + // It will have numData data points and inputData->count parameters. + s = gsl_multifit_fdfsolver_alloc(T, numData, inputData->count); + + // Initialize the GSL minimizer to use function defined by the data + // structure "f" and x.vector as an initial guess for the parameters. gsl_multifit_fdfsolver_set(s, &f, &x.vector); - // print_state(iter, s); + + // Each iteration of the following loop will perform one step in an + // attempt to minimized chi-squared for the function. The loop exits + // either when the change in parameters is small enough, or when the + // maximum number of iterations is reached. do { iter++; @@ -280,7 +304,4 @@ status = gsl_multifit_fdfsolver_iterate(s); - // printf("status = %s\n", gsl_strerror(status)); - // print_state(iter, s); - // If there was a problem, abort. if (status) { @@ -288,15 +309,12 @@ } - // Checks whether the (L2-norm) computed derivative and the difference - // between the real/actual for that test x-vector. If were close - // enough exit loop. - // Test if the parameters changed by a small enough amount. status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4); - } while (status == GSL_CONTINUE && iter < 500); - - - - + } while (status == GSL_CONTINUE && iter < MAX_LMM_NUM_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; @@ -314,13 +332,35 @@ } + // Calculate the chi-squared for the derived solution. *chiSq = gsl_blas_dnrm2(s->f); + + // Free all allocated memory + // NOTE: Free x. gsl_multifit_fdfsolver_free(s); - // Free all allocated memory psFree(xInit); psFree(inputData); + // Bye bye. return(initialGuess); } + + + + + + + + + + + + + + +// NOTE: rewrite so there is no maximum order for the polynomials. +#define MAX_POLY_ORDER 10 +#define MAX_POLYNOMIAL_TERMS (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2) +int MyInfoLevel = 0; /** @brief This procedure calculates various combinations of powers of x and y * and stores them in the data structure sums[][]. After it completes: Index: /trunk/psLib/src/dataManip/psMinimize.h =================================================================== --- /trunk/psLib/src/dataManip/psMinimize.h (revision 1092) +++ /trunk/psLib/src/dataManip/psMinimize.h (revision 1093) @@ -12,29 +12,33 @@ */ -/** This routine must ninimize a particular non-linear function */ -psFloatArray * -psMinimize(float (*myFunction)(const psFloatArray *restrict), ///< Function to minimize - psFloatArray *restrict initialGuess, ///< Initial guess - psIntArray *restrict paramMask //!< 1 = fit for parameter, 0 = hold parameter constant +/** This routine must minimize a non-linear function */ +psVector * +psMinimize(float (*myFunction)(const psVector *restrict), ///< Function to minimize + psVector *restrict initialGuess, ///< Initial guess + psVector *restrict paramMask //!< 1 = fit for parameter, 0 = hold parameter constant ); /** Minimize chi^2 for input data */ -psFloatArray * -psMinimizeChi2(float (*evalModel)(const psFloatArray *restrict, - const psFloatArray *restrict), ///< Model to fit; (domain and params) - const psFloatArray *restrict domain, ///< The domain values for the corresponding measurements - const psFloatArray *restrict data, ///< Data to fit - const psFloatArray *restrict errors, ///< Errors in the data - psFloatArray *restrict initialGuess, ///< Initial guess - const psIntArray *restrict paramMask ///< 1 = fit for parameter, 0 = hold parameter constant +psVector * +psMinimizeChi2(float (*evalModel)(const psVector *restrict, + const psVector *restrict), ///< Model to fit; (domain and params) + float (*DevalModel)(const psVector *restrict, + const psVector *restrict, + int), ///< Derivative of model to fit; (domain and params) + const psImage *restrict domain, ///< The domain values for the corresponding measurements + const psVector *restrict data, ///< Data to fit + const psVector *restrict errors, ///< Errors in the data + psVector *restrict initialGuess, ///< Initial guess + const psVector *restrict paramMask, ///< 1 = fit for parameter, 0 = hold parameter constant + float *chiSq ); -/** Derive a polynomial fit by chi^2 minimisation --- can be done analytically */ +/** Derive a polynomial fit by chi^2 minimisation (analytically) */ psPolynomial1D * psGetArrayPolynomial(psPolynomial1D myPoly, ///< Polynomial to fit - const psFloatArray *restrict x, ///< Ordinates (or NULL to just use the indices) - const psFloatArray *restrict y, ///< Coordinates - const psFloatArray *restrict yErr ///< Errors in coordinates, or NULL + const psVector *restrict x, ///< Ordinates (or NULL to just use the indices) + const psVector *restrict y, ///< Coordinates + const psVector *restrict yErr ///< Errors in coordinates, or NULL ); Index: /trunk/psLib/src/math/psMinimize.c =================================================================== --- /trunk/psLib/src/math/psMinimize.c (revision 1092) +++ /trunk/psLib/src/math/psMinimize.c (revision 1093) @@ -19,30 +19,28 @@ #include "psFunctions.h" #include "psSort.h" +#include "psMinimize.h" #include "float.h" #include -// NOTE: rewrite so there is no maximum order for the polynomials. -#define MAX_POLY_ORDER 10 -#define MAX_POLYNOMIAL_TERMS (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2) -int MyInfoLevel = 0; - + +#define MAX_LMM_NUM_ITERATIONS 500 typedef struct { size_t n; // Number of data points points in domain. int count; // Number of non-masked parameters. - psVector *restrict initialGuess; - psImage *domain; // - psVector *data; - psVector *errors; - psVector *paramMask; - // The first argument to evalModel() and - // d_evalModel() specifies the data point. - // It must have the same size as the second - // dimension of *domain. - // The second argument must have the same size - // as *initialGuess and *paramMask. + psVector *restrict initialGuess; + const psImage *restrict domain; + const psVector *restrict data; + const psVector *restrict errors; + const psVector *restrict paramMask; float (*evalModel) (psVector *, psVector *); - float (*d_evalModel) (psVector *, psVector *); + float (*d_evalModel) (psVector *, psVector *, int); } psModelData; +// The first argument to evalModel() and +// d_evalModel() specifies the data point. +// It must have the same size as the second +// dimension of *domain. +// The second argument must have the same size +// as *initialGuess and *paramMask. @@ -65,10 +63,10 @@ for the user supplied function which is to be minimized. The GSL minimization routines have no knowledge of the user-supplied function. They -only clal this routine, which then calls the user-supplied function. The +only call this routine, which then calls the user-supplied function. The arguments are: x: These are the parameter which are to be varied by GSL in order to minimized chi2 over the data set. params: - f + f: @@ -79,16 +77,20 @@ gsl_vector *f) { - psImage *domain = ((psModelData *)params)->domain; - psVector *data = ((psModelData *)params)->data; - psVector *errors = ((psModelData *) params)->errors; - psVector *mask = ((psModelData *) params)->paramMask; - psVector *initialGuess = ((psModelData *)params)->initialGuess; + int i; // Loop index variable. + int j; // Loop index variable. + float tmpf; // Temporary floating point variable. + const psImage *restrict domain = ((psModelData *)params)->domain; + const psVector *restrict data = ((psModelData *)params)->data; + const psVector *restrict errors = ((psModelData *) params)->errors; + const psVector *restrict mask = ((psModelData *) params)->paramMask; + psVector *restrict initialGuess = ((psModelData *)params)->initialGuess; + float (*evalModel)(psVector *, psVector *) = ((psModelData *) params)->evalModel; psVector *inputParameterList = NULL; - // float (*evalModel)(psVector *, psVector *) = ((psModelData *) params)->evalModel; - size_t i; - int j; - float tmpf; - - // The GSL routines will call this functions with the masked parameters + psVector *tmpVecPtr = NULL; + + + tmpVecPtr = psVectorAlloc(domain->numCols, PS_TYPE_F32); + + // The GSL routines will call this function with the masked parameters // removed. However, the user-supplied function (to be modified) does not // have those parameters removed. Here will create a new parameter list @@ -112,8 +114,13 @@ for (i=0;inumRows;i++) { - //GUS tmpf = evalModel(inputParameterList, domain->data.F32[i]); + for (j=0;jn;j++) { + tmpVecPtr->data.F32[j] = domain->data.F32[i][j]; + } + tmpf = evalModel(inputParameterList, tmpVecPtr); gsl_vector_set(f, i, (tmpf - data->data.F32[i])/errors->data.F32[i]); } + psFree(inputParameterList); + psFree(tmpVecPtr); return GSL_SUCCESS; } @@ -123,16 +130,17 @@ gsl_matrix *J) { - // size_t n = ((psModelData *)params)->n; - psImage *domain = ((psModelData *)params)->domain; - // psVector *data = ((psModelData *)params)->data; - psVector *errors = ((psModelData *) params)->errors; - psVector *mask = ((psModelData *) params)->paramMask; - psVector *initialGuess = ((psModelData *)params)->initialGuess; + const psImage *restrict domain = ((psModelData *)params)->domain; + const psVector *restrict errors = ((psModelData *) params)->errors; + const psVector *restrict mask = ((psModelData *) params)->paramMask; + psVector *restrict initialGuess = ((psModelData *)params)->initialGuess; psVector *inputParameterList = NULL; - // int count = 0; - // float (*d_evalModel) = ((psModelData *) params)->evalModel; + psVector *tmpVecPtr = NULL; + float (*d_evalModel)(psVector *, psVector *, int) = ((psModelData *) params)->d_evalModel; + size_t i; int j; float tmpf; + + tmpVecPtr = psVectorAlloc(domain->numCols, PS_TYPE_F32); // The GSL routines will call this functions with the masked parameters @@ -157,11 +165,16 @@ for (i=0;inumRows;i++) { + for (j=0;jn;j++) { + tmpVecPtr->data.F32[j] = domain->data.F32[i][j]; + } + for (j=0;jn;j++) { - //GUS tmpf = d_evalModel(inputParameterList, domain->data.F32[i], j); - + tmpf = d_evalModel(inputParameterList, tmpVecPtr, j); gsl_matrix_set(J, i, j, (tmpf/errors->data.F32[i])); } } + psFree(inputParameterList); + psFree(tmpVecPtr); return GSL_SUCCESS; } @@ -185,6 +198,6 @@ psVector * psMinimizeChi2(float (*evalModel)(const psVector *restrict, const psVector *restrict), - float (*DevalModel)(psVector *, psVector *), - const psVector *restrict domain, + float (*DevalModel)(const psVector *restrict, const psVector *restrict, int), + const psImage *restrict domain, const psVector *restrict data, const psVector *restrict errors, @@ -193,20 +206,27 @@ float *chiSq) { - size_t n = domain->n; + + int numData = domain->numCols; // Number of data points + gsl_multifit_fdfsolver *s = NULL; // GSL data structure. + int status; // Return status for the GSL solver. + int i = 0; // Loop index variable. + int j = 0; // Loop index variable. + int iter = 0; // Iteration counter. + psModelData *inputData = NULL; // Contains data that the user-supplied + // function/derivate need. + gsl_multifit_function_fdf f; // GSL structure that contains the + // functions/derivative to be solved. + double *xInit = NULL; // The initial guess at the parameters + // with masked parameters removed. const gsl_multifit_fdfsolver_type *T = gsl_multifit_fdfsolver_lmsder; - gsl_multifit_fdfsolver *s = NULL; - int status; - size_t i = 0; - int j = 0; - size_t iter = 0; - psModelData *inputData = NULL; - gsl_multifit_function_fdf f; - double *xInit = NULL; + // This tells GSL to use the Levenberg- + // Marquardt algorithm for chi2 + // minimization. inputData = (psModelData *) psAlloc(sizeof(psModelData)); - inputData->n = n; + inputData->n = numData; inputData->count = 0; inputData->initialGuess = initialGuess; - // inputData->domain = domain; + inputData->domain = domain; // inputData->data = data; // inputData->errors = errors; @@ -226,10 +246,7 @@ } - - // The initial guess at the parameters for the function are written into // the vector inputParameterList. If the paramMask is not NULL, then those - // parameters are masked out. How can this possibly work? The user- - // supplied function will require a fixed number of parameters. + // parameters are masked out. xInit = (double *) psAlloc(inputData->count * sizeof(double)); @@ -257,22 +274,29 @@ r = gsl_rng_alloc(type); + // Initialize the main data structure used by the GSL solver. This will + // contain pointers to the function to be minimized, it's derivative + // function, the number of data points, the number of free parameters, + // and the data structures those functions use. + f.f = &gsl_function_f; f.df = &gsl_function_df; f.fdf = &gsl_my_function_fdf; - f.n = n; + f.n = numData; f.p = inputData->count; f.params = &inputData; - // This tells GSL to use the Levenberg-Marquardt algorithm for chi2 - // minimization. - // T = gsl_multifit_fdfsolver_lmsder; - - // Creates an instance of the GSL solver that we will be iterating on. - s = gsl_multifit_fdfsolver_alloc(T, n, inputData->count); - - // Initialize the GSL minimizer. + // Create an instance of the GSL solver that we will be iterating on. + // It will have numData data points and inputData->count parameters. + s = gsl_multifit_fdfsolver_alloc(T, numData, inputData->count); + + // Initialize the GSL minimizer to use function defined by the data + // structure "f" and x.vector as an initial guess for the parameters. gsl_multifit_fdfsolver_set(s, &f, &x.vector); - // print_state(iter, s); + + // Each iteration of the following loop will perform one step in an + // attempt to minimized chi-squared for the function. The loop exits + // either when the change in parameters is small enough, or when the + // maximum number of iterations is reached. do { iter++; @@ -280,7 +304,4 @@ status = gsl_multifit_fdfsolver_iterate(s); - // printf("status = %s\n", gsl_strerror(status)); - // print_state(iter, s); - // If there was a problem, abort. if (status) { @@ -288,15 +309,12 @@ } - // Checks whether the (L2-norm) computed derivative and the difference - // between the real/actual for that test x-vector. If were close - // enough exit loop. - // Test if the parameters changed by a small enough amount. status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4); - } while (status == GSL_CONTINUE && iter < 500); - - - - + } while (status == GSL_CONTINUE && iter < MAX_LMM_NUM_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; @@ -314,13 +332,35 @@ } + // Calculate the chi-squared for the derived solution. *chiSq = gsl_blas_dnrm2(s->f); + + // Free all allocated memory + // NOTE: Free x. gsl_multifit_fdfsolver_free(s); - // Free all allocated memory psFree(xInit); psFree(inputData); + // Bye bye. return(initialGuess); } + + + + + + + + + + + + + + +// NOTE: rewrite so there is no maximum order for the polynomials. +#define MAX_POLY_ORDER 10 +#define MAX_POLYNOMIAL_TERMS (((MAX_POLY_ORDER+1) * (MAX_POLY_ORDER + 2)) / 2) +int MyInfoLevel = 0; /** @brief This procedure calculates various combinations of powers of x and y * and stores them in the data structure sums[][]. After it completes: Index: /trunk/psLib/src/math/psMinimize.h =================================================================== --- /trunk/psLib/src/math/psMinimize.h (revision 1092) +++ /trunk/psLib/src/math/psMinimize.h (revision 1093) @@ -12,29 +12,33 @@ */ -/** This routine must ninimize a particular non-linear function */ -psFloatArray * -psMinimize(float (*myFunction)(const psFloatArray *restrict), ///< Function to minimize - psFloatArray *restrict initialGuess, ///< Initial guess - psIntArray *restrict paramMask //!< 1 = fit for parameter, 0 = hold parameter constant +/** This routine must minimize a non-linear function */ +psVector * +psMinimize(float (*myFunction)(const psVector *restrict), ///< Function to minimize + psVector *restrict initialGuess, ///< Initial guess + psVector *restrict paramMask //!< 1 = fit for parameter, 0 = hold parameter constant ); /** Minimize chi^2 for input data */ -psFloatArray * -psMinimizeChi2(float (*evalModel)(const psFloatArray *restrict, - const psFloatArray *restrict), ///< Model to fit; (domain and params) - const psFloatArray *restrict domain, ///< The domain values for the corresponding measurements - const psFloatArray *restrict data, ///< Data to fit - const psFloatArray *restrict errors, ///< Errors in the data - psFloatArray *restrict initialGuess, ///< Initial guess - const psIntArray *restrict paramMask ///< 1 = fit for parameter, 0 = hold parameter constant +psVector * +psMinimizeChi2(float (*evalModel)(const psVector *restrict, + const psVector *restrict), ///< Model to fit; (domain and params) + float (*DevalModel)(const psVector *restrict, + const psVector *restrict, + int), ///< Derivative of model to fit; (domain and params) + const psImage *restrict domain, ///< The domain values for the corresponding measurements + const psVector *restrict data, ///< Data to fit + const psVector *restrict errors, ///< Errors in the data + psVector *restrict initialGuess, ///< Initial guess + const psVector *restrict paramMask, ///< 1 = fit for parameter, 0 = hold parameter constant + float *chiSq ); -/** Derive a polynomial fit by chi^2 minimisation --- can be done analytically */ +/** Derive a polynomial fit by chi^2 minimisation (analytically) */ psPolynomial1D * psGetArrayPolynomial(psPolynomial1D myPoly, ///< Polynomial to fit - const psFloatArray *restrict x, ///< Ordinates (or NULL to just use the indices) - const psFloatArray *restrict y, ///< Coordinates - const psFloatArray *restrict yErr ///< Errors in coordinates, or NULL + const psVector *restrict x, ///< Ordinates (or NULL to just use the indices) + const psVector *restrict y, ///< Coordinates + const psVector *restrict yErr ///< Errors in coordinates, or NULL );