Index: trunk/psLib/src/math/psMinimizeLMM.c
===================================================================
--- trunk/psLib/src/math/psMinimizeLMM.c	(revision 14320)
+++ trunk/psLib/src/math/psMinimizeLMM.c	(revision 15046)
@@ -10,6 +10,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-07-20 00:22:24 $
+ *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-28 00:35:17 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -54,4 +54,5 @@
 /*****************************************************************************/
 
+// Alpha & Beta only represent unmasked values
 bool psMinLM_GuessABP(
     psImage  *Alpha,
@@ -76,87 +77,65 @@
     }
 
-    // XXX the LU Decomposition code requires F64.  this is now incompatible
-    // with the rest of this code. drop this segment when the changes are tested.
-    # define USE_LU_DECOMP 0
-    # if (USE_LU_DECOMP)
-        psVector *LUv = NULL;
-    psImage  *LUm = NULL;
-    psImage  *A   = NULL;
-    psF32    det;
-    // LU decomposition version
-    psTrace("psLib.math", 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);
-    if (LUm == NULL) {
-        psError(PS_ERR_UNKNOWN, false, "psMatrixLUD() returned NULL\n");
-    }
-    Beta  = psMatrixLUSolve(Beta, LUm, beta, LUv);
-    if (Beta == NULL) {
-        psError(PS_ERR_UNKNOWN, false, "psMatrixLUSolve() returned NULL\n");
-    }
-    Alpha = psMatrixInvert(Alpha, A, &det);
-    psFree(LUm);
-    psFree(LUv);
-    psFree(A);
-
-    if (Alpha == NULL) {
-        psError(PS_ERR_UNKNOWN, false, "psMatrixInvert() returned NULL\n");
-        return(false);
-    }
-    # else
-        // gauss-jordan version
-        psTrace("psLib.math", 5, "using Gauss-J version");
-
-    // set new guess values (creates matrix A)
+    assert (alpha->numCols == beta->n);
+    assert (alpha->numCols == alpha->numRows);
+
+    // set new guess values, applying (1+lambda) scaling to pivots
     Beta = psVectorCopy(Beta, beta, PS_TYPE_F32);
     Alpha = psImageCopy(Alpha, alpha, PS_TYPE_F32);
-    for (int j = 0; j < params->n; j++) {
-        if ((paramMask != NULL) && (paramMask->data.U8[j])) {
-            continue;
-        }
+    for (int j = 0; j < Alpha->numCols; j++) {
         Alpha->data.F32[j][j] = alpha->data.F32[j][j] * (1.0 + lambda);
     }
 
+    // error and clear above if kept?
     if (false == psMatrixGJSolveF32(Alpha, Beta)) {
-        // psError(PS_ERR_UNKNOWN, false, "singular matrix in Guess ABP\n");
         psTrace ("psLib.math", 4, "singular matrix in Guess ABP\n");
         return(false);
     }
-    # endif
 
     // measure linear model prediction
     // (we must do this before truncating Beta below)
     if (dLinear) {
-	*dLinear = psMinLM_dLinear(Beta, beta, paramMask, lambda);
+	*dLinear = psMinLM_dLinear(Beta, beta, lambda);
+    }
+
+    // full-length Beta for checkLimits functions 
+    psVector *tmpBeta = psVectorAlloc(params->n, PS_TYPE_F32);
+    psVectorInit (tmpBeta, 0.0);
+
+    // set tmpBeta values which are not masked
+    for (int j = 0, n = 0; j < params->n; j++) {
+	if (paramMask && (paramMask->data.U8[j])) continue;
+	tmpBeta->data.F32[j] = Beta->data.F32[n];
+	n++;
     }
 
     // apply Beta to get new Params values
     for (int j = 0; j < params->n; j++) {
-        if ((paramMask != NULL) && (paramMask->data.U8[j])) {
+        if (paramMask && (paramMask->data.U8[j])) {
             Params->data.F32[j] = params->data.F32[j];
             continue;
         }
         // apply beta limits
-        if (checkLimits)
-            checkLimits (PS_MINIMIZE_BETA_LIMIT, j, Params->data.F32, Beta->data.F32);
-        Params->data.F32[j] = params->data.F32[j] - Beta->data.F32[j];
+        if (checkLimits) {
+            checkLimits (PS_MINIMIZE_BETA_LIMIT, j, Params->data.F32, tmpBeta->data.F32);
+	}
+
+        Params->data.F32[j] = params->data.F32[j] - tmpBeta->data.F32[j];
 
         // compare new params to param limits
-        if (checkLimits)
-            checkLimits (PS_MINIMIZE_PARAM_MIN,  j, Params->data.F32, Beta->data.F32);
-        if (checkLimits)
-            checkLimits (PS_MINIMIZE_PARAM_MAX,  j, Params->data.F32, Beta->data.F32);
-    }
+        if (checkLimits) {
+            checkLimits (PS_MINIMIZE_PARAM_MIN,  j, Params->data.F32, tmpBeta->data.F32);
+            checkLimits (PS_MINIMIZE_PARAM_MAX,  j, Params->data.F32, tmpBeta->data.F32);
+	}
+    }
+
+    // apply tmpBeta after limits have been checked
+    for (int j = 0, n = 0; j < params->n; j++) {
+	if (paramMask && (paramMask->data.U8[j])) continue;
+	Beta->data.F32[n] = tmpBeta->data.F32[j];
+	n++;
+    }
+
+    psFree (tmpBeta);
     return(true);
 }
@@ -172,11 +151,22 @@
 {
     psTrace("psLib.math", 3, "---- begin ----\n");
+
     // allocate internal arrays (current vs Guess)
-    psImage  *alpha  = psImageAlloc (params->n, params->n, PS_TYPE_F32);
-    psImage  *Alpha  = psImageAlloc (params->n, params->n, PS_TYPE_F32);
-    psVector *beta   = psVectorAlloc(params->n, PS_TYPE_F32);
+    psImage *Alpha = NULL;
+    psVector *Beta = NULL;
+
+    psVectorInit (delta, 0.0);
+
+    // Alpha & Beta only contain elements to represent the unmasked parameters
+    // if none are available, return false
+    if (!psMinLM_AllocAB (&Alpha, &Beta, params, paramMask)) {
+	return false;
+    }
+
+    psImage *alpha   = psImageAlloc(Alpha->numCols, Alpha->numRows, PS_TYPE_F32);
     psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32);
+
     psVector *dy     = NULL;
-    bool rc = true;
+    bool retValue = true;
 
     // the user provides the error or NULL.  we need to convert
@@ -190,21 +180,22 @@
 
     // XXX should we give up if chisq is nan?
-    psF32 rcF32 = psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func);
-    if (isnan(rcF32)) {
-        psTrace ("psLib.math", 5, "psMinLM_SetABX() returned a NAN.\n");
-        rc = false;
-    }
+    psF32 chisq = psMinLM_SetABX(alpha, Beta, params, paramMask, x, y, dy, func);
+    if (isnan(chisq)) {
+        psTrace ("psLib.math", 5, "psMinLM_SetABX() returned a NAN chisq.\n");
+        retValue = false;
+    }
+
     psTrace("psLib.math", 5, "psMinLM_SetABX() was succesful\n");
     // dump some useful info if trace is defined
     if (psTraceGetLevel("psLib.math") >= 6) {
         p_psImagePrint(psTraceGetDestination(), alpha, "alpha guess (0)");
-        p_psVectorPrint(psTraceGetDestination(), beta, "beta guess (0)");
+        p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (0)");
         p_psVectorPrint(psTraceGetDestination(), params, "params guess (0)");
     }
 
-    bool rcBool = psMinLM_GuessABP(Alpha, delta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL);
-    if (rcBool == false) {
+    bool status = psMinLM_GuessABP(Alpha, delta, Params, alpha, Beta, params, paramMask, NULL, 0.0, NULL);
+    if (!status) {
         psTrace ("psLib.math", 5, "psMinLM_GuessABP() returned FALSE.\n");
-        rc = false;
+        retValue = false;
     }
     psTrace("psLib.math", 5, "psMinLM_GuessABP() was succesful\n");
@@ -217,5 +208,5 @@
     psFree(alpha);
     psFree(Alpha);
-    psFree(beta);
+    psFree(Beta);
     psFree(Params);
     if (yWt == NULL) {
@@ -223,5 +214,5 @@
     }
     psTrace("psLib.math", 3, "---- end ----\n");
-    return(rc);
+    return(retValue);
 }
 
@@ -230,5 +221,4 @@
     const psVector *Beta,
     const psVector *beta,
-    const psVector *paramMask,
     psF32 lambda)
 {
@@ -240,10 +230,10 @@
 
     float dh = 0.0, sh = 0.0, Sh = 0.0;
+
+    // beta only counts unmasked parameters
     for (int i = 0; i < beta->n; i++) {
-        if ((paramMask != NULL) && (paramMask->data.U8[i])) continue;
 	dh = lambda*B[i] + b[i];
 	sh = 0.5*B[i]*dh;
 	Sh += sh;
-	// fprintf (stderr, "%f, %f  * %f %f %f * :  %f : %f + %f\n", B[i], b[i], dh, sh, Sh, lambda*PS_SQR(B[i]) + B[i]*b[i], lambda*PS_SQR(B[i]), B[i]*b[i]);
         dLinear += lambda*PS_SQR(B[i]) + B[i]*b[i];
     }
@@ -262,5 +252,4 @@
     psMinimizeLMChi2Func func)
 {
-    // XXX: Check vector sizes.
     PS_ASSERT_IMAGE_NON_NULL(alpha, NAN);
     PS_ASSERT_VECTOR_NON_NULL(beta, NAN);
@@ -281,14 +270,11 @@
     psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32);
 
-    // zero alpha and beta for summing below
-    for (psS32 j = 0; j < params->n; j++) {
-        for (psS32 k = 0; k < params->n; k++) {
-            alpha->data.F32[j][k] = 0;
-        }
-        beta->data.F32[j] = 0;
-    }
+    // zero alpha, beta, and chisq for summing below
+    psImageInit (alpha, 0.0);
+    psVectorInit (beta, 0.0);
     chisq = 0.0;
 
-    // calculate chisq, alpha, beta
+    // calculate chisq, alpha, beta. alpha & beta only represent unmasked parameters; skip
+    // masked ones
     for (psS32 i = 0; i < y->n; i++) {
         ymodel = func(deriv, params, (psVector *) x->data[i]);
@@ -296,40 +282,28 @@
         delta = ymodel - y->data.F32[i];
         chisq += PS_SQR(delta) * dy->data.F32[i];
-        if (isnan(dy->data.F32[i]))
-            psAbort("nan in weights");
-        if (isnan(delta))
-            psAbort("nan in delta");
-        if (isnan(chisq))
-            psAbort("nan in chisq");
-
-        for (psS32 j = 0; j < params->n; j++) {
-            if ((paramMask != NULL) && (paramMask->data.U8[j])) {
-                continue;
+        assert (!isnan(dy->data.F32[i]));
+        assert (!isnan(delta));
+        assert (!isnan(chisq));
+
+	// we track alpha,beta and params,deriv separately
+        for (int j = 0, J = 0; j < params->n; j++) {
+            if (paramMask && (paramMask->data.U8[j])) continue;
+
+            weight = deriv->data.F32[j] * dy->data.F32[i];
+
+            for (int k = 0, K = 0; k <= j; k++) {
+                if (paramMask && (paramMask->data.U8[k])) continue;
+                alpha->data.F32[J][K] += weight * deriv->data.F32[k];
+		K++;
             }
-            weight = deriv->data.F32[j] * dy->data.F32[i];
-            for (psS32 k = 0; k <= j; k++) {
-                if ((paramMask != NULL) && (paramMask->data.U8[k])) {
-                    continue;
-                }
-                alpha->data.F32[j][k] += weight * deriv->data.F32[k];
-            }
-            beta->data.F32[j] += weight * delta;
+            beta->data.F32[J] += weight * delta;
+	    J++;
         }
     }
 
     // calculate lower-left half of alpha
-    for (psS32 j = 1; j < params->n; j++) {
-        for (psS32 k = 0; k < j; k++) {
+    for (int j = 1; j < alpha->numCols; j++) {
+        for (int k = 0; k < j; k++) {
             alpha->data.F32[k][j] = alpha->data.F32[j][k];
-        }
-    }
-
-    // fill in pivots if we apply a mask
-    if (paramMask != NULL) {
-        for (psS32 j = 0; j < params->n; j++) {
-            if (paramMask->data.U8[j]) {
-                alpha->data.F32[j][j] = 1;
-                beta->data.F32[j] = 1;
-            }
         }
     }
@@ -398,9 +372,16 @@
 
     // allocate internal arrays (current vs Guess)
-    psImage *alpha   = psImageAlloc(params->n, params->n, PS_TYPE_F32);
-    psImage *Alpha   = psImageAlloc(params->n, params->n, PS_TYPE_F32);
-    psVector *beta   = psVectorAlloc(params->n, PS_TYPE_F32);
-    psVector *Beta   = psVectorAlloc(params->n, PS_TYPE_F32);
+    psImage *Alpha = NULL;
+    psVector *Beta = NULL;
+
+    // Alpha & Beta only contain elements to represent the unmasked parameters
+    if (!psMinLM_AllocAB (&Alpha, &Beta, params, paramMask)) {
+	psAbort ("programming error: no unmasked parameters to be fit\n");
+    }
+    
+    psImage *alpha   = psImageAlloc(Alpha->numCols, Alpha->numRows, PS_TYPE_F32);
+    psVector *beta   = psVectorAlloc(Beta->n, PS_TYPE_F32);
     psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32);
+
     psVector *dy     = NULL;
     psF32 Chisq = 0.0;
@@ -494,7 +475,21 @@
     // construct & return the covariance matrix (if requested)
     if (covar != NULL) {
-        if (!psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) {
+        if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) {
             psTrace ("psLib.math", 5, "failure to calculate covariance matrix\n");
         }
+	// set covar values which are not masked
+	psImageInit (covar, 0.0);
+	for (int j = 0, J = 0; j < params->n; j++) {
+	    if (paramMask && (paramMask->data.U8[j])) {
+		covar->data.F32[j][j] = 1.0;
+		continue;
+	    }
+	    for (int k = 0, K = 0; k < params->n; k++) {
+		if (paramMask && (paramMask->data.U8[k])) continue;
+		covar->data.F32[j][k] = Alpha->data.F32[J][K];
+		K++;
+	    }
+	    J++;
+	}
     }
 
@@ -516,4 +511,30 @@
 }
 
+bool psMinLM_AllocAB (psImage **Alpha, psVector **Beta, const psVector *params, const psVector *paramMask) {
+
+    assert (Alpha);
+    assert (Beta);
+    assert (params);
+
+    int nParams = params->n;
+
+    // count unmasked parameters
+    if (paramMask) {
+	nParams = 0;
+	for (int i = 0; i < paramMask->n; i++) {
+	    if (paramMask->data.U8[i]) continue;
+	    nParams ++;
+	}
+    }
+
+    if (nParams == 0) { 
+	return false;
+    }
+
+    *Alpha = psImageAlloc(nParams, nParams, PS_TYPE_F32);
+    *Beta  = psVectorAlloc(nParams, PS_TYPE_F32);
+    return true;
+}
+
 static void minimizationFree(psMinimization *min)
 {