Changeset 1407 for trunk/psLib/src/dataManip/psMinimize.c

Timestamp:

Aug 6, 2004, 2:06:06 PM (22 years ago)

Author:

desonia

Message:

another attempt to get astyle to get it right.

File:

• trunk/psLib/src/dataManip/psMinimize.c (modified) (63 diffs)

trunk/psLib/src/dataManip/psMinimize.c

@@ +1 @@
+
 /** @file  psMinimize.c
  *  \brief basic minimization functions
@@ -9 +10 @@
  *  @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 <stdlib.h>
@@ -40 +44 @@
 #include "psMinimize.h"
 #include "psMatrix.h"
-/*****************************************************************************/
+
+/*****************************************************************************/
+
 /* DEFINE STATEMENTS                                                         */
+
 /*****************************************************************************/
 #define MAX_LMM_ITERATIONS 100
@@ -83 +90 @@
 }
 
-
-/*****************************************************************************/
+/*****************************************************************************/
+
 /* 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 +111 @@
 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 +121 @@
 
 /*****************************************************************************/
+
 /* GLOBAL VARIABLES                                                          */
+
 /*****************************************************************************/
 
@@ -119 +129 @@
 
 /*****************************************************************************/
+
 /* FILE STATIC VARIABLES                                                     */
+
 /*****************************************************************************/
 
@@ -125 +137 @@
 
 /*****************************************************************************/
+
 /* FUNCTION IMPLEMENTATION - LOCAL                                           */
+
 /*****************************************************************************/
 
@@ -143 +157 @@
  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 +175 @@
     if (mask != NULL) {
         j = 0;
-        for (i=0;i<mask->n;i++) {
+        for (i = 0; i < mask->n; i++) {
             if (mask->data.U8[i] != 0) {
                 inputParameterList->data.F32[i] = initialGuess->data.F32[i];
@@ -171 +183 @@
         }
     } else {
-        for (i=0;i<initialGuess->n;i++) {
+        for (i = 0; i < initialGuess->n; i++) {
             inputParameterList->data.F32[i] = gsl_vector_get(params, i);
         }
@@ -181 +193 @@
     // Free allocated memory and return the value of the function.
     psFree(inputParameterList);
-    return(tmpf);
+    return (tmpf);
 }
 
@@ -199 +211 @@
  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 +230 @@
     if (mask != NULL) {
         j = 0;
-        for (i=0;i<mask->n;i++) {
+        for (i = 0; i < mask->n; i++) {
             if (mask->data.U8[i] != 0) {
                 inputParameterList->data.F32[i] = initialGuess->data.F32[i];
@@ -228 +238 @@
         }
     } else {
-        for (i=0;i<initialGuess->n;i++) {
+        for (i = 0; i < initialGuess->n; i++) {
             inputParameterList->data.F32[i] = gsl_vector_get(params, i);
         }
@@ -235 +245 @@
     // Evaluate the derivative w.r.t. each parameter.
     // NOTE: we can probably remove the calls for masked parameters.
-    for (i=0;i<initialGuess->n;i++) {
+    for (i = 0; i < initialGuess->n; i++) {
         tmpf = d_evalModel(inputParameterList, coord, i);
         gsl_vector_set(df, i, tmpf);
@@ -247 +257 @@
     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 +288 @@
   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 +312 @@
     if (mask != NULL) {
         j = 0;
-        for (i=0;i<mask->n;i++) {
+        for (i = 0; i < mask->n; i++) {
             if (mask->data.U8[i] != 0) {
                 inputParameterList->data.F32[i] = initialGuess->data.F32[i];
@@ -315 +320 @@
         }
     } else {
-        for (i=0;i<initialGuess->n;i++) {
+        for (i = 0; i < initialGuess->n; i++) {
             inputParameterList->data.F32[i] = gsl_vector_get(params, i);
         }
@@ -321 +326 @@
 
     // Evaluate the function at each data point.
-    for (i=0;i<domain->numRows;i++) {
-        for (j=0;j<domain->numCols;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 +358 @@
  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 +383 @@
     if (mask != NULL) {
         j = 0;
-        for (i=0;i<mask->n;i++) {
+        for (i = 0; i < mask->n; i++) {
             if (mask->data.U8[i] != 0) {
                 inputParameterList->data.F32[i] = initialGuess->data.F32[i];
@@ -388 +391 @@
         }
     } else {
-        for (i=0;i<initialGuess->n;i++) {
+        for (i = 0; i < initialGuess->n; i++) {
             inputParameterList->data.F32[i] = gsl_vector_get(params, i);
         }
@@ -394 +397 @@
 
     // Evaluate the derivtaive at each data point, and w.r.t. each parameter.
-    for (i=0;i<domain->numRows;i++) {
-        for (j=0;j<tmpVecPtr->n;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;j<inputParameterList->n;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 +413 @@
 }
 
-
-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 +420 @@
     return GSL_SUCCESS;
 }
-
 
 /******************************************************************************
@@ -428 +426 @@
 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 +443 @@
  *****************************************************************************/
 psVector *psBuildImageScalingFactors(int x)
-
 {
     int i = 0;                  // loop index variable.
     psVector *imageScalingFactors = NULL;
 
-
     imageScalingFactors = psVectorAlloc(x, PS_TYPE_F32);
 
-    for (i=0;i<x;i++) {
-        imageScalingFactors->data.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 +474 @@
     *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 +518 @@
     // 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 +530 @@
         #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 +538 @@
     psFree(error);
     psFree(col);
-    for(j=1;j<=N;j++) {
+    for (j = 1; j <= N; j++) {
         psFree(y[j]);
     }
@@ -554 +544 @@
 }
 
-
-
-
-/*****************************************************************************/
+/*****************************************************************************/
+
 /* FUNCTION IMPLEMENTATION - PUBLIC                                          */
-/*****************************************************************************/
-
-
+
+/*****************************************************************************/
 
 /******************************************************************************
@@ -569 +556 @@
 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 +592 @@
     // for the parameters.
     if (paramMask != NULL) {
-        for (i=0;i<paramMask->n;i++) {
+        for (i = 0; i < paramMask->n; i++) {
             if (paramMask->data.U8[i] != 0) {
                 inputData.paramCount++;
@@ -613 +598 @@
         }
     } else {
-        inputData.paramCount= initialGuess->n;
+        inputData.paramCount = initialGuess->n;
     }
 
@@ -622 +607 @@
     if (paramMask != NULL) {
         j = 0;
-        for (i=0;i<initialGuess->n;i++) {
+        for (i = 0; i < initialGuess->n; i++) {
             if (paramMask->data.U8[i] == 0) {
                 gsl_vector_set(x, j++, initialGuess->data.F32[i]);
@@ -628 +613 @@
         }
     } else {
-        for (i=0;i<initialGuess->n;i++) {
+        for (i = 0; i < initialGuess->n; i++) {
             gsl_vector_set(x, i, initialGuess->data.F32[i]);
         }
@@ -651 +636 @@
 
         if (status == GSL_SUCCESS)
-            printf ("Minimum found at:\n");
+            printf("Minimum found at:\n");
 
     } while (status == GSL_CONTINUE && iter < MAX_MINIMIZE_ITERATIONS);
@@ -660 +645 @@
     if (paramMask != NULL) {
         j = 0;
-        for (i=0;i<initialGuess->n;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 +653 @@
         }
     } else {
-        for (i=0;i<initialGuess->n;i++) {
+        for (i = 0; i < initialGuess->n; i++) {
             initialGuess->data.F32[i] = gsl_vector_get(s->x, i);
         }
     }
-    return(initialGuess);
-}
-
-
-
-
-
+    return (initialGuess);
+}
 
 /******************************************************************************
@@ -684 +664 @@
     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 +701 @@
     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 +723 @@
     // for the parameters.
     if (paramMask != NULL) {
-        for (i=0;i<paramMask->n;i++) {
+        for (i = 0; i < paramMask->n; i++) {
             if (paramMask->data.U8[i] != 0) {
                 inputData.paramCount++;
@@ -749 +729 @@
         }
     } else {
-        inputData.paramCount= initialGuess->n;
+        inputData.paramCount = initialGuess->n;
     }
 
@@ -755 +735 @@
     // 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;i<initialGuess->n;i++) {
+        for (i = 0; i < initialGuess->n; i++) {
             if (paramMask->data.U8[i] == 0) {
                 xInit[j++] = initialGuess->data.F32[i];
@@ -764 +744 @@
         }
     } else {
-        for (i=0;i<initialGuess->n;i++) {
+        for (i = 0; i < initialGuess->n; i++) {
             xInit[i] = initialGuess->data.F32[i];
         }
@@ -771 +751 @@
     const gsl_rng_type *type;
     gsl_rng *r;
+
     gsl_rng_env_setup();
 
@@ -789 +770 @@
 
     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 +778 @@
     do {
         iter++;
-        for (i=0;i<initialGuess->n;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 +784 @@
         status = gsl_multifit_fdfsolver_iterate(s);
         printf("gsl_multifit_fdfsolver_iterate() status is %s\n", gsl_strerror(status));
-        for (i=0;i<initialGuess->n;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 +792 @@
             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 +810 @@
     } 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 +815 @@
     if (paramMask != NULL) {
         j = 0;
-        for (i=0;i<initialGuess->n;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 +823 @@
         }
     } else {
-        for (i=0;i<initialGuess->n;i++) {
+        for (i = 0; i < initialGuess->n; i++) {
             initialGuess->data.F32[i] = gsl_vector_get(s->x, i);
         }
@@ -857 +837 @@
 
     // Bye bye.
-    return(initialGuess);
-}
-
+    return (initialGuess);
+}
 
 /******************************************************************************
@@ -868 +847 @@
 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 +856 @@
     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 +863 @@
     psVector *xSums = NULL;
 
-    //    printf("psVectorFitPolynomial1D()\n");
-    //    for (i=0;i<x->n;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;i<x->n;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 +878 @@
     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;i<X->n;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;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];
+    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 +918 @@
     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;i<x->n;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;i<x->n;i++) {
+    // printf("HMMM: psEvalPolynomial1D(%f) is %f\n", x->data.F64[i], psEvalPolynomial1D(x->data.F64[i],
+    // myPoly));
+    // }
 
     psFree(A);
@@ -959 +936 @@
     psFree(xSums);
 
-    return(myPoly);
-}
+    return (myPoly);
+}