update to dynamicalsystems v3 (#130)

* use Dt and not diffeq

* rename pinteg to pds

* update trajectory animator

* repalce integ with ds

* fix orbit diagram

* fix psos scan

* fix cobweb

* give up on Poincare map for now

* `bump version

Author: Datseris

Project.toml

@@ -1,7 +1,7 @@
 name = "InteractiveDynamics"
 uuid = "ec714cd0-5f51-11eb-0b6e-452e7367ff84"
 repo = "https://github.com/JuliaDynamics/InteractiveDynamics.jl.git"
-version = "0.21.12"
+version = "0.22.0"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"

examples/dynamicalsystems/brainscan.jl

@@ -3,12 +3,14 @@ using OrdinaryDiffEq
 
 ds = Systems.henonheiles()
 diffeq = (alg = Vern9(), abstol = 1e-9, reltol = 1e-9)
+ds = CoupledODEs(ds, diffeq)
+
 u0s = [
     [0.0, -0.25, 0.42081, 0.0],
     [0.0, 0.1, 0.5, 0.0],
-    [0.0, -0.31596, 0.354461, 0.0591255]
+    [0.0, -0.31596, 0.354461, 0.0591255],
 ]
-trs = [trajectory(ds, 10000, u0; diffeq)[:, SVector(1,2,3)] for u0 ∈ u0s]
+trs = [trajectory(ds, 10000, u0)[1][:, SVector(1,2,3)] for u0 ∈ u0s]
 j = 2 # the dimension of the plane
 
 brainscan_poincaresos(trs, j; linekw = (transparency = true,))

examples/dynamicalsystems/cobweb.jl

@@ -1,9 +1,10 @@
-using InteractiveDynamics, GLMakie, DynamicalSystems
+using DynamicalSystems
+using InteractiveDynamics, GLMakie
 
 # the second range is a convenience for intermittency example of logistic
 rrange = 1:0.001:4.0
 # rrange = (rc = 1 + sqrt(8); [rc, rc - 1e-5, rc - 1e-3])
 
 lo = Systems.logistic(0.4; r=rrange[1])
 
-interactive_cobweb(lo, rrange, 5)
+interactive_cobweb(lo, rrange, 5);

examples/dynamicalsystems/orbitdiagram.jl

@@ -1,14 +1,21 @@
-using InteractiveDynamics, GLMakie
 using DynamicalSystems
+using InteractiveDynamics, GLMakie
 
 i = 1
-p_index = 1
 
-systems = [(Systems.logistic(), 2.5, 4.0, "r", "logistic map"),
-           (Systems.henon(), 0.8, 1.4, "a", "Hénon map"),
-           (Systems.standardmap(), 0.0, 1.2, "k", "standard map")]
+ds = Systems.lorenz()
+pmap = PoincareMap(ds, (2, 0.0); rootkw = (xrtol = 1e-8, atol = 1e-9))
+
+systems = [(Systems.logistic(), 2.5, 4.0, "r", "logistic map", 1),
+           (Systems.henon(), 0.8, 1.4, "a", "Hénon map", 1),
+           (Systems.standardmap(), 0.0, 1.2, "k", "standard map", 1),
+           # Poincare map
+           (pmap, 100.0, 200.0, "ρ", "Lorenz Poincare map", 2),
+           # Stroboscopic map of duffing
+
+]
 
-ds, p_min, p_max, parname, t = systems[2]
+ds, p_min, p_max, parname, t, p_index = systems[4]
 t = "orbit diagram for the "*t
 
 fig, oddata = interactive_orbitdiagram(ds, p_index, p_min, p_max, i;

examples/dynamicalsystems/psos.jl

@@ -1,7 +1,8 @@
 using InteractiveDynamics, GLMakie, OrdinaryDiffEq, DynamicalSystems
-diffeq = (alg = Vern9(), abstol = 1e-9, reltol = 1e-9)
 
 hh = Systems.henonheiles()
+diffeq = (alg = Vern9(), abstol = 1e-9, reltol = 1e-9)
+hh = CoupledODEs(hh, diffeq)
 
 potential(x, y) = 0.5(x^2 + y^2) + (x^2*y - (y^3)/3)
 energy(x,y,px,py) = 0.5(px^2 + py^2) + potential(x,y)
@@ -21,7 +22,7 @@ plane = (1, 0.0) # first variable crossing 0
 # Interactive PSOS with random colors
 scene, state = interactive_poincaresos(
     hh, plane, (2, 4), complete;
-    labels = ("q₂" , "p₂"), diffeq
+    labels = ("q₂" , "p₂")
 );
 
 # %% Coloring points using a custom function

examples/dynamicalsystems/trajectory_animator.jl

@@ -1,10 +1,8 @@
-using InteractiveDynamics
 using DynamicalSystems, GLMakie
-using OrdinaryDiffEq
+using InteractiveDynamics
 
 F, G, a, b = 6.886, 1.347, 0.255, 4.0
-ds = Systems.lorenz84(; F, G, a, b)
-diffeq = (alg = Tsit5(), adaptive = false, dt = 0.01)
+ds = PredefinedDynamicalSystems.lorenz84(; F, G, a, b)
 
 u1 = [0.1, 0.1, 0.1] # periodic
 u2 = u1 .+ 1e-3     # fixed point
@@ -13,13 +11,12 @@ u4 = [-1.5, 1.2, 1.3] .+ 21e-9 # chaotic 2
 u0s = [u1, u2, u3, u4]
 
 interactive_evolution(
-    ds, u0s; tail = 1000, diffeq, fade = true,
+    ds, u0s; tail = 1000, fade = true,
     tsidxs = [1,2],
     # tsidxs = nothing, # comment/uncomment this line to remove timeseries
 )
 
 # %% Lorenz63 with parameters and additional plotted elements
-diffeq = (alg = Tsit5(), adaptive = false, dt = 0.01)
 ps = Dict(
     1 => 1:0.1:30,
     2 => 10:0.1:50,
@@ -40,10 +37,9 @@ u3 = [20,10,40.0]
 u0s = [u1, u3]
 
 idxs = [1,2,3]
-diffeq = (alg = Tsit5(), dt = 0.01, adaptive = false)
 
 figure, obs, step, slidervals = interactive_evolution(
-    ds, u0s; ps, idxs, tail = 1000, diffeq, pnames, lims
+    ds, u0s; ps, idxs, tail = 1000, pnames, lims
 )
 
 # Use the `slidervals` observable to plot fixed points
@@ -58,16 +54,14 @@ scatter!(ax, fpobs; markersize = 15, marker = :diamond, color = :black)
 
 # %% Custom animation
 using DynamicalSystems, InteractiveDynamics, GLMakie
-using OrdinaryDiffEq: Tsit5
 using LinearAlgebra: dot, norm
-
 ds = Systems.thomas_cyclical(b = 0.2)
 u0s = ([3, 1, 1.], [1, 3, 1.], [1, 1, 3.])
-diffeq = (alg = Tsit5(), adaptive = false, dt = 0.05)
+Δt = 0.05
 
 fig, obs, step, = interactive_evolution(
-    ds, u0s; tail = 1000, diffeq, add_controls = false, tsidxs = nothing,
-    idxs = [1, 2, 3],
+    ds, u0s; tail = 1000, add_controls = false, tsidxs = nothing,
+    idxs = [1, 2, 3], Δt,
     figure = (resolution = (1200, 600),),
 )
 axss = content(fig[1,1][1,1])

src/chaos/brainscan.jl

@@ -55,14 +55,14 @@ function brainscan_poincaresos(
     for i in 1:length(As)
         A = As[i]
         # plot 3D trajectories
-        lines!(ax, A.data; color = colors[i], transparency = false, linekw...)
+        lines!(ax, vec(A); color = colors[i], transparency = false, linekw...)
 
         # Poincare sos
         psos = lift(y) do y
             DynamicalSystems.poincaresos(A, (j, y); direction, warning = false)
         end
-        psos2d = lift(p -> p[:, otheridxs].data, psos)
-        psos3d = lift(p -> p.data, psos)
+        psos2d = lift(p -> vec(p[:, otheridxs]), psos)
+        psos3d = lift(p -> vec(p), psos)
         Makie.scatter!(axp, psos2d; color = colors[i], scatterkw...)
         Makie.scatter!(ax, psos3d; color = colors[i], markersize = 5, scatterkw...)
     end

src/chaos/cobweb.jl

@@ -36,6 +36,7 @@ function interactive_cobweb(
 
 @assert O ≥ 1
 @assert DynamicalSystems.dimension(ds) == 1
+@assert !DynamicalSystems.isinplace(ds)
 
 figure = Figure(resolution = (1000, 800))
 axts = figure[1, :] = Axis(figure)
@@ -44,7 +45,7 @@ axmap = figure[2, :] = Axis(figure)
 sg = SliderGrid(figure[3, :],
     (label = pname, range = prange),
     (label = "n", range = 1:Tmax, startvalue = 20),
-    (label = "x₀", range = x0s, startvalue = ds.u0),
+    (label = "x₀", range = x0s, startvalue = DynamicalSystems.initial_state(ds)[1]),
 )
 r_observable, L, x0 = [sg.sliders[i].value for i in 1:3]
 
@@ -54,7 +55,14 @@ function seriespoints(x)
     return [Point2f(n[i], x[i]) for i in 1:length(x)]
 end
 
-x = Observable(DynamicalSystems.trajectory(ds, L[], x0[]))
+function get_float_trajectory(ds, L, x0)
+    y = x0 isa Real ? DynamicalSystems.SVector(x0) : x0
+    z = DynamicalSystems.trajectory(ds, L, y)[1]
+    w = reinterpret(Float64, vec(z))
+    return Float32.(w)
+end
+
+x = Observable(get_float_trajectory(ds, L[], x0[]))
 xn = lift(a -> seriespoints(a), x)
 lines!(axts, xn; color = trajcolor, lw = 2.0)
 scatter!(axts, xn; color = trajcolor, markersize = 10)
@@ -64,9 +72,14 @@ ylims!(axts, xmin, xmax)
 # Cobweb diagram
 DynamicalSystems.set_parameter!(ds, pindex, prange[1])
 
-fobs = Any[Observable(ds.f.(xs, Ref(ds.p), 0))]
+fobs = Any[Observable(
+    map(x -> DynamicalSystems.dynamic_rule(ds)(x, DynamicalSystems.current_parameters(ds), 0)[1], xs)
+)]
 for order in 2:O
-    push!(fobs, Observable(ds.f.(fobs[end][], Ref(ds.p), 0)))
+    push!(fobs, Observable(
+        map(x -> DynamicalSystems.dynamic_rule(ds)(x, DynamicalSystems.current_parameters(ds), 0)[1], fobs[end][])
+        # DynamicalSystems.dynamic_rule(ds).(fobs[end][], Ref(DynamicalSystems.current_parameters(ds)), 0)
+    ))
 end
 
 # plot diagonal and fⁿ
@@ -87,20 +100,23 @@ ylims!(axmap, xmin, xmax)
 # On trigger r-slider update all plots:
 on(r_observable) do r
     DynamicalSystems.set_parameter!(ds, pindex, r)
-    x[] = DynamicalSystems.trajectory(ds, L[], x0[])
-    fobs[1][] = ds.f.(xs, Ref(ds.p), 0)
+    x[] = get_float_trajectory(ds, L[], x0[])
+    fobs[1][] = map(x -> DynamicalSystems.dynamic_rule(ds)(x, DynamicalSystems.current_parameters(ds), 0)[1], xs)
     for order in 2:O
-        fobs[order][] = ds.f.(fobs[order-1][], Ref(ds.p), 0)
+        fobs[order][] = map(
+            x -> DynamicalSystems.dynamic_rule(ds)(x, DynamicalSystems.current_parameters(ds), 0)[1], fobs[order-1][]
+        )
     end
 end
 
 on(L) do l
-    x[] = DynamicalSystems.trajectory(ds, l, x0[])
+    res = get_float_trajectory(ds, l, x0[])
+    x[] = res
     xlims!(axts, 0, l) # this is better than autolimits
 end
 
 on(x0) do u
-    x[] = DynamicalSystems.trajectory(ds, L[], u)
+    x[] = get_float_trajectory(ds, L[], u)
 end
 
 # Finally add buttons to hide/show elements of the plot