Skip to content

Commit

Permalink
New algorithm for detection of UPOs (#328)
Browse files Browse the repository at this point in the history
* separating lambda matrix functions for use in DL algo

* extend constant multiple range

* fixing exports

* docstrings only on exported functions

* adding include

* davidchacklai algorithm

* tests

* documentation (haven't run yet)

* docs fix and citing fix

* adding DocumenterCitations.jl

---------

Co-authored-by: Datseris <datseris.george@gmail.com>
  • Loading branch information
JonasKoziorek and Datseris committed Mar 26, 2024
1 parent 14adea9 commit 1c424f7
Show file tree
Hide file tree
Showing 12 changed files with 675 additions and 95 deletions.
5 changes: 3 additions & 2 deletions docs/Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,11 @@ ChaosTools = "608a59af-f2a3-5ad4-90b4-758bdf3122a7"
ComplexityMeasures = "ab4b797d-85ee-42ba-b621-05d793b346a2"
DelayEmbeddings = "5732040d-69e3-5649-938a-b6b4f237613f"
Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244"
DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8"
DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
Neighborhood = "645ca80c-8b79-4109-87ea-e1f58159d116"
OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
TimeseriesSurrogates = "c804724b-8c18-5caa-8579-6025a0767c70"
TimeseriesSurrogates = "c804724b-8c18-5caa-8579-6025a0767c70"
14 changes: 7 additions & 7 deletions docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,7 @@ pages = [
"dimreduction.md",
"periodicity.md",
"rareevents.md",
"references.md",
]

import Downloads
Expand All @@ -20,16 +21,15 @@ Downloads.download(
)
include("build_docs_with_style.jl")

# TODO: Port all citations to use this:
# using DocumenterCitations
using DocumenterCitations

# bib = CitationBibliography(
# joinpath(@__DIR__, "refs.bib");
# style=:authoryear
# )
bib = CitationBibliography(
joinpath(@__DIR__, "refs.bib");
style=:authoryear
)

build_docs_with_style(pages, ChaosTools, DynamicalSystemsBase, Neighborhood;
# bib, # TODO: Enable bib
bib,
# TODO: Fix warnings so that instead of:
warnonly = true,
# we can have:
Expand Down
171 changes: 171 additions & 0 deletions docs/refs.bib
Original file line number Diff line number Diff line change
@@ -0,0 +1,171 @@
@string{aamop = "Adv. At. Mol. Opt. Phys."}
@string{aarc = "Autom. Rem. Contr."}
@string{ac = "Anal. Chem."}
@string{acie = "Angew. Chem. Int. Ed."}
@string{aipa = "AIP Advances"}
@string{ajp = "Am. J. Phys."}
@string{algo = "Algorithmica"}
@string{ao = "Appl. Opt."}
@string{ap = "Adv. Phys."}
@string{apb = "Appl. Phys. B"}
@string{apl = "Appl. Phys. Lett."}
@string{apx = "Adv. Phys. X"}
@string{aqt = "Adv. Quantum Tech."}
@string{arcmp = "Annu. Rev. Condens. Matter Phys."}
@string{arpc = "Annu. Rev. Phys. Chem."}
@string{astrocomp = "Astron. Comput."}
@string{atms = "ACM Trans. Math. Softw."}
@string{avsqs = "AVS Quantum Sci."}
@string{bit = "BIT"}
@string{bstj = "Bell System Tech. J."}
@string{cacm = "Commun. ACM"}
@string{ccyb = "Control Cybern."}
@string{cmp = "Commun. Math. Phys."}
@string{computj = "Comput. J."}
@string{contp = "Contemp. Phys."}
@string{cp = "Chem. Phys."}
@string{cpam = "Commun. Pur. Appl. Math."}
@string{cpc = "Comput. Phys. Commun."}
@string{cpl = "Chem. Phys. Lett."}
@string{cse = "Comput. Sci. Eng."}
@string{csj = "CEAS Space J."}
@string{ecyb = "Engrg. Cybernetics"}
@string{ejc = "Eur. J. Control"}
@string{ejp = "Eur. J. Phys."}
@string{electr = "Electronics"}
@string{entr = "Entropy"}
@string{epjb = "Eur. Phys. J. B"}
@string{epjd = "Eur. Phys. J. D"}
@string{epjp = "Eur. Phys. J. Plus"}
@string{epjqt = "EPJ Quantum Technol."}
@string{epl = "Europhys. Lett."}
@string{farad = "Faraday Disc."}
@string{foundphys = "Found. Phys."}
@string{fp = "Fortschr. Phys."}
@string{icta = "IET Control Theory Appl."}
@string{ijqe = "IEEE J. Quantum Electron."}
@string{ijqi = "Int. J. Quantum Inform."}
@string{ijtp = "Int. J. Theor. Phys."}
@string{imajam = "IMA J. Appl. Math."}
@string{ip = "Inverse Problems"}
@string{itac = "IEEE Trans. Automat. Contr."}
@string{itas = "IEEE Trans. on Appl. Superc."}
@string{jap = "J. Appl. Phys."}
@string{jcam = "J. Comput. Appl. Math"}
@string{jcmpp = "J. Comput. Phys."}
@string{jcp = "J. Chem. Phys."}
@string{jcpm = "J. Phys. Condens. Matter"}
@string{jcss = "J. Comput. System Sci."}
@string{jctn = "J. Comput. Theor. Nanos."}
@string{jlum = "J. Lumin."}
@string{jmo = "J. Mod. Opt."}
@string{jmp = "J. Math. Phys."}
@string{jmr = "J. Magnet. Res."}
@string{jmra = "J. Magnet. Res. A"}
@string{job = "J. Optics B"}
@string{jors = "J. Open Res. Softw."}
@string{josab = "J. Opt. Soc. Am. B"}
@string{joss = "J. Open Source Softw."}
@string{jota = "J. Optim. Theor. Appl."}
@string{jpa = "J. Phys. A"}
@string{jpamt = "J. Phys. A: Math. Theor."}
@string{jpb = "J. Phys. B"}
@string{jpc = "J. Phys. Chem."}
@string{jpca = "J. Phys. Chem. A"}
@string{jpcm = "J. Phys.: Condens. Matter"}
@string{jsp = "J .Stat. Phys."}
@string{mc = "Math. Comput."}
@string{mlst = "Mach. Learn.: Sci. Technol."}
@string{nams = "Notices Amer. Math. Soc."}
@string{nat = "Nature"}
@string{natcom = "Nat. Commun."}
@string{natmeth = "Nat. Methods"}
@string{natnano = "Nat. Nano."}
@string{natphot = "Nat. Photon."}
@string{natphys = "Nat. Phys."}
@string{njp = "New J. Phys."}
@string{npjqi = "npj Quantum Inf"}
@string{oc = "Opt. Comm."}
@string{oe = "Opt. Express"}
@string{os = "Opt. Spectr."}
@string{physd = "Physica D"}
@string{physrep = "Phys. Rep."}
@string{pire = "Proc. IRE"}
@string{pl = "Phys. Lett."}
@string{pla = "Phys. Lett. A"}
@string{plms = "Proc. Lond. Math. Soc."}
@string{pnas = "Proc. Natl. Acad. Sci. U.S.A"}
@string{pr = "Phys. Rev."}
@string{pra = "Phys. Rev. A"}
@string{prapl = "Phys. Rev. Applied"}
@string{prb = "Phys. Rev. B"}
@string{prc = "Phys. Rev. C"}
@string{prd = "Phys. Rev. D"}
@string{pre = "Phys. Rev. E"}
@string{prl = "Phys. Rev. Lett."}
@string{prr = "Phys. Rev. Research"}
@string{prsa = "Proc. R. Soc. A"}
@string{prx = "Phys. Rev. X"}
@string{prxq = "PRX Quantum"}
@string{ps = "Phys. Scripta"}
@string{pt = "Phys. Today"}
@string{ptrsa = "Phil. Trans. R. Soc. A"}
@string{qam = "Q. Appl. Math."}
@string{qic = "Quantum Info. Comput."}
@string{qip = "Quantum Inf. Process."}
@string{qso = "Quantum Semiclass. Opt."}
@string{qst = "Quantum Sci. Technol."}
@string{quant = "Quantum"}
@string{rmp = "Rev. Mod. Phys."}
@string{rms = "Russ. Math. Surv."}
@string{rpp = "Rep. Prog. Phys."}
@string{rsi = "Rev. Sci. Instr."}
@string{sb = "Sci. Bull."}
@string{sci = "Science"}
@string{sciam = "Sci. Am."}
@string{scis = "Sci. China Inf. Sci."}
@string{siamjc = "SIAM J. Comput."}
@string{siamjsc = "SIAM J. Sci. Comput."}
@string{siamrev = "SIAM Rev."}
@string{sp = "Sig. Process."}
@string{spp = "SciPost Phys."}
@string{sr = "Sci. Rep."}
@string{sst = "Supercond. Sci. Technol."}
@string{widm = "WIREs Data Mining Knowl Discov."}
@string{zp = "Z. Phys."}

% https://github.com/Humans-of-Julia/BibParser.jl/issues/32
@string{XXXclearparser = ""}




@article{Davidchack1999,
title = {Efficient algorithm for detecting unstable periodic orbits in chaotic systems},
volume = {60},
ISSN = {1095-3787},
url = {http://dx.doi.org/10.1103/PhysRevE.60.6172},
DOI = {10.1103/physreve.60.6172},
number = {5},
journal = {Physical Review E},
publisher = {American Physical Society (APS)},
author = {Davidchack, Ruslan L. and Lai, Ying-Cheng},
year = {1999},
month = nov,
pages = {6172–6175}
}

@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}
}
124 changes: 120 additions & 4 deletions docs/src/periodicity.md
Original file line number Diff line number Diff line change
Expand Up @@ -58,9 +58,13 @@ of low dimensional dynamical systems is affected by the position and the stabili

Finding unstable (or stable) periodic orbits of a discrete mapping analytically
rapidly becomes impossible for higher orders of fixed points.
Fortunately there is a numeric algorithm due to
Schmelcher & Diakonos which allows such a computation. Notice that even though
the algorithm can find stable fixed points, it is mainly aimed at *unstable* ones.
Fortunately there are numerical algorithms that allow their detection.

### Schmelcher & Diakonos

First of the algorithms was proposed by Schmelcher & Diakonos [Schmelcher1997](@cite).
Notice that even though the algorithm can find stable fixed points, it is mainly
aimed at *unstable* ones.

The functions `periodicorbits` and `lambdamatrix` implement the algorithm:
```@docs
Expand All @@ -69,7 +73,7 @@ lambdamatrix
lambdaperms
```

### Standard Map example
#### Standard Map example
For example, let's find the fixed points of the Standard map of order 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 21×21 density of initial conditions.
Expand Down Expand Up @@ -152,6 +156,118 @@ Okay, this output is great, and we can tell that it is correct because:
order $n$ come in (possible multiples of) $2n$-sized pairs (see e.g. order 5).
This is a direct consequence of the Poincaré–Birkhoff theorem.

### Davidchack & Lai

An extension of the previous algorithm was proposed by Davidchack & Lai [Davidchack1999](@cite).
It works similarly, but it uses smarter seeding and an improved transformation rule.

The functions `davidchacklai` implements the algorithm:
```@docs
davidchacklai
```

#### Logistic Map example

The idea of periodic orbits can be illustrated easily on 1D maps. Finding all periodic orbits of period
$n$ is equivalent to finding all points $x$ such that $f^{n}(x)=x$, where $f^{n}$ is $n$-th composition of $f$. Hence, solving $f^{n}(x)-x=0$ yields such points. However, this is impossible analytically.
Let's see how `davidchacklai` deals with it:

First let's start with finding first $9$ periodic orbits of the logistic map for parameter $3.72$.

```@example MAIN
using ChaosTools
using CairoMakie
logistic_rule(x, p, n) = @inbounds SVector(p[1]*x[1]*(1 - x[1]))
ds = DeterministicIteratedMap(logistic_rule, SVector(0.4), [3.72])
seeds = [SVector(i) for i in LinRange(0.0, 1.0, 10)]
output = davidchacklai(ds, 9, seeds, 6; abstol=1e-6, disttol=1e-12)
```

Now to check the results, let's plot the periodic orbits of period $6$, $7$, $8$ and $9$.

```@example MAIN
function ydata(ds, order, xdata)
ydata = typeof(current_state(ds)[1])[]
for x in xdata
reinit!(ds, x)
step!(ds, order)
push!(ydata, current_state(ds)[1])
end
return ydata
end
fig = Figure()
x = LinRange(0.0, 1.0, 1000)
for (order, position) in zip([6,7,8,9], [(1,1), (1,2), (2,1), (2,2)])
fpsx = output[order]
y = ydata(ds, order, [SVector(x0) for x0 in x])
fpsy = ydata(ds, order, fpsx)
axis = Axis(fig[position...])
axis.title = "Order $order"
lines!(axis, x, x, color=:black, linewidth=0.5)
lines!(axis, x, y, color = :blue, linewidth=0.7)
scatter!(axis, [i[1] for i in fpsx], fpsy, color = :red, markersize=5)
end
fig
```
Points $x$ which fulfill $f^{n}(x)=x$ can be interpreted as an intersection of the function $f^{n}(x)$ and the identity $x$. Our result is correct because all the points of the intersection between the identity and the logistic map were found.

#### Henon Map example

Let's try to use `davidchacklai` in higher dimension. We will try to detect
all periodic points of Henon map of period `1` to `14`.

```@example MAIN
using ChaosTools, CairoMakie
function henon(u0=zeros(2); a = 1.4, b = 0.3)
return DeterministicIteratedMap(henon_rule, u0, [a,b])
end
henon_rule(x, p, n) = SVector{2}(1.0 - p[1]*x[1]^2 + x[2], p[2]*x[1])
ds = henon()
xs = LinRange(-3.0, 3.0, 10)
ys = LinRange(-10.0, 10.0, 10)
seeds = [SVector{2}(x,y) for x in xs for y in ys]
n = 14
m = 6
output = davidchacklai(ds, n, seeds, m; abstol=1e-7, disttol=1e-10)
fig = Figure()
ax = Axis(fig[1,1])
for result in output
scatter!(ax, [x[1] for x in result], [x[2] for x in result], markersize=8, color=:blue)
end
fig
```

The theory of periodic orbits states that UPOs form sort of a skeleton of the chaotic attractor. Our results supports this claim since it closely resembles the Henon attractor.

Note that in this case parameter `m` has to be set to at least `6`. Otherwise, the algorithm
fails to detect orbits of higher periods correctly.

To check if the detected points are indeed periodic, we can do the following test:

```@example MAIN
orbit14 = output[end]
c = 0
for x in orbit14
set_state!(ds, x)
step!(ds, 14)
xn = current_state(ds)
if ChaosTools.norm(xn - x) > 1e-10
c += 1
end
end
"$c non-periodic points found in the 14th order orbit."
```

The same test can be applied to orbits of lower periods.




## Estimating the Period

The function [`estimate_period`](@ref) offers ways for estimating the period (either exact for periodic timeseries, or approximate for near-periodic ones) of a given timeseries.
Expand Down
4 changes: 4 additions & 0 deletions docs/src/references.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,4 @@
# References

```@bibliography
```
2 changes: 2 additions & 0 deletions src/ChaosTools.jl
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +19,9 @@ include("stability/fixedpoints.jl")
include("periodicity/period.jl")
include("periodicity/yin.jl")
include("periodicity/custombintree.jl")
include("periodicity/lambdamatrix.jl")
include("periodicity/periodicorbits.jl")
include("periodicity/davidchacklai.jl")

include("rareevents/mean_return_times/mrt_api.jl")

Expand Down
6 changes: 5 additions & 1 deletion src/periodicity/custombintree.jl
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
using DataStructures: RBTree
using DataStructures: RBTree, search_node
import Base: <, ==

# This structure is used to overload the behavior of < and == for the use in a binary tree.
Expand Down Expand Up @@ -50,4 +50,8 @@ end

function storefp!(collection, fp, abstol)
push!(collection, VectorWithEpsRadius{typeof(fp)}(fp, abstol))
end

function previously_detected(tree, fp, abstol)
return !isnothing(search_node(tree, VectorWithEpsRadius{typeof(fp)}(fp, abstol)).data)
end
Loading

0 comments on commit 1c424f7

Please sign in to comment.