Correct behavior of default u at trajectory (#226)

* Correct behavior of default u at trajectory

* fix test suite (hopefully)

* fucking poincare map doesn't work i've spent 4 hours on it already

* holy fucking shit i finally fixed it oh my gopd

* changelog

* finally i truly fixed it fucking hell

Author: Datseris

CHANGELOG.md

@@ -1,3 +1,15 @@
+# v3.12.0
+
+- Crucial bugfix of the `trajectory` function. While it is documented that the default value
+of `u` is the `current_state(ds)`, the source code was incorrectly using `initial_state(ds)`. This is now fixed and `current_state` is used as was stated. Note that `t0` remains as `initial_time` as this is how it was documented.
+
+- The `initial_state(PoincareMap)` used to return the initial state of the parent system
+  which was almost never on the Poincare plane. Now it properly returns the first state
+  on the Poincare plane, which coincides with `current_state(PoincareMap)` right after
+  constructing the map.
+
+- `PoincareMap` had extra clarifications in the docstring about what's going on.
+
 # v3.11.0
 
 Brand new dynamical system `CoupledSDEs` that represents stochastic differential equations. It also comes with a dedicated documentation page.

Project.toml

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

src/core/trajectory.jl

@@ -33,7 +33,7 @@ trajectory (but remember to convert the output of `solve` to a `StateSpaceSet`).
   Note: if you mix integer and symbolic indexing be sure to initialize the array
   as `Any` so that integers `1, 2, ...` are not converted to symbolic expressions.
 """
-function trajectory(ds::DynamicalSystem, T, u0 = initial_state(ds);
+function trajectory(ds::DynamicalSystem, T, u0 = current_state(ds);
         save_idxs = nothing, t0 = initial_time(ds), kwargs...
     )
     accessor = svector_access(save_idxs)

src/derived_systems/poincare/poincaremap.jl

@@ -28,7 +28,8 @@ See also [`StroboscopicMap`](@ref), [`poincaresos`](@ref).
 * `direction = -1`: Only crossings with `sign(direction)` are considered to belong to
   the surface of section. Negative direction means going from less than ``b``
   to greater than ``b``.
-* `u0 = nothing`: Specify an initial state.
+* `u0 = nothing`: Specify an initial state. If `nothing` it is the
+  `current_state(ds)`.
 * `rootkw = (xrtol = 1e-6, atol = 1e-8)`: A `NamedTuple` of keyword arguments
   passed to `find_zero` from [Roots.jl](https://github.com/JuliaMath/Roots.jl).
 * `Tmax = 1e3`: The argument `Tmax` exists so that the integrator can terminate instead
@@ -74,11 +75,16 @@ and root finding from Roots.jl, to create a high accuracy estimate of the sectio
 4. For the special case of `plane` being a `Tuple{Int, <:Real}`, a special `reinit!` method
    is allowed with input state of length `D-1` instead of `D`, i.e., a reduced state already
    on the hyperplane that is then converted into the `D` dimensional state.
+5. The `initial_state(pmap)` returns the state initial state of the map. This is not `u0`
+   because `u0` is evolved forwards until it resides on the Poincaré plane.
+6. In the `reinit!` function, the `t0` keyword denotes the starting time of the
+   continuous time dynamical system, as the starting time of the `PoincareMap` is by
+   definition always 0.
 
 ## Example
 
 ```julia
-using DynamicalSystemsBase
+using DynamicalSystemsBase, PredefinedDynamicalSystems
 ds = Systems.rikitake(zeros(3); μ = 0.47, α = 1.0)
 pmap = poincaremap(ds, (3, 0.0))
 step!(pmap)
@@ -102,6 +108,7 @@ mutable struct PoincareMap{I<:ContinuousTimeDynamicalSystem, F, P, R, V} <: Disc
     # (i.e., given a D-1 dimensional state, create the full D-dimensional state)
     dummy::Vector{Float64}
     diffidxs::Vector{Int}
+    state_initial::V
 end
 Base.parent(pmap::PoincareMap) = pmap.ds
 
@@ -122,9 +129,10 @@ function PoincareMap(
     diffidxs = _indices_on_poincare_plane(plane, D)
 	pmap = PoincareMap(
         ds, plane_distance, planecrossing, Tmax, rootkw,
-        v, current_time(ds), 0, dummy, diffidxs
+        v, current_time(ds), 0, dummy, diffidxs, recursivecopy(v),
     )
-    step!(pmap)
+    step!(pmap) # this ensures that the state is on the hyperplane
+    pmap.state_initial = recursivecopy(current_state(pmap))
     pmap.t = 0 # first step is 0
     return pmap
 end
@@ -147,6 +155,7 @@ end
 initial_time(pmap::PoincareMap) = 0
 current_time(pmap::PoincareMap) = pmap.t
 current_state(pmap::PoincareMap) = pmap.state_on_plane
+initial_state(pmap::PoincareMap) = pmap.state_initial
 
 """
     current_crossing_time(pmap::PoincareMap) → tcross
@@ -169,8 +178,9 @@ end
 SciMLBase.step!(pmap::PoincareMap, n::Int, s = true) = (for _ ∈ 1:n; step!(pmap); end; pmap)
 
 function SciMLBase.reinit!(pmap::PoincareMap, u0::AbstractArray = initial_state(pmap);
-        t0 = initial_time(pmap), p = current_parameters(pmap)
+        t0 = initial_time(pmap.ds), p = current_parameters(pmap)
     )
+    uc = current_state(pmap)
     if length(u0) == dimension(pmap)
 	    u = u0
     elseif length(u0) == dimension(pmap) - 1
@@ -179,8 +189,13 @@ function SciMLBase.reinit!(pmap::PoincareMap, u0::AbstractArray = initial_state(
         error("Dimension of state for `PoincareMap` reinit is inappropriate.")
     end
     reinit!(pmap.ds, u; t0, p)
-    step!(pmap) # always step once to reach the PSOS
-    pmap.t = 0 # first step is 0
+    # Check if we have to step the map or not
+    if u0 == initial_state(pmap)
+        pmap.state_on_plane = copy(u0)
+    elseif u0 ≠ uc # if we were on current state (like in `trajectory`) then we don't step
+        step!(pmap) # step once to reach the PSOS
+    end
+    pmap.t = 0 # first step is always 0
     pmap
 end
 
@@ -223,20 +238,18 @@ function poincare_step!(ds, plane_distance, planecrossing, Tmax, rootkw)
     # Otherwise evolve until juuuuuust crossing the plane
     tprev = t0
     while side < 0
-        (current_time(ds) - t0) > Tmax && break
+        (current_time(ds) - t0) > Tmax && return (nothing, nothing)
         tprev = current_time(ds)
         step!(ds)
         side = planecrossing(current_state(ds))
     end
-    while side ≥ 0
-        (current_time(ds) - t0) > Tmax && break
+    while side ≥ 0 # when side becomes negative again we've just crossed the plane!
+        (current_time(ds) - t0) > Tmax && return (nothing, nothing)
         tprev = current_time(ds)
         step!(ds)
         side = planecrossing(current_state(ds))
     end
-    # we evolved too long and no crossing, return nothing
-    (current_time(ds) - t0) > Tmax && return (nothing, nothing)
-    # Else, we're guaranteed to have time window before and after the plane
+    # We're now guaranteed to have time window before and after the plane
     time_window = (tprev, current_time(ds))
     tcross = Roots.find_zero(plane_distance, time_window, ROOTS_ALG; rootkw...)
     ucross = ds(tcross)

test/poincare.jl

@@ -21,9 +21,6 @@ u0 = [
 Γ = 0.9
 p = [μ, ν, Γ]
 
-gissinger_oop = CoupledODEs(gissinger_rule, u0, p)
-gissinger_iip = CoupledODEs(gissinger_rule_iip, recursivecopy(u0), p)
-
 # Define appropriate hyperplane for gissinger system
 plane1 = (1, 0.0)
 # I want hyperperplane defined by these two points:
@@ -34,7 +31,7 @@ gis_plane(μ) = [cross(Np(μ), Nm(μ))..., 0]
 plane2 = gis_plane(μ)
 
 function poincare_tests(ds, pmap, plane)
-    P, t = trajectory(pmap, 10)
+    P, t = trajectory(pmap, 10, initial_state(pmap))
     reinit!(ds)
     P2 = poincaresos(ds, plane, 100; rootkw = (xrtol = 1e-12, atol = 1e-12))
     @test length(P) == 11
@@ -54,10 +51,10 @@ end
 
 @testset "IIP=$(IIP) plane=$(P)" for IIP in (false, true), P in (1, 2)
     rule = !IIP ? gissinger_rule : gissinger_rule_iip
-    ds = CoupledODEs(rule, recursivecopy(u0), p)
+    ds = CoupledODEs(rule, u0, p)
     plane = P == 1 ? plane1 : plane2
     pmap = PoincareMap(ds, plane; rootkw = (xrtol = 1e-12, atol = 1e-12))
-    u0pmap = recursivecopy(current_state(pmap))
+    u0pmap = deepcopy(current_state(pmap))
     test_dynamical_system(pmap, u0pmap, p;
     idt=true, iip=IIP, test_trajectory = true,
     test_init_state_equiv=false, u0init = initial_state(pmap))
@@ -77,7 +74,7 @@ end
     u0 = fill(10.0, 3)
     p = [10, 28, 8/3]
     plane = (1, 0.0)
-    ds = CoupledODEs(lorenz_rule, recursivecopy(u0), p)
+    ds = CoupledODEs(lorenz_rule, deepcopy(u0), p)
 
     pmap = PoincareMap(ds, plane; rootkw = (xrtol = 1e-12, atol = 1e-12))
     P, t = trajectory(pmap, 10)
@@ -91,6 +88,7 @@ end
 end
 
 @testset "poincare of dataset" begin
+    gissinger_oop = CoupledODEs(gissinger_rule, u0, p)
     X, t = trajectory(gissinger_oop, 1000.0)
     A = poincaresos(X, plane1)
     @test dimension(A) == 3

test/test_system_function.jl

@@ -114,10 +114,10 @@ function test_dynamical_system(ds, u0, p0; idt, iip,
                 @test Y[1] == X[end]
 
                 # obtain only first variable
-                Z, t = trajectory(ds, 10; save_idxs = [1])
+                Z, t = trajectory(ds, 10, u0init; save_idxs = [1])
                 @test length(Z) == 11
                 @test dimension(Z) == 1
-                @test Z[1][1] == u0[1]
+                @test Z[1][1] == u0init[1]
             else
                 reinit!(ds)
                 @test current_state(ds) == u0
@@ -133,7 +133,7 @@ function test_dynamical_system(ds, u0, p0; idt, iip,
                 @test t2[1] > t[end]
 
                 # obtain only first variable
-                Z, t = trajectory(ds, 3; save_idxs = [1], Δt = 1)
+                Z, t = trajectory(ds, 3, u0; save_idxs = [1], Δt = 1)
                 @test length(Z) == 4
                 @test dimension(Z) == 1
                 @test Z[1][1] == u0[1]