<p dir="auto">I am having an issue with what seems to be a very infrequent infinite loop in <code class="notranslate">RungeKutta._step_impl</code> in <code class="notranslate">scipy/integrate/_ivp/rk.py</code>. Running a bunch of experiments in parallel on large number of cores, I maybe see it once a week, and not in any easily reproducible fashion, as far as I can tell.</p>
<p dir="auto">I caught it after, as far as I can tell, <code class="notranslate">scipy.integrate</code> was running for a couple of days without completing. I am guessing <code class="notranslate">RungeKutta._step_impl</code> is the one not completing, but perhaps it's higher in the call stack.</p>
<p dir="auto">I examined the running process using gdb with python extensions. At the point where <code class="notranslate">rk_step</code> is called inside the <code class="notranslate">while</code> loop, the following are some of the relevant values (more details are in the gdb session log below):</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="self.n == 6
self.t == t == 0
self.y == y == numpy.ndarray([2.19183852e-33, 0, 1.40129846e-45, 0, 1.39532281e-34, 0])
self.h_abs == -NaN
h_abs == NaN
self.max_step == max_step == Inf
min_step == 4.9406564584124654e-323
self.direction == 1
order == 4
h == t_new == -NaN
self.t_bound = 0.02"><pre class="notranslate"><code class="notranslate">self.n == 6
self.t == t == 0
self.y == y == numpy.ndarray([2.19183852e-33, 0, 1.40129846e-45, 0, 1.39532281e-34, 0])
self.h_abs == -NaN
h_abs == NaN
self.max_step == max_step == Inf
min_step == 4.9406564584124654e-323
self.direction == 1
order == 4
h == t_new == -NaN
self.t_bound = 0.02
</code></pre></div>
<p dir="auto">I am running:</p>
<ul dir="auto">
<li>Ubuntu 16.04</li>
<li>python 3.5.2</li>
<li>numpy 1.15.4</li>
<li>scipy 1.3.3</li>
</ul>
<p dir="auto">More complete <code class="notranslate">gdb</code> output:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="(gdb) py-bt
Traceback (most recent call first):
[...]
  File &quot;[...]/envs/model/model.py&quot;, line 215, in _tire_dynamics_rear
    * (np.tan(alpha)/(1+K))**2 )
  File &quot;[...]/envs/model/model.py&quot;, line 83, in _dynamics
    F_xr, F_yr = self._tire_dynamics_rear(v_x, cmd_vx, alpha_r)
  File &quot;[...]/envs/model/model.py&quot;, line 55, in &lt;lambda&gt;
    fun=(lambda t, X: self._dynamics(X, t, U)),
  File &quot;/usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/base.py&quot;, line 21, in fun_wrapped
    return np.asarray(fun(t, y), dtype=dtype)
  File &quot;/usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/base.py&quot;, line 139, in fun
    return self.fun_single(t, y)
  File &quot;/usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/rk.py&quot;, line 70, in rk_step
    K[s + 1] = fun(t + c * h, y + dy)
  File &quot;/usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/rk.py&quot;, line 143, in _step_impl
    self.B, self.C, self.E, self.K)
  File &quot;/usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/base.py&quot;, line 182, in step
    success, message = self._step_impl()
  File &quot;/usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/ivp.py&quot;, line 502, in solve_ivp
    message = solver.step()
  File &quot;[...]/envs/model/model.py&quot;, line 56, in state_transition
    t_span=t, y0=X, atol=1e-5)
  [...]
(gdb) py-up
[...]
(gdb) py-up
#24 Frame 0x61ffa18, for file /usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/rk.py, line 143, in _step_impl (self=&lt;RK45(t=&lt;float at remote 0x7f80dcc71510&gt;, direction=&lt;numpy.float64 at remote 0x7f80dccc5ee8&gt;, K=&lt;numpy.ndarray at remote 0x7f80dcc53d00&gt;, vectorized=False, _fun=&lt;function at remote 0x7f80dcc51f28&gt;, y_old=None, status='running', f=&lt;numpy.ndarray at remote 0x7f80dcc53a80&gt;, y=&lt;numpy.ndarray at remote 0x7f80dcc53350&gt;, fun_vectorized=&lt;function at remote 0x7f80dcc51950&gt;, t_old=None, fun=&lt;function at remote 0x7f80dcc51d90&gt;, t_bound=&lt;float at remote 0x7f80dcc71090&gt;, njev=0, fun_single=&lt;function at remote 0x7f80dcc51f28&gt;, nfev=2149842846, nlu=0, atol=&lt;numpy.ndarray at remote 0x7f80dcc537b0&gt;, h_abs=&lt;numpy.float64 at remote 0x7f80dccdb2d0&gt;, rtol=&lt;float at remote 0x7f809c0ca6a8&gt;, n=6, max_step=&lt;float at remote 0x7f8184eb3420&gt;) at remote 0x7f80a0c9bb00&gt;, t=&lt;float at remote 0x7f80dcc71510&gt;, y=&lt;numpy.ndarray at remote 0x7f80dcc53350&gt;, max_step=&lt;float at remote 0x7f8184eb3420&gt;, rtol=&lt;float at remote 0x7f...(truncated)
    self.B, self.C, self.E, self.K)
(gdb) py-print self
local 'self' = &lt;RK45(t=&lt;float at remote 0x7f80dcc71510&gt;, direction=&lt;numpy.float64 at remote 0x7f80dccc5ee8&gt;, K=&lt;numpy.ndarray at remote 0x7f80dcc53d00&gt;, vectorized=False, _fun=&lt;function at remote 0x7f80dcc51f28&gt;, y_old=None, status='running', f=&lt;numpy.ndarray at remote 0x7f80dcc53a80&gt;, y=&lt;numpy.ndarray at remote 0x7f80dcc53350&gt;, fun_vectorized=&lt;function at remote 0x7f80dcc51950&gt;, t_old=None, fun=&lt;function at remote 0x7f80dcc51d90&gt;, t_bound=&lt;float at remote 0x7f80dcc71090&gt;, njev=0, fun_single=&lt;function at remote 0x7f80dcc51f28&gt;, nfev=2149842846, nlu=0, atol=&lt;numpy.ndarray at remote 0x7f80dcc537b0&gt;, h_abs=&lt;numpy.float64 at remote 0x7f80dccdb2d0&gt;, rtol=&lt;float at remote 0x7f809c0ca6a8&gt;, n=6, max_step=&lt;float at remote 0x7f8184eb3420&gt;) at remote 0x7f80a0c9bb00&gt;
(gdb) py-print t
local 't' = &lt;float at remote 0x7f80dcc71510&gt;
(gdb) p * (PyFloatObject*) 0x7f80dcc71510
$36 = {ob_base = {ob_refcnt = 7, ob_type = 0xa3c2c0 &lt;PyFloat_Type&gt;}, ob_fval = 0}
(gdb) py-print y
local 'y' = &lt;numpy.ndarray at remote 0x7f80dcc53350&gt;
(gdb) p *(PyArrayObject_fields*) 0x7f80dcc53350
$37 = {ob_base = {ob_refcnt = 12, ob_type = 0x7f81838f7180},
  data = 0x9361740 &quot;4D\350\372^\a\005\300\320:\240v&lt;M\374?\221\361\240\266#\343\347\277\250\037\335d\323\271\376?O\210g\377\221\331\307?\332\v+1\227\315\301?\320\246\235\006&quot;, nd = 1, dimensions = 0x7397860,
  strides = 0x7397868, base = &lt;numpy.ndarray at remote 0x7f80dcc53710&gt;, descr = 0x7f81838f7ec0, flags = 1283,
  weakreflist = 0x0}
(gdb) p ((PyArrayObject_fields*) 0x7f80dcc53350)-&gt;dimensions
$38 = (npy_intp *) 0x7397860
(gdb)
$39 = (npy_intp *) 0x7397860
(gdb) p *((PyArrayObject_fields*) 0x7f80dcc53350)-&gt;dimensions
$40 = 6
(gdb) p (float[6])(((PyArrayObject_fields*) 0x7f80dcc53350)-&gt;data)
$41 = {2.19183852e-33, 0, 1.40129846e-45, 0, 1.39532281e-34, 0}
(gdb) py-print h_abs
local 'h_abs' = &lt;numpy.float64 at remote 0x7f80dccdbca8&gt;
(gdb) p * (PyFloatObject*) 0x7f80dccdb2d0
$42 = {ob_base = {ob_refcnt = 1, ob_type = 0x7f8183909600}, ob_fval = -nan(0x8000000000000)}
(gdb) p * (PyFloatObject*) 0x7f80dccdbca8
$43 = {ob_base = {ob_refcnt = 1, ob_type = 0x7f8183909600}, ob_fval = nan(0x8000000000000)}
(gdb) py-print max_step
local 'max_step' = &lt;float at remote 0x7f8184eb3420&gt;
(gdb) p * (PyFloatObject*) 0x7f8184eb3420
$44 = {ob_base = {ob_refcnt = 246, ob_type = 0xa3c2c0 &lt;PyFloat_Type&gt;}, ob_fval = inf}
(gdb) py-print min_step
local 'min_step' = &lt;numpy.float64 at remote 0x7f80dcd0f360&gt;
(gdb) p * (PyFloatObject*) 0x7f80dcd0f360
$45 = {ob_base = {ob_refcnt = 1, ob_type = 0x7f8183909600}, ob_fval = 4.9406564584124654e-323}
(gdb) p * (PyFloatObject*) 0x7f80dccc5ee8
$46 = {ob_base = {ob_refcnt = 1, ob_type = 0x7f8183909600}, ob_fval = 1}
(gdb) py-print order
local 'order' = 4
(gdb) py-print step_accepted
local 'step_accepted' = False
(gdb) py-print h
local 'h' = &lt;numpy.float64 at remote 0x7f80dccd4270&gt;
(gdb) p * (PyFloatObject*) 0x7f80dccd4270
$47 = {ob_base = {ob_refcnt = 3, ob_type = 0x7f8183909600}, ob_fval = -nan(0x8000000000000)}
(gdb) p * (PyFloatObject*) 0x7f80dcc71090
$48 = {ob_base = {ob_refcnt = 2, ob_type = 0xa3c2c0 &lt;PyFloat_Type&gt;}, ob_fval = 0.02}
(gdb) py-print t_new
local 't_new' = &lt;numpy.float64 at remote 0x7f80dcdacf18&gt;
(gdb) p * (PyFloatObject*) 0x7f80dcdacf18
$49 = {ob_base = {ob_refcnt = 1, ob_type = 0x7f8183909600}, ob_fval = -nan(0x8000000000000)}"><pre class="notranslate"><code class="notranslate">(gdb) py-bt
Traceback (most recent call first):
[...]
  File "[...]/envs/model/model.py", line 215, in _tire_dynamics_rear
    * (np.tan(alpha)/(1+K))**2 )
  File "[...]/envs/model/model.py", line 83, in _dynamics
    F_xr, F_yr = self._tire_dynamics_rear(v_x, cmd_vx, alpha_r)
  File "[...]/envs/model/model.py", line 55, in &lt;lambda&gt;
    fun=(lambda t, X: self._dynamics(X, t, U)),
  File "/usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/base.py", line 21, in fun_wrapped
    return np.asarray(fun(t, y), dtype=dtype)
  File "/usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/base.py", line 139, in fun
    return self.fun_single(t, y)
  File "/usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/rk.py", line 70, in rk_step
    K[s + 1] = fun(t + c * h, y + dy)
  File "/usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/rk.py", line 143, in _step_impl
    self.B, self.C, self.E, self.K)
  File "/usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/base.py", line 182, in step
    success, message = self._step_impl()
  File "/usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/ivp.py", line 502, in solve_ivp
    message = solver.step()
  File "[...]/envs/model/model.py", line 56, in state_transition
    t_span=t, y0=X, atol=1e-5)
  [...]
(gdb) py-up
[...]
(gdb) py-up
#24 Frame 0x61ffa18, for file /usr/local/lib/python3.5/dist-packages/scipy/integrate/_ivp/rk.py, line 143, in _step_impl (self=&lt;RK45(t=&lt;float at remote 0x7f80dcc71510&gt;, direction=&lt;numpy.float64 at remote 0x7f80dccc5ee8&gt;, K=&lt;numpy.ndarray at remote 0x7f80dcc53d00&gt;, vectorized=False, _fun=&lt;function at remote 0x7f80dcc51f28&gt;, y_old=None, status='running', f=&lt;numpy.ndarray at remote 0x7f80dcc53a80&gt;, y=&lt;numpy.ndarray at remote 0x7f80dcc53350&gt;, fun_vectorized=&lt;function at remote 0x7f80dcc51950&gt;, t_old=None, fun=&lt;function at remote 0x7f80dcc51d90&gt;, t_bound=&lt;float at remote 0x7f80dcc71090&gt;, njev=0, fun_single=&lt;function at remote 0x7f80dcc51f28&gt;, nfev=2149842846, nlu=0, atol=&lt;numpy.ndarray at remote 0x7f80dcc537b0&gt;, h_abs=&lt;numpy.float64 at remote 0x7f80dccdb2d0&gt;, rtol=&lt;float at remote 0x7f809c0ca6a8&gt;, n=6, max_step=&lt;float at remote 0x7f8184eb3420&gt;) at remote 0x7f80a0c9bb00&gt;, t=&lt;float at remote 0x7f80dcc71510&gt;, y=&lt;numpy.ndarray at remote 0x7f80dcc53350&gt;, max_step=&lt;float at remote 0x7f8184eb3420&gt;, rtol=&lt;float at remote 0x7f...(truncated)
    self.B, self.C, self.E, self.K)
(gdb) py-print self
local 'self' = &lt;RK45(t=&lt;float at remote 0x7f80dcc71510&gt;, direction=&lt;numpy.float64 at remote 0x7f80dccc5ee8&gt;, K=&lt;numpy.ndarray at remote 0x7f80dcc53d00&gt;, vectorized=False, _fun=&lt;function at remote 0x7f80dcc51f28&gt;, y_old=None, status='running', f=&lt;numpy.ndarray at remote 0x7f80dcc53a80&gt;, y=&lt;numpy.ndarray at remote 0x7f80dcc53350&gt;, fun_vectorized=&lt;function at remote 0x7f80dcc51950&gt;, t_old=None, fun=&lt;function at remote 0x7f80dcc51d90&gt;, t_bound=&lt;float at remote 0x7f80dcc71090&gt;, njev=0, fun_single=&lt;function at remote 0x7f80dcc51f28&gt;, nfev=2149842846, nlu=0, atol=&lt;numpy.ndarray at remote 0x7f80dcc537b0&gt;, h_abs=&lt;numpy.float64 at remote 0x7f80dccdb2d0&gt;, rtol=&lt;float at remote 0x7f809c0ca6a8&gt;, n=6, max_step=&lt;float at remote 0x7f8184eb3420&gt;) at remote 0x7f80a0c9bb00&gt;
(gdb) py-print t
local 't' = &lt;float at remote 0x7f80dcc71510&gt;
(gdb) p * (PyFloatObject*) 0x7f80dcc71510
$36 = {ob_base = {ob_refcnt = 7, ob_type = 0xa3c2c0 &lt;PyFloat_Type&gt;}, ob_fval = 0}
(gdb) py-print y
local 'y' = &lt;numpy.ndarray at remote 0x7f80dcc53350&gt;
(gdb) p *(PyArrayObject_fields*) 0x7f80dcc53350
$37 = {ob_base = {ob_refcnt = 12, ob_type = 0x7f81838f7180},
  data = 0x9361740 "4D\350\372^\a\005\300\320:\240v&lt;M\374?\221\361\240\266#\343\347\277\250\037\335d\323\271\376?O\210g\377\221\331\307?\332\v+1\227\315\301?\320\246\235\006", nd = 1, dimensions = 0x7397860,
  strides = 0x7397868, base = &lt;numpy.ndarray at remote 0x7f80dcc53710&gt;, descr = 0x7f81838f7ec0, flags = 1283,
  weakreflist = 0x0}
(gdb) p ((PyArrayObject_fields*) 0x7f80dcc53350)-&gt;dimensions
$38 = (npy_intp *) 0x7397860
(gdb)
$39 = (npy_intp *) 0x7397860
(gdb) p *((PyArrayObject_fields*) 0x7f80dcc53350)-&gt;dimensions
$40 = 6
(gdb) p (float[6])(((PyArrayObject_fields*) 0x7f80dcc53350)-&gt;data)
$41 = {2.19183852e-33, 0, 1.40129846e-45, 0, 1.39532281e-34, 0}
(gdb) py-print h_abs
local 'h_abs' = &lt;numpy.float64 at remote 0x7f80dccdbca8&gt;
(gdb) p * (PyFloatObject*) 0x7f80dccdb2d0
$42 = {ob_base = {ob_refcnt = 1, ob_type = 0x7f8183909600}, ob_fval = -nan(0x8000000000000)}
(gdb) p * (PyFloatObject*) 0x7f80dccdbca8
$43 = {ob_base = {ob_refcnt = 1, ob_type = 0x7f8183909600}, ob_fval = nan(0x8000000000000)}
(gdb) py-print max_step
local 'max_step' = &lt;float at remote 0x7f8184eb3420&gt;
(gdb) p * (PyFloatObject*) 0x7f8184eb3420
$44 = {ob_base = {ob_refcnt = 246, ob_type = 0xa3c2c0 &lt;PyFloat_Type&gt;}, ob_fval = inf}
(gdb) py-print min_step
local 'min_step' = &lt;numpy.float64 at remote 0x7f80dcd0f360&gt;
(gdb) p * (PyFloatObject*) 0x7f80dcd0f360
$45 = {ob_base = {ob_refcnt = 1, ob_type = 0x7f8183909600}, ob_fval = 4.9406564584124654e-323}
(gdb) p * (PyFloatObject*) 0x7f80dccc5ee8
$46 = {ob_base = {ob_refcnt = 1, ob_type = 0x7f8183909600}, ob_fval = 1}
(gdb) py-print order
local 'order' = 4
(gdb) py-print step_accepted
local 'step_accepted' = False
(gdb) py-print h
local 'h' = &lt;numpy.float64 at remote 0x7f80dccd4270&gt;
(gdb) p * (PyFloatObject*) 0x7f80dccd4270
$47 = {ob_base = {ob_refcnt = 3, ob_type = 0x7f8183909600}, ob_fval = -nan(0x8000000000000)}
(gdb) p * (PyFloatObject*) 0x7f80dcc71090
$48 = {ob_base = {ob_refcnt = 2, ob_type = 0xa3c2c0 &lt;PyFloat_Type&gt;}, ob_fval = 0.02}
(gdb) py-print t_new
local 't_new' = &lt;numpy.float64 at remote 0x7f80dcdacf18&gt;
(gdb) p * (PyFloatObject*) 0x7f80dcdacf18
$49 = {ob_base = {ob_refcnt = 1, ob_type = 0x7f8183909600}, ob_fval = -nan(0x8000000000000)}
</code></pre></div>
<p dir="auto">I am not able to share all the code, but the <code class="notranslate">model.py</code> file is essentially the same as <a href="https://github.com/r-pad/aa_simulation/blob/990f548047/envs/model/model.py">https://github.com/r-pad/aa_simulation/blob/990f548047/envs/model/model.py</a></p>
<h3 dir="auto">Describe your issue.</h3>
<p dir="auto"><code class="notranslate">scipy.integrate._ivp.common.select_initial_step()</code> returns <code class="notranslate">numpy.nan</code> if <code class="notranslate">max(abs(y0)) == numpy.inf</code></p>
<p dir="auto">This causes programs to lock, when calling <code class="notranslate">solve_ivp()</code> with such an initial value.</p>
<p dir="auto">Probably, <code class="notranslate">select_initial_step()</code> (or its callers) should raise an Error in that situation.</p>
<h3 dir="auto">Reproducing Code Example</h3>
<div class="highlight highlight-source-python notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="Call `solve_ivp()` with at least one component of `y0` being infinite (probably a NaN component would cause the same problem)."><pre class="notranslate"><span class="pl-v">Call</span> <span class="pl-s">`solve_ivp()`</span> <span class="pl-k">with</span> <span class="pl-s1">at</span> <span class="pl-s1">least</span> <span class="pl-s1">one</span> <span class="pl-s1">component</span> <span class="pl-s1">of</span> <span class="pl-s">`y0`</span> <span class="pl-s1">being</span> <span class="pl-en">infinite</span> (<span class="pl-s1">probably</span> <span class="pl-s1">a</span> <span class="pl-v">NaN</span> <span class="pl-s1">component</span> <span class="pl-s1">would</span> <span class="pl-s1">cause</span> <span class="pl-s1">the</span> <span class="pl-s1">same</span> <span class="pl-s1">problem</span>).</pre></div>
<h3 dir="auto">Error message</h3>
<div class="highlight highlight-source-shell notranslate position-relative overflow-auto" dir="auto" data-snippet-clipboard-copy-content="None, program goes into infinite loop, when calling `solve_ivp()`"><pre class="notranslate">None, program goes into infinite loop, when calling <span class="pl-en">`solve_ivp</span>()<span class="pl-s"><span class="pl-pds">`</span></span></pre></div>
<h3 dir="auto">SciPy/NumPy/Python version information</h3>
<p dir="auto">1.9.3 1.23.5 sys.version_info(major=3, minor=10, micro=6, releaselevel='final', serial=0)</p>