Index: trunk/psLib/src/dataManip/psMinimize.c =================================================================== --- trunk/psLib/src/dataManip/psMinimize.c (revision 1406) +++ trunk/psLib/src/dataManip/psMinimize.c (revision 1407) @@ -1,2 +1,3 @@ + /** @file psMinimize.c * \brief basic minimization functions @@ -9,11 +10,14 @@ * @author George Gusciora, MHPCC * - * @version $Revision: 1.27 $ $Name: not supported by cvs2svn $ - * @date $Date: 2004-08-06 22:34:05 $ + * @version $Revision: 1.28 $ $Name: not supported by cvs2svn $ + * @date $Date: 2004-08-07 00:06:06 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii */ -/*****************************************************************************/ + +/*****************************************************************************/ + /* INCLUDE FILES */ + /*****************************************************************************/ #include @@ -40,6 +44,9 @@ #include "psMinimize.h" #include "psMatrix.h" -/*****************************************************************************/ + +/*****************************************************************************/ + /* DEFINE STATEMENTS */ + /*****************************************************************************/ #define MAX_LMM_ITERATIONS 100 @@ -83,16 +90,17 @@ } - -/*****************************************************************************/ +/*****************************************************************************/ + /* TYPE DEFINITIONS */ + /*****************************************************************************/ typedef struct { - size_t n; // Number of data points points in domain. - int paramCount; // Number of non-masked parameters. - psVector *restrict initialGuess; - const psImage *restrict domain; + size_t n; // Number of data points points in domain. + int paramCount; // Number of non-masked parameters. + psVector *restrict initialGuess; + const psImage *restrict domain; const psVector *restrict data; - const psVector *restrict errors; + const psVector *restrict errors; const psVector *restrict paramMask; float (*evalModel) (const psVector *, const psVector *); @@ -103,7 +111,7 @@ typedef struct { - int paramCount; // Number of non-masked parameters. - psVector *restrict initialGuess; - const psVector *restrict coord; + int paramCount; // Number of non-masked parameters. + psVector *restrict initialGuess; + const psVector *restrict coord; const psVector *restrict paramMask; float (*evalModel) (const psVector *, const psVector *); @@ -113,5 +121,7 @@ /*****************************************************************************/ + /* GLOBAL VARIABLES */ + /*****************************************************************************/ @@ -119,5 +129,7 @@ /*****************************************************************************/ + /* FILE STATIC VARIABLES */ + /*****************************************************************************/ @@ -125,5 +137,7 @@ /*****************************************************************************/ + /* FUNCTION IMPLEMENTATION - LOCAL */ + /*****************************************************************************/ @@ -143,15 +157,13 @@ guess at the parameters, an option parameter mask, etc. *****************************************************************************/ -double p_psMinFunc(const gsl_vector *params, - void *funcData) -{ - int i; // Loop index variable. - int j; // Loop index variable. - float tmpf; // Temporary floating point variable. - const psVector *restrict coord = ((psMinimizeData *) funcData)->coord; - const psVector *restrict mask = ((psMinimizeData *) funcData)->paramMask; - psVector *restrict initialGuess = ((psMinimizeData *)funcData)->initialGuess; - float (*evalModel)(const psVector *, const psVector *) = - ((psMinimizeData *) funcData)->evalModel; +double p_psMinFunc(const gsl_vector * params, void *funcData) +{ + int i; // Loop index variable. + int j; // Loop index variable. + float tmpf; // Temporary floating point variable. + const psVector *restrict coord = ((psMinimizeData *) funcData)->coord; + const psVector *restrict mask = ((psMinimizeData *) funcData)->paramMask; + psVector *restrict initialGuess = ((psMinimizeData *) funcData)->initialGuess; + float (*evalModel) (const psVector *, const psVector *) = ((psMinimizeData *) funcData)->evalModel; psVector *inputParameterList = NULL; @@ -163,5 +175,5 @@ if (mask != NULL) { j = 0; - for (i=0;in;i++) { + for (i = 0; i < mask->n; i++) { if (mask->data.U8[i] != 0) { inputParameterList->data.F32[i] = initialGuess->data.F32[i]; @@ -171,5 +183,5 @@ } } else { - for (i=0;in;i++) { + for (i = 0; i < initialGuess->n; i++) { inputParameterList->data.F32[i] = gsl_vector_get(params, i); } @@ -181,5 +193,5 @@ // Free allocated memory and return the value of the function. psFree(inputParameterList); - return(tmpf); + return (tmpf); } @@ -199,15 +211,13 @@ in "params" and return those derivatives in this psVector. *****************************************************************************/ -void p_psMinFuncDeriv(const gsl_vector *params, - void *funcData, - gsl_vector *df) -{ - int i; // Loop index variable. - int j; // Loop index variable. - float tmpf; // Temporary floating point variable. - const psVector *restrict coord = ((psMinimizeData *) funcData)->coord; - const psVector *restrict mask = ((psMinimizeData *) funcData)->paramMask; - psVector *restrict initialGuess = ((psMinimizeData *)funcData)->initialGuess; - float (*d_evalModel)(const psVector *, const psVector *, int) = +void p_psMinFuncDeriv(const gsl_vector * params, void *funcData, gsl_vector * df) +{ + int i; // Loop index variable. + int j; // Loop index variable. + float tmpf; // Temporary floating point variable. + const psVector *restrict coord = ((psMinimizeData *) funcData)->coord; + const psVector *restrict mask = ((psMinimizeData *) funcData)->paramMask; + psVector *restrict initialGuess = ((psMinimizeData *) funcData)->initialGuess; + float (*d_evalModel) (const psVector *, const psVector *, int) = ((psMinimizeData *) funcData)->d_evalModel; psVector *inputParameterList = NULL; @@ -220,5 +230,5 @@ if (mask != NULL) { j = 0; - for (i=0;in;i++) { + for (i = 0; i < mask->n; i++) { if (mask->data.U8[i] != 0) { inputParameterList->data.F32[i] = initialGuess->data.F32[i]; @@ -228,5 +238,5 @@ } } else { - for (i=0;in;i++) { + for (i = 0; i < initialGuess->n; i++) { inputParameterList->data.F32[i] = gsl_vector_get(params, i); } @@ -235,5 +245,5 @@ // Evaluate the derivative w.r.t. each parameter. // NOTE: we can probably remove the calls for masked parameters. - for (i=0;in;i++) { + for (i = 0; i < initialGuess->n; i++) { tmpf = d_evalModel(inputParameterList, coord, i); gsl_vector_set(df, i, tmpf); @@ -247,8 +257,5 @@ Compute both p_psMinFunc and p_psMinFuncDeriv together. *****************************************************************************/ -void p_psMinFuncFuncDeriv(const gsl_vector *params, - void *funcData, - double *f, - gsl_vector *df) +void p_psMinFuncFuncDeriv(const gsl_vector * params, void *funcData, double *f, gsl_vector * df) { *f = p_psMinFunc(params, funcData); @@ -281,17 +288,15 @@ expected value and divide by the error. *****************************************************************************/ -int p_psMinChi2Func(const gsl_vector *params, - void *funcData, - gsl_vector *outData) -{ - int i; // Loop index variable. - int j; // Loop index variable. - float tmpf; // Temporary floating point variable. - const psImage *restrict domain = ((psMinChi2Data *)funcData)->domain; - const psVector *restrict data = ((psMinChi2Data *)funcData)->data; - const psVector *restrict errors = ((psMinChi2Data *) funcData)->errors; - const psVector *restrict mask = ((psMinChi2Data *) funcData)->paramMask; - psVector *restrict initialGuess = ((psMinChi2Data *)funcData)->initialGuess; - float (*evalModel)(const psVector *, const psVector *) = ((psMinChi2Data *) funcData)->evalModel; +int p_psMinChi2Func(const gsl_vector * params, void *funcData, gsl_vector * outData) +{ + int i; // Loop index variable. + int j; // Loop index variable. + float tmpf; // Temporary floating point variable. + const psImage *restrict domain = ((psMinChi2Data *) funcData)->domain; + const psVector *restrict data = ((psMinChi2Data *) funcData)->data; + const psVector *restrict errors = ((psMinChi2Data *) funcData)->errors; + const psVector *restrict mask = ((psMinChi2Data *) funcData)->paramMask; + psVector *restrict initialGuess = ((psMinChi2Data *) funcData)->initialGuess; + float (*evalModel) (const psVector *, const psVector *) = ((psMinChi2Data *) funcData)->evalModel; psVector *inputParameterList = NULL; psVector *tmpVecPtr = NULL; @@ -307,5 +312,5 @@ if (mask != NULL) { j = 0; - for (i=0;in;i++) { + for (i = 0; i < mask->n; i++) { if (mask->data.U8[i] != 0) { inputParameterList->data.F32[i] = initialGuess->data.F32[i]; @@ -315,5 +320,5 @@ } } else { - for (i=0;in;i++) { + for (i = 0; i < initialGuess->n; i++) { inputParameterList->data.F32[i] = gsl_vector_get(params, i); } @@ -321,12 +326,11 @@ // Evaluate the function at each data point. - for (i=0;inumRows;i++) { - for (j=0;jnumCols;j++) { + for (i = 0; i < domain->numRows; i++) { + for (j = 0; j < domain->numCols; j++) { tmpVecPtr->data.F32[j] = domain->data.F32[i][j]; } tmpf = evalModel(tmpVecPtr, inputParameterList); - gsl_vector_set(outData, i, (tmpf - data->data.F32[i])/ - errors->data.F32[i]); + gsl_vector_set(outData, i, (tmpf - data->data.F32[i]) / errors->data.F32[i]); } @@ -354,15 +358,14 @@ and returned in this data structure. *****************************************************************************/ -int p_psMinChi2FuncDeriv(const gsl_vector *params, - void *funcData, - gsl_matrix *J) -{ - const psImage *restrict domain = ((psMinChi2Data *)funcData)->domain; - const psVector *restrict errors = ((psMinChi2Data *) funcData)->errors; - const psVector *restrict mask = ((psMinChi2Data *) funcData)->paramMask; - psVector *restrict initialGuess = ((psMinChi2Data *)funcData)->initialGuess; +int p_psMinChi2FuncDeriv(const gsl_vector * params, void *funcData, gsl_matrix * J) +{ + const psImage *restrict domain = ((psMinChi2Data *) funcData)->domain; + const psVector *restrict errors = ((psMinChi2Data *) funcData)->errors; + const psVector *restrict mask = ((psMinChi2Data *) funcData)->paramMask; + psVector *restrict initialGuess = ((psMinChi2Data *) funcData)->initialGuess; psVector *inputParameterList = NULL; psVector *tmpVecPtr = NULL; - float (*d_evalModel)(const psVector *, const psVector *, int) = ((psMinChi2Data *) funcData)->d_evalModel; + float (*d_evalModel) (const psVector *, const psVector *, int) = + ((psMinChi2Data *) funcData)->d_evalModel; size_t i; @@ -380,5 +383,5 @@ if (mask != NULL) { j = 0; - for (i=0;in;i++) { + for (i = 0; i < mask->n; i++) { if (mask->data.U8[i] != 0) { inputParameterList->data.F32[i] = initialGuess->data.F32[i]; @@ -388,5 +391,5 @@ } } else { - for (i=0;in;i++) { + for (i = 0; i < initialGuess->n; i++) { inputParameterList->data.F32[i] = gsl_vector_get(params, i); } @@ -394,12 +397,12 @@ // Evaluate the derivtaive at each data point, and w.r.t. each parameter. - for (i=0;inumRows;i++) { - for (j=0;jn;j++) { + for (i = 0; i < domain->numRows; i++) { + for (j = 0; j < tmpVecPtr->n; j++) { tmpVecPtr->data.F32[j] = domain->data.F32[i][j]; } - for (j=0;jn;j++) { + for (j = 0; j < inputParameterList->n; j++) { tmpf = d_evalModel(tmpVecPtr, inputParameterList, j); - gsl_matrix_set(J, i, j, (tmpf/errors->data.F32[i])); + gsl_matrix_set(J, i, j, (tmpf / errors->data.F32[i])); } } @@ -410,9 +413,5 @@ } - -int p_psMinChi2FuncFuncDeriv(const gsl_vector *params, - void *funcData, - gsl_vector *f, - gsl_matrix *J) +int p_psMinChi2FuncFuncDeriv(const gsl_vector * params, void *funcData, gsl_vector * f, gsl_matrix * J) { p_psMinChi2Func(params, funcData, f); @@ -421,5 +420,4 @@ return GSL_SUCCESS; } - /****************************************************************************** @@ -428,18 +426,15 @@ returned as a psVector sums. *****************************************************************************/ -void p_psBuildSums1D(double x, - int polyOrder, - psVector *sums) -{ - int i = 0; - double xSum = 0.0; +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++) { + for (i = 0; i <= polyOrder; i++) { sums->data.F64[i] = xSum; - xSum*= x; - } -} - + xSum *= x; + } +} /****************************************************************************** @@ -448,17 +443,15 @@ *****************************************************************************/ 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); + for (i = 0; i < x; i++) { + imageScalingFactors->data.F32[i] = (((float)2 * i) / ((float)x)) - 1.0; + } + + return (imageScalingFactors); } @@ -481,41 +474,39 @@ *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 +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 + 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 + 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. + 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)); + 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++) { + 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++) { + // lubksb(A, N, indx, col); + for (i = 1; i <= N; i++) { y[i][j] = col[i]; } @@ -527,9 +518,9 @@ // terms and check if they are consistent with zero. - numPolyTerms = (((polyOrder+1) * (polyOrder + 2)) / 2); + numPolyTerms = (((polyOrder + 1) * (polyOrder + 2)) / 2); lastTerm = numPolyTerms + 1; firstTerm = lastTerm - polyOrder; *flag = 1; - for (i=firstTerm; i<=lastTerm; i++) { + for (i = firstTerm; i <= lastTerm; i++) { #ifdef DARWIN error[i] = (float)sqrt(y[i][i]); @@ -539,6 +530,5 @@ #endif - if (!((B[i] <= (2.0f * error[i])) && - ((-2.0f * error[i]) <= B[i]))) { + if (!((B[i] <= (2.0f * error[i])) && ((-2.0f * error[i]) <= B[i]))) { *flag = 0; } @@ -548,5 +538,5 @@ psFree(error); psFree(col); - for(j=1;j<=N;j++) { + for (j = 1; j <= N; j++) { psFree(y[j]); } @@ -554,12 +544,9 @@ } - - - -/*****************************************************************************/ +/*****************************************************************************/ + /* FUNCTION IMPLEMENTATION - PUBLIC */ -/*****************************************************************************/ - - + +/*****************************************************************************/ /****************************************************************************** @@ -569,10 +556,8 @@ 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) +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; @@ -607,5 +592,5 @@ // for the parameters. if (paramMask != NULL) { - for (i=0;in;i++) { + for (i = 0; i < paramMask->n; i++) { if (paramMask->data.U8[i] != 0) { inputData.paramCount++; @@ -613,5 +598,5 @@ } } else { - inputData.paramCount= initialGuess->n; + inputData.paramCount = initialGuess->n; } @@ -622,5 +607,5 @@ if (paramMask != NULL) { j = 0; - for (i=0;in;i++) { + for (i = 0; i < initialGuess->n; i++) { if (paramMask->data.U8[i] == 0) { gsl_vector_set(x, j++, initialGuess->data.F32[i]); @@ -628,5 +613,5 @@ } } else { - for (i=0;in;i++) { + for (i = 0; i < initialGuess->n; i++) { gsl_vector_set(x, i, initialGuess->data.F32[i]); } @@ -651,5 +636,5 @@ if (status == GSL_SUCCESS) - printf ("Minimum found at:\n"); + printf("Minimum found at:\n"); } while (status == GSL_CONTINUE && iter < MAX_MINIMIZE_ITERATIONS); @@ -660,5 +645,5 @@ if (paramMask != NULL) { j = 0; - for (i=0;in;i++) { + for (i = 0; i < initialGuess->n; i++) { if (paramMask->data.U8[i] == 0) { initialGuess->data.F32[i] = gsl_vector_get(s->x, j++); @@ -668,15 +653,10 @@ } } else { - for (i=0;in;i++) { + for (i = 0; i < initialGuess->n; i++) { initialGuess->data.F32[i] = gsl_vector_get(s->x, i); } } - return(initialGuess); -} - - - - - + return (initialGuess); +} /****************************************************************************** @@ -684,28 +664,28 @@ such that they best fit the supplied data points. *****************************************************************************/ -psVector * -psMinimizeChi2(float (*evalModel)(const psVector *restrict, const psVector *restrict), - float (*DevalModel)(const psVector *restrict, const psVector *restrict, int), - const psImage *restrict domain, - const psVector *restrict data, - const psVector *restrict errors, - psVector *restrict initialGuess, - const psVector *restrict paramMask, - float *chiSq) -{ - int numData = domain->numRows; // Number of data points - 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. - gsl_multifit_function_fdf f; // GSL structure that contains the +psVector *psMinimizeChi2(float (*evalModel) (const psVector * restrict, const psVector * restrict), + float (*DevalModel) (const psVector * restrict, const psVector * restrict, int), + const psImage * restrict domain, + const psVector * restrict data, + const psVector * restrict errors, + psVector * restrict initialGuess, const psVector * restrict paramMask, float *chiSq) +{ + int numData = domain->numRows; // Number of data points + 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. + 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; + // This tells GSL to use the Levenberg- // Marquardt algorithm for chi2 // minimization. - gsl_multifit_fdfsolver *s; // GSL data structure. + gsl_multifit_fdfsolver *s; // GSL data structure. psMinChi2Data inputData; float chiSqOld = 0.0; @@ -721,5 +701,5 @@ PS_CHECK_VECTOR_SIZE_EQUAL(data, errors); if (domain->numRows != data->n) { - psAbort(__func__,"Number of data points and data values not equal."); + psAbort(__func__, "Number of data points and data values not equal."); } if (paramMask != NULL) { @@ -743,5 +723,5 @@ // for the parameters. if (paramMask != NULL) { - for (i=0;in;i++) { + for (i = 0; i < paramMask->n; i++) { if (paramMask->data.U8[i] != 0) { inputData.paramCount++; @@ -749,5 +729,5 @@ } } else { - inputData.paramCount= initialGuess->n; + inputData.paramCount = initialGuess->n; } @@ -755,8 +735,8 @@ // the vector inputParameterList. If the paramMask is not NULL, then those // parameters are masked out. - xInit = (double *) psAlloc(inputData.paramCount * sizeof(double)); + xInit = (double *)psAlloc(inputData.paramCount * sizeof(double)); if (paramMask != NULL) { j = 0; - for (i=0;in;i++) { + for (i = 0; i < initialGuess->n; i++) { if (paramMask->data.U8[i] == 0) { xInit[j++] = initialGuess->data.F32[i]; @@ -764,5 +744,5 @@ } } else { - for (i=0;in;i++) { + for (i = 0; i < initialGuess->n; i++) { xInit[i] = initialGuess->data.F32[i]; } @@ -771,4 +751,5 @@ const gsl_rng_type *type; gsl_rng *r; + gsl_rng_env_setup(); @@ -789,4 +770,5 @@ gsl_vector_view x = gsl_vector_view_array(xInit, inputData.paramCount); + T = gsl_multifit_fdfsolver_lmsder; s = gsl_multifit_fdfsolver_alloc(T, numData, inputData.paramCount); @@ -796,5 +778,5 @@ do { iter++; - for (i=0;in;i++) { + for (i = 0; i < initialGuess->n; i++) { printf("Iteration %d: parameter %d is %.3f\n", iter, i, gsl_vector_get(s->x, i)); } @@ -802,5 +784,5 @@ status = gsl_multifit_fdfsolver_iterate(s); printf("gsl_multifit_fdfsolver_iterate() status is %s\n", gsl_strerror(status)); - for (i=0;in;i++) { + for (i = 0; i < initialGuess->n; i++) { printf("Iteration %d: parameter %d is %.3f\n", iter, i, gsl_vector_get(s->x, i)); } @@ -810,9 +792,8 @@ psAbort(__func__, "gsl_multifit_fdfsolver_iterate(%s)\n", gsl_strerror(status)); } - // Test if the parameters changed by a small enough amount. // NOTE: This wasn't working right when the parameters fit exactly. // Figure out why. - // status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4); + // status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4); // We test for convergence if chiSquared changes by less than 1.0 @@ -829,5 +810,4 @@ } while (status == GSL_CONTINUE && iter < MAX_LMM_ITERATIONS); - // In the above steps we had removed the masked elements from the // the solver. This next code blocks puts those masked elements @@ -835,5 +815,5 @@ if (paramMask != NULL) { j = 0; - for (i=0;in;i++) { + for (i = 0; i < initialGuess->n; i++) { if (paramMask->data.U8[i] == 0) { initialGuess->data.F32[i] = gsl_vector_get(s->x, j++); @@ -843,5 +823,5 @@ } } else { - for (i=0;in;i++) { + for (i = 0; i < initialGuess->n; i++) { initialGuess->data.F32[i] = gsl_vector_get(s->x, i); } @@ -857,7 +837,6 @@ // Bye bye. - return(initialGuess); -} - + return (initialGuess); +} /****************************************************************************** @@ -868,9 +847,7 @@ NOTE: yErr is currently ignored. *****************************************************************************/ -psPolynomial1D * -psVectorFitPolynomial1D(psPolynomial1D *myPoly, - const psVector *restrict x, - const psVector *restrict y, - const psVector *restrict yErr) +psPolynomial1D *psVectorFitPolynomial1D(psPolynomial1D * myPoly, + const psVector * restrict x, + const psVector * restrict y, const psVector * restrict yErr) { int polyOrder = myPoly->n; @@ -879,5 +856,5 @@ psVector *B = NULL; psVector *outPerm = NULL; - psVector *X = NULL; // NOTE: do we need this? + psVector *X = NULL; // NOTE: do we need this? psVector *coeffs = NULL; int i = 0; @@ -886,8 +863,8 @@ psVector *xSums = NULL; - // printf("psVectorFitPolynomial1D()\n"); - // for (i=0;in;i++) { - // printf("(x, y, yErr) is (%f, %f, %f)\n", x->data.F64[i], y->data.F64[i], yErr->data.F64[i]); - // } + // printf("psVectorFitPolynomial1D()\n"); + // for (i=0;in;i++) { + // printf("(x, y, yErr) is (%f, %f, %f)\n", x->data.F64[i], y->data.F64[i], yErr->data.F64[i]); + // } PS_CHECK_NULL_1DPOLY(myPoly); @@ -901,37 +878,37 @@ PS_CHECK_VECTOR_SIZE_EQUAL(y, yErr); - A = psImageAlloc(polyOrder, polyOrder, PS_TYPE_F64); - ALUD = psImageAlloc(polyOrder, polyOrder, PS_TYPE_F64); - B = psVectorAlloc(polyOrder, PS_TYPE_F64); - coeffs = psVectorAlloc(polyOrder, PS_TYPE_F64); - X = psVectorAlloc(x->n, PS_TYPE_F64); + A = psImageAlloc(polyOrder, polyOrder, PS_TYPE_F64); + ALUD = psImageAlloc(polyOrder, polyOrder, PS_TYPE_F64); + B = psVectorAlloc(polyOrder, PS_TYPE_F64); + coeffs = psVectorAlloc(polyOrder, PS_TYPE_F64); + X = psVectorAlloc(x->n, PS_TYPE_F64); outPerm = psVectorAlloc(polyOrder, PS_TYPE_F64); - xSums = psVectorAlloc(1+2*polyOrder, PS_TYPE_F64); + xSums = psVectorAlloc(1 + 2 * polyOrder, PS_TYPE_F64); // Initialize data structures. - for(i=0;i<(polyOrder);i++) { + for (i = 0; i < (polyOrder); i++) { B->data.F64[i] = 0.0; coeffs->data.F64[i] = 0.0; outPerm->data.F64[i] = 0.0; - for(j=0;j<(polyOrder);j++) { + for (j = 0; j < (polyOrder); j++) { A->data.F64[i][j] = 0.0; ALUD->data.F64[i][j] = 0.0; } } - for (i=0;in;i++) { + for (i = 0; i < X->n; i++) { X->data.F64[i] = x->data.F64[i]; } // Build the B and A data structs. - for (i=0;in;i++) { - p_psBuildSums1D(X->data.F64[i], 2*polyOrder, xSums); - - for(k=0;k<(polyOrder);k++) { - B->data.F64[k]+= y->data.F64[i] * xSums->data.F64[k]; - } - - for(k=0;k<(polyOrder);k++) { - for(j=0;j<(polyOrder);j++) { - A->data.F64[k][j]+= xSums->data.F64[k+j]; + for (i = 0; i < X->n; i++) { + p_psBuildSums1D(X->data.F64[i], 2 * polyOrder, xSums); + + for (k = 0; k < (polyOrder); k++) { + B->data.F64[k] += y->data.F64[i] * xSums->data.F64[k]; + } + + for (k = 0; k < (polyOrder); k++) { + for (j = 0; j < (polyOrder); j++) { + A->data.F64[k][j] += xSums->data.F64[k + j]; } } @@ -941,13 +918,13 @@ coeffs = psMatrixLUSolve(coeffs, ALUD, B, outPerm); - for(k=0;k<(polyOrder);k++) { + for (k = 0; k < (polyOrder); k++) { myPoly->coeff[k] = coeffs->data.F64[k]; - // printf("myPoly->coeff[%d] is %f\n", k, myPoly->coeff[k]); - } - - - // for (i=0;in;i++) { - // printf("HMMM: psEvalPolynomial1D(%f) is %f\n", x->data.F64[i], psEvalPolynomial1D(x->data.F64[i], myPoly)); - // } + // printf("myPoly->coeff[%d] is %f\n", k, myPoly->coeff[k]); + } + + // for (i=0;in;i++) { + // printf("HMMM: psEvalPolynomial1D(%f) is %f\n", x->data.F64[i], psEvalPolynomial1D(x->data.F64[i], + // myPoly)); + // } psFree(A); @@ -959,4 +936,4 @@ psFree(xSums); - return(myPoly); -} + return (myPoly); +}