CompatHelper: bump compat for StateSpaceSets to 2, (keep existing compat) (#217)

* CompatHelper: bump compat for StateSpaceSets to 2, (keep existing compat)

* allow `container` keyword

* changelog

* fix incorrect usage of size

* use exclusively only v2 of SSEts

* use dimension for poincaresos

* simplify field t of PoincareMap

* and fix incorrect AA type annotation

* final typo fix

* revert type annotation

* correct reinit

* changelog

---------

Co-authored-by: CompatHelper Julia <[email protected]>
Co-authored-by: Datseris <[email protected]>

Author: 3 people

CHANGELOG.md

@@ -1,6 +1,12 @@
+# v3.10.0
+
+- OrdinaryDiffEq.jl dependency reduced to OrdinaryDiffEqTsit5.jl.
+- Updated to StateSpaceSets.jl v2: 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

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"
@@ -21,7 +21,7 @@ OrdinaryDiffEqTsit5 = "1.1"
 Reexport = "1"
 Roots = "1, 2"
 SciMLBase = "1.19.5, 2"
-StateSpaceSets = "1"
+StateSpaceSets = "2"
 Statistics = "1"
 SymbolicIndexingInterface = "0.3.4"
 julia = "1.9"

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

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)

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,10 +122,10 @@ 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
+    pmap.t = 0 # first step is 0
     return pmap
 end
 
@@ -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
@@ -169,7 +171,6 @@ SciMLBase.step!(pmap::PoincareMap, n::Int, s = true) = (for _ ∈ 1:n; step!(pma
 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
@@ -179,7 +180,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
 

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