Add examples to the docs (#166)

* add two examples for construction

* add complexe xample

* better badges

* pate docs

Author: Datseris

.gitignore

@@ -4,3 +4,6 @@
 docs/build
 deps/build.log
 Manifest.toml
+*style.jl
+*.scss
+*.css

CHANGELOG.md

@@ -1,6 +1,11 @@
 # v3.0
 
-Complete rewrite of the package. The DynamicalSystems.jl v3 changelog summarizes the highlights. Here we will try to list all changes, but it will be difficult to provide a changelog given that everything was changed.
+Complete rewrite of the package. The DynamicalSystems.jl v3 changelog summarizes the highlights. Here we will list all changes to _this specific package_.
+
+## Majorly Breaking
+- The `Discrete/ContinuousDynamicalSystem` constructors no longer accept a Jacobian. Use the dedicated `TangentDynamicalSystem` for something that represents the tangent space and can be given to downstream functions such as `lyapunovspectrum`. As a result, none of the predefined systems come with a hand coded Jacobian. The function is still available for manual use nevertheless.
+- The keyword `diffeq` does not exist anymore and is not given to any downstream functions such as `lyapunovspectrum`. The only struct that cares about DifferentialEquations.jl arguments is `CoupledODEs` so it is the only one that accepts `diffeq` keyword.
+- `trajectory` now returns the actual trajectory _and_ the time vector: `X, t = trajectory(ds, ...)`.
 
 ## Enhancements
 - `DynamicalSystem` now defines a proper, well-thought-out interface that is implemented from all its concrete subtypes. See its docstring for details.
@@ -16,10 +21,6 @@ Complete rewrite of the package. The DynamicalSystems.jl v3 changelog summarizes
 - `stroboscopicmap -> StroboscopicMap`
 - `poincaremap -> PoincareMap`
 
-## Majorly Breaking
-- The `Discrete/ContinuousDynamicalSystem` constructors no longer accept a Jacobian. Use the dedicated `TangentDynamicalSystem` for something that represents the tangent space and can be given to downstream functions such as `lyapunovspectrum`. As a result, none of the predefined systems come with a hand coded Jacobian. The function is still available for manual use nevertheless.
-- The keyword `diffeq` does not exist anymore and is not given to any downstream functions such as `lyapunovspectrum`. The only struct that cares about DifferentialEquations.jl arguments is `CoupledODEs` so it is the only one that accepts `diffeq` keyword.
-- `trajectory` now returns the actual trajectory _and_ the time vector: `X, t = trajectory(ds, ...)`
 
 # v2.9
 - All code related to the poincare map integrator have been moved here from ChaosTools.jl.

README.md

@@ -1,7 +1,8 @@
 # DynamicalSystemsBase.jl
 
-[![](https://img.shields.io/badge/docs-dev-blue.svg)](https://JuliaDynamics.github.io/DynamicalSystemsBase.jl/dev)
-[![](https://img.shields.io/badge/DOI-10.1007/978-3-030-91032-7-purple)](https://link.springer.com/book/10.1007/978-3-030-91032-7)
+[![](https://img.shields.io/badge/docs-dev-lightblue.svg)](https://JuliaDynamics.github.io/DynamicalSystemsBase.jl/dev)
+[![](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaDynamics.github.io/DynamicalSystemsBase.jl/stable)
+[![](https://img.shields.io/badge/DOI-10.1007%2F978--3--030--91032--7-purple)](https://link.springer.com/book/10.1007/978-3-030-91032-7)
 [![CI](https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/workflows/CI/badge.svg)](https://github.com/JuliaDynamics/DynamicalSystemsBase.jl/actions?query=workflow%3ACI)
 [![codecov](https://codecov.io/gh/JuliaDynamics/DynamicalSystemsBase.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/JuliaDynamics/DynamicalSystemsBase.jl)
 [![Package Downloads](https://shields.io/endpoint?url=https://pkgs.genieframework.com/api/v1/badge/DynamicalSystemsBase)](https://pkgs.genieframework.com?packages=DynamicalSystemsBase)

docs/Project.toml

@@ -3,4 +3,5 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"
 DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StateSpaceSets = "40b095a5-5852-4c12-98c7-d43bf788e795"

docs/make.jl

@@ -1,50 +1,26 @@
 cd(@__DIR__)
-CI = get(ENV, "CI", nothing) == "true" || get(ENV, "GITHUB_TOKEN", nothing) !== nothing
-using DynamicalSystemsBase # comes from global environment in CI
-using Documenter
-using DocumenterTools: Themes
-ENV["JULIA_DEBUG"] = "Documenter"
-using CairoMakie
 
-# %% JuliaDynamics theme
-# It includes themeing for the HTML build
-# and themeing for the Makie plotting
-
-using DocumenterTools: Themes
-for file in ("juliadynamics-lightdefs.scss", "juliadynamics-darkdefs.scss", "juliadynamics-style.scss")
-    filepath = joinpath(@__DIR__, file)
-    if !isfile(filepath)
-        download("https://raw.githubusercontent.com/JuliaDynamics/doctheme/master/$file", joinpath(@__DIR__, file))
-    end
-end
-# create the themes
-for w in ("light", "dark")
-    header = read(joinpath(@__DIR__, "juliadynamics-style.scss"), String)
-    theme = read(joinpath(@__DIR__, "juliadynamics-$(w)defs.scss"), String)
-    write(joinpath(@__DIR__, "juliadynamics-$(w).scss"), header*"\n"*theme)
-end
-# compile the themes
-Themes.compile(joinpath(@__DIR__, "juliadynamics-light.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-light.css"))
-Themes.compile(joinpath(@__DIR__, "juliadynamics-dark.scss"), joinpath(@__DIR__, "src/assets/themes/documenter-dark.css"))
+import Downloads
+Downloads.download(
+    "https://raw.githubusercontent.com/JuliaDynamics/doctheme/master/apply_style.jl",
+    joinpath(@__DIR__, "apply_style.jl")
+)
+include("apply_style.jl")
 
-# %% Build docs
-ENV["JULIA_DEBUG"] = "Documenter"
+using DynamicalSystemsBase
 
 DYNAMICALSYSTEMSBASE_PAGES = [
     "index.md",
 ]
 using DynamicalSystemsBase.SciMLBase
 
-include("style.jl")
-
 makedocs(
     modules = [DynamicalSystemsBase, SciMLBase, StateSpaceSets],
     format = Documenter.HTML(
         prettyurls = CI,
         assets = [
             asset("https://fonts.googleapis.com/css?family=Montserrat|Source+Code+Pro&display=swap", class=:css),
         ],
-        collapselevel = 3,
     ),
     sitename = "DynamicalSystemsBase.jl",
     authors = "George Datseris",

docs/src/index.md

@@ -59,6 +59,15 @@ current_crossing_time
 poincaresos
 ```
 
+## `TangentDynamicalSystem`
+```@docs
+CoreDynamicalSystem
+TangentDynamicalSystem
+current_deviations
+set_deviations!
+orthonormal
+```
+
 ## `ProjectedDynamicalSystem`
 ```@docs
 ProjectedDynamicalSystem
@@ -71,21 +80,17 @@ initial_states
 current_states
 ```
 
-## `TangentDynamicalSystem`
+## `ArbitrarySteppable`
 ```@docs
-CoreDynamicalSystem
-TangentDynamicalSystem
-current_deviations
-set_deviations!
-orthonormal
+ArbitrarySteppable
 ```
 
 ## Parallelization
 
-Since `DynamicalSystem`s are mutable, one needs to copy them before parallelizing, to avoid having to deal with complicated race conditions etc. The simplest way is with `deepcopy`. Here is an example block that shows how to parallelize calling some expensive function (e.g., calculating the Lyapunov exponent) over a parameter range:
+Since `DynamicalSystem`s are mutable, one needs to copy them before parallelizing, to avoid having to deal with complicated race conditions etc. The simplest way is with `deepcopy`. Here is an example block that shows how to parallelize calling some expensive function (e.g., calculating the Lyapunov exponent) over a parameter range using `Threads`:
 
 ```julia
-ds = DynamicalSystem(f, u, p)
+ds = DynamicalSystem(f, u, p) # some concrete implementation
 parameters = 0:0.01:1
 outputs = zeros(length(parameters))
 
@@ -99,3 +104,186 @@ Threads.@threads for i in eachindex(parameters)
     outputs[i] = expensive_function(system, args...)
 end
 ```
+
+## Examples
+
+### Iterated map, out of place
+
+Let's make the [Hénon map](https://en.wikipedia.org/wiki/H%C3%A9non_map) as an example.
+
+```@example MAIN
+using DynamicalSystemsBase
+
+henon_rule(x, p, n) = SVector(1.0 - p[1]*x[1]^2 + x[2], p[2]*x[1])
+u0 = zeros(2)
+p0 = [1.4, 0.3]
+
+henon = DeterministicIteratedMap(henon_rule, u0, p0)
+```
+
+and get a trajectory
+
+```@example MAIN
+X, t = trajectory(henon, 10000; Ttr = 100)
+X
+```
+
+### Coupled ODEs, in place
+
+Let's make the Lorenz system
+[Hénon map](https://en.wikipedia.org/wiki/H%C3%A9non_map) as an example.
+The system is small, and therefore should utilize the out of place syntax, but for the case of example, we will use the in-place syntax.
+We'll also use a high accuracy solver from OrdinaryDiffEq.jl.
+
+```@example MAIN
+using DynamicalSystemsBase
+using OrdinaryDiffEq: Vern9
+
+@inbounds function lorenz_rule!(du, u, p, t)
+    σ = p[1]; ρ = p[2]; β = p[3]
+    du[1] = σ*(u[2]-u[1])
+    du[2] = u[1]*(ρ-u[3]) - u[2]
+    du[3] = u[1]*u[2] - β*u[3]
+    return nothing
+end
+
+u0 = [0, 10.0, 0]
+p0 = [10, 28, 8/3]
+diffeq = (alg = Vern9(), abstol = 1e-9, reltol = 1e-9)
+
+lorenz = CoupledODEs(lorenz_rule!, u0, p0; diffeq)
+```
+
+and get a trajectory
+
+```@example MAIN
+X, t = trajectory(lorenz, 1000; Δt = 0.05, Ttr = 10)
+X
+```
+
+
+### Advanced example
+This is an advanced example of making an in-place implementation of coupled [standard maps](https://en.wikipedia.org/wiki/Standard_map). It will utilize a handcoded Jacobian, a sparse matrix for the Jacobinan, a default initial Jacobian matrix, as well as function-like-objects as the dynamic rule.
+
+Coupled standard maps is a deterministic iterated map that can have arbitrary number of equations of motion, since you can couple `N` standard maps which are 2D maps, like so:
+
+```math
+\theta_{i}' = \theta_i + p_{i}' \\
+p_{i}' = p_i + k_i\sin(\theta_i) - \Gamma \left[\sin(\theta_{i+1} - \theta_{i}) + \sin(\theta_{i-1} - \theta_{i}) \right]
+```
+
+To model this, we will make a dedicated `struct`, which is parameterized on the
+number of coupled maps:
+
+```@example MAIN
+struct CoupledStandardMaps{N}
+    idxs::SVector{N, Int}
+    idxsm1::SVector{N, Int}
+    idxsp1::SVector{N, Int}
+end
+```
+
+(what these fields are will become apparent later)
+
+We initialize the struct with the amount of standard maps we want to couple,
+and we also define appropriate parameters:
+
+```@example MAIN
+M = 5  # couple number
+u0 = 0.001rand(2M) #initial state
+ks = 0.9ones(M) # nonlinearity parameters
+Γ = 1.0 # coupling strength
+p = (ks, Γ) # parameter container
+
+# Create struct:
+SV = SVector{M, Int}
+idxs = SV(1:M...) # indexes of thetas
+idxsm1 = SV(circshift(idxs, +1)...)  #indexes of thetas - 1
+idxsp1 = SV(circshift(idxs, -1)...)  #indexes of thetas + 1
+# So that:
+# x[i] ≡ θᵢ
+# x[[idxsp1[i]]] ≡ θᵢ+₁
+# x[[idxsm1[i]]] ≡ θᵢ-₁
+csm = CoupledStandardMaps{M}(idxs, idxsm1, idxsp1)
+```
+
+We will now use this struct to define a [function-like-object](https://docs.julialang.org/en/v1/manual/methods/#Function-like-objects), a Type that also acts as a function
+
+```@example MAIN
+function (f::CoupledStandardMaps{N})(xnew::AbstractVector, x, p, n) where {N}
+    ks, Γ = p
+    @inbounds for i in f.idxs
+
+        xnew[i+N] = mod2pi(
+            x[i+N] + ks[i]*sin(x[i]) -
+            Γ*(sin(x[f.idxsp1[i]] - x[i]) + sin(x[f.idxsm1[i]] - x[i]))
+        )
+
+        xnew[i] = mod2pi(x[i] + xnew[i+N])
+    end
+    return nothing
+end
+```
+
+We will use *the same* `struct` to create a function for the Jacobian:
+
+```@example MAIN
+function (f::CoupledStandardMaps{M})(
+    J::AbstractMatrix, x, p, n) where {M}
+
+    ks, Γ = p
+    # x[i] ≡ θᵢ
+    # x[[idxsp1[i]]] ≡ θᵢ+₁
+    # x[[idxsm1[i]]] ≡ θᵢ-₁
+    @inbounds for i in f.idxs
+        cosθ = cos(x[i])
+        cosθp= cos(x[f.idxsp1[i]] - x[i])
+        cosθm= cos(x[f.idxsm1[i]] - x[i])
+        J[i+M, i] = ks[i]*cosθ + Γ*(cosθp + cosθm)
+        J[i+M, f.idxsm1[i]] = - Γ*cosθm
+        J[i+M, f.idxsp1[i]] = - Γ*cosθp
+        J[i, i] = 1 + J[i+M, i]
+        J[i, f.idxsm1[i]] = J[i+M, f.idxsm1[i]]
+        J[i, f.idxsp1[i]] = J[i+M, f.idxsp1[i]]
+    end
+    return nothing
+end
+```
+
+This is possible because the system state is a `Vector` while the Jacobian is a `Matrix`, so multiple dispatch can differentiate between the two.
+
+Notice in addition, that the Jacobian function accesses *only half the elements of the matrix*. This is intentional, and takes advantage of the fact that the
+other half is constant. We can leverage this further, by making the Jacobian a sparse matrix. Because the `DynamicalSystem` constructors allow us to give in a pre-initialized Jacobian matrix, we take advantage of that and create:
+```@example MAIN
+using SparseArrays
+J = zeros(eltype(u0), 2M, 2M)
+# Set ∂/∂p entries (they are eye(M,M))
+# And they dont change they are constants
+for i in idxs
+    J[i, i+M] = 1
+    J[i+M, i+M] = 1
+end
+sparseJ = sparse(J)
+
+csm(sparseJ, u0, p, 0) # apply Jacobian to initial state
+sparseJ
+```
+
+Now we are ready to create our dynamical system
+
+```@example MAIN
+ds = DeterministicIteratedMap(csm, u0, p)
+```
+
+Of course, the reason we went through all this trouble was to make a [`TangentDynamicalSystem`](@ref), that can actually use the Jacobian function.
+
+```@example MAIN
+tands = TangentDynamicalSystem(ds; J = csm, J0 = sparseJ, k = M)
+```
+
+```@example MAIN
+step!(tands, 5)
+current_deviations(tands)
+```
+
+(the deviation vectors will increase in magnitude rapidly because the dynamical system is chaotic)

docs/style.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 `AbstractStateSpaceSet` and are finite.
+are represented using `StateSpaceSet` 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.
 

src/core/dynamicalsystem_interface.jl

@@ -17,7 +17,7 @@ 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.
 
-See also [`StroboscopicMap`](@ref), [`poincaresos`](@ref), [`produce_orbitdiagram`](@ref).
+See also [`StroboscopicMap`](@ref), [`poincaresos`](@ref).
 
 ## Keyword arguments