Changeset 6322 for trunk/psLib/src/math/psMinimizeLMM.c

Timestamp:

Feb 3, 2006, 12:05:22 PM (20 years ago)

Author:

gusciora

Message:

Misc code cleaning

File:

• trunk/psLib/src/math/psMinimizeLMM.c (modified) (15 diffs)

trunk/psLib/src/math/psMinimizeLMM.c

@@ -10 +10 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.3 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-01-27 20:08:58 $
+ *  @version $Revision: 1.4 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-03 22:05:22 $
  *
  *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -48 +48 @@
 
 // XXX EAM : can we use static copies of LUv, LUm, A?
-// XXX: Add trace messages, check return codes.
 psBool p_psMinLM_GuessABP(
     psImage  *Alpha,
@@ -75 +74 @@
     A = psImageCopy(NULL, alpha, PS_TYPE_F64);
     for (int j = 0; j < params->n; j++) {
-        if ((paramMask != NULL) && (paramMask->data.U8[j]))
+        if ((paramMask != NULL) && (paramMask->data.U8[j])) {
             continue;
+        }
         A->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
     }
@@ -83 +83 @@
     // 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
@@ -94 +107 @@
     Alpha = psImageCopy(Alpha, alpha, PS_TYPE_F64);
     for (int j = 0; j < params->n; j++) {
-        if ((paramMask != NULL) && (paramMask->data.U8[j]))
+        if ((paramMask != NULL) && (paramMask->data.U8[j])) {
             continue;
+        }
         Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
     }
 
-    // XXX: Check error codes!
-    psGaussJordan(Alpha, Beta);
-    psFree(A);
-    psFree(LUm);
-    psFree(LUv);
+    if (false == psGaussJordan(Alpha, Beta)) {
+        psError(PS_ERR_UNKNOWN, false, "psMatrixInvert() returned NULL\n");
+        return(false);
+    }
     # endif
 
     // apply Beta to get new Params values
     for (int j = 0; j < params->n; j++) {
-        if ((paramMask != NULL) && (paramMask->data.U8[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
@@ -126 +140 @@
         }
     }
-    # if (USE_LU_DECOMP)
-        psFree(A);
-    psFree(LUm);
-    psFree(LUv);
-    # endif
 
     return(true);
@@ -136 +145 @@
 
 
-// XXX: Add trace messages, check return codes.
 bool psMinimizeGaussNewtonDelta(
     psVector *delta,
@@ -146 +154 @@
     psMinimizeLMChi2Func func)
 {
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
     // allocate internal arrays (current vs Guess)
     psImage  *alpha  = psImageAlloc (params->n, params->n, PS_TYPE_F64);
@@ -152 +161 @@
     psVector *Params = psVectorAlloc(params->n, PS_TYPE_F64);
     psVector *dy     = NULL;
+    psBool rc = true;
 
     // the user provides the error or NULL.  we need to convert
@@ -162 +172 @@
     }
 
-    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);
+    psF64 rcF64 = p_psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func);
+    if (isnan(rcF64)) {
+        psError(PS_ERR_UNKNOWN, false, "p_psMinLM_SetABX() retruned a NAN.\n");
+        rc = false;
+    }
+    psTrace(__func__, 5, "p_psMinLM_SetABX() was succesful\n", __func__);
+
+    psBool rcBool = p_psMinLM_GuessABP(Alpha, delta, Params, alpha, beta, params, paramMask, NULL, NULL, NULL, 0.0);
+    if (rcBool == false) {
+        psError(PS_ERR_UNKNOWN, false, "p_psMinLM_GuessABP() retruned FALSE.\n");
+        rc = false;
+    }
+    psTrace(__func__, 5, "p_psMinLM_GuessABP() was succesful\n", __func__);
 
     psFree(alpha);
@@ -172 +193 @@
         psFree(dy);
     }
-    return (true);
+    psTrace(__func__, 3, "---- %s() end ----\n", __func__);
+    return(rc);
 }
 
@@ -235 +257 @@
 
         for (psS32 j = 0; j < params->n; j++) {
-            if ((paramMask != NULL) && (paramMask->data.U8[j]))
+            if ((paramMask != NULL) && (paramMask->data.U8[j])) {
                 continue;
+            }
             weight = deriv->data.F32[j] * dy->data.F32[i];
             for (psS32 k = 0; k <= j; k++) {
-                if ((paramMask != NULL) && (paramMask->data.U8[k]))
+                if ((paramMask != NULL) && (paramMask->data.U8[k])) {
                     continue;
+                }
                 alpha->data.F64[j][k] += weight * deriv->data.F32[k];
             }
@@ -394 +418 @@
     }
     if (psTraceGetLevel (__func__) >= 6) {
-        //XXX:  p_psVectorPrintRow(psTraceGetDestination(), Params, "params guess");
+        psTrace(__func__, 6, "The current Param vector: \n");
+        for (psS32 i = 0 ; i < Params->n ; i++) {
+            psTrace(__func__, 6, "Params[%d] is %f\n", Params->data.F32[i]);
+        }
     }
 
@@ -416 +443 @@
         }
         if (psTraceGetLevel(__func__) >= 6) {
-            //XXX: p_psVectorPrintRow(psTraceGetDestination(), Params, "params guess");
+            if (psTraceGetLevel (__func__) >= 6) {
+                psTrace(__func__, 6, "The current Param vector: \n");
+                for (psS32 i = 0 ; i < Params->n ; i++) {
+                    psTrace(__func__, 6, "Params[%d] is %f\n", Params->data.F32[i]);
+                }
+            }
         }
 
@@ -628 +660 @@
 bool psMemCheckConstrain(psPtr tmp)
 {
-    return( psMemGetDeallocator(tmp) == (psFreeFunc)constrainFree );
-}
+    return(psMemGetDeallocator(tmp) == (psFreeFunc) constrainFree);
+}