From 7162bd5bb6eca0eab172e1f96a51638b294b78f4 Mon Sep 17 00:00:00 2001 From: George Datseris Date: Wed, 13 Nov 2024 15:57:42 +0000 Subject: [PATCH] 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 --- CHANGELOG.md | 12 +++++++ Project.toml | 2 +- src/core/trajectory.jl | 2 +- src/derived_systems/poincare/poincaremap.jl | 39 ++++++++++++++------- test/poincare.jl | 12 +++---- test/test_system_function.jl | 6 ++-- 6 files changed, 48 insertions(+), 25 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index fb7d18f1..d2ed5395 100644 --- a/CHANGELOG.md +++ b/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. diff --git a/Project.toml b/Project.toml index e9371dc9..0a9ad46d 100644 --- a/Project.toml +++ b/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" diff --git a/src/core/trajectory.jl b/src/core/trajectory.jl index 0455eb41..acf59c9c 100644 --- a/src/core/trajectory.jl +++ b/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) diff --git a/src/derived_systems/poincare/poincaremap.jl b/src/derived_systems/poincare/poincaremap.jl index 9d811615..b119f7a0 100644 --- a/src/derived_systems/poincare/poincaremap.jl +++ b/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) diff --git a/test/poincare.jl b/test/poincare.jl index 3e894326..2f3bcaf7 100644 --- a/test/poincare.jl +++ b/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 diff --git a/test/test_system_function.jl b/test/test_system_function.jl index 556cb96a..51b46dd4 100644 --- a/test/test_system_function.jl +++ b/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]