Bug report and fix for adaptive step size control in ODE solving

Bug report and fix for adaptive step size control in ODE solving

[Frank Reininghaus]

[-- Attachment #1: Type: text/plain, Size: 3309 bytes --]

Hi,

I've found some problems with adaptive step size control which can
lead to an infinite loop under special circumstances. I've attached
three test case programs where the rhs of the ODE to be solved is a
step function. All programs lead to infinite loops either in the the
loop in my code (which should be equivalent to the loop in the example 
program
in the documentation) or in gsl_odeiv_evolve_apply () (in the file
ode-initval/evolve.c) for different reasons.

The bug ocurred on an AMD Athlon 64 in 32-bit mode with the current
CVS version of GSL. 'gcc --version' yields
gcc (GCC) 4.1.3 20070929 (prerelease) (Ubuntu 4.1.2-16ubuntu2).
I've also checked the testcases on a Pentium III with
gcc (GCC) 4.1.2 20061115 (prerelease) (SUSE Linux).

In testcase1.c, the function rhs (t) yields 1 for t>=1.0 and 0
otherwise. Start and end points for the ODE solver are t=0.0 and
t=2.0, respectively, and the initial value is y=0.0. The problem is
that for very small epsabs, the step size h0 in gsl_odeiv_evolve_apply
() is reduced so much [by calls to gsl_odeiv_control_hadjust ()] that
t0 + h0 seems to be equal  to t0 due to the limited precision (h0/t0
is smaller than GSL_DBL_EPSILON, defined in gsl_machine.h). Therefore,
the ODE solving process gets stuck at the point of discontinuity
(t=1.0).

This behaviour occurs not only with gsl_odeiv_step_rk2, but also with
other solvers (the value of epsabs necessary to get into the infinite
loop might vary), except for the implicit solvers, but some of these
yield very inaccurate results.

In testcase2.c, the rhs is 1e300 for t>=0 and 0 otherwise. Integrating
from t=-1.0 to t=1.0 causes step size control to reduce h0 successively
at the point of discontinuity until it becomes zero, and this makes
the ODE solver get into the same infinite loop.

In testcase3.c, the problem is different. The only difference from
testcase2.c is the choice of a different solver (rkf45 instead of
rk2). Now the step size h0 is not reduced to zero, but only to
4.94066e-324, which seems to be the smallest possible double on my
machine. The problem here is that the function std_control_hadjust ()
multiplies the step size in line 113 in ode-initval/cstd.c with
r=0.535431. This is rounded to 4.94066e-324 by the CPU, i.e., the step
size is unchanged. Nevertheless, the function returns
GSL_ODEIV_HADJ_DEC, and this makes gsl_odeiv_evolve_apply () believe
that the step size was reduced, and therefore it jumps back to the label
'try_step' and does exactly the same thing as before with the unchanged
step size, so we end up in an infinite loop in gsl_odeiv_evolve_apply ().

I've attached a suggestion for a fix. In the current CVS version,
gsl_odeiv_evolve_apply () jumps back to the label 'try_step' if
gsl_odeiv_control_hadjust () returns GSL_ODEIV_HADJ_DEC, but in my
version, it first makes sure that h0 does not get smaller than the
maximum of GSL_DBL_MIN and GSL_DBL_EPSILON * t0 to avoid the problems
described above. Before it jumps back to 'try_step', it checks if h0
is really smaller than in the previous step (we would get into an
infinite loop otherwise).

Note that std_control_hadjust () still returns GSL_ODEIV_HADJ_DEC
although the step size is unchanged in testcase3.c, but the new
version of gsl_odeiv_evolve_apply () detects this.

Regards,
Frank

[-- Attachment #2: testcase1.c --]
[-- Type: text/x-csrc, Size: 1067 bytes --]

/* compile and link with
gcc -o testcase1 testcase1.c -lgsl -lgslcblas -lm
*/

#include <stdio.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_errno.h>

int rhs (double t, const double * y, double * dydt, void * params) {
  if (t >= 1.0)
    dydt [0] = 1;
  else
    dydt [0] = 0;

  return 0;
}

int main (int argc, char ** argv) {
  /* infinite loop for epsabs = 1e-18, but not for 1e-17 */
  double epsabs = 1e-18;
  double epsrel = 1e-6;

  const gsl_odeiv_step_type * T = gsl_odeiv_step_rk2;
  gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 1);
  gsl_odeiv_control * c = gsl_odeiv_control_y_new (epsabs, epsrel);
  gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (1);
  gsl_odeiv_system sys = {rhs, 0, 1, 0};
     
  double t = 0.0;
  double h = 1e-6;
  double y = 0.0;
     
  while (t < 2.0) {
    int status = gsl_odeiv_evolve_apply (e, c, s,
					 &sys, 
					 &t, 2,
					 &h, &y);
    if (status != GSL_SUCCESS)
      break;
  }

  printf ("%g\n", y);
     
  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);
  return 0;
}


[-- Attachment #3: testcase2.c --]
[-- Type: text/x-csrc, Size: 1074 bytes --]

/* compile and link with
gcc -o testcase2 testcase2.c -lgsl -lgslcblas -lm
*/

#include <stdio.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_errno.h>

int rhs (double t, const double * y, double * dydt, void * params) {
  if (t >= 0.0)
    dydt [0] = 1e300;
  else
    dydt [0] = 0;

  return 0;
}

int main (int argc, char ** argv) {
  /* infinite loop for epsabs = 1e-25, but not for 1e-24 */
  double epsabs = 1e-25;
  double epsrel = 1e-6;

  const gsl_odeiv_step_type * T = gsl_odeiv_step_rk2;
  gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 1);
  gsl_odeiv_control * c = gsl_odeiv_control_y_new (epsabs, epsrel);
  gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (1);
  gsl_odeiv_system sys = {rhs, 0, 1, 0};
     
  double t = -1.0;
  double h = 1e-6;
  double y = 0.0;
     
  while (t < 1.0) {
    int status = gsl_odeiv_evolve_apply (e, c, s,
					 &sys, 
					 &t, 1.0,
					 &h, &y);
    if (status != GSL_SUCCESS)
      break;
  }

  printf ("%g\n", y);
     
  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);
  return 0;
}


[-- Attachment #4: testcase3.c --]
[-- Type: text/x-csrc, Size: 1076 bytes --]

/* compile and link with
gcc -o testcase3 testcase3.c -lgsl -lgslcblas -lm
*/

#include <stdio.h>
#include <gsl/gsl_odeiv.h>
#include <gsl/gsl_errno.h>

int rhs (double t, const double * y, double * dydt, void * params) {
  if (t >= 0.0)
    dydt [0] = 1e300;
  else
    dydt [0] = 0;

  return 0;
}

int main (int argc, char ** argv) {
  /* infinite loop for epsabs = 1e-26, but not for 1e-25 */
  double epsabs = 1e-26;
  double epsrel = 1e-6;

  const gsl_odeiv_step_type * T = gsl_odeiv_step_rkf45;
  gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 1);
  gsl_odeiv_control * c = gsl_odeiv_control_y_new (epsabs, epsrel);
  gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (1);
  gsl_odeiv_system sys = {rhs, 0, 1, 0};
     
  double t = -1.0;
  double h = 1e-6;
  double y = 0.0;
     
  while (t < 1.0) {
    int status = gsl_odeiv_evolve_apply (e, c, s,
					 &sys, 
					 &t, 1.0,
					 &h, &y);
    if (status != GSL_SUCCESS)
      break;
  }

  printf ("%g\n", y);
     
  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);
  return 0;
}


[-- Attachment #5: evolve.c.patch --]
[-- Type: text/x-diff, Size: 1236 bytes --]

--- cvs/gsl/ode-initval/evolve.c	2007-12-14 23:53:49.000000000 +0100
+++ cvs-no-ode-bug/gsl/ode-initval/evolve.c	2007-12-26 22:51:34.000000000 +0100
@@ -199,15 +199,28 @@
   if (con != NULL)
     {
       /* Check error and attempt to adjust the step. */
+      const double h0_old = h0;
       const int hadjust_status 
         = gsl_odeiv_control_hadjust (con, step, y, e->yerr, e->dydt_out, &h0);
 
       if (hadjust_status == GSL_ODEIV_HADJ_DEC)
         {
-          /* Step was decreased. Undo and go back to try again. */
-          DBL_MEMCPY (y, e->y0, dydt->dimension);
-          e->failed_steps++;
-          goto try_step;
+	  /* Step size h0 was decreased. Make sure it does not get too small:
+	     It should be positive (i.e., h0 >= GSL_DBL_MIN) and not 
+	     smaller than GSL_DBL_EPSILON * t0 (because t0 + h0 has to be 
+	     distinguishable from t0) */
+
+	  h0 = GSL_MAX_DBL (h0, GSL_MAX_DBL (GSL_DBL_MIN, GSL_DBL_EPSILON * t0));
+
+          /* If the step size h0 is still smaller than in the previous step, 
+	     undo and go back to try again. */
+
+	  if (h0 < h0_old) 
+	    {
+	      DBL_MEMCPY (y, e->y0, dydt->dimension);
+	      e->failed_steps++;
+	      goto try_step;
+	    }
         }
     }
 

Re: Bug report and fix for adaptive step size control in ODE solving

[Kevin H. Hobbs]

[-- Attachment #1: Type: text/plain, Size: 827 bytes --]

On Wed, 2007-12-26 at 23:14 +0100, Frank Reininghaus wrote:
> Hi,
> 
> I've found some problems with adaptive step size control which can
> lead to an infinite loop under special circumstances. 
<snip>
> In testcase1.c, the function rhs (t) yields 1 for t>=1.0 and 0
> otherwise. 

I may have encountered this when using step functions for current
injection into neuron models. I'd like to be included in the replies. 

My workaround was to adjust the step size so that the step ends at the
time of the current injection change. I had to tell the model to change
the current through the comm pointer so that the end point of the step
before the change and the start point of the step after the change have
the same current level as the rest of the step. This workaround is
probably too complex for beginners.

[-- Attachment #2: This is a digitally signed message part --]
[-- Type: application/pgp-signature, Size: 189 bytes --]

Re: Bug report and fix for adaptive step size control in ODE solving

[Jochen Küpper]

[-- Attachment #1: Type: text/plain, Size: 2107 bytes --]

Dear Frank, All,


On 26.12.2007, at 23:14, Frank Reininghaus wrote:

> I've found some problems with adaptive step size control which can  
> lead to an infinite loop under special circumstances.


Nice write-up of this problem!

Although I have not looked through your cases in detail, I can  
conform that we have had similar infinite-loop problems, for example,  
on AMD Opteron systems running Linux and operated in 64bit mode  
(gcc-4.x with -m64).

Because I did not have time to look into it in detail, and because  
the q&d-hack below works for us, so far we use the following code  
snippet to work around it:

                 while(t < tend) {
                     if(GSL_SUCCESS != gsl_odeiv_evolve_apply(evolve,  
control, step,
                                                               
&system, &t, tend, &h, y)) {
                         throw Error("Flight::operator(): integration  
problem in GSL");
                     }
                     if(h < hmin) {
                         h = hinit;
                     }
                 }

with (for example)
   hinit = tstep * 0.1
   hmin  = tstep * 0.01

We use this with gsl_odeiv_step_rkf45, gsl_odeiv_step_rkck, or  
gsl_odeiv_step_rk8pd.


We would be quite interested in seeing this fixed in GSL, obviously.


I'll try to find some time to test your patch. However, why don't you  
test what you really describe your comment? This would look like the  
following (untested):


   /* Step size h0 was decreased. Make sure it does not get too small:
      It should be positive (i.e., h0 >= GSL_DBL_MIN) and large enough
      to make t0 + h0 distinguishable from t0, i.e., at least of size
      GSL_DBL_EPSILON*t0 */

   h0 = GSL_MAX_DBL(GSL_DBL_MIN, (t0 != t0+h0) ? h0 :  
GSL_DBL_EPSILON*t0);


This also avoids the multiplication in most cases.

Greetings,
Jochen
-- 
Einigkeit und Recht und Freiheit                http://www.Jochen- 
Kuepper.de
     Liberté, Égalité, Fraternité                GnuPG key: CC1B0B4D
         Sex, drugs and rock-n-roll



[-- Attachment #2: This is a digitally signed message part --]
[-- Type: application/pgp-signature, Size: 186 bytes --]

Re: Bug report and fix for adaptive step size control in ODE solving

[Frank Reininghaus]

Hi Jochen,

Jochen Küpper wrote:
> However, why don't you test what you really describe your comment? 
> This would look like the following (untested):
>
>   h0 = GSL_MAX_DBL(GSL_DBL_MIN, (t0 != t0+h0) ? h0 : GSL_DBL_EPSILON*t0);
Actually, I had been thinking of something like that first, but then I 
was afraid that an optimising compiler might transform 't0 != t0+h0' 
into '0 != h0', making the patch useless for the 'h0 < 
GSL_DBL_EPSILON*t0' problem.

Regards,
Frank

RE: Bug report and fix for adaptive step size control in ODE solving

[Michael Stauffer]

Hi,

Thanks for the patch. I've had a problem with an infinite, or at least
near-infinite loop, using the ODE solver, so I'm hoping this will help me.
I'm new to GSL and ODE solvers, so I appreciate someone who understands the
complexities taking a whack at it.

Cheers,
Michael

> -----Original Message-----
> From: Frank Reininghaus [mailto:frank78ac@googlemail.com] 
> Sent: Wednesday, December 26, 2007 5:14 PM
> To: gsl-discuss@sourceware.org
> Subject: Bug report and fix for adaptive step size control in 
> ODE solving
> 
> 
> Hi,
> 
> I've found some problems with adaptive step size control 
> which can lead to an infinite loop under special 
> circumstances. I've attached three test case programs where 
> the rhs of the ODE to be solved is a step function. All 
> programs lead to infinite loops either in the the loop in my 
> code (which should be equivalent to the loop in the example 
> program
> in the documentation) or in gsl_odeiv_evolve_apply () (in the file
> ode-initval/evolve.c) for different reasons.
> 
> The bug ocurred on an AMD Athlon 64 in 32-bit mode with the 
> current CVS version of GSL. 'gcc --version' yields gcc (GCC) 
> 4.1.3 20070929 (prerelease) (Ubuntu 4.1.2-16ubuntu2). I've 
> also checked the testcases on a Pentium III with gcc (GCC) 
> 4.1.2 20061115 (prerelease) (SUSE Linux).
> 
> In testcase1.c, the function rhs (t) yields 1 for t>=1.0 and 
> 0 otherwise. Start and end points for the ODE solver are 
> t=0.0 and t=2.0, respectively, and the initial value is 
> y=0.0. The problem is that for very small epsabs, the step 
> size h0 in gsl_odeiv_evolve_apply
> () is reduced so much [by calls to gsl_odeiv_control_hadjust 
> ()] that t0 + h0 seems to be equal  to t0 due to the limited 
> precision (h0/t0 is smaller than GSL_DBL_EPSILON, defined in 
> gsl_machine.h). Therefore, the ODE solving process gets stuck 
> at the point of discontinuity (t=1.0).
> 
> This behaviour occurs not only with gsl_odeiv_step_rk2, but 
> also with other solvers (the value of epsabs necessary to get 
> into the infinite loop might vary), except for the implicit 
> solvers, but some of these yield very inaccurate results.
> 
> In testcase2.c, the rhs is 1e300 for t>=0 and 0 otherwise. 
> Integrating from t=-1.0 to t=1.0 causes step size control to 
> reduce h0 successively at the point of discontinuity until it 
> becomes zero, and this makes the ODE solver get into the same 
> infinite loop.
> 
> In testcase3.c, the problem is different. The only difference 
> from testcase2.c is the choice of a different solver (rkf45 
> instead of rk2). Now the step size h0 is not reduced to zero, 
> but only to 4.94066e-324, which seems to be the smallest 
> possible double on my machine. The problem here is that the 
> function std_control_hadjust () multiplies the step size in 
> line 113 in ode-initval/cstd.c with r=0.535431. This is 
> rounded to 4.94066e-324 by the CPU, i.e., the step size is 
> unchanged. Nevertheless, the function returns 
> GSL_ODEIV_HADJ_DEC, and this makes gsl_odeiv_evolve_apply () 
> believe that the step size was reduced, and therefore it 
> jumps back to the label 'try_step' and does exactly the same 
> thing as before with the unchanged step size, so we end up in 
> an infinite loop in gsl_odeiv_evolve_apply ().
> 
> I've attached a suggestion for a fix. In the current CVS 
> version, gsl_odeiv_evolve_apply () jumps back to the label 
> 'try_step' if gsl_odeiv_control_hadjust () returns 
> GSL_ODEIV_HADJ_DEC, but in my version, it first makes sure 
> that h0 does not get smaller than the maximum of GSL_DBL_MIN 
> and GSL_DBL_EPSILON * t0 to avoid the problems described 
> above. Before it jumps back to 'try_step', it checks if h0 is 
> really smaller than in the previous step (we would get into 
> an infinite loop otherwise).
> 
> Note that std_control_hadjust () still returns 
> GSL_ODEIV_HADJ_DEC although the step size is unchanged in 
> testcase3.c, but the new version of gsl_odeiv_evolve_apply () 
> detects this.
> 
> Regards,
> Frank
> 

Unwanted compiler optimization (was: Bug report and fix for adaptive step size control in ODE solving)

[Jochen Küpper]

[-- Attachment #1: Type: text/plain, Size: 974 bytes --]

Dear Optimization specialists,

On 02.01.2008, at 19:46, Frank Reininghaus wrote:

> Jochen Küpper wrote:
>
>> However, why don't you test what you really describe your comment?  
>> This would look like the following (untested):
>>
>>  h0 = GSL_MAX_DBL(GSL_DBL_MIN, (t0 != t0+h0) ? h0 :  
>> GSL_DBL_EPSILON*t0);
>
> Actually, I had been thinking of something like that first, but then  
> I was afraid that an optimising compiler might transform 't0 !=  
> t0+h0' into '0 != h0', making the patch useless for the 'h0 <  
> GSL_DBL_EPSILON*t0' problem.


Is this an optimization that is considered valid even for floating- 
point calculations in compiler building?

Maybe we should ask our friends from the GCC development about it...

Greetings,
Jochen
-- 
Einigkeit und Recht und Freiheit                http://www.Jochen-Kuepper.de
     Liberté, Égalité, Fraternité                GnuPG key: CC1B0B4D
         Sex, drugs and rock-n-roll



[-- Attachment #2: This is a digitally signed message part --]
[-- Type: application/pgp-signature, Size: 186 bytes --]

Re: Bug report and fix for adaptive step size control in ODE solving

[Brian Gough]

At Wed, 26 Dec 2007 23:14:17 +0100,
Frank Reininghaus wrote:
> I've found some problems with adaptive step size control which can
> lead to an infinite loop under special circumstances. I've attached
> three test case programs where the rhs of the ODE to be solved is a
> step function. All programs lead to infinite loops either in the the
> loop in my code (which should be equivalent to the loop in the example 
> program
> in the documentation) or in gsl_odeiv_evolve_apply () (in the file
> ode-initval/evolve.c) for different reasons.

Thanks for the bug report. This is really useful.
I am adding the examples to the test suite right now.
Then I'll see what the best solution is.

-- 
Brian Gough

Re: Bug report and fix for adaptive step size control in ODE solving

[Brian Gough]

At Wed, 26 Dec 2007 23:14:17 +0100,
Frank Reininghaus wrote:
> I've attached a suggestion for a fix. In the current CVS version,
> gsl_odeiv_evolve_apply () jumps back to the label 'try_step' if
> gsl_odeiv_control_hadjust () returns GSL_ODEIV_HADJ_DEC, but in my
> version, it first makes sure that h0 does not get smaller than the
> maximum of GSL_DBL_MIN and GSL_DBL_EPSILON * t0 to avoid the problems
> described above. Before it jumps back to 'try_step', it checks if h0
> is really smaller than in the previous step (we would get into an
> infinite loop otherwise).
> 
> Note that std_control_hadjust () still returns GSL_ODEIV_HADJ_DEC
> although the step size is unchanged in testcase3.c, but the new
> version of gsl_odeiv_evolve_apply () detects this.

I've implemented it as follows, using GSL_COERCE_DBL to avoid any
optimisation or excess precision problems.  This is checked in along
with the test cases.

-- 
Brian Gough

 if (con != NULL)
    {
      /* Check error and attempt to adjust the step. */

      double h_old = h0;

      const int hadjust_status 
        = gsl_odeiv_control_hadjust (con, step, y, e->yerr, e->dydt_out, &h0);

      if (hadjust_status == GSL_ODEIV_HADJ_DEC)
        {
          /* Check that the reported status is correct (i.e. an actual
             decrease in h0 occured) and the suggested h0 will change
             the time by at least 1 ulp */

          double t_curr = GSL_COERCE_DBL(*t);
          double t_next = GSL_COERCE_DBL((*t) + h0);

          if (fabs(h0) < fabs(h_old) && t_next != t_curr) 
            {
              /* Step was decreased. Undo step, and try again with new h0. */
              DBL_MEMCPY (y, e->y0, dydt->dimension);
              e->failed_steps++;
              goto try_step;
            }
          else
            {
              h0 = h_old; /* keep current step size */
            }
        }
    }

  *h = h0;  /* suggest step size for next time-step */

Re: Unwanted compiler optimization (was: Bug report and fix for adaptive step size control in ODE solving)

[Brian Gough]

At Thu, 3 Jan 2008 10:00:54 +0100,
Jochen Küpper wrote:
> Is this an optimization that is considered valid even for floating- 
> point calculations in compiler building?

For general variables no, I don't think so - it only is used if there are
constants on both sides and -ffast-math/-funsafe-math-optimizations is
enabled.  e.g. x+1==1 is converted to x==0 and x+2=3 is converted to x==1.
If -ffast-math is not used any comparisons should not be modified.