New algorithm for detection of UPOs

docs/Project.toml

@@ -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"

docs/make.jl

@@ -11,6 +11,7 @@ pages = [
     "dimreduction.md",
     "periodicity.md",
     "rareevents.md",
+    "references.md",
 ]
 
 import Downloads
@@ -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:

docs/refs.bib

@@ -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}
+}

docs/src/periodicity.md

@@ -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
@@ -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.
@@ -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.

docs/src/references.md

@@ -0,0 +1,4 @@
+# References
+
+```@bibliography
+```

src/ChaosTools.jl

@@ -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")
 

src/periodicity/custombintree.jl

@@ -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.
@@ -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