BUG: integrate: make select_initial_step aware of integration interval

[valentineap]

Reference issue

Closes gh-17341
Closes gh-8848
Closes gh-9198

What does this implement/fix?

#17341 reports that solve_ivp(func, (t0, t_bound), f0) would evaluate func() outside the range (t0, t_bound) if |t_bound-t0| is small.

This issue was traced to select_initial_step, which was not aware of the far limit, t_bound, and evaluates func(t0+h0) for some internally-determined value of h0.

Fix:

• add optional argument t_bound to select_initial_step and ensure that that h0<=|t_bound-t0|
• modify all solvers within _ivp that use select_initial_step to take advantage of this argument

Additional information

If t0 == t_bound, select_initial_step will now return 0. This is a logical result but has potential for divide-by-zero issues.

[Patol75]

It is not clear to me if what you are dealing with here is an actual bug or an implementation choice, given that code execution does not raise. I can see that @nmayorov contributed most of the code you are modifying here, so he probably knows. In any case, why would not just changing the return statement in select_initial_step to return min(100 * h0, h1, abs(t_bound - t0)), in addition to incorporating the t_bound argument, be sufficient?

[valentineap]

Hi @Patol75, in my view it is a clear bug, as there is no requirement that the function func is well-defined (mathematically or computationally) outside the domain of integration. The integrator should never attempt to evaluate func outside the domain, as doing so may result in errors or other unexpected behaviour.

Thus the need for a more complex solution than you propose -- the point is to avoid ever calling func for these values, not just to avoid using them.

[Patol75]

That is a fair point indeed. The current implementation would require users to handle function calls outside of the expected domain in the case of exotic functions, which does not seem ideal. Nonetheless, I am still wondering if this is a known issue with regards to the implemented scheme, in which case applying a patch seems reasonable, or if there is a mistake in the implementation of the scheme.

@mdhaber Do you mind sharing your opinion on these changes?

[valentineap]

Just to highlight: this is a fundamental issue for many physically-motivated problems, rather than just exotic/esoteric applications. Often one wishes to integrate a function that depends on the physical properties of some body, and these are (obviously!) not defined outside the boundary of that body.

[laurent90git]

Maybe t_bound should not be optional, since this function is only used by solve_ivp solvers which are all affected by the issue you want to resolve.
By the way, the python-wrapped Fortran solver LSODA os solve_ivp also ensures that t_bound is not exceeded.

[valentineap]

I'd tend to agree -- I just wasn't sure whether there were any external uses that would be broken by an interface change.

[mdhaber]

Despite not having an underscore, select_initial_step is not public, so if it doesn't break any uses of the function within SciPy, it is arguably OK to change this without a default. Are all the uses of select_initial_step inside SciPy addressed here?

[laurent90git]

Not sure if returning 0 will not affect the calling solver. Would it be best to simply raise an error, if that is not already done before ?

[valentineap]

In principle, yes -- but I think that would then necessitate refactoring some of the integration routines to handle tbound == t0 as a special case; at present they call select_initial_step and do various other things before determining that they should return 0.

[valentineap]

I think the current version works ok with the solvers as they stand -- the existing test-suite includes the case tbound==t0.

[Patol75]

I have just had another look at this PR. I believe there are two main issues highlighted here:

• select_initial_step can return values that exceed the integration span
• Null integration spans are left unchecked

For the second item, it is true that integrate solvers can handle a null span and return without errors claiming that they have successfully reached the end of the integration interval. I am not quite sure if this is a meaningful way to deal with such a situation. The immediate alternative I am seeing is to add a check in check_arguments to interrupt execution if there is basically nothing to be done.

For the first item, select_initial_step is only called by RK-family solvers, Radau, and BDF — LSODA has its own first-step selection algorithm, which does not suffer from the issue highlighted here, as reported by @laurent90git. The five mentioned solvers all execute the following code:

        if first_step is None:
            self.h_abs = select_initial_step(
                self.fun, self.t, self.y, self.f, self.direction,
                self.error_estimator_order, self.rtol, self.atol)
        else:
            self.h_abs = validate_first_step(first_step, t0, t_bound)

We could change the above to

        if first_step is None:
            first_step = select_initial_step(
                self.fun, self.t, self.y, self.f, self.direction,
                self.error_estimator_order, self.rtol, self.atol)
        self.h_abs = validate_first_step(first_step, t0, t_bound)

Within validate_first_step, we have

    if first_step <= 0:
        raise ValueError("`first_step` must be positive.")
    if first_step > np.abs(t_bound - t0):
        raise ValueError("`first_step` exceeds bounds.")
    return first_step

This could change to

    if first_step <= 0:
        raise ValueError("`first_step` must be positive.")
    return min(first_step, np.abs(t_bound - t0))

What would you think of such a change? If you do not believe it is a good idea, then any objections to the changes proposed by @valentineap?

[valentineap]

As a matter of general principle: it seems reasonable to handle the case where t_0 == t_bound at the earliest possible opportunity -- it is axiomatic that the result should be zero and no work needs to be done to discover that.

As for @Patol75's suggested fix to the validate_first_step function: this is a straightforward fix to the problem I identified. However, the downside is that if the user has chosen to explicitly specify their own value for first_step but has provided an inappropriate value, this will silently be overridden. I am not sure this is desirable -- if a user seeks control over this option they would probably prefer to be alerted that their chosen value is unsuitable.

[mdhaber]

Do you mind sharing your opinion on these changes?

Oops, I missed the ping.

Does this address gh-9198 and/or gh-8848? If so, I'd support such a change.

That said, I haven't taken a close look at how this works, and I'd need to ask some basic questions to facilitate review. LMK if you'd like me to take a closer look.

[mdhaber]

Nit: the logic h = min(h, (t_bound-t0)/direction) is used below. Consider using the same here idea here.

[valentineap]

Fixed.

[mdhaber]

Despite not having an underscore, select_initial_step is not public, so if it doesn't break any uses of the function within SciPy, it is arguably OK to change this without a default. Are all the uses of select_initial_step inside SciPy addressed here?

[Patol75]

Do you mind sharing your opinion on these changes?

Oops, I missed the ping.

Does this address gh-9198 and/or gh-8848? If so, I'd support such a change.

That said, I haven't taken a close look at how this works, and I'd need to ask some basic questions to facilitate review. LMK if you'd like me to take a closer look.

Yes, the changes here would prevent any function evaluation outside the integration range specified by the user. It is relatively straightforward to do so, either in the way the PR proposes or the way my message above suggests. But, as I stressed in my first reply here, it is not necessarily obvious for such behaviour to be a bug — consider, for example, @WarrenWeckesser reply in gh-9198.

[mdhaber]

I'll limit my reply to the original concern: "solve_ivp(func, (t0, t_bound), f0) would evaluate func() outside the range (t0, t_bound)". You mentioned two separate issues in #17348 (comment), but I haven't considered them.

But, as I stressed in my first reply here, it is not necessarily obvious for such behaviour to be a bug

@drhagen was involved in creating the solvers and addressed that in #9198 (comment), stating that it was a bug. After more carefully reading the code and related issues, I think it is something to address, whether it is called a bug or not, but I lean toward "bug" because it seems to have been inadvertent and does not seem to be a useful "feature".

Looking at the original discussion (e.g. here), t_bound was initially known as t_crit, and it was supposed to be the time beyond which the integrand would _not be evaluated because the integrate might not be well defined beyond this time. I agree that this intent is not 100% explicit in the description "Boundary time — the integration won’t continue beyond it.", but it seems to be clarified in the code. There are several places where it is used to prevent function evaluation beyond the specified time.

scipy/scipy/integrate/_ivp/rk.py

Lines 138 to 139 in c1ed5ec

	if self.direction * (t_new - self.t_bound) > 0:
	t_new = self.t_bound

In any case, SUNDIALS apparently has an option to turn this behavior on and off. We don't, so something should change.

If there is a demonstrable benefit to evaluating the integrand beyond t_bound, let's reconsider, but right now, it seems like the intent was to prevent evaluation of the integrand beyond t_bound (not just to prevent steps after it is exceeded). This has come up in three issues, and this PR sounds like a reasonable solution to the problem, so I'd support it. In addition, I would prevent function evaluations beyond t0 + max_step.

(One could argue that only preventing function evaluations beyond max_step from the initial/current time would be sufficient. I agree that it would be enough to allow the user to prevent select_initial_step from evaluating the integrand beyond t_bound, but I still think the default should be to prevent this without requiring max_step.)

[valentineap]

I completely agree with @mdhaber. I think it is important to highlight that this issue (at least #17348; I haven't fully investigated the various other linked issues) is solely associated with the select_initial_step routine. If a user specifies their own step this gets passed into validate_initial_step which raises an error if this is larger than the integration domain.

[mdhaber]

@valentineap Great. In that case, please consider:

• adjusting this PR to ensure that max_step is respected by select_initial_step and demonstrating this with a unit test
• investigating whether you believe this closes solve_ivp RK45 can evaluate the function at times later than t_span[1] #9198 and/or integrate.solve_ivp try to evaluate the function with much larger value than t_span #8848.
• Add unit tests, preferably one for each issue that this closes unless the issues are essentially identical. For instance, integrate.solve_ivp try to evaluate the function with much larger value than t_span #8848 has some simple code that could be adapted to become a unit test. You might check that the code runs without emitting a warning and that the final values of y are close to those produced by odeint.

[valentineap]

Closes #8848

Closes #9198

[valentineap]

@mdhaber OK, finally getting around to dealing with this..

This PR modifies scipy.integrate._ivp.common.select_initial_step() to ensure that:

• The inital step is no larger than the requested integration interval, and
• The initial step is no larger than the specified maximum integration step.

It will close #8848, #9198, and #17341. Example code from each of these issues is now included in the test suite: test_ivp.test_tbound_larger_interval(), test_ivp.test_tbound_oscillator() and test_ivp.test_tbound_small_interval() respectively. Additionally test_ivp.test_inital_maxstep() verifies that max_step is honoured.

Enforcing max_step within select_initial_step broke one existing test: test_ivp.test_max_step(). This includes a test that integration fails if max_step is too small. It assumed failure on the first step. With the change, the first step succeeds and failure occurs on the second step. It seems to me that this is still acceptable from the point of view of 'passing' the test, and I have modified it accordingly.

I think this is ready to be considered for merge into main.

[valentineap]

@nmayorov Thanks for the comments -- I believe I have now addressed all of these. Let me know if you think I've missed anything.

@j-bowhay I have updated the assert statements in the tests that are affected by this PR; I haven't looked at any of the pre-existing tests.

[nmayorov]

Looks much cleaner now.

What about a test for t0 == t_bound, wouldn't it be too hard to implement it?

[valentineap]

@nmayorov I had added test_zero_interval which checks the case t0 == t_bound. Perhaps I'm misunderstanding what you wanted?

I've fixed the 0. -> 0.0 issue.

[nmayorov]

@nmayorov I had added test_zero_interval which checks the case t0 == t_bound. Perhaps I'm misunderstanding what you wanted?

I've fixed the 0. -> 0.0 issue.

Just didn't spot it quickly, excellent.

[valentineap]

Great -- so unless there's anything I've missed, I think that's all of @nmayorov's comments addressed.

I don't think the 'failing' check(ci/circleci: build_scipy) is anything to do with the changes made in this PR.

@tylerjereddy Perhaps this is now suitable for merge?

[nmayorov]

Some very minor comments.

[valentineap]

Thanks @nmayorov @j-bowhay -- I've incorporated your formatting fixes.

Note that this PR should close gh-8848 and gh-9198 in addition to gh-17341 for which it is currently marked.

[j-bowhay]

@valentineap can you edit the description linking those issues too please

[valentineap]

@j-bowhay Done.

[nmayorov]

To me everything looks good.

[nmayorov]

I'll go ahead and merge it.

I believe it's a great fix closing a lot of issues. True, there is always a risk of breaking something in numerical algorithms, but this change looks rather safe and makes the t_bound and max_step arguments properly defined.

Thank you @valentineap for the work and insisting on finishing it.

@tylerjereddy can be included in 0.14.0 then.