Changeset 1775 for trunk/psLib/src/math/psMinimize.c

Timestamp:

Sep 9, 2004, 4:55:59 PM (22 years ago)

Author:

gusciora

Message:

...

File:

• trunk/psLib/src/math/psMinimize.c (modified) (3 diffs)

trunk/psLib/src/math/psMinimize.c

@@ -9 +9 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.35 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-09-08 23:32:03 $
+ *  @version $Revision: 1.36 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-09-10 02:52:02 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -846 +846 @@
 /******************************************************************************
 GUS
- 
-XXX: bug report for the psMatrix covar argument.
  *****************************************************************************/
 bool psMinimizeLM(psMinimization *min,
@@ -856 +854 @@
                   psMinimizeLMFunc func)
 {
-    /*
-        psVector *beta = psVectorAlloc(params->n, PS_TYPE_F32);
-        psImage *alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32);
-        psVector *deriv;
-        int i;
-        int k;
-        float tmp;
-     
-        tmp = func(deriv, params, coords);
-        for (k=0;k<params->n;k++) {
-            beta->data.F32[k] = 0.0;
+    psVector *beta = psVectorAlloc(params->n, PS_TYPE_F32);
+    psVector *perm = psVectorAlloc(params->n, PS_TYPE_F32);
+    psVector *newParams = psVectorAlloc(params->n, PS_TYPE_F32);
+    psImage *alpha = psImageAlloc(params->n, params->n, PS_TYPE_F32);
+    psImage *A = psImageAlloc(params->n, params->n, PS_TYPE_F32);
+    psImage *aOut = psImageAlloc(params->n, params->n, PS_TYPE_F32);
+    psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32);
+    psVector *newDeriv = psVectorAlloc(params->n, PS_TYPE_F32);
+    int i;
+    int j;
+    int k;
+    float newValue;
+    float oldValue;
+    float lamda = 1.0;
+
+    min->lastDelta = -HUGE;
+    min->iter = 0;
+    min->value = func(deriv, params, coords);
+
+    while ((min->lastDelta > min->tol) &&
+            (min->iter < min->maxIter)) {
+
+        for (i=0;i<params->n;i++) {
+            beta->data.F32[i] = deriv->data.F32[i];
+        }
+
+        for (j=0;j<params->n;j++) {
+            for (k=0;k<params->n;k++) {
+                alpha->data.F32[j][k] = deriv->data.F32[j] * deriv->data.F32[j];
+            }
+        }
+        for (j=0;j<params->n;j++) {
+            for (k=0;k<params->n;k++) {
+                if (j != k) {
+                    A->data.F32[j][k] = alpha->data.F32[j][k];
+                } else {
+                    A->data.F32[j][k] = (1.0 + lamda) * alpha->data.F32[j][k];
+                }
+            }
+        }
+        aOut = psMatrixLUD(aOut, perm, A);
+        newParams = psMatrixLUSolve(newParams, aOut, beta, perm);
+
+        oldValue = min->value;
+        newValue = func(newDeriv, newParams, coords);
+        if (oldValue > newValue) {
+            min->lastDelta = oldValue - newValue;
+            min->value = newValue;
+
             for (i=0;i<params->n;i++) {
-                beta->data.F32[k]+= deriv
-            }
-    */
+                params->data.F32[i] = newParams->data.F32[i];
+                deriv->data.F32[i] = newDeriv->data.F32[i];
+            }
+            min->value = func(deriv, params, coords);
+
+            lamda*= 0.1;
+        } else {
+            lamda*= 10.0;
+        }
+        min->iter++;
+    }
+    psFree(beta);
+    psFree(perm);
+    psFree(newParams);
+    psFree(alpha);
+    psFree(A);
+    psFree(aOut);
+    psFree(deriv);
+    psFree(newDeriv);
+
+    if ((min->iter < min->maxIter) ||
+            (min->lastDelta <= min->tol)) {
+        return(true);
+    }
+
     return(false);
 }