Index: trunk/psLib/src/math/psMinimize.c
===================================================================
--- trunk/psLib/src/math/psMinimize.c	(revision 5113)
+++ trunk/psLib/src/math/psMinimize.c	(revision 5175)
@@ -10,11 +10,12 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.140 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2005-09-23 23:01:30 $
+ *  @version $Revision: 1.141 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-09-29 02:13:54 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
  *
  *  XXX: must follow coding name standards on local functions.
- *
+ *  XXX: put local functions in front.
+ * 
  */
 /*****************************************************************************/
@@ -33,9 +34,6 @@
 /*****************************************************************************/
 /* DEFINE STATEMENTS                                                         */
-/* What are the following macros for?                                        */
 /*****************************************************************************/
-#define PS_SEG psLib
-#define PS_PWD dataManip
-#define PS_FILE psMinimize
+
 /*****************************************************************************/
 /* TYPE DEFINITIONS                                                          */
@@ -44,9 +42,9 @@
 /*****************************************************************************/
 /* GLOBAL VARIABLES                                                          */
+/* XXX: Do these conform to code standard?         */
 /*****************************************************************************/
 static psMinimizeChi2PowellFunc Chi2PowellFunc = NULL;
 static psVector *myValue;
 static psVector *myError;
-
 /*****************************************************************************/
 /* FILE STATIC VARIABLES                                                     */
@@ -64,48 +62,135 @@
  ******************************************************************************
  *****************************************************************************/
-// measure the distance to the minimum assuming a linear model
-// XXX: Move this to the pub area.
-bool psMinimizeGaussNewtonDelta (psVector *delta,
-                                 const psVector *params,
-                                 const psVector *paramMask,
-                                 const psArray  *x,
-                                 const psVector *y,
-                                 const psVector *yErr,
-                                 psMinimizeLMChi2Func func)
+// XXX EAM : can we use static copies of LUv, LUm, A?
+psBool p_psMinLM_GuessABP(
+    psImage  *Alpha,
+    psVector *Beta,
+    psVector *Params,
+    const psImage  *alpha,
+    const psVector *beta,
+    const psVector *params,
+    const psVector *paramMask,
+    const psVector *beta_lim,
+    const psVector *params_min,
+    const psVector *params_max,
+    psF64 lambda)
+{
+    # define USE_LU_DECOMP 1
+    # if (USE_LU_DECOMP)
+        psVector *LUv = NULL;
+    psImage  *LUm = NULL;
+    psImage  *A   = NULL;
+    psF32    det;
+
+    // LU decomposition version
+    psTrace(__func__, 5, "using LUD version\n");
+
+    // set new guess values (creates matrix A)
+    A = psImageCopy(NULL, alpha, PS_TYPE_F64);
+    for (int j = 0; j < params->n; j++) {
+        if ((paramMask != NULL) && (paramMask->data.U8[j]))
+            continue;
+        A->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
+    }
+
+    // solve A*beta = Beta (Alpha = 1/A)
+    // these operations do not modify the input values (creates LUm, LUv)
+    LUm   = psMatrixLUD(NULL, &LUv, A);
+    Beta  = psMatrixLUSolve(Beta, LUm, beta, LUv);
+    Alpha = psMatrixInvert(Alpha, A, &det);
+
+    # else
+        // gauss-jordan version
+        psTrace(__func__, 5, "using Gauss-J version");
+
+    // set new guess values (creates matrix A)
+    Beta = psVectorCopy(Beta, beta, PS_TYPE_F64);
+    Alpha = psImageCopy(Alpha, alpha, PS_TYPE_F64);
+    for (int j = 0; j < params->n; j++) {
+        if ((paramMask != NULL) && (paramMask->data.U8[j]))
+            continue;
+        Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
+    }
+
+    psGaussJordan(Alpha, Beta);
+    psFree(A);
+    psFree(LUm);
+    psFree(LUv);
+    # endif
+
+    // apply Beta to get new Params values
+    for (int j = 0; j < params->n; j++) {
+        if ((paramMask != NULL) && (paramMask->data.U8[j]))
+            continue;
+        // Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j];
+        // compare Beta to beta limits
+        if (beta_lim != NULL) {
+            if (fabs(Beta->data.F64[j]) > fabs(beta_lim->data.F32[j])) {
+                Beta->data.F64[j] = (Beta->data.F64[j] > 0) ? fabs(beta_lim->data.F32[j]) : -fabs(beta_lim->data.F32[j]);
+            }
+        }
+        Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j];
+        // compare new params to param limits
+        if (params_max != NULL) {
+            Params->data.F32[j] = PS_MIN (Params->data.F32[j], params_max->data.F32[j]);
+        }
+        if (params_min != NULL) {
+            Params->data.F32[j] = PS_MAX (Params->data.F32[j], params_min->data.F32[j]);
+        }
+    }
+    # if (USE_LU_DECOMP)
+        psFree(A);
+    psFree(LUm);
+    psFree(LUv);
+    # endif
+
+    return(true);
+}
+
+
+bool psMinimizeGaussNewtonDelta(
+    psVector *delta,
+    const psVector *params,
+    const psVector *paramMask,
+    const psArray  *x,
+    const psVector *y,
+    const psVector *yWt,
+    psMinimizeLMChi2Func func)
 {
 
     // allocate internal arrays (current vs Guess)
-    psImage  *alpha  = psImageAlloc  (params->n, params->n, PS_TYPE_F64);
-    psImage  *Alpha  = psImageAlloc  (params->n, params->n, PS_TYPE_F64);
-    psVector *beta   = psVectorAlloc (params->n, PS_TYPE_F64);
-    psVector *Params = psVectorAlloc (params->n, PS_TYPE_F64);
-    psVector *dy     = psVectorAlloc (y->n, PS_TYPE_F32);
+    psImage  *alpha  = psImageAlloc (params->n, params->n, PS_TYPE_F64);
+    psImage  *Alpha  = psImageAlloc (params->n, params->n, PS_TYPE_F64);
+    psVector *beta   = psVectorAlloc(params->n, PS_TYPE_F64);
+    psVector *Params = psVectorAlloc(params->n, PS_TYPE_F64);
+    psVector *dy     = NULL;
 
     // the user provides the error or NULL.  we need to convert
     // to appropriate weights
-    if (yErr != NULL) {
-        for (int i = 0; i < dy->n; i++) {
-            dy->data.F32[i] = 1.0 / PS_SQR (yErr->data.F32[i]);
-        }
+    if (yWt != NULL) {
+        dy = (psVector *) yWt;
     } else {
-        for (int i = 0; i < dy->n; i++) {
-            dy->data.F32[i] = 1.0;
-        }
-    }
-
-    p_psMinLM_SetABX (alpha, beta, params, paramMask, x, y, dy, func);
-    p_psMinLM_GuessABP (Alpha, delta, Params, alpha, beta, params, paramMask, 0.0);
-
-    psFree (alpha);
-    psFree (Alpha);
-    psFree (beta);
-    psFree (Params);
-    psFree (dy);
+        dy = psVectorAlloc(y->n, PS_TYPE_F32);
+        psVectorInit(dy, 1.0);
+    }
+
+    p_psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func);
+    p_psMinLM_GuessABP(Alpha, delta, Params, alpha, beta, params, paramMask, NULL, NULL, NULL, 0.0);
+
+    psFree(alpha);
+    psFree(Alpha);
+    psFree(beta);
+    psFree(Params);
+    if (yWt == NULL) {
+        psFree(dy);
+    }
     return (true);
 }
 
-
 // measure linear model prediction
-psF64 p_psMinLM_dLinear (const psVector *Beta, const psVector *beta, psF64 lambda)
+psF64 p_psMinLM_dLinear(
+    const psVector *Beta,
+    const psVector *beta,
+    psF64 lambda)
 {
 
@@ -117,12 +202,90 @@
         dLinear += lambda*PS_SQR(B[i]) + B[i]*b[i];
     }
-    return (0.5*dLinear);
-}
+    return(0.5*dLinear);
+}
+
+// XXX EAM: this needs to respect the mask on params
+// alpha, beta, params are already allocated
+psF64 p_psMinLM_SetABX(
+    psImage  *alpha,
+    psVector *beta,
+    const psVector *params,
+    const psVector *paramMask,
+    const psArray  *x,
+    const psVector *y,
+    const psVector *dy,
+    psMinimizeLMChi2Func func)
+{
+    PS_ASSERT_IMAGE_NON_NULL(alpha, NAN);
+    PS_ASSERT_VECTOR_NON_NULL(beta, NAN);
+    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
+    PS_ASSERT_PTR_NON_NULL(x, NAN);
+    PS_ASSERT_VECTOR_NON_NULL(y, NAN);
+    PS_ASSERT_VECTOR_NON_NULL(dy, NAN);
+
+    psF64 chisq;
+    psF64 delta;
+    psF64 weight;
+    psF64 ymodel;
+    psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32);
+
+    // zero alpha and beta for summing below
+    for (int j = 0; j < params->n; j++) {
+        for (int k = 0; k < params->n; k++) {
+            alpha->data.F64[j][k] = 0;
+        }
+        beta->data.F64[j] = 0;
+    }
+    chisq = 0.0;
+
+    // calculate chisq, alpha, beta
+    for (int i = 0; i < y->n; i++) {
+        ymodel = func(deriv, params, (psVector *) x->data[i]);
+
+        delta = ymodel - y->data.F32[i];
+        chisq += PS_SQR(delta) * dy->data.F32[i];
+
+        for (int j = 0; j < params->n; j++) {
+            if ((paramMask != NULL) && (paramMask->data.U8[j]))
+                continue;
+            weight = deriv->data.F32[j] * dy->data.F32[i];
+            for (int k = 0; k <= j; k++) {
+                if ((paramMask != NULL) && (paramMask->data.U8[k]))
+                    continue;
+                alpha->data.F64[j][k] += weight * deriv->data.F32[k];
+            }
+            beta->data.F64[j] += weight * delta;
+        }
+    }
+
+    // calculate lower-left half of alpha
+    for (int j = 1; j < params->n; j++) {
+        for (int k = 0; k < j; k++) {
+            alpha->data.F64[k][j] = alpha->data.F64[j][k];
+        }
+    }
+
+    // fill in pivots if we apply a mask
+    if (paramMask != NULL) {
+        for (int j = 0; j < params->n; j++) {
+            if (paramMask->data.U8[j]) {
+                alpha->data.F64[j][j] = 1;
+                beta->data.F64[j] = 1;
+            }
+        }
+    }
+
+    psFree(deriv);
+    return(chisq);
+}
+
 
 /******************************************************************************
-psMinimizeLMChi2():  This routine will take an procedure which calculates
-an arbitrary function and it's derivative and minimize the chi-squared match
-between that function at the specified coords and the specified value at
-those coords.
+psMinimizeLMChi2():  This routine will take an procedure which calculates an
+arbitrary function and it's derivative and minimize the chi-squared match
+between that function at the specified coords and the specified value at those
+coords.
+ 
+XXX: Put the ASSERTS in.
  
 XXX EAM this is my re-implementation of MinLM
@@ -135,14 +298,17 @@
 XXX: Change the whole thing to F64, if input data is F32, convert it.
  *****************************************************************************/
-psBool psMinimizeLMChi2(psMinimization *min,
-                        psImage *covar,
-                        psVector *params,
-                        const psVector *paramMask,
-                        const psArray *x,
-                        const psVector *y,
-                        const psVector *yErr,
-                        psMinimizeLMChi2Func func)
-{
+psBool psMinimizeLMChi2(
+    psMinimization *min,
+    psImage *covar,
+    psVector *params,
+    const psVector *paramMask,
+    const psArray *x,
+    const psVector *y,
+    const psVector *yWt,
+    psMinimizeLMChi2Func func)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     PS_ASSERT_PTR_NON_NULL(min, false);
+    // XXX: If covar not NULL, do asserts...
     PS_ASSERT_VECTOR_NON_NULL(params, false);
     PS_ASSERT_VECTOR_NON_EMPTY(params, false);
@@ -162,7 +328,7 @@
     PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false);
     PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, false);
-    if (yErr != NULL) {
-        PS_ASSERT_VECTOR_TYPE(yErr, PS_TYPE_F32, false);
-        PS_ASSERT_VECTORS_SIZE_EQUAL(y, yErr, false);
+    if (yWt != NULL) {
+        PS_ASSERT_VECTOR_TYPE(yWt, PS_TYPE_F32, false);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(y, yWt, false);
     }
     PS_ASSERT_PTR_NON_NULL(func, false);
@@ -173,74 +339,84 @@
 
     // allocate internal arrays (current vs Guess)
-    psImage *alpha   = psImageAlloc (params->n, params->n, PS_TYPE_F64);
-    psImage *Alpha   = psImageAlloc (params->n, params->n, PS_TYPE_F64);
-    psVector *beta   = psVectorAlloc (params->n, PS_TYPE_F64);
-    psVector *Beta   = psVectorAlloc (params->n, PS_TYPE_F64);
-    psVector *Params = psVectorAlloc (params->n, PS_TYPE_F32);
+    psImage *alpha   = psImageAlloc(params->n, params->n, PS_TYPE_F64);
+    psImage *Alpha   = psImageAlloc(params->n, params->n, PS_TYPE_F64);
+    psVector *beta   = psVectorAlloc(params->n, PS_TYPE_F64);
+    psVector *Beta   = psVectorAlloc(params->n, PS_TYPE_F64);
+    psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32);
     psVector *dy     = NULL;
     psF64 Chisq = 0.0;
     psF64 lambda = 0.001;
-
-    // XXX EAM: why is this needed here? the value is not used, and the memory
-    // is allocated above.  However, if I drop it, I get weird answers or
-    // crashes.
-    Params = psVectorCopy (Params, params, PS_TYPE_F32);
+    // XXX: Code this properly.  Don't use mustFree00.
+    psBool mustFree00 = false;
+    psVector *beta_lim = NULL;
+    psVector *param_min = NULL;
+    psVector *param_max = NULL;
+
+    // if we are provided a covar image, we expect to find these three vectors in first three rows
+    if (covar != NULL) {
+        mustFree00 = true;
+        beta_lim  = psVectorAlloc(params->n, PS_TYPE_F32);
+        param_min = psVectorAlloc(params->n, PS_TYPE_F32);
+        param_max = psVectorAlloc(params->n, PS_TYPE_F32);
+        for (int i = 0; i < params->n; i++) {
+            beta_lim->data.F32[i] = covar->data.F64[0][i];
+            param_min->data.F32[i] = covar->data.F64[1][i];
+            param_max->data.F32[i] = covar->data.F64[2][i];
+        }
+        psImageRecycle(covar, params->n, params->n, PS_TYPE_F64);
+    }
+
+    // why is this needed here??? the initial guess on params is provided by the user
+    Params = psVectorCopy(Params, params, PS_TYPE_F32);
 
     // the user provides the error or NULL.  we need to convert
     // to appropriate weights
-    dy = psVectorAlloc (y->n, PS_TYPE_F32);
-    if (yErr != NULL) {
-        for (int i = 0; i < dy->n; i++) {
-            dy->data.F32[i] = 1.0 / PS_SQR (yErr->data.F32[i]);
-        }
+    if (yWt != NULL) {
+        dy = (psVector *) yWt;
     } else {
-        for (int i = 0; i < dy->n; i++) {
-            dy->data.F32[i] = 1.0;
-        }
+        dy = psVectorAlloc(y->n, PS_TYPE_F32);
+        psVectorInit(dy, 1.0);
     }
 
     // calculate initial alpha and beta, set chisq (min->value)
-    min->value = p_psMinLM_SetABX (alpha, beta, params, paramMask, x, y, dy, func);
-    # ifndef PS_NO_TRACE
+    min->value = p_psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func);
+    if (isnan(min->value)) {
+        min->iter = min->maxIter;
+        return(false);
+    }
     // dump some useful info if trace is defined
-    /* XXX: This code is seg faulting.
-        if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) {
-            p_psImagePrint  (psTraceGetDestination(), alpha, "alpha guess");
-            p_psVectorPrint (psTraceGetDestination(), beta, "beta guess");
-            p_psVectorPrint (psTraceGetDestination(), params, "params guess");
-        }
-        if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 4) {
-            // XXX: Where is this?
-            // p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess");
-        }
-    */
-    # endif /* PS_NO_TRACE */
-
-    min->lastDelta = min->tol + 1.0;
+    if (psTraceGetLevel(__func__) >= 6) {
+        p_psImagePrint(psTraceGetDestination(), alpha, "alpha guess (0)");
+        p_psVectorPrint(psTraceGetDestination(), beta, "beta guess (0)");
+        p_psVectorPrint(psTraceGetDestination(), params, "params guess (0)");
+    }
+    if (psTraceGetLevel (__func__) >= 6) {
+        //XXX:  p_psVectorPrintRow(psTraceGetDestination(), Params, "params guess");
+    }
+
     // iterate until the tolerance is reached, or give up
-    while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {
+    while ((min->iter < min->maxIter) && ((min->lastDelta > min->tol) || !isfinite(min->lastDelta))) {
+        psTrace(__func__, 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
+        psTrace(__func__, 5, "Last delta is %f.  Min->tol is %f.\n", min->lastDelta, min->tol);
+
         // set a new guess for Alpha, Beta, Params
-        p_psMinLM_GuessABP (Alpha, Beta, Params, alpha, beta, params, paramMask, lambda);
+        p_psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask,
+                           beta_lim, param_min, param_max, lambda);
 
         // measure linear model prediction
-        psF64 dLinear = p_psMinLM_dLinear (Beta, beta, lambda);
-
-        # ifndef PS_NO_TRACE
+        psF64 dLinear = p_psMinLM_dLinear(Beta, beta, lambda);
+
         // dump some useful info if trace is defined
-        /* XXX: This code is seg faulting.
-                if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) {
-                    p_psImagePrint  (psTraceGetDestination(), Alpha, "alpha guess");
-                    p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess");
-                    p_psVectorPrint (psTraceGetDestination(), Params, "params guess");
-                }
-                if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 4) {
-                    // XXX: Where is this?
-                    // p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess");
-                }
-        */
-        # endif /* PS_NO_TRACE */
+        if (psTraceGetLevel(__func__) >= 6) {
+            p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (1)");
+            p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (1)");
+            p_psVectorPrint(psTraceGetDestination(), Params, "params guess (1)");
+        }
+        if (psTraceGetLevel(__func__) >= 6) {
+            //XXX: p_psVectorPrintRow(psTraceGetDestination(), Params, "params guess");
+        }
 
         // calculate Chisq for new guess, update Alpha & Beta
-        Chisq = p_psMinLM_SetABX (Alpha, Beta, Params, paramMask, x, y, dy, func);
+        Chisq = p_psMinLM_SetABX(Alpha, Beta, Params, paramMask, x, y, dy, func);
 
         // XXX EAM alternate convergence criterion:
@@ -250,15 +426,13 @@
         psF64 rho = (min->value - Chisq) / dLinear;
 
-        psTrace (".psLib.dataManip.psMinimizeLMChi2", 4, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value, Chisq, min->lastDelta, rho);
-        # ifndef PS_NO_TRACE
+        psTrace(__func__, 5, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value,
+                Chisq, min->lastDelta, rho);
+
         // dump some useful info if trace is defined
-        /* XXX: This code is seg faulting.
-                if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) {
-                    p_psImagePrint  (psTraceGetDestination(), Alpha, "alpha guess");
-                    p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess");
-                    p_psVectorPrint (psTraceGetDestination(), Params, "params guess");
-                }
-        */
-        # endif /* PS_NO_TRACE */
+        if (psTraceGetLevel(__func__) >= 6) {
+            p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (2)");
+            p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (2)");
+            p_psVectorPrint(psTraceGetDestination(), Params, "params guess (2)");
+        }
 
         /* if (Chisq < min->value) {  */
@@ -266,180 +440,48 @@
             min->lastDelta = (min->value - Chisq) / (dy->n - params->n);
             min->value = Chisq;
-            alpha  = psImageCopy (alpha, Alpha, PS_TYPE_F64);
-            beta   = psVectorCopy (beta, Beta, PS_TYPE_F64);
-            params = psVectorCopy (params, Params, PS_TYPE_F32);
+            alpha  = psImageCopy(alpha, Alpha, PS_TYPE_F64);
+            beta   = psVectorCopy(beta, Beta, PS_TYPE_F64);
+            params = psVectorCopy(params, Params, PS_TYPE_F32);
             lambda *= 0.1;
         } else {
             lambda *= 10.0;
         }
-        min->iter ++;
-
-        //        printf("CONDITIONS: (%f > %f) && (%d < %d)\n", min->lastDelta, min->tol, min->iter, min->maxIter);
-    }
-    psTrace (".psLib.dataManip.psMinimizeLMChi2", 4, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter);
+        min->iter++;
+    }
+    psTrace(__func__, 5, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter);
 
     // construct & return the covariance matrix (if requested)
     if (covar != NULL) {
-        p_psMinLM_GuessABP (covar, Beta, Params, alpha, beta, params, paramMask, 0.0);
+        p_psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask,
+                           beta_lim, param_min, param_max, 0.0);
     }
 
     // free the internal temporary data
-    psFree (alpha);
-    psFree (Alpha);
-    psFree (beta);
-    psFree (Beta);
-    psFree (Params);
-    psFree (dy);
-
+    psFree(alpha);
+    psFree(Alpha);
+    psFree(beta);
+    psFree(Beta);
+    psFree(Params);
+    if (yWt == NULL) {
+        psFree(dy);
+    }
+    if (mustFree00 == true) {
+        psFree(beta_lim);
+        psFree(param_min);
+        psFree(param_max);
+    }
     if (min->iter == min->maxIter) {
-        return (false);
-    }
-    return (true);
-}
-
-// XXX EAM: this needs to respect the mask on params
-// alpha, beta, params are already allocated
-psF64 p_psMinLM_SetABX (psImage  *alpha,
-                        psVector *beta,
-                        const psVector *params,
-                        const psVector *paramMask,
-                        const psArray  *x,
-                        const psVector *y,
-                        const psVector *dy,
-                        psMinimizeLMChi2Func func)
-{
-    PS_ASSERT_IMAGE_NON_NULL(alpha, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(beta, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(params, NAN);
-    PS_ASSERT_PTR_NON_NULL(x, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(y, NAN);
-    PS_ASSERT_VECTOR_NON_NULL(dy, NAN);
-
-    psF64 chisq;
-    psF64 delta;
-    psF64 weight;
-    psF64 ymodel;
-    psVector *deriv = psVectorAlloc (params->n, PS_TYPE_F32);
-
-    // zero alpha and beta for summing below
-    for (int j = 0; j < params->n; j++) {
-        for (int k = 0; k < params->n; k++) {
-            alpha->data.F64[j][k] = 0;
-        }
-        beta->data.F64[j] = 0;
-    }
-    chisq = 0.0;
-
-    // calculate chisq, alpha, beta
-    for (int i = 0; i < y->n; i++) {
-        ymodel = func (deriv, params, (psVector *) x->data[i]);
-
-        delta = ymodel - y->data.F32[i];
-        chisq += PS_SQR (delta) * dy->data.F32[i];
-
-        for (int j = 0; j < params->n; j++) {
-            if ((paramMask != NULL) && (paramMask->data.U8[j]))
-                continue;
-            weight = deriv->data.F32[j] * dy->data.F32[i];
-            for (int k = 0; k <= j; k++) {
-                if ((paramMask != NULL) && (paramMask->data.U8[k]))
-                    continue;
-                alpha->data.F64[j][k] += weight * deriv->data.F32[k];
-            }
-            beta->data.F64[j] += weight * delta;
-        }
-    }
-
-    // calculate lower-left half of alpha
-    for (int j = 1; j < params->n; j++) {
-        for (int k = 0; k < j; k++) {
-            alpha->data.F64[k][j] = alpha->data.F64[j][k];
-        }
-    }
-
-    // fill in pivots if we apply a mask
-    if (paramMask != NULL) {
-        for (int j = 0; j < params->n; j++) {
-            if (paramMask->data.U8[j]) {
-                alpha->data.F64[j][j] = 1;
-                beta->data.F64[j] = 1;
-            }
-        }
-    }
-
-    psFree (deriv);
-    return (chisq);
-}
-
-// XXX EAM : can we use static copies of LUv, LUm, A?
-psBool p_psMinLM_GuessABP (psImage  *Alpha,
-                           psVector *Beta,
-                           psVector *Params,
-                           const psImage  *alpha,
-                           const psVector *beta,
-                           const psVector *params,
-                           const psVector *paramMask,
-                           psF64 lambda)
-{
-
-    # define USE_LU_DECOMP 1
-    # if (USE_LU_DECOMP)
-        psVector *LUv = NULL;
-    psImage  *LUm = NULL;
-    psImage  *A   = NULL;
-    psF32    det;
-
-    // LU decomposition version
-    psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using LUD version\n");
-
-    // set new guess values (creates matrix A)
-    A = psImageCopy (NULL, alpha, PS_TYPE_F64);
-    for (int j = 0; j < params->n; j++) {
-        if ((paramMask != NULL) && (paramMask->data.U8[j]))
-            continue;
-        A->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
-    }
-
-    // solve A*beta = Beta (Alpha = 1/A)
-    // these operations do not modify the input values (creates LUm, LUv)
-    LUm   = psMatrixLUD (NULL, &LUv, A);
-    Beta  = psMatrixLUSolve (Beta, LUm, beta, LUv);
-    Alpha = psMatrixInvert (Alpha, A, &det);
-
-    # else
-        // gauss-jordan version
-        psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using Gauss-J version");
-
-    // set new guess values (creates matrix A)
-    Beta = psVectorCopy (Beta, beta, PS_TYPE_F64);
-    Alpha = psImageCopy (Alpha, alpha, PS_TYPE_F64);
-    for (int j = 0; j < params->n; j++) {
-        if ((paramMask != NULL) && (paramMask->data.U8[j]))
-            continue;
-        Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
-    }
-
-    psGaussJordan (Alpha, Beta);
-    # endif
-
-    // apply Beta to get new Params values
-    for (int j = 0; j < params->n; j++) {
-        if ((paramMask != NULL) && (paramMask->data.U8[j]))
-            continue;
-        Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j];
-    }
-
-    # if (USE_LU_DECOMP)
-        psFree (A);
-    psFree (LUm);
-    psFree (LUv);
-    # endif
-
-    return true;
+        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        return(false);
+    }
+    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
+    return(true);
 }
 
 // XXX EAM : temporary gauss-jordan solver based on gene's
 // version based on the Numerical Recipes version
-bool psGaussJordan (psImage *a, psVector *b)
+bool psGaussJordan(
+    psImage *a,
+    psVector *b)
 {
 
@@ -455,7 +497,7 @@
     vector = b->data.F64;
 
-    indxc = psAlloc (Nx*sizeof(int));
-    indxr = psAlloc (Nx*sizeof(int));
-    ipiv  = psAlloc (Nx*sizeof(int));
+    indxc = psAlloc(Nx*sizeof(int));
+    indxr = psAlloc(Nx*sizeof(int));
+    ipiv  = psAlloc(Nx*sizeof(int));
     for (j = 0; j < Nx; j++) {
         ipiv[j] = 0;
@@ -476,5 +518,5 @@
                     if (ipiv[k] == 0) {
                         if (fabs (matrix[j][k]) >= big) {
-                            big  = fabs (matrix[j][k]);
+                            big  = fabs(matrix[j][k]);
                             irow = j;
                             icol = k;
@@ -492,7 +534,7 @@
         if (irow != icol) {
             for (l = 0; l < Nx; l++) {
-                PS_SWAP (matrix[irow][l], matrix[icol][l]);
-            }
-            PS_SWAP (vector[irow], vector[icol]);
+                PS_SWAP(matrix[irow][l], matrix[icol][l]);
+            }
+            PS_SWAP(vector[irow], vector[icol]);
         }
         indxr[i] = irow;
@@ -524,18 +566,18 @@
         if (indxr[l] != indxc[l]) {
             for (k = 0; k < Nx; k++) {
-                PS_SWAP (matrix[k][indxr[l]], matrix[k][indxc[l]]);
-            }
-        }
-    }
-    psFree (ipiv);
-    psFree (indxr);
-    psFree (indxc);
-    return (true);
+                PS_SWAP(matrix[k][indxr[l]], matrix[k][indxc[l]]);
+            }
+        }
+    }
+    psFree(ipiv);
+    psFree(indxr);
+    psFree(indxc);
+    return(true);
 
 fescape:
-    psFree (ipiv);
-    psFree (indxr);
-    psFree (indxc);
-    return (false);
+    psFree(ipiv);
+    psFree(indxr);
+    psFree(indxc);
+    return(false);
 }
 
@@ -565,5 +607,5 @@
 bool psMemCheckMinimization(psPtr ptr)
 {
-    return ( psMemGetDeallocator(ptr) == (psFreeFunc)minimizationFree );
+    return( psMemGetDeallocator(ptr) == (psFreeFunc)minimizationFree );
 }