Skip to content

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

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
8 changes: 7 additions & 1 deletion CHANGELOG.md
Original file line number Diff line number Diff line change
@@ -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
Expand Down
4 changes: 2 additions & 2 deletions Project.toml
Original file line number Diff line number Diff line change
@@ -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"
Expand All @@ -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"
20 changes: 12 additions & 8 deletions src/core/trajectory.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand All @@ -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).
Expand All @@ -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)
Expand All @@ -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)
Expand All @@ -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
Expand Down
2 changes: 1 addition & 1 deletion src/derived_systems/poincare/hyperplane.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
19 changes: 10 additions & 9 deletions src/derived_systems/poincare/poincaremap.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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.
Expand All @@ -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}
Expand All @@ -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

Expand All @@ -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

"""
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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

Expand Down
6 changes: 4 additions & 2 deletions test/test_system_function.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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
Expand Down
Loading