let nmsimplex handle non-finite function values properly

["Alxneit-Kamber Ivo" <ivo dot alxneit at psi dot ch>]

hi all

finally, i managed to teach 'gsl_fminimizer_nmsimplex' how to handle
non-finite function values (finite.patch). sorry, it took far too long
(since the original post of march 9th this year).

the logic behind the patch is that we can generally treat non-finite
function values (POSINF, NEGINF, NAN) as "no improvment reached" in
'gsl_fminimizer_iterate()'. only when we contract the whole simplex
with respect to its lowest point we have to abort in case we encounter
non-finite function values. the same is true in 'gsl_fminimizer_set()'.
this ensures that the simplex never contains non-finite values.

this also allows to cleanly handle any failures of the user supplied
function: the function should return a non-finite value (e.g. 'NAN')
in such a case. this should probably be documented.

also included are two minor clean ups.
leak_plug.patch: plugs memory leak if 'gsl_fminimizer_alloc()' fails.
dim_test.patch: checks for a correct number of parameters in
'gsl_fminimizer_alloc()' and for the correct sizes of the various vector
arguments passed around.

-- 
Dr. Ivo Alxneit
Laboratory for Solar Technology   phone: +41 56 310 4092
Paul Scherrer Institute             fax: +41 56 310 2688
CH-5232 Villigen                   http://solar.web.psi.ch
Switzerland                   gnupg key: 0x515E30C7

--- gsl-1.9_orig/multimin/simplex.c	2007-04-25 08:19:12.132001235 +0200
+++ gsl-1.9/multimin/simplex.c	2007-04-25 08:52:23.222445769 +0200
@@ -61,11 +61,6 @@
   size_t i, j;
   double newval, mp;
 
-  if (x1->size1 < 2)
-    {
-      GSL_ERROR ("simplex cannot have less than two corners!", GSL_EFAILED);
-    }
-
   for (j = 0; j < x1->size2; j++)
     {
       mp = 0.0;
@@ -184,6 +179,11 @@
 {
   nmsimplex_state_t *state = (nmsimplex_state_t *) vstate;
 
+  if (n == 0)
+    {
+      GSL_ERROR("invalid number of parameters specified", GSL_EINVAL);
+    }
+
   state->x1 = gsl_matrix_alloc (n + 1, n);
 
   if (state->x1 == NULL)
@@ -228,6 +228,16 @@
 
   gsl_vector *xtemp = state->ws1;
 
+  if (xtemp->size != x->size)
+    {
+      GSL_ERROR("incompatible size of x", GSL_EINVAL);
+    }
+
+  if (xtemp->size != step_size->size)
+    {
+      GSL_ERROR("incompatible size of step_size", GSL_EINVAL);
+    }
+
   /* first point is the original x0 */
 
   val = GSL_MULTIMIN_FN_EVAL (f, x);
@@ -294,6 +304,12 @@
   int status;
   double val, val2;
 
+
+  if (xc->size != x->size)
+    {
+      GSL_ERROR("incompatible size of x", GSL_EINVAL);
+    }
+
   /* get index of highest, second highest and lowest point */
 
   dhi = ds_hi = dlo = gsl_vector_get (y1, 0);

--- gsl-1.9_orig/multimin/simplex.c	2007-04-25 08:19:12.132001235 +0200
+++ gsl-1.9/multimin/simplex.c	2007-04-25 08:22:16.306020239 +0200
@@ -103,6 +103,8 @@
   size_t i, j;
   double newval;
 
+  int status = GSL_SUCCESS;
+
   for (i = 0; i < x1->size1; i++)
     {
       if (i != best)
@@ -119,10 +121,19 @@
           gsl_matrix_get_row (xc, x1, i);
           newval = GSL_MULTIMIN_FN_EVAL (f, xc);
           gsl_vector_set (y1, i, newval);
+
+	  /* notify caller that we found at least one bad function value.
+	     we finish the contraction (and do not abort) to allow the user
+	     to handle the situation */
+
+          if(!gsl_finite(newval))
+	    {
+	      status = GSL_EBADFUNC;
+	    }
         }
     }
 
-  return GSL_SUCCESS;
+  return status;
 }
 
 static int
@@ -231,6 +242,12 @@
   /* first point is the original x0 */
 
   val = GSL_MULTIMIN_FN_EVAL (f, x);
+  
+  if (!gsl_finite(val))
+    {
+      GSL_ERROR("non-finite function value encountered", GSL_EBADFUNC);
+    }
+
   gsl_matrix_set_row (state->x1, 0, x);
   gsl_vector_set (state->y1, 0, val);
 
@@ -248,6 +265,12 @@
       val = gsl_vector_get (xtemp, i) + gsl_vector_get (step_size, i);
       gsl_vector_set (xtemp, i, val);
       val = GSL_MULTIMIN_FN_EVAL (f, xtemp);
+  
+      if (!gsl_finite(val))
+        {
+          GSL_ERROR("non-finite function value encountered", GSL_EBADFUNC);
+        }
+
       gsl_matrix_set_row (state->x1, i + 1, xtemp);
       gsl_vector_set (state->y1, i + 1, val);
     }
@@ -324,14 +347,14 @@
 
   val = nmsimplex_move_corner (-1.0, state, hi, xc, f);
 
-  if (val < gsl_vector_get (y1, lo))
+  if (gsl_finite(val) && val < gsl_vector_get (y1, lo))
     {
 
       /* reflected point becomes lowest point, try expansion */
 
       val2 = nmsimplex_move_corner (-2.0, state, hi, xc2, f);
 
-      if (val2 < gsl_vector_get (y1, lo))
+      if (gsl_finite(val2) && val2 < gsl_vector_get (y1, lo))
         {
           gsl_matrix_set_row (x1, hi, xc2);
           gsl_vector_set (y1, hi, val2);
@@ -343,11 +366,13 @@
         }
     }
 
-  /* reflection does not improve things enough */
+  /* reflection does not improve things enough
+     or
+     we got a non-finite (illegal) function value */
 
-  else if (val > gsl_vector_get (y1, s_hi))
+  else if (!gsl_finite(val) || val > gsl_vector_get (y1, s_hi))
     {
-      if (val <= gsl_vector_get (y1, hi))
+      if (gsl_finite(val) && val <= gsl_vector_get (y1, hi))
         {
 
           /* if trial point is better than highest point, replace 
@@ -361,7 +386,7 @@
 
       val2 = nmsimplex_move_corner (0.5, state, hi, xc2, f);
 
-      if (val2 <= gsl_vector_get (y1, hi))
+      if (gsl_finite(val2) && val2 <= gsl_vector_get (y1, hi))
         {
           gsl_matrix_set_row (state->x1, hi, xc2);
           gsl_vector_set (y1, hi, val2);
@@ -373,7 +398,7 @@
           /* contract the whole simplex in respect to the best point */
 
           status = nmsimplex_contract_by_best (state, lo, xc, f);
-          if (status != 0)
+          if (status != GSL_SUCCESS)
             {
               GSL_ERROR ("nmsimplex_contract_by_best failed", GSL_EFAILED);
             }

--- gsl-1.9_orig/multimin/simplex.c	2005-06-26 15:25:35.123496646 +0200
+++ gsl-1.9/multimin/simplex.c	2007-04-20 15:59:51.609295741 +0200
@@ -195,6 +195,7 @@
 
   if (state->y1 == NULL)
     {
+      gsl_matrix_free(state->x1);
       GSL_ERROR ("failed to allocate space for y", GSL_ENOMEM);
     }
 
@@ -202,6 +203,8 @@
 
   if (state->ws1 == NULL)
     {
+      gsl_matrix_free(state->x1);
+      gsl_vector_free(state->y1);
       GSL_ERROR ("failed to allocate space for ws1", GSL_ENOMEM);
     }
 
@@ -209,6 +212,9 @@
 
   if (state->ws2 == NULL)
     {
+      gsl_matrix_free(state->x1);
+      gsl_vector_free(state->y1);
+      gsl_vector_free(state->ws1);
       GSL_ERROR ("failed to allocate space for ws2", GSL_ENOMEM);
     }

Attachment: signature.asc
Description: This is a digitally signed message part