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

Timestamp:

Aug 26, 2004, 10:10:25 AM (22 years ago)

Author:

gusciora

Message:

...

File:

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

trunk/psLib/src/math/psMinimize.c

@@ -9 +9 @@
  *  @author George Gusciora, MHPCC
  *
- *  @version $Revision: 1.32 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2004-08-24 19:59:05 $
+ *  @version $Revision: 1.33 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2004-08-26 20:10:25 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -939 +939 @@
 }
 
+#define STEP_SIZE 0.10
+/******************************************************************************
+    This routine takes as input an arbitrary function of a single parameter.
+    This function produces as output a bracket [a, b, c] such that f(b) is
+    less than f(a) and f(b).
+ 
+    Algorithm: XXX completely ad hoc: start with the user-supplied starting
+    parameter and call that b.  Calculate a/c as a fractional amount
+    smaller/larger than b.  Repeat this process until a local minimum is
+    found.
+ *****************************************************************************/
+psVector *p_psDetermineBracket(psVector *params,
+                               const psArray *coords,
+                               psMinimizePowellFunc func)
+{
+    float a = 0.0;
+    float b = 0.0;
+    float c = 0.0;
+    float fa = 0.0;
+    float fb = 0.0;
+    float fc = 0.0;
+    int iter = 100;
+    float aDir = 0.0;
+    float cDir = 0.0;
+    float new_aDir = 0.0;
+    float new_cDir = 0.0;
+    psVector *bracket = psVectorAlloc(3, PS_TYPE_F32);
+    psVector *dummyParams = psVectorAlloc(1, PS_TYPE_F32);
+
+    a = b = c = params->data.F32[0];
+    a = a * (1.0 - STEP_SIZE);
+    c = c * (1.0 + STEP_SIZE);
+
+    dummyParams->data.F32[0] = a;
+    fa = func(dummyParams, coords);
+    dummyParams->data.F32[0] = b;
+    fb = func(dummyParams, coords);
+    dummyParams->data.F32[0] = c;
+    fc = func(dummyParams, coords);
+    if (fa < fb) {
+        aDir = -1;
+    } else {
+        aDir = 1;
+    }
+
+    if (fc < fb) {
+        cDir = -1;
+    } else {
+        cDir = 1;
+    }
+
+    while (iter > 0) {
+        if ((b < a) && (b < c)) {
+            bracket->data.F32[0] = a;
+            bracket->data.F32[1] = b;
+            bracket->data.F32[2] = c;
+            return(bracket);
+        }
+        a = a * (1.0 - STEP_SIZE);
+        c = c * (1.0 + STEP_SIZE);
+        dummyParams->data.F32[0] = a;
+        fa = func(dummyParams, coords);
+        dummyParams->data.F32[0] = c;
+        fc = func(dummyParams, coords);
+
+        if (fa < fb) {
+            new_aDir = -1;
+        } else {
+            new_aDir = 1;
+        }
+
+        if (fc < fb) {
+            new_cDir = -1;
+        } else {
+            new_cDir = 1;
+        }
+        if ((new_aDir == 1) && (aDir == -1)) {
+            bracket->data.F32[0] = a;
+            bracket->data.F32[1] = bracket->data.F32[0] / (1.0 - STEP_SIZE);
+            bracket->data.F32[2] = bracket->data.F32[1] / (1.0 - STEP_SIZE);
+            return(bracket);
+        }
+
+        if ((new_cDir == 1) && (cDir == -1)) {
+            bracket->data.F32[2] = c;
+            bracket->data.F32[1] = bracket->data.F32[0] / (1.0 + STEP_SIZE);
+            bracket->data.F32[0] = bracket->data.F32[1] / (1.0 + STEP_SIZE);
+            return(bracket);
+        }
+        aDir = new_aDir;
+        cDir = new_cDir;
+    }
+    psFree(dummyParams);
+    psFree(bracket);
+    return(NULL);
+}
+
+
+
+/******************************************************************************
+    This routine must minimize a function of exactly one parameter.
+ *****************************************************************************/
+bool psMinimize1DFunc(psMinimization *min,
+                      psVector *params,
+                      const psArray *coords,
+                      psMinimizePowellFunc func)
+{
+    psVector *bracket;
+    psVector *dummyParams = psVectorAlloc(1, PS_TYPE_F32);
+    float a = 0.0;
+    float b = 0.0;
+    float c = 0.0;
+    float n = 0.0;
+    float fb = 0.0;
+    float fn = 0.0;
+
+    if (params->n != 1) {
+        psAbort(__func__, "params is larger than a single dimension.");
+    }
+    bracket = p_psDetermineBracket(params, coords, func);
+
+
+    min->iter = 0;
+    while (1) {
+        min->iter++;
+
+        a = bracket->data.F32[0];
+        b = bracket->data.F32[1];
+        c = bracket->data.F32[2];
+        dummyParams->data.F32[0] = b;
+        fb = func(dummyParams, coords);
+
+        // We determine which is the biggest segment in [a,b,c] then split
+        // that with the point n.
+        if ((b-a) > (c-b)) {
+            // This is the golden section formula
+            dummyParams->data.F32[0] = n = a + (0.69 * (b-a));
+            fn = func(dummyParams, coords);
+            if (fn > fb) {
+                // a = n, b = b, c = c
+                bracket->data.F32[0] = n;
+            } else {
+                // a = a, b = n, c = b
+                bracket->data.F32[1] = n;
+                bracket->data.F32[2] = b;
+            }
+        } else {
+            dummyParams->data.F32[0] = b + (0.69 * (c-b));
+            fn = func(dummyParams, coords);
+            if (fn > fb) {
+                // a = a, b = b, c = n
+                bracket->data.F32[2] = n;
+            } else {
+                // a = b, b = n, c = c
+                bracket->data.F32[0] = b;
+                bracket->data.F32[1] = n;
+            }
+        }
+    }
+
+    psFree(dummyParams);
+}
+
 bool psMinimizePowell(psMinimization *min,
                       psVector *params,