finalize v3 (#164)

* first draft of arbitrary steppable

* different multiline printing

* print ODE solve stuff

* increase default poincare tolerance

* pmap improvements

* more duffing tests

* test refactor

* require statespacesets v1

* rename for SSSets v1

* finish source of arbitrary steppable

* finish arbitrary steppable as well!

* better docstring

* overall safety increase

* use svector/smatrix instead of sarray

Author: Datseris

CHANGELOG.md

@@ -96,7 +96,7 @@ In addition, be sure that you have version `SimpleDiffEq 0.3.0` or greater (whic
 
 # v1.1
 
-* All functionality related to `neighborhood`, `reconstruct` and `Dataset` has moved to a new package: [`DelayEmbeddings`](https://github.com/JuliaDynamics/DelayEmbeddings.jl). It is reexported by `DynamicalSystemsBase`, which is now (as the name suggests) the basis package for defining `DynamicalSystem`.
+* All functionality related to `neighborhood`, `reconstruct` and `StateSpaceSet` has moved to a new package: [`DelayEmbeddings`](https://github.com/JuliaDynamics/DelayEmbeddings.jl). It is reexported by `DynamicalSystemsBase`, which is now (as the name suggests) the basis package for defining `DynamicalSystem`.
 
 # v1.0
 
@@ -115,7 +115,7 @@ In addition, be sure that you have version `SimpleDiffEq 0.3.0` or greater (whic
 ## Breaking
 * Dropped support of Julia versions < 0.7
 * `Reconstruction` does not exist anymore. Instead, `reconstruct` is used in it's
-  place. The function now also always returns a `Dataset`.
+  place. The function now also always returns a `StateSpaceSet`.
 * In the `reconstruct` function, `D` now stands for the number of temporal neighbors
   which is **one less** than the dimensionality of the reconstructed space.
     *  This change was done because it is more intuitive and more scalable when considering multiple timeseries or spatio temporal timeseries.
@@ -157,7 +157,7 @@ very happy to say that these changes reduced *tremendously* the source code!
 * Multi-time Reconstruction
 * Multi-dimensional, multi-time Reconstruction
 * Methods for estimating Reconstruction parameters are now here.
-* Reconstruction from Dataset
+* Reconstruction from StateSpaceSet
 ## Breaking
 * Corrected name `helies -> heiles`
 * Reconstructions now have field `delay` which is the delay time/times and are
@@ -193,7 +193,7 @@ very happy to say that these changes reduced *tremendously* the source code!
   parameters are passed directly into the function!!!
 
 # Non-breaking
-* Improved the algorithm that converts a Dataset to a Matrix. It is now not only
+* Improved the algorithm that converts a StateSpaceSet to a Matrix. It is now not only
   faster, but also more clear!
 
 # v0.4.1
@@ -218,15 +218,15 @@ very happy to say that these changes reduced *tremendously* the source code!
 ## Non-breaking
 * Bugfix of `eltype` of `Reconstruction`.
 * Added `circlemap` to `Systems`.
-* Bugfix on `Dataset` that incorrect methods were being called due to missing `<:`.
-* Method `Base.getindex(d::AbstractDataset{D,T}, ::Colon, j<:AbstractVector{Int}) ` now exists!
+* Bugfix on `StateSpaceSet` that incorrect methods were being called due to missing `<:`.
+* Method `Base.getindex(d::AbstractStateSpaceSet{D,T}, ::Colon, j<:AbstractVector{Int}) ` now exists!
 * Added function `evolve!`.
 * Clearly state that we do not support matrices in the equations of motion.
 * Bugfixes regarding ODEProblem.
 
 # v0.3.2
 ## Non-breaking
-* Orders of magnitude speed-up in conversions between `Matrix` and `Dataset`,
+* Orders of magnitude speed-up in conversions between `Matrix` and `StateSpaceSet`,
   because now methods use `transpose` internally and only `reinterpret`.
 
 # v0.3.1

Project.toml

@@ -20,5 +20,5 @@ OrdinaryDiffEq = "6"
 Reexport = "1"
 Roots = "1, 2"
 SciMLBase = "1.19.5"
-StateSpaceSets = "0.1"
+StateSpaceSets = "1"
 julia = "1.8"

docs/src/index.md

@@ -34,7 +34,7 @@ set_parameters!
 ```@docs
 step!(::DynamicalSystem, args...; kwargs...)
 trajectory
-Dataset
+StateSpaceSet
 ```
 
 ## `DeterministicIteratedMap`

src/core/dynamicalsystem_interface.jl

@@ -43,7 +43,7 @@ Concrete subtypes typically also contain more information than the above 3 items
 
 In this scope dynamical systems have a known dynamic rule `f` defined as a
 standard Julia function. _Observed_ or _measured_ data from a dynamical system
-are represented using `AbstractDataset` and are finite.
+are represented using `AbstractStateSpaceSet` and are finite.
 Such data are obtained from the [`trajectory`](@ref) function or
 from an experimental measurement of a dynamical system with an unknown dynamic rule.
 
@@ -294,15 +294,15 @@ end
 # API - step and reset
 ###########################################################################################
 """
-    step!(ds::DiscreteTimeDynamicalSystem [, dt::Integer]) → ds
+    step!(ds::DiscreteTimeDynamicalSystem [, n::Integer]) → ds
 
-Evolve the discrete time dynamical system for 1 or `dt` steps.
+Evolve the discrete time dynamical system for 1 or `n` steps.
 
     step!(ds::ContinuousTimeDynamicalSystem, [, dt::Real [, stop_at_tdt]]) → ds
 
 Evolve the continuous time dynamical system for one integration step.
 
-Alternative, if a `dt` is given, then progress the integration until
+Alternatively, if a `dt` is given, then progress the integration until
 there is a temporal difference `≥ dt` (so, step _at least_ for `dt` time).
 
 When `true` is passed to the optional third argument,

src/core/pretty_printing.jl

@@ -6,10 +6,14 @@
 Base.summary(ds::DynamicalSystem) =
 "$(dimension(ds))-dimensional $(nameof(typeof(ds)))"
 
+function Base.show(io::IO, ds::DynamicalSystem)
+    print(io, summary(ds))
+end
+
 # Extend this function to return a vector of `Pair`s of `"description" => value`
 additional_details(::DynamicalSystem) = []
 
-function Base.show(io::IO, ds::DynamicalSystem)
+function Base.show(io::IO, ::MIME"text/plain", ds::DynamicalSystem)
     descriptors = [
         "deterministic" => isdeterministic(ds),
         "discrete time" => isdiscretetime(ds),

src/core/trajectory.jl

@@ -8,7 +8,7 @@ after evolving it for total time `T`.
 `u0` is the state given given to [`reinit!`](@ref) prior to time evolution
 and defaults to [`initial_state(ds)`](@ref).
 
-See [`Dataset`](@ref) for info on how to use `X`.
+See [`StateSpaceSet`](@ref) for info on how to use `X`.
 The returned time vector is `t = (t0+Ttr):Δt:(t0+Ttr+T)`.
 
 ## Keyword arguments
@@ -44,7 +44,7 @@ function trajectory_discrete(ds, T; Δt::Integer = 1, Ttr::Integer = 0, accessor
         step!(ds, Δt)
         data[i] = SVector{X, ET}(obtain_access(current_state(ds), accessor))
     end
-    return Dataset(data), tvec
+    return StateSpaceSet(data), tvec
 end
 
 function trajectory_continuous(ds, T; Δt = 0.1, Ttr = 0.0, accessor = nothing)
@@ -61,7 +61,7 @@ function trajectory_continuous(ds, T; Δt = 0.1, Ttr = 0.0, accessor = nothing)
         end
         sol[i] = SVector{X, ET}(obtain_access(ds(t), accessor))
     end
-    return Dataset(sol), tvec
+    return StateSpaceSet(sol), tvec
 end
 
 # Util functions for `trajectory`

src/core_systems/arbitrary_steppable.jl

@@ -0,0 +1,96 @@
+export ArbitrarySteppable
+
+"""
+    ArbitrarySteppable <: DiscreteTimeDynamicalSystem
+    ArbitrarySteppable(
+        model, step!, extract_state, extract_parameters, reset_model!;
+        isdeterministic = true, set_state = reinit!,
+    )
+
+A dynamical system generated by an arbitrary "model" that can be stepped _in-place_
+with some function `step!(model)` for 1 step.
+The state of the model is extracted by the `extract_state(model) -> u` function
+The parameters of the model are extracted by the `extract_parameters(model) -> p` function.
+The system may be re-initialized, via [`reinit!`](@ref), with the `reset_model!`
+user-provided function that must have the call signature
+```julia
+reset_model!(model, u, p)
+```
+given a (potentially new) state `u` and parameter container `p`, both of which will
+default to the initial ones in the [`reinit!`](@ref) call.
+
+`ArbitrarySteppable` exists to provide the DynamicalSystems.jl interface to models from
+other packages that could be used within the DynamicalSystems.jl library.
+`ArbitrarySteppable` follows the [`DynamicalSystem`](@ref) interface with
+the following adjustments:
+
+- [`initial_time`](@ref) is always 0, as time counts the steps the model has
+  taken since creation or last [`reinit!`](@ref) call.
+- [`set_state!`](@ref) is the same as [`reinit!`](@ref) by default.
+  If not, the keyword argument `set_state` is a function `set_state(model, u)`
+  that sets the state of the model to `u`.
+- The keyword `isdeterministic` should be set properly, as it decides whether
+  downstream algorithms should error or not.
+"""
+struct ArbitrarySteppable{M, F, R, ES, EP, U, P, S} <: DiscreteTimeDynamicalSystem
+    model::M
+    step::F
+    reinit::R
+    extract_state::ES
+    extract_parameters::EP
+	t::Base.RefValue{Int}
+    u0::U
+    p0::P
+    isdeterministic::Bool
+    set_state::S
+end
+
+function ArbitrarySteppable(model, step, extract_state, extract_parameters, reinit;
+        isdeterministic = true, set_state = (ds, u) -> reinit(ds, u, current_parameters(ds)),
+    )
+    p0 = deepcopy(extract_parameters(model))
+    u0 = deepcopy(extract_state(model))
+    if !(u0 isa AbstractArray{<:Real})
+        @warn """
+        The state isn't an `AbstractArray{<:Real}` and that may cause problems
+        in downstream functions.
+        """
+    end
+    return ArbitrarySteppable(model, step, reinit, extract_state, extract_parameters,
+    Ref(0), u0, p0, isdeterministic, set_state)
+end
+
+# Extend DynamicalSystems.jl interface
+current_state(ds::ArbitrarySteppable) = ds.extract_state(ds.model)
+current_parameters(ds::ArbitrarySteppable) = ds.extract_parameters(ds.model)
+isdeterministic(ds::ArbitrarySteppable) = ds.isdeterministic
+dynamic_rule(ds::ArbitrarySteppable) = ds.step
+current_time(ds::ArbitrarySteppable) = ds.t[]
+initial_time(ds::ArbitrarySteppable) = 0
+SciMLBase.isinplace(ds::ArbitrarySteppable) = true
+
+function SciMLBase.step!(ds::ArbitrarySteppable, n::Int = 1, stop::Bool = true)
+    for _ in 1:n
+        ds.step(ds.model)
+    end
+    ds.t[] = ds.t[] + n
+    return ds
+end
+
+function SciMLBase.reinit!(ds::ArbitrarySteppable, u = initial_state(ds);
+        p = current_parameters(ds), t0 = 0, # t0 is not used but required for downstream.
+    )
+    isnothing(u) && return
+    ds.reinit(ds.model, u, p)
+    ds.t[] = 0
+    return ds
+end
+
+function set_state!(ds::ArbitrarySteppable, u)
+    ds.set_state(ds.model, u)
+    return ds
+end
+
+additional_details(ds::ArbitrarySteppable) = [
+    "model type" => nameof(typeof(ds.model)),
+]

src/core_systems/continuous_time_ode.jl

@@ -103,6 +103,14 @@ function SciMLBase.ODEProblem(ds::CoupledODEs{IIP}, tspan = (initial_time(ds), I
     ODEProblem{IIP}(dynamic_rule(ds), initial_state(ds), tspan, current_parameters(ds))
 end
 
+# Pretty print
+function additional_details(ds::CoupledODEs)
+    solver, remaining = _decompose_into_solver_and_remaining(ds.diffeq)
+    return ["ODE solver" => string(nameof(typeof(solver))),
+        "ODE kwargs" => remaining,
+    ]
+end
+
 ###########################################################################################
 # Extend interface and extend for `DEIntegrator`
 ###########################################################################################
@@ -135,7 +143,13 @@ current_time(integ::DEIntegrator) = integ.t
 initial_time(integ::DEIntegrator) = integ.sol.prob.tspan[1]
 
 function set_state!(integ::DEIntegrator, u)
-    integ.u = recursivecopy(u)
+    if integ.u isa Array{<:Real}
+        integ.u .= u
+    elseif integ.u isa Union{SVector, SMatrix}
+        integ.u = u
+    else
+        integ.u = recursivecopy(u)
+    end
     u_modified!(integ, true)
     return
 end

src/derived_systems/parallel_systems.jl

@@ -182,7 +182,18 @@ initial_states(pdtds::PDTDS) = [initial_state(ds) for ds in pdtds.systems]
 
 # Set stuff
 set_parameter!(pdtds::PDTDS) = for ds in pdtds.systems; set_parameter!(ds, args...); end
-set_state!(pdtds::PDTDS, u, i::Int = 1) = set_state!(pdtds.systems[i], u)
+function set_state!(pdtds::PDTDS, u, i::Int = 1)
+    # We need to set state in all systems, in case this does
+    # some kind of resetting, e.g., the `u_modified!` stuff.
+    for (k, ds) in enumerate(pdtds.systems)
+        if k == i
+            set_state!(ds, u)
+        else
+            set_state!(ds, current_state(ds))
+        end
+    end
+end
+
 function SciMLBase.reinit!(pdtds::PDTDS, states = initial_states(pdtds); kwargs...)
     for (ds, s) in zip(pdtds.systems, states); reinit!(ds, s; kwargs...); end
 end

src/derived_systems/poincare/hyperplane.jl

@@ -75,25 +75,25 @@ end
 # TODO: Nice improvement would be to use cubic interpolation instead of linear,
 # using points i-2, i-1, i, i+1
 """
-    poincaresos(A::AbstractDataset, plane; kwargs...) → P::Dataset
+    poincaresos(A::AbstractStateSpaceSet, plane; kwargs...) → P::StateSpaceSet
 
 Calculate the Poincaré surface of section of the given dataset with the given `plane`
 by performing linear interpolation betweeen points that sandwich the hyperplane.
 
 Argument `plane` and keywords `direction, warning, save_idxs`
 are the same as in the method below.
 """
-function poincaresos(A::AbstractDataset, plane;
+function poincaresos(A::AbstractStateSpaceSet, plane;
         direction = -1, warning = true, save_idxs = dimension(A)
     )
     check_hyperplane_match(plane, size(A, 2))
     i = typeof(save_idxs) <: Int ? save_idxs : SVector{length(save_idxs), Int}(save_idxs...)
     planecrossing = PlaneCrossing(plane, direction > 0)
     data = poincaresos(A, planecrossing, i)
     warning && length(data) == 0 && @warn PSOS_ERROR
-    return Dataset(data)
+    return StateSpaceSet(data)
 end
-function poincaresos(A::Dataset, planecrossing::PlaneCrossing, j)
+function poincaresos(A::StateSpaceSet, planecrossing::PlaneCrossing, j)
     i, L = 1, length(A)
     data = _initialize_output(A[1], j)
     # Check if initial condition is already on the plane