From 1c424f7bb461b9abef6ea7e24a2f62a23a7cd41e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jon=C3=A1=C5=A1=20Koziorek?= <73384756+JonasKoziorek@users.noreply.github.com> Date: Tue, 26 Mar 2024 12:15:47 +0100 Subject: [PATCH] New algorithm for detection of UPOs (#328) * 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 --- docs/Project.toml | 5 +- docs/make.jl | 14 +- docs/refs.bib | 171 +++++++++++++++++++++++ docs/src/periodicity.md | 124 ++++++++++++++++- docs/src/references.md | 4 + src/ChaosTools.jl | 2 + src/periodicity/custombintree.jl | 6 +- src/periodicity/davidchacklai.jl | 218 ++++++++++++++++++++++++++++++ src/periodicity/lambdamatrix.jl | 85 ++++++++++++ src/periodicity/periodicorbits.jl | 83 +----------- test/periodicity/davidchacklai.jl | 57 ++++++++ test/runtests.jl | 1 + 12 files changed, 675 insertions(+), 95 deletions(-) create mode 100644 docs/refs.bib create mode 100644 docs/src/references.md create mode 100644 src/periodicity/davidchacklai.jl create mode 100644 src/periodicity/lambdamatrix.jl create mode 100644 test/periodicity/davidchacklai.jl diff --git a/docs/Project.toml b/docs/Project.toml index aa643f5a..0a703bb1 100644 --- a/docs/Project.toml +++ b/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" \ No newline at end of file +TimeseriesSurrogates = "c804724b-8c18-5caa-8579-6025a0767c70" diff --git a/docs/make.jl b/docs/make.jl index 0e9e6d13..6e25d510 100644 --- a/docs/make.jl +++ b/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: diff --git a/docs/refs.bib b/docs/refs.bib new file mode 100644 index 00000000..8f2084e9 --- /dev/null +++ b/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} +} \ No newline at end of file diff --git a/docs/src/periodicity.md b/docs/src/periodicity.md index 4f95acd9..bd9416b8 100644 --- a/docs/src/periodicity.md +++ b/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. diff --git a/docs/src/references.md b/docs/src/references.md new file mode 100644 index 00000000..4b476771 --- /dev/null +++ b/docs/src/references.md @@ -0,0 +1,4 @@ +# References + +```@bibliography +``` \ No newline at end of file diff --git a/src/ChaosTools.jl b/src/ChaosTools.jl index e70a8271..a59bb4e1 100644 --- a/src/ChaosTools.jl +++ b/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") diff --git a/src/periodicity/custombintree.jl b/src/periodicity/custombintree.jl index 75aaf819..d47a8d46 100644 --- a/src/periodicity/custombintree.jl +++ b/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 \ No newline at end of file diff --git a/src/periodicity/davidchacklai.jl b/src/periodicity/davidchacklai.jl new file mode 100644 index 00000000..a69bc475 --- /dev/null +++ b/src/periodicity/davidchacklai.jl @@ -0,0 +1,218 @@ +export davidchacklai + +""" + davidchacklai(ds::DeterministicIteratedMap, n::Int, ics, m::Int; kwargs...) -> fps + +Find periodic orbits `fps` of orders `1` to `n` for the map `ds` +using the algorithm propesed by Davidchack & Lai[Davidchack1999](@cite). +`ics` is a collection of initial conditions (container of vectors) to be evolved. +`ics` will be used to detect periodic orbits of orders `1` to `m`. These `m` +periodic orbits will be used to detect periodic orbits of order `m+1` to `n`. +`fps` is a vector with `n` elements. `i`-th element is a periodic orbit of order `i`. + +## Keyword arguments + +* `β = nothing`: If it is nothing, then `β(n) = 10*1.2^n`. Otherwise can be a + function that takes period `n` and return a number. It is a parameter mentioned + in the paper[Davidchack1999](@cite). +* `maxiters = nothing`: If it is nothing, then initial condition will be iterated + `max(100, 4*β(p))` times (where `p` is the order of the periodic orbit) + before claiming it has not converged. If an integer, then it is the maximum + amount of iterations an initial condition will be iterated before claiming + it has not converged. +* `disttol = 1e-10`: Distance tolerance. If `norm(f^{n}(x)-x) < disttol` + where `f^{n}` is the `n`-th iterate of the dynamic rule `f`, then `x` + is an `n`-periodic point. +* `abstol = 1e-8`: A detected periodic point isn't stored if it is in `abstol` + neighborhood of some previously detected point. Distance is measured by + euclidian norm. If you are getting duplicate periodic points, decrease this value. + +## Description + +The algorithm is an extension of Schmelcher & Diakonos[Schmelcher1997](@cite) +implemented in [`periodicorbits`](@ref). + +The algorithm 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} + +[\\beta |g(\\mathbf{x}_{n}| C^{T} - J(\\mathbf{x}_{n})]^{-1} g(\\mathbf{x}_{n}) +``` +where +```math +g(\\mathbf{x}_{n}) = f^{n}(\\mathbf{x}_{n}) - \\mathbf{x}_{n} +``` +and +```math +J(\\mathbf{x}_{n}) = \\frac{\\partial g(\\mathbf{x}_{n})}{\\partial \\mathbf{x}_{n}} +```` + +The main difference between Schmelcher & Diakonos[Schmelcher1997](@cite) and +Davidchack & Lai[Davidchack1999](@cite) is that the latter uses periodic points of +previous period as seeds to detect periodic points of the next period. +Additionally, [`periodicorbits`](@ref) only detects periodic points of a given order, +while `davidchacklai` detects periodic points of all orders up to `n`. + + +## Important note + +For low periods `n` circa less than 6, you should select `m = n` otherwise the algorithm +detect periodic orbits correctly. For higher periods, you can select `m` as 6. +You can use initial grid of points for `ics`. Increase `m` in case the orbits are +not being detected correctly. + +""" +function davidchacklai( + ds::DeterministicIteratedMap, + n::Int, + ics, + m::Int = 1; + kwargs... + ) + if isinplace(ds) + throw(ArgumentError("`davidchacklai` currently supports only out of place systems.")) + end + + if (n < 1) + throw(ArgumentError("`n` must be a positive integer.")) + end + + if !(1 <= m <= n) + throw(ArgumentError("`m` must be in [1, `n`=$(n)]")) + end + + type = typeof(current_state(ds)) + fps = storage(type, n) + detection!(fps, ds, n, ics, m; kwargs...) + return output(fps, type, n) +end + +function detection!(fps, ds, n, ics, m; β=nothing, kwargs...) + betagen=betagenerator(β) + indss, signss = lambdaperms(dimension(ds)) + C_matrices = [cmatrix(inds,signs) for inds in indss, signs in signss] + + initial_detection!(fps, ds, ics, m, betagen, C_matrices; kwargs...) + main_detection!(fps, ds, n, betagen, C_matrices; kwargs...) +end + +function betagenerator(β) + if isnothing(β) + return n-> 10*1.2^(n) + else + return β + end +end + +function initial_detection!(fps, ds, ics, m, betagenerator, C_matrices; kwargs...) + for i in 1:m + detect_orbits(fps[i], ds, i, ics, betagenerator(i), C_matrices; kwargs...) + end +end + +function main_detection!(fps, ds, n, betagenerator, C_matrices; kwargs...) + for period = 2:n + β = betagenerator(period) + previousfps = fps[period-1] + currentfps = fps[period] + nextfps = fps[period+1] + for (container, seed, order) in [ + (currentfps, previousfps, period), + (nextfps, currentfps, period+1), + (currentfps, nextfps, period) + ] + detect_orbits(container, ds, order, seed, β, C_matrices; kwargs...) + end + end +end + +function _detect_orbits!(fps, ds, n, seed, C, β; disttol::Float64=1e-10, abstol::Float64=1e-8, maxiters=nothing) + for x in seed + for _ in 1:(isnothing(maxiters) ? max(100, 4*β) : maxiters) + xn = DL_rule(x, β, C, ds, n) + if converged(ds, xn, n, disttol) + if previously_detected(fps, xn, abstol) == false + completeorbit!(fps, ds, xn, n, disttol, abstol) + end + break + end + x = xn + end + end +end + +function completeorbit!(fps, ds, xn, n, disttol, abstol) + traj = trajectory(ds, n, xn)[1] + for t in traj + if converged(ds, t, n, disttol) + storefp!(fps, t, abstol) + end + end +end + +function converged(ds, xn, n, disttol) + return norm(g(ds, xn, n)) < disttol +end + + +function detect_orbits( + fps::RBTree{VectorWithEpsRadius{T}}, + ds::DeterministicIteratedMap, + n::Integer, + seed::AbstractVector{D}, + β, + C_matrices; + kwargs... + ) where {T, D} + for C in C_matrices + _detect_orbits!(fps, ds, n, seed, C, β; kwargs...) + end +end + +function detect_orbits( + fps::RBTree{VectorWithEpsRadius{T}}, + ds::DeterministicIteratedMap, + n::Integer, + seed::RBTree{VectorWithEpsRadius{T}}, + β, + C_matrices; + kwargs... + ) where {T} + for C in C_matrices + _detect_orbits!(fps, ds, n, tovector(seed), C, β; kwargs...) + end +end + +function DL_rule(x, β, C, ds, n) + Jx = DynamicalSystemsBase.ForwardDiff.jacobian(x0 -> g(ds, x0, n), x) + gx = g(ds, x, n) + xn = x + inv(β*norm(gx)*C' - Jx) * gx + return xn +end + +function g(ds, state, n) + p = current_parameters(ds) + newst = state + # not using step!(ds, n) to allow automatic jacobian + for _ = 1:n + newst = ds.f(newst, p, 1.0) + end + return newst - state +end + +function output(fps, type, n) + output = Vector{Vector{type}}(undef, n) + for i in 1:n # not including periodic orbit n+1 because it may be incomplete + output[i] = tovector(fps[i]) + end + return output +end + +function storage(type, n) + storage = Vector{RBTree{VectorWithEpsRadius{type}}}(undef, n+1) + for i in 1:n+1 + storage[i] = fpcollection(type) + end + return storage +end \ No newline at end of file diff --git a/src/periodicity/lambdamatrix.jl b/src/periodicity/lambdamatrix.jl new file mode 100644 index 00000000..c184153b --- /dev/null +++ b/src/periodicity/lambdamatrix.jl @@ -0,0 +1,85 @@ +export lambdamatrix, lambdaperms + +########################################################################################## +# 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 function +[`periodicorbits`](@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] 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] for details. + +[^Pingel2000]: D. Pingel *et al.*, Phys. Rev. E **62**, pp 2119 (2000) + +[^Diakonos1998]: F. K. Diakonos *et al.*, Phys. Rev. Lett. **81**, pp 4349 (1998) +""" +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.")) + 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) diff --git a/src/periodicity/periodicorbits.jl b/src/periodicity/periodicorbits.jl index 54bdeae9..66dfc4e1 100644 --- a/src/periodicity/periodicorbits.jl +++ b/src/periodicity/periodicorbits.jl @@ -1,7 +1,7 @@ using Combinatorics: permutations, multiset_permutations using Random: randperm -export lambdamatrix, lambdaperms, periodicorbits +export periodicorbits ########################################################################################## # Docstring and core function @@ -124,83 +124,4 @@ function Sk(ds::DeterministicIteratedMap, prevst, o::Int, Λ) step!(ds, o) # TODO: For IIP systems optimizations can be done here to not re-allocate vectors... return prevst, prevst + Λ*(current_state(ds) .- prevst) -end - -########################################################################################## -# 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 function -[`periodicorbits`](@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] 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] for details. - -[^Pingel2000]: D. Pingel *et al.*, Phys. Rev. E **62**, pp 2119 (2000) - -[^Diakonos1998]: F. K. Diakonos *et al.*, Phys. Rev. Lett. **81**, pp 4349 (1998) -""" -function lambdamatrix(λ::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.")) - 0 < λ < 1 || throw(ArgumentError("λ must be in (0,1)")) - 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(λ), (D,D)) - for i in 1:D - a[(i-1)*D + inds[i]] = λ*(sings[i] > 0 ? +1 : -1) - end - return SMatrix{D,D}(a) -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 - - +end \ No newline at end of file diff --git a/test/periodicity/davidchacklai.jl b/test/periodicity/davidchacklai.jl new file mode 100644 index 00000000..290ae908 --- /dev/null +++ b/test/periodicity/davidchacklai.jl @@ -0,0 +1,57 @@ +using ChaosTools +using Test + +# test if high period orbits are indeed periodic +@testset "Henon map" begin + henon_rule(x, p, n) = SVector{2}(1.0 - p[1]*x[1]^2 + x[2], p[2]*x[1]) + ds = DeterministicIteratedMap(henon_rule, zeros(2), [1.4, 0.3]) + xs = range(0, stop = 2π, length = 5); ys = copy(xs) + ics = [SVector(x,y) for x in xs for y in ys] + o = 10 + m = 6 + fps = davidchacklai(ds, o, ics, m; abstol=1e-8, disttol=1e-12) + tol = 1e-12 + for x0 in fps[end] + set_state!(ds, x0) + step!(ds, o) + xn = current_state(ds) + @test isapprox(x0, xn; atol = tol) + end +end + + +# test analytical values and correct length +@testset "Logistic map" begin + r = 1+sqrt(8) + logistic_rule(x, p, n) = @inbounds SVector(p[1]*x[1]*(1 - x[1])) + ds = DeterministicIteratedMap(logistic_rule, SVector(0.5), [r]) + xs = LinRange(0, 1, 5) + ics = [SVector(x) for x in xs] + o = 10 + m = 5 + fps = davidchacklai(ds, o, ics, m; abstol=1e-6, disttol=1e-13) + tol = 1e-13 + + ## known analytically + @test length(fps[1]) == 2 + @test isapprox(SVector(0.0), minimum(fps[1]); atol = tol) + @test isapprox(SVector((r-1)/r), maximum(fps[1]); atol = tol) + + sort!(fps[3]) + fps[3] = [round.(x, digits=7) for x in fps[3]] + @test length(fps[3]) == 5 + @test isapprox(SVector(0.0), fps[1][1]; atol = tol) + @test isapprox(SVector(0.1599288), fps[3][2]; atol = tol) + @test isapprox(SVector(0.5143553), fps[3][3]; atol = tol) + @test isapprox(SVector(0.7387961), fps[3][4]; atol = tol) + @test isapprox(SVector(0.9563178), fps[3][5]; atol = tol) + ## + + ## check if last periodic orbit indeed contains periodic points + for x0 in fps[end] + set_state!(ds, x0) + step!(ds, o) + xn = current_state(ds) + @test isapprox(x0, xn; atol = tol) + end +end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 45365ddc..cd667a06 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -28,6 +28,7 @@ testfile("chaosdetection/expansionentropy.jl") include("stability/fixedpoints.jl") include("periodicity/periodicorbits.jl") +include("periodicity/davidchacklai.jl") include("periodicity/period.jl") # testfile("rareevents/return_time_tests.jl", "Return times")