[Bug-gsl] gsl multimin/simplex.c

[Steve M. Robbins]

Hello,

It looks to me like there's a bug in simplex.c, function
nmsimplex_iterate().  I'm using the CVS code as of today, 2003-04-15.

The initialization code has a loop intended to locate the indices
of the highest, second highest, and lowest points of the simplex.
It will fail if the second highest point is *after* the highest
point in the sequence examined.  

Here is the original code.

  /* get index of highest, second highest and lowest point */
 
  dhi = gsl_vector_get (y1, 0);
  dlo = dhi;
 
  for (i = 1; i < n; i++)
    {
      val = (gsl_vector_get (y1, i));
      if (val < dlo)
        {
          dlo = val;
          lo = i;
        }
      else if (val > dhi)
        {
          s_hi = hi;
          dhi = val;
          hi = i;
        }
    }

Here's a fix.

Index: simplex.c
===================================================================
RCS file: /cvs/gsl/gsl/multimin/simplex.c,v
retrieving revision 1.3
diff -u -b -B -r1.3 simplex.c
--- simplex.c   5 Oct 2002 21:35:27 -0000       1.3
+++ simplex.c   15 Apr 2003 15:29:58 -0000
@@ -289,14 +289,13 @@
   size_t n = y1->size;
   size_t i;
   size_t hi = 0, s_hi = 0, lo = 0;
-  double dhi, dlo;
+  double dhi, ds_hi, dlo;
   int status;
   double val, val2;
 
   /* get index of highest, second highest and lowest point */
 
-  dhi = gsl_vector_get (y1, 0);
-  dlo = dhi;
+  dhi = ds_hi = dlo = gsl_vector_get (y1, 0);
 
   for (i = 1; i < n; i++)
     {
@@ -308,9 +307,15 @@
        }
       else if (val > dhi)
        {
+         ds_hi = dhi;
          s_hi = hi;
          dhi = val;
          hi = i;
+       }
+      else if (val > ds_hi)
+       {
+         ds_hi = val;
+         s_hi = i;
        }
     }
 

Cheers,
-Steve