Porting Schmelcher-Diakonos

docs/refs.bib

@@ -8,4 +8,49 @@ @article{Dednam2014
   year = {2014},
   copyright = {arXiv.org perpetual,  non-exclusive license},
   journal = {}
+}
+
+@article{Schmelcher1997,
+  title = {Detecting Unstable Periodic Orbits of Chaotic Dynamical Systems},
+  volume = {78},
+  ISSN = {1079-7114},
+  url = {http://dx.doi.org/10.1103/PhysRevLett.78.4733},
+  DOI = {10.1103/physrevlett.78.4733},
+  number = {25},
+  journal = {Physical Review Letters},
+  publisher = {American Physical Society (APS)},
+  author = {Schmelcher,  P. and Diakonos,  F. K.},
+  year = {1997},
+  month = jun,
+  pages = {4733–4736}
+}
+
+@article{Pingel2000,
+  title = {Theory and applications of the systematic detection of unstable periodic orbits in dynamical systems},
+  volume = {62},
+  ISSN = {1095-3787},
+  url = {http://dx.doi.org/10.1103/physreve.62.2119},
+  DOI = {10.1103/physreve.62.2119},
+  number = {2},
+  journal = {Physical Review E},
+  publisher = {American Physical Society (APS)},
+  author = {Pingel,  Detlef and Schmelcher,  Peter and Diakonos,  Fotis K. and Biham,  Ofer},
+  year = {2000},
+  month = aug,
+  pages = {2119–2134}
+}
+
+@article{Diakonos1998,
+  title = {Systematic Computation of the Least Unstable Periodic Orbits in Chaotic Attractors},
+  volume = {81},
+  ISSN = {1079-7114},
+  url = {http://dx.doi.org/10.1103/PhysRevLett.81.4349},
+  DOI = {10.1103/physrevlett.81.4349},
+  number = {20},
+  journal = {Physical Review Letters},
+  publisher = {American Physical Society (APS)},
+  author = {Diakonos,  Fotis K. and Schmelcher,  Peter and Biham,  Ofer},
+  year = {1998},
+  month = nov,
+  pages = {4349–4352}
 }

docs/src/algorithms.md

@@ -28,8 +28,105 @@ function plot_result(res, ds; azimuth = 1.3 * pi, elevation = 0.3 * pi)
 end
 
 ig = InitialGuess(SVector(2.0, 5.0, 10.0), 10.2)
-alg = OptimizedShooting(Δt=0.01, n=3)
+OSalg = OptimizedShooting(Δt=0.01, n=3)
 ds = CoupledODEs(roessler_rule, [1.0, -2.0, 0.1], [0.15, 0.2, 3.5])
-res = periodic_orbit(ds, alg, ig)
+res = periodic_orbit(ds, OSalg, ig)
 plot_result(res, ds; azimuth = 1.3pi, elevation=0.1pi)
-```
+```
+## Schmelcher & Diakonos
+
+```@docs
+SchmelcherDiakonos
+lambdamatrix
+lambdaperms
+```
+
+### Standard Map example
+For example, let's find the fixed points of the Standard map of period 2, 3, 4, 5, 6
+and 8. We will use all permutations for the `signs` but only one for the `inds`.
+We will also only use one `λ` value, and a 11×11 density of initial conditions.
+
+First, initialize everything
+```@example MAIN
+using PeriodicOrbits
+
+function standardmap_rule(x, k, n)
+    theta = x[1]; p = x[2]
+    p += k[1]*sin(theta)
+    theta += p
+    return SVector(mod2pi(theta), mod2pi(p))
+end
+
+standardmap = DeterministicIteratedMap(standardmap_rule, rand(2), [1.0])
+xs = range(0, stop = 2π, length = 11); ys = copy(xs)
+ics = InitialGuess[InitialGuess(SVector{2}(x,y), nothing) for x in xs for y in ys]
+
+# All permutations of [±1, ±1]:
+signss = lambdaperms(2)[2] # second entry are the signs
+
+# I know from personal research I only need this `inds`:
+indss = [[1,2]] # <- must be container of vectors!
+
+λs = [0.005] # <- vector of numbers
+
+periods = [2, 3, 4, 5, 6, 8]
+ALLFP = Vector{PeriodicOrbit}[]
+
+standardmap
+```
+Then, do the necessary computations for all periods
+
+```@example MAIN
+for o in periods
+    SDalg = SchmelcherDiakonos(o=o, λs=λs, indss=indss, signss=signss, maxiters=30000)
+    FP = periodic_orbits(standardmap, SDalg, ics)
+    FP = uniquepos(FP; atol=1e-5)
+    push!(ALLFP, FP)
+end
+```
+
+Plot the phase space of the standard map
+```@example MAIN
+using CairoMakie
+iters = 1000
+dataset = trajectory(standardmap, iters)[1]
+for x in xs
+    for y in ys
+        append!(dataset, trajectory(standardmap, iters, [x, y])[1])
+    end
+end
+
+fig = Figure()
+ax = Axis(fig[1,1]; xlabel = L"\theta", ylabel = L"p",
+    limits = ((xs[1],xs[end]), (xs[1],xs[end]))
+)
+scatter!(ax, dataset[:, 1], dataset[:, 2]; markersize = 1, color = "black")
+fig
+```
+
+and finally, plot the fixed points
+```@example MAIN
+markers = [:diamond, :utriangle, :rect, :pentagon, :hexagon, :circle]
+
+for i in eachindex(ALLFP)
+    FP = ALLFP[i]
+    o = periods[i]
+    points = Tuple{Float64, Float64}[]
+    for po in FP
+        append!(points, [Tuple(x) for x in po.points])
+    end
+    println(points)
+    scatter!(ax, points; marker=markers[i], color = Cycled(i),
+        markersize = 30 - 2i, strokecolor = "grey", strokewidth = 1, label = "period $o"
+    )
+end
+axislegend(ax)
+fig
+```
+
+Okay, this output is great, and we can tell that it is correct because:
+
+1. Fixed points of period $n$ are also fixed points of period $2n, 3n, 4n, ...$
+2. Besides fixed points of previous periods, *original* fixed points of
+   period $n$ come in (possible multiples of) $2n$-sized pairs (see e.g. period 5).
+   This is a direct consequence of the Poincaré–Birkhoff theorem.

docs/src/api.md

@@ -9,6 +9,7 @@ periodic_orbits
 ```
 ## Algorithms for Discrete-Time Systems
 
+- [`SchmelcherDiakonos`](@ref)
 
 ## Algorithms for Continuous-Time Systems
 

src/PeriodicOrbits.jl

@@ -17,6 +17,8 @@ include("api.jl")
 include("stability.jl")
 include("minimal_period.jl")
 include("pretty_printing.jl")
+include("lambdamatrix.jl")
+include("algorithms/schmelcher_diakonos.jl")
 include("algorithms/optimized_shooting.jl")
 
 end

src/algorithms/schmelcher_diakonos.jl

@@ -0,0 +1,88 @@
+export SchmelcherDiakonos, periodic_orbits
+
+using LinearAlgebra: norm
+
+
+"""
+    SchmelcherDiakonos(; kwargs...)
+
+Detect periodic orbits of `ds <: DiscreteTimeDynamicalSystem` using algorithm
+proposed by Schmelcher & Diakonos [Schmelcher1997](@cite).
+
+## Keyword arguments
+
+* `o` : period of the periodic orbit
+* `λs` : vector of λ parameters, see [Schmelcher1997](@cite) for details
+* `indss` : vector of vectors of indices for the permutation matrix
+* `signss` : vector of vectors of signs for the permutation matrix
+* `maxiters=10000` : maximum amount of iterations an initial guess will be iterated before 
+  claiming it has not converged
+* `inftol=10.0` : if a state reaches `norm(state) ≥ inftol` it is assumed that it has 
+  escaped to infinity (and is thus abandoned)
+* `disttol=1e-10` : distance tolerance. If the 2-norm of a previous state with the next one 
+  is `≤ disttol` then it has converged to a fixed point
+
+## Description
+
+The algorithm used can detect periodic orbits
+by turning fixed points of the original
+map `ds` to stable ones, through the transformation
+```math
+\\mathbf{x}_{n+1} = \\mathbf{x}_n +
+\\mathbf{\\Lambda}_k\\left(f^{(o)}(\\mathbf{x}_n) - \\mathbf{x}_n\\right)
+```
+The index ``k`` counts the various
+possible ``\\mathbf{\\Lambda}_k``.
+
+## Performance notes
+
+*All* initial guesses are
+evolved for *all* ``\\mathbf{\\Lambda}_k`` which can very quickly lead to
+long computation times.
+"""
+@kwdef struct SchmelcherDiakonos <: PeriodicOrbitFinder
+    o::Int64
+    λs::Vector{Float64}
+    indss::Vector{Vector{Int64}}
+    signss::Vector{Vector{Int64}}
+    maxiters::Int64 = 10000
+    disttol::Float64 = 1e-10
+    inftol::Float64 = 10.0
+end
+
+function periodic_orbits(ds::DiscreteTimeDynamicalSystem, alg::SchmelcherDiakonos, igs::Vector{InitialGuess})
+    type = typeof(current_state(ds))
+    POs = type[]
+    for λ in alg.λs, inds in alg.indss, sings in alg.signss
+        Λ = lambdamatrix(λ, inds, sings)
+        _periodic_orbits!(POs, ds, alg, igs, Λ)
+    end
+    po = PeriodicOrbit[PeriodicOrbit(ds, fp, alg.o) for fp in POs]
+    return po
+end
+
+
+function _periodic_orbits!(POs, ds, alg, igs, Λ)
+    igs = [ig.u0 for ig in igs]
+    for ig in igs
+        reinit!(ds, ig)
+        previg = ig
+        for _ in 1:alg.maxiters
+            previg, ig = Sk(ds, previg, alg.o, Λ)
+            norm(ig) > alg.inftol && break
+
+            if norm(previg - ig) < alg.disttol
+                push!(POs, ig)
+                break
+            end
+            previg = ig
+        end
+    end
+end
+
+function Sk(ds, state, o, Λ)
+    reinit!(ds, state)
+    step!(ds, o)
+    # TODO: For IIP systems optimizations can be done here to not re-allocate vectors...
+    return state, state + Λ*(current_state(ds) .- state)
+end

src/api.jl

@@ -199,13 +199,18 @@ end
 
 
 """
-    uniquepos(pos::Vector{PeriodicOrbit}, atol=1e-6) → Vector{PeriodicOrbit}
+    uniquepos(pos::Vector{PeriodicOrbit}; atol=1e-6) → Vector{PeriodicOrbit}
 
 Return a vector of unique periodic orbits from the vector `pos` of periodic orbits.
 By unique we mean that the distance between any two periodic orbits in the vector is 
 greater than `atol`. To see details about the distance function, see [`podistance`](@ref).
+
+## Keyword arguments
+
+* `atol` : minimal distance between two periodic orbits for them to be considered unique.
+
 """
-function uniquepos(pos::Vector{PeriodicOrbit{D,B,R}}, atol::Real=1e-6) where {D,B,R}
+function uniquepos(pos::Vector{PeriodicOrbit}; atol::Real=1e-6)
     length(pos) == 0 && return pos
     unique_pos = typeof(pos[end])[]
     pos = copy(pos)

src/lambdamatrix.jl

@@ -0,0 +1,84 @@
+export lambdamatrix, lambdaperms
+
+using Combinatorics: permutations, multiset_permutations
+using Random: randperm
+
+##########################################################################################
+# Lambda matrix stuff
+##########################################################################################
+"""
+    lambdamatrix(λ, inds::Vector{Int}, sings) -> Λk
+
+Return the matrix ``\\mathbf{\\Lambda}_k`` used to create a new
+dynamical system with some unstable fixed points turned to stable in the algorithm
+[`SchmelcherDiakonos`](@ref).
+
+## Arguments
+
+1. `λ<:Real` : the multiplier of the ``C_k`` matrix, with `0 < λ < 1`.
+2. `inds::Vector{Int}` :
+   The `i`-th entry of this vector gives the *row* of the nonzero element of the `i`th
+   column of ``C_k``.
+3. `sings::Vector{<:Real}` : The element of the `i`-th column of ``C_k`` is +1
+   if `signs[i] > 0` and -1 otherwise (`sings` can also be `Bool` vector).
+
+Calling `lambdamatrix(λ, D::Int)`
+creates a random ``\\mathbf{\\Lambda}_k`` by randomly generating
+an `inds` and a `signs` from all possible combinations. The *collections*
+of all these combinations can be obtained from the function [`lambdaperms`](@ref).
+
+## Description
+
+Each element
+of `inds` *must be unique* such that the resulting matrix is orthogonal
+and represents the group of special reflections and permutations.
+
+Deciding the appropriate values for `λ, inds, sings` is not trivial. However, in
+ref.[Pingel2000](@cite) there is a lot of information that can help with that decision. Also,
+by appropriately choosing various values for `λ`, one can sort periodic
+orbits from e.g. least unstable to most unstable, see [Diakonos1998](@cite) for details.
+"""
+function lambdamatrix(λ::Real, inds::AbstractVector{<:Integer},
+    sings::AbstractVector{<:Real})
+    0 < λ < 1 || throw(ArgumentError("λ must be in (0,1)"))
+    return cmatrix(λ, inds, sings)
+end
+
+function lambdamatrix(λ::T, D::Integer) where {T<:Real}
+    positions = randperm(D)
+    signs = rand(Bool, D)
+    lambdamatrix(λ, positions, signs)
+end
+
+"""
+    lambdaperms(D) -> indperms, singperms
+
+Return two collections that each contain all possible combinations of indices (total of
+``D!``) and signs (total of ``2^D``) for dimension `D` (see [`lambdamatrix`](@ref)).
+"""
+function lambdaperms(D::Integer)
+    indperms = collect(permutations([1:D;], D))
+    p = trues(D)
+    singperms = [p[:]] #need p[:] because p is mutated afterwards
+    for i = 1:D
+        p[i] = false; append!(singperms, multiset_permutations(p, D))
+    end
+    return indperms, singperms
+end
+
+
+function cmatrix(constant::Real, inds::AbstractVector{<:Integer},
+    sings::AbstractVector{<:Real})
+    # this function seems to be super inefficient
+    D = length(inds)
+    D != length(sings) && throw(ArgumentError("inds and sings must have equal size."))
+    Base.unique(inds)!=inds && throw(ArgumentError("All elements of inds must be unique."))
+    # This has to be improved to not create intermediate arrays!!!
+    a = zeros(typeof(constant), (D,D))
+    for i in 1:D
+        a[(i-1)*D + inds[i]] = constant*(sings[i] > 0 ? +1 : -1)
+    end
+    return SMatrix{D,D}(a)
+end
+
+cmatrix(inds::AbstractVector{<:Integer}, signs::AbstractVector{<:Real}) = cmatrix(1.0, inds, signs)