From 18d002af114341f0cffc7930927fe44222face3d Mon Sep 17 00:00:00 2001 From: CompatHelper Julia Date: Sat, 21 Sep 2024 12:39:09 +0000 Subject: [PATCH 01/12] CompatHelper: bump compat for StateSpaceSets to 2, (keep existing compat) --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index f0221b8c..22ad3bba 100644 --- a/Project.toml +++ b/Project.toml @@ -21,7 +21,7 @@ OrdinaryDiffEqTsit5 = "1.1" Reexport = "1" Roots = "1, 2" SciMLBase = "1.19.5, 2" -StateSpaceSets = "1" +StateSpaceSets = "1, 2" Statistics = "1" SymbolicIndexingInterface = "0.3.4" julia = "1.9" From 7438b5c2ce64fdd944063f8c76fa5723325820a3 Mon Sep 17 00:00:00 2001 From: Datseris Date: Sun, 22 Sep 2024 11:12:54 +0100 Subject: [PATCH 02/12] allow `container` keyword --- src/core/trajectory.jl | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/src/core/trajectory.jl b/src/core/trajectory.jl index a4e3b9ca..0455eb41 100644 --- a/src/core/trajectory.jl +++ b/src/core/trajectory.jl @@ -4,18 +4,19 @@ export trajectory trajectory(ds::DynamicalSystem, T [, u0]; kwargs...) → X, t Evolve `ds` for a total time of `T` and return its trajectory `X`, sampled at equal time -intervals, and corresponding time vector. +intervals, and corresponding time vector. `X` is a `StateSpaceSet`. Optionally provide a starting state `u0` which is `current_state(ds)` by default. 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. +If time evolution diverged, or in general failed, before `T`, +the remaining of the trajectory is set to the last valid point. `trajectory` is a very simple function provided for convenience. For continuous time systems, it doesn't play well with callbacks, use `DifferentialEquations.solve` if you want a trajectory/timeseries -that works with callbacks. +that works with callbacks, or in general you want more flexibility in the generated +trajectory (but remember to convert the output of `solve` to a `StateSpaceSet`). ## Keyword arguments @@ -24,6 +25,8 @@ that works with callbacks. have access to unicode, the keyword `Dt` can be used instead. * `Ttr = 0`: Transient time to evolve the initial state before starting saving states. * `t0 = initial_time(ds)`: Starting time. +* `container = SVector`: Type of vector that will represent the state space points + that will be included in the `StateSpaceSet` output. See `StateSpaceSet` for valid options. * `save_idxs::AbstractVector`: Which variables to output in `X`. It can be any type of index that can be given to [`observe_state`](@ref). Defaults to `1:dimension(ds)` (all dynamic variables). @@ -43,7 +46,8 @@ function trajectory(ds::DynamicalSystem, T, u0 = initial_state(ds); end function trajectory_discrete(ds, T; - Dt::Integer = 1, Δt::Integer = Dt, Ttr::Integer = 0, accessor = nothing + Dt::Integer = 1, Δt::Integer = Dt, Ttr::Integer = 0, accessor = nothing, + kw... # container keyword ) ET = eltype(current_state(ds)) t0 = current_time(ds) @@ -65,10 +69,10 @@ function trajectory_discrete(ds, T; break end end - return StateSpaceSet(data), tvec + return StateSpaceSet(data; kw...), tvec end -function trajectory_continuous(ds, T; Dt = 0.1, Δt = Dt, Ttr = 0.0, accessor = nothing) +function trajectory_continuous(ds, T; Dt = 0.1, Δt = Dt, Ttr = 0.0, accessor = nothing, kw...) D = dimension(ds) t0 = current_time(ds) tvec = (t0+Ttr):Δt:(t0+T+Ttr) @@ -92,7 +96,7 @@ function trajectory_continuous(ds, T; Dt = 0.1, Δt = Dt, Ttr = 0.0, accessor = end end - return StateSpaceSet(sol), tvec + return StateSpaceSet(sol; kw...), tvec end # Util functions for `trajectory` to make indexing static vectors From d458646881be72ded9391f407a0170b65194dfbb Mon Sep 17 00:00:00 2001 From: Datseris Date: Sun, 22 Sep 2024 11:12:57 +0100 Subject: [PATCH 03/12] changelog --- CHANGELOG.md | 8 +++++++- Project.toml | 2 +- 2 files changed, 8 insertions(+), 2 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index ec799380..48d01e96 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,6 +1,12 @@ +# v3.10.0 + +- OrdinaryDiffEq.jl dependency reduced to OrdinaryDiffEqTsit5.jl. +- Updated to StateSpaceSets.jl: now the keyword `container` can be given to + the `trajectory` function. + # v3.9.0 -A new function `jacobian` that generates Jacobian rule for any `ds<:CoreDynamicalSystem` via +A new function `jacobian` that generates Jacobian rule for any `ds<:CoreDynamicalSystem` via automatic differentiation with `ForwardDiff.jl`. # v3.8.0 diff --git a/Project.toml b/Project.toml index 22ad3bba..925c0889 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.10.0" +version = "3.10.1" [deps] ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" From 814bbc5cc33e0c5cd5343c4c999202a059531874 Mon Sep 17 00:00:00 2001 From: Datseris Date: Sun, 22 Sep 2024 11:25:03 +0100 Subject: [PATCH 04/12] fix incorrect usage of size --- test/test_system_function.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/test/test_system_function.jl b/test/test_system_function.jl index 1c1017a0..f4906e1c 100644 --- a/test/test_system_function.jl +++ b/test/test_system_function.jl @@ -111,7 +111,8 @@ function test_dynamical_system(ds, u0, p0; idt, iip, # obtain only first variable Z, t = trajectory(ds, 10; save_idxs = [1]) - @test size(Z) == (11, 1) + @test length(Z) == 11 + @test dimension(Z) == 1 @test Z[1][1] == u0[1] else reinit!(ds) @@ -129,7 +130,8 @@ function test_dynamical_system(ds, u0, p0; idt, iip, # obtain only first variable Z, t = trajectory(ds, 3; save_idxs = [1], Δt = 1) - @test size(Z) == (4, 1) + @test length(Z) == 4 + @test dimension(Z) == 1 @test Z[1][1] == u0[1] end end From 50c4d88558aeb2d787b1956465cf1728e054c724 Mon Sep 17 00:00:00 2001 From: Datseris Date: Sun, 22 Sep 2024 11:30:44 +0100 Subject: [PATCH 05/12] use exclusively only v2 of SSEts --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 925c0889..9d1bd230 100644 --- a/Project.toml +++ b/Project.toml @@ -21,7 +21,7 @@ OrdinaryDiffEqTsit5 = "1.1" Reexport = "1" Roots = "1, 2" SciMLBase = "1.19.5, 2" -StateSpaceSets = "1, 2" +StateSpaceSets = "2" Statistics = "1" SymbolicIndexingInterface = "0.3.4" julia = "1.9" From 83c3805d08306cef0da94d81e4b9510785c8a573 Mon Sep 17 00:00:00 2001 From: Datseris Date: Sun, 22 Sep 2024 13:33:22 +0100 Subject: [PATCH 06/12] use dimension for poincaresos --- src/derived_systems/poincare/hyperplane.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/derived_systems/poincare/hyperplane.jl b/src/derived_systems/poincare/hyperplane.jl index 9c1f9b40..573b2c3e 100644 --- a/src/derived_systems/poincare/hyperplane.jl +++ b/src/derived_systems/poincare/hyperplane.jl @@ -84,7 +84,7 @@ are the same as in [`PoincareMap`](@ref). function poincaresos(A::AbstractStateSpaceSet, plane; direction = -1, warning = true, save_idxs = 1:dimension(A) ) - check_hyperplane_match(plane, size(A, 2)) + check_hyperplane_match(plane, dimension(A)) i = typeof(save_idxs) <: Int ? save_idxs : SVector{length(save_idxs), Int}(save_idxs...) planecrossing = PlaneCrossing(plane, direction > 0) data = poincaresos(A, planecrossing, i) From 55d3ff8ad796c2b4f5e8a39d8ea6b8393c0080ab Mon Sep 17 00:00:00 2001 From: Datseris Date: Sun, 22 Sep 2024 13:42:27 +0100 Subject: [PATCH 07/12] simplify field t of PoincareMap --- src/derived_systems/poincare/poincaremap.jl | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/src/derived_systems/poincare/poincaremap.jl b/src/derived_systems/poincare/poincaremap.jl index 604983eb..9ccd470b 100644 --- a/src/derived_systems/poincare/poincaremap.jl +++ b/src/derived_systems/poincare/poincaremap.jl @@ -19,6 +19,8 @@ A discrete time dynamical system that produces iterations over the Poincaré map of the given continuous time `ds`. This map is defined as the sequence of points on the Poincaré surface of section, which is defined by the `plane` argument. +Iterating `pmap` also mutates `ds` which is referrenced in `pmap`. + See also [`StroboscopicMap`](@ref), [`poincaresos`](@ref). ## Keyword arguments @@ -67,8 +69,8 @@ and root finding from Roots.jl, to create a high accuracy estimate of the sectio surface of section. [`initial_time`](@ref) is always `0` and [`current_time`](@ref) is current iteration number. 3. A new function [`current_crossing_time`](@ref) returns the real time corresponding - to the latest crossing of the hyperplane, which is what the [`current_state(ds)`](@ref) - corresponds to as well. + to the latest crossing of the hyperplane. The corresponding state on the hyperplane + is `current_state(pmap)` as expected. 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. @@ -95,7 +97,7 @@ mutable struct PoincareMap{I<:ContinuousTimeDynamicalSystem, F, P, R, V} <: Disc rootkw::R state_on_plane::V tcross::Float64 - t::Base.RefValue{Int} + t::Int # These two fields are for setting the state of the pmap from the plane # (i.e., given a D-1 dimensional state, create the full D-dimensional state) dummy::Vector{Float64} @@ -120,7 +122,7 @@ function PoincareMap( diffidxs = _indices_on_poincare_plane(plane, D) pmap = PoincareMap( ds, plane_distance, planecrossing, Tmax, rootkw, - v, current_time(ds), Ref(0), dummy, diffidxs + v, current_time(ds), 0, dummy, diffidxs ) step!(pmap) pmap.t[] = 0 # first step is 0 @@ -143,7 +145,7 @@ for f in (:initial_state, :current_parameters, :initial_parameters, :referrenced @eval $(f)(pmap::PoincareMap, args...) = $(f)(pmap.ds, args...) end initial_time(pmap::PoincareMap) = 0 -current_time(pmap::PoincareMap) = pmap.t[] +current_time(pmap::PoincareMap) = pmap.t current_state(pmap::PoincareMap) = pmap.state_on_plane """ @@ -160,7 +162,7 @@ function SciMLBase.step!(pmap::PoincareMap) else pmap.state_on_plane = u # this is always a brand new vector pmap.tcross = t - pmap.t[] = pmap.t[] + 1 + pmap.t += 1 return pmap end end @@ -179,7 +181,7 @@ function SciMLBase.reinit!(pmap::PoincareMap, u0::AbstractArray = initial_state( end reinit!(pmap.ds, u; t0, p) step!(pmap) # always step once to reach the PSOS - pmap.t[] = 0 # first step is 0 + pmap.t = 0 # first step is 0 pmap end From 50a5b801cb6681a35acc35de8166fd7c415f9422 Mon Sep 17 00:00:00 2001 From: Datseris Date: Sun, 22 Sep 2024 13:43:18 +0100 Subject: [PATCH 08/12] and fix incorrect AA type annotation --- src/derived_systems/poincare/poincaremap.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/derived_systems/poincare/poincaremap.jl b/src/derived_systems/poincare/poincaremap.jl index 9ccd470b..36d53b8e 100644 --- a/src/derived_systems/poincare/poincaremap.jl +++ b/src/derived_systems/poincare/poincaremap.jl @@ -168,7 +168,7 @@ function SciMLBase.step!(pmap::PoincareMap) 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); +function SciMLBase.reinit!(pmap::PoincareMap, u0 = initial_state(pmap); t0 = initial_time(pmap), p = current_parameters(pmap) ) isnothing(u0) && return pmap From eb4b076331d319bbc11274f814ff78bbf5a7815f Mon Sep 17 00:00:00 2001 From: Datseris Date: Sun, 22 Sep 2024 14:35:30 +0100 Subject: [PATCH 09/12] final typo fix --- src/derived_systems/poincare/poincaremap.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/derived_systems/poincare/poincaremap.jl b/src/derived_systems/poincare/poincaremap.jl index 36d53b8e..12e95f3c 100644 --- a/src/derived_systems/poincare/poincaremap.jl +++ b/src/derived_systems/poincare/poincaremap.jl @@ -125,7 +125,7 @@ function PoincareMap( v, current_time(ds), 0, dummy, diffidxs ) step!(pmap) - pmap.t[] = 0 # first step is 0 + pmap.t = 0 # first step is 0 return pmap end From 5e91611317faf339c6fb7b725f80d76aebb1b41d Mon Sep 17 00:00:00 2001 From: Datseris Date: Sun, 22 Sep 2024 14:57:47 +0100 Subject: [PATCH 10/12] revert type annotation --- src/derived_systems/poincare/poincaremap.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/derived_systems/poincare/poincaremap.jl b/src/derived_systems/poincare/poincaremap.jl index 12e95f3c..68c57898 100644 --- a/src/derived_systems/poincare/poincaremap.jl +++ b/src/derived_systems/poincare/poincaremap.jl @@ -168,7 +168,7 @@ function SciMLBase.step!(pmap::PoincareMap) end SciMLBase.step!(pmap::PoincareMap, n::Int, s = true) = (for _ ∈ 1:n; step!(pmap); end; pmap) -function SciMLBase.reinit!(pmap::PoincareMap, u0 = initial_state(pmap); +function SciMLBase.reinit!(pmap::PoincareMap, u0::Union{AbstractArray, Nothing} = initial_state(pmap); t0 = initial_time(pmap), p = current_parameters(pmap) ) isnothing(u0) && return pmap From 8f419d68e7f73b1aaaaf3b95fbb81e0c8b892e18 Mon Sep 17 00:00:00 2001 From: Datseris Date: Sun, 22 Sep 2024 14:58:21 +0100 Subject: [PATCH 11/12] correct reinit --- src/derived_systems/poincare/poincaremap.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/derived_systems/poincare/poincaremap.jl b/src/derived_systems/poincare/poincaremap.jl index 68c57898..9d811615 100644 --- a/src/derived_systems/poincare/poincaremap.jl +++ b/src/derived_systems/poincare/poincaremap.jl @@ -168,10 +168,9 @@ function SciMLBase.step!(pmap::PoincareMap) end SciMLBase.step!(pmap::PoincareMap, n::Int, s = true) = (for _ ∈ 1:n; step!(pmap); end; pmap) -function SciMLBase.reinit!(pmap::PoincareMap, u0::Union{AbstractArray, Nothing} = initial_state(pmap); +function SciMLBase.reinit!(pmap::PoincareMap, u0::AbstractArray = initial_state(pmap); t0 = initial_time(pmap), p = current_parameters(pmap) ) - isnothing(u0) && return pmap if length(u0) == dimension(pmap) u = u0 elseif length(u0) == dimension(pmap) - 1 From 68551f1794b2b8d1cd73833724a9020630baa801 Mon Sep 17 00:00:00 2001 From: Datseris Date: Sun, 22 Sep 2024 15:11:17 +0100 Subject: [PATCH 12/12] changelog --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 48d01e96..6f336e3c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,7 +1,7 @@ # v3.10.0 - OrdinaryDiffEq.jl dependency reduced to OrdinaryDiffEqTsit5.jl. -- Updated to StateSpaceSets.jl: now the keyword `container` can be given to +- Updated to StateSpaceSets.jl v2: now the keyword `container` can be given to the `trajectory` function. # v3.9.0