fix infinite while loop problem in trajectory (#194)

* fix infinite while loop problem in trajectory

* only reference ODEintegrator

* AbstractODEIntegrator is tehe correct name

Author: Datseris

Project.toml

@@ -1,7 +1,7 @@
 name = "DynamicalSystemsBase"
 uuid = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
 repo = "https://github.com/JuliaDynamics/DynamicalSystemsBase.jl.git"
-version = "3.5.1"
+version = "3.5.2"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"

src/core/trajectory.jl

@@ -3,14 +3,15 @@ export trajectory
 """
     trajectory(ds::DynamicalSystem, T [, u0]; kwargs...) → X, t
 
-Return a dataset `X` that will contain the trajectory of the system `ds`,
-after evolving it for total time `T`.
-`u0` is the state given given to [`reinit!`](@ref) prior to time evolution
-and defaults to [`initial_state(ds)`](@ref).
+Evolve `ds` for a total time of `T` and return its trajectory `X`, sampled at equal time
+intervals, and corresponding time vector.
+Optionally provide a starting state `u0` which is `current_state(ds)` by default.
 
-See [`StateSpaceSet`](@ref) for info on how to use `X`.
 The returned time vector is `t = (t0+Ttr):Δt:(t0+Ttr+T)`.
 
+If time evolution diverged before `T`, the remaining of the trajectory is set
+to the last valid point.
+
 ## Keyword arguments
 
 * `Δt`:  Time step of value output. For discrete time systems it must be an integer.
@@ -69,6 +70,7 @@ function trajectory_continuous(ds, T; Dt = 0.1, Δt = Dt, Ttr = 0.0, accessor =
     for (i, t) in enumerate(tvec)
         while t > current_time(ds)
             step!(ds)
+            successful_step(ds) || break
         end
         sol[i] = SVector{X, ET}(obtain_access(ds(t), accessor))
         if !successful_step(ds)

src/core_systems/continuous_time_ode.jl

@@ -86,8 +86,8 @@ function CoupledODEs(f, u0, p = SciMLBase.NullParameters(); t0 = 0, diffeq = DEF
     prob = ODEProblem{IIP}(f, s, (T(t0), T(Inf)), p)
     return CoupledODEs(prob, diffeq)
 end
-# This preserves the referrenced MTK system
-CoupledODEs(ds::CoupledODEs, diffeq) = CoupledODEs(ODEProblem(ds), diffeq)
+# This preserves the referrenced MTK system and the originally passed diffeq kwargs
+CoupledODEs(ds::CoupledODEs, diffeq) = CoupledODEs(ODEProblem(ds), merge(ds.diffeq, diffeq))
 # Below `special_kwargs` is undocumented internal option for passing `internalnorm`
 function CoupledODEs(prob::ODEProblem, diffeq = DEFAULT_DIFFEQ; special_kwargs...)
     if haskey(special_kwargs, :diffeq)
@@ -159,8 +159,17 @@ current_state(integ::DEIntegrator) = integ.u
 current_time(integ::DEIntegrator) = integ.t
 initial_time(integ::DEIntegrator) = integ.sol.prob.tspan[1]
 
-using SciMLBase: successful_retcode, check_error
-successful_step(integ::DEIntegrator) = successful_retcode(check_error(integ))
+# For checking successful step, the `SciMLBase.step!` function checks
+# `integ.sol.retcode in (ReturnCode.Default, ReturnCode.Success) || break`.
+# But the actual API call would be `successful_retcode(check_error(integ))`.
+# The latter however is already used in `step!(integ)` so there is no reason to re-do it.
+# Besides, within DynamicalSystems.jl the integration is never expected to terminate.
+# Nevertheless here we extend explicitly only for ODE stuff because it may be that for
+# other type of DEIntegrators a different step interruption is possible.
+function successful_step(integ::SciMLBase.AbstractODEIntegrator)
+    rcode = integ.sol.retcode
+    return rcode == SciMLBase.ReturnCode.Default || rcode == SciMLBase.ReturnCode.Success
+end
 
 function set_state!(integ::DEIntegrator, u)
     if integ.u isa Array{<:Real}

src/core_systems/discrete_time_map.jl

@@ -88,6 +88,8 @@ function SciMLBase.step!(ds::DeterministicIteratedMap{true})
     return ds
 end
 function SciMLBase.step!(ds::DeterministicIteratedMap{true}, N, stop_at_tdt = true)
+    # We do not check for errors here for performance (of not computing equals to NaN
+    # or Inf all the time), because in the end detecting the error early doesn't help
     for _ in 1:N
         ds.dummy, ds.u = ds.u, ds.dummy
         ds.f(ds.u, ds.dummy, ds.p, ds.t)