diff --git a/CHANGELOG.md b/CHANGELOG.md index 817e12f..35dda8e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,3 +1,7 @@ + +# 0.9.0 +* Most functions have moved to ChaosTools.jl + # 0.8.0 * Added fixed resolution uncertainty exponent estimation * improved uncertainty exponent estimation diff --git a/Project.toml b/Project.toml index 108690c..997d175 100644 --- a/Project.toml +++ b/Project.toml @@ -1,27 +1,21 @@ name = "Basins" uuid = "7e5b28ef-eee2-49be-86fd-c343c6921134" authors = ["Alexandre Wagemakers"] -version = "0.8.0" +version = "0.9.0" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" -DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" DynamicalSystems = "61744808-ddfa-5f27-97ff-6e42cc95d634" ImageFiltering = "6a3955dd-da59-5b1f-98d4-e7296123deb5" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -LsqFit = "2fda8390-95c7-5789-9bda-21331edee243" NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce" -RegionTrees = "dee08c22-ab7f-5625-9660-a9af2021b33f" Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665" [compat] Combinatorics = "0.7, 1" -DataStructures = "0.18" -DynamicalSystems = "1.6.0" +DynamicalSystems = "2.0" ImageFiltering = "0.6" -LsqFit = "0.12.0" NearestNeighbors = "0.4" -RegionTrees = "0.3" Roots = "0.7, 0.8, 1.0" julia = "1.5" diff --git a/README.md b/README.md index 39d0264..e926739 100644 --- a/README.md +++ b/README.md @@ -1,264 +1,20 @@ Basins.jl ========= -This Julia package computes basins of attraction of dynamical systems in the phase plane and also -several metrics over the basins. The algorithm computes the basin without prior knowledge of the attractors. +This Julia package is now deprecated in favor of [ChaosTools.jl](https://github.com/JuliaDynamics/ChaosTools.jl). However some experimental stuff remains in the package: Wada detection methods and Saddles computation. -This package depends heavily on the package DynamicalSystems and DifferentialEquations for the definitions of the dynamical systems. However it is possible to define custom integrators. The package provides the following metrics: -1. Basins of attraction -2. Basin entropy -3. Uncertainty exponent -4. Wada detection -5. Basin saddles -6. Basin stability -7. Examples +1. Wada detection +2. Basin saddles +3. Examples -## 1 - Computing the basins of attraction -The technique used to compute the basin of attraction is described in ref. [1]. It consists in tracking the trajectory on the plane and coloring the points of according to the attractor it leads to. This technique is very efficient for 2D basins. -The algorithm gives back a matrix with the attractor numbered from 1 to N. If an attractor exists outside the defined grid of if the trajectory escapes, this initial condition is labelled -1. It may happens for example if there is a fixed point not on the Poincaré map. +## 1 - Detection of the property of Wada -### 1.1 - Stroboscopic Maps - -First define a dynamical system on the plane, for example with a *stroboscopic* map or Poincaré section. For example we can set up an dynamical system with a stroboscopic map defined: - -```jl -using Basins, DynamicalSystems, DifferentialEquations -ω=1.; F=0.2 -ds =Systems.duffing([0.1, 0.25]; ω = ω, f = F, d = 0.15, β = -1) -integ = integrator(ds; alg=Tsit5(), reltol=1e-8, save_everystep=false) -``` - -Now we define the grid of ICs that we want to analyze and launch the procedure: - -```jl -xg = range(-2.2,2.2,length=200) -yg = range(-2.2,2.2,length=200) -bsn=basins_map2D(xg, yg, integ; T=2π/ω) -``` - -The keyword arguments are: -* `T` : the period of the stroboscopic map. -* `idxs` : the indices of the variable to track on the plane. By default the initial conditions of other variables are set to zero. - -The function returns a structure `bsn` with several fields of interests: -* `bsn.basin` is a matrix that contains the information of the basins of attraction. The attractors are numbered from 1 to N and each element -correspond to an initial condition on the grid. -* `bsn.xg` and `bsn.yg` are the grid vectors. -* `bsn.attractors` is a collection of vectors with the location of the attractors found. - -Now we can plot the nice result of the computation: - -```jl -using Plots -plot(xg,yg,bsn.basin', seriestype=:heatmap) - -``` - -![image](https://i.imgur.com/R2veb5tl.png) - -### 1.2 - Poincaré Maps - -Another example with a Poincaré map: -```jl -using Plots -using DynamicalSystems -using Basins - -ds = Systems.rikitake(μ = 0.47, α = 1.0) -integ=integrator(ds) -``` - -Once the integrator has been set, the Poincaré map can defined on a plane: - -```jl -xg=range(-6.,6.,length=200) -yg=range(-6.,6.,length=200) -pmap = poincaremap(ds, (3, 0.), Tmax=1e6; idxs = 1:2, rootkw = (xrtol = 1e-8, atol = 1e-8), reltol=1e-9) - -@time bsn = basin_poincare_map(xg, yg, pmap) - -plot(xg,yg,bsn.basin',seriestype=:heatmap) -``` - -The arguments are: -* `pmap` : A Poincaré map as defined in [ChaosTools.jl](https://github.com/JuliaDynamics/ChaosTools.jl) - - -![image](https://i.imgur.com/hKcOiwTl.png) - - -### 1.3 - Discrete Maps - -The process to compute the basin of a discrete map is very similar: - -```jl -function newton_map(dz,z, p, n) - f(x) = x^p[1]-1 - df(x)= p[1]*x^(p[1]-1) - z1 = z[1] + im*z[2] - dz1 = f(z1)/df(z1) - z1 = z1 - dz1 - dz[1]=real(z1) - dz[2]=imag(z1) - return -end - -# dummy Jacobian function to keep the initializator happy -function newton_map_J(J,z0, p, n) - return -end - -ds = DiscreteDynamicalSystem(newton_map,[0.1, 0.2], [3] , newton_map_J) -integ = integrator(ds) - -xg=range(-1.5,1.5,length=200) -yg=range(-1.5,1.5,length=200) - -bsn=basin_discrete_map(xg, yg, integ) -``` - -![image](https://i.imgur.com/ppHlGPbl.png) - - -### 1.4 - Custom differential equations and low level functions. - -Supose we want to define a custom ODE and compute the basin of attraction on a defined -Poincaré map: - -```jl -using DifferentialEquations -using Basins - -@inline @inbounds function duffing(u, p, t) - d = p[1]; F = p[2]; omega = p[3] - du1 = u[2] - du2 = -d*u[2] + u[1] - u[1]^3 + F*sin(omega*t) - return SVector{2}(du1, du2) -end - -d=0.15; F=0.2; ω = 0.5848 -ds = ContinuousDynamicalSystem(duffing, rand(2), [d, F, ω]) -integ = integrator(ds; alg=Tsit5(), reltol=1e-8, save_everystep=false) -xg = range(-2.2,2.2,length=200) -yg = range(-2.2,2.2,length=200) - -iter_f! = (integ) -> step!(integ, 2π/ω, true) -reinit_f! = (integ,y) -> reinit!(integ, [y...]) -get_u = (integ) -> integ.u[1:2] - -bsn = draw_basin(xg, yg, integ, iter_f!, reinit_f!, get_u) -``` - -The following anonymous functions are important: -* iter_f! : defines a function that iterates the system one step on the map. -* reinit_f! : sets the initial conditions on the map. Remember that only the -initial conditions on the map must be set. -* get_u : it is a custom function to get the state of the integrator only for the variables -defined on the plane - -### 1.6 Basins in Higher Dimensions - -When you cannot define a Stroboscopic map or a well defined Poincaré map you can always try -the general method for higher dimensions. It is slower and may requires some tuning. The algorithm -looks for atractors on a 2D grid. The initial conditions are set on this grid and all others variables -are set to zero by default. - -### Usage - -```jl -ds = Systems.magnetic_pendulum(γ=1, d=0.2, α=0.2, ω=0.8, N=3) -integ = integrator(ds, u0=[0,0,0,0], reltol=1e-9) -xg=range(-4,4,length=150) -yg=range(-4,4,length=150) -@time bsn = basins_general(xg, yg, integ; dt=1., idxs=1:2) -``` - -Keyword parameters are: -* `dt` : this is the time step. It is recomended to use a value above 1. The result may vary a little -depending on this time step. -* `idxs` : Indices of the variables defined on the plane. - - -![image](https://imgur.com/qgBHZ8Ml.png) - -### 1.5 - Notes about the method - -This method identifies the attractors and their basins of attraction on the grid without prior knowledge about the -system. At the end of a successfull computation the function returns a structure BasinInfo with usefull information -on the basin defined by the grid (`xg`,`yg`). There is an important member named `basin` that contains the estimation -of the basins and also of the attractors. For its content see the following section `Structure of the basin`. - -From now on we will refer to the final attractor or an initial condition to its *number*, *odd numbers* are assigned -to basins and *even numbers* are assigned to attractors. The method starts by picking the first available initial -condition not yet numbered. The dynamical system is then iterated until one of the following condition happens: -* The trajectory hits a known attractor already numbered: the initial condition is collored with corresponding odd number. -* The trajectory diverges or hits an attractor outside the defined grid: the initial condition is set to -1 -* The trajectory hits a known basins 10 times in a row: the initial condition belongs to that basin and is numbered accordingly. -* The trajectory hits 60 times in a row an unnumbered cell: it is considered an attractor and is labelled with a even number. - -Regarding performace, this method is at worst as fast as tracking the attractors. In most cases there is a signicative improvement -in speed. - -### Structure of the basin: - -The basin of attraction is organized in the followin way: -* The atractors points are *even numbers* in the matrix. For example, 2 and 4 refer to distinct attractors. -* The basins are collored with *odd numbers*, `2n+1` corresponding the attractor `2n`. -* If the trajectory diverges or converge to an atractor outside the defined grid it is numbered -1 - -## 2 - Computataion of the Basin Entropy - -The [Basin Entropy](https://doi.org/10.1007/978-3-319-68109-2_2) is a measure of the impredictability of the basin of attraction of a dynamical system. An important feature of the basins of attraction is that for a value above log(2) we can say that the basin is fractalized. - -### Usage - -Once the basin of attraction has been computed, the computing the Basin Entropy is easy: - -```jl -using Basins, DynamicalSystems, DifferentialEquations -ω=1.; F = 0.2 -ds =Systems.duffing([0.1, 0.25]; ω = ω, f = F, d = 0.15, β = -1) -integ_df = integrator(ds; alg=Tsit5(), reltol=1e-8, save_everystep=false) -xg = range(-2.2,2.2,length=200); yg = range(-2.2,2.2,length=200) -bsn = basins_map2D(xg, yg, integ_df; T=2*pi/ω) - -Sb,Sbb = basin_entropy(bsn; eps_x=20, eps_y=20) -``` -The arguments of `basin_entropy` are: -* `basin` : The basin computed on a grid. -* `eps_x`, `eps_y` : size of the window that samples the basin to compute the entropy. - - -## 3 - Computation of the uncertainty exponent of a basin of attraction - -The [uncertainty exponent](https://en.wikipedia.org/wiki/Uncertainty_exponent) is conected to the [box-counting dimension](https://en.wikipedia.org/wiki/Box-counting_dimension). For a given resolution of the original basin, a sampling of the basin is done until the the fraction of uncertain boxes converges. The process is repeated for different box sizes and then the exponent is estimated. - - -### Usage - -```jl -using Basins, DynamicalSystems, DifferentialEquations -ω=1.; F = 0.2 -ds =Systems.duffing([0.1, 0.25]; ω = ω, f = F, d = 0.15, β = -1) -integ_df = integrator(ds; alg=Tsit5(), reltol=1e-8, save_everystep=false) -xg = range(-2.2,2.2,length=200); yg = range(-2.2,2.2,length=200) -bsn = basins_map2D(xg, yg, integ_df; T=2*pi/ω) - -bd = box_counting_dim(xg, yg, bsn) - -# uncertainty exponent is the dimension of the plane minus the box-couting dimension -ue = 2-bd -``` - - -## 4 - Detection of the property of Wada - -### 4.1 - Merge Method +### 1.1 - Merge Method The [Wada property](https://en.wikipedia.org/wiki/Lakes_of_Wada) in basins of attraction is an amazing feature of some basins. It is not trivial at all to demonstrate rigurously this property. There are however computational approaches that gives hints about the presence of this property in a basin of attraction. One of the fastest approach is the [Merging Method](https://doi.org/10.1038/s41598-018-28119-0). The algorithm gives the maximum and minimum Haussdorff distances between merged basins. A good rule of thumb to discard the Wada property is to check if the maximum distance is large in comparison to the resolution of the basin, i.e., if the number of pixel is large. @@ -300,7 +56,7 @@ epsilon = xg[2]-xg[1] @show dmin = min_dist/epsilon ``` -### 4.2 - Grid Method +### 1.2 - Grid Method Another method available and much more accurate is the [Grid Method](https://doi.org/10.1038/srep16579). It divides the grid and scrutinize the boundary to test if all the attractors are present in every point of the boundary. It may be very long to get an answer since the number of points to test duplicates at each step. The algorithm returns a vector with the proportion of boxes with 1 to N attractor. For example if the vector W[N] is above 0.95 we have all the initial boxes in the boundary on the grid with N attractors. It is therefore a strong evidence that we have a Wada boundary. @@ -340,7 +96,7 @@ The algorithm returns: * `W` contains a vector with the proportion of boxes in the boundary of `k` attractor. A good criterion to decide if the boundary is Wada is to look at `W[N]` with N the number of attractors. If this number is above 0.95 we can conclude that the boundary is Wada. -## 5 - Computation of the saddle embedded in the boundary [EXPERIMENTAL FEATURE] +## 2 - Computation of the saddle embedded in the boundary [EXPERIMENTAL FEATURE] There is an invariant subset of the boundary which is invariant under the forward iteration of the dynamical system. This set is called the chaotic set, chaotic saddle or simply saddle set. It is possible to compute an approximation arbitrarily close to the saddle with the saddle straddle method. For a detailed description of the method see [1]. This method requires two `generalized basins` such that the algorithm focus on the boundary between these two sets. We divide the basins in two class such that `bas_A ∪ bas_B = [1:N]` and `bas_A ∩ bas_B = ∅` with `N` the number of attractors. @@ -394,30 +150,9 @@ Keyword arguments are: ![image](https://i.imgur.com/pQLDO0Ol.png) -## 6 - Computation of the Basin Stability - -The Basin Stability [6] measures the relative sizes of the basin. Larger basin are considered more stable since a small perturbation or error in the initial conditions is less likely to change the attractor. - -### Usage - -```jl -using Basins, DynamicalSystems, DifferentialEquations -ω=1.; F = 0.2 -ds =Systems.duffing([0.1, 0.25]; ω = ω, f = F, d = 0.15, β = -1) -integ_df = integrator(ds; alg=Tsit5(), reltol=1e-8, save_everystep=false) -xg = range(-2.2,2.2,length=200); yg = range(-2.2,2.2,length=200) -bsn = basins_map2D(xg, yg, integ_df; T=2*pi/ω) - -@show basin_stability(bsn) -``` - -## 7 - More examples - -You can find more examples in `src/examples` - ## References -[1] H. E. Nusse and J. A. Yorke, Dynamics: numerical explorations, Springer, New York, 2012 +[1] H. E. Nusse and J. A. Yorke, Dynamics: numerical explorations, Springer, New York, 1997 [2] A. Daza, A. Wagemakers, B. Georgeot, D. Guéry-Odelin and M. A. F. Sanjuán, Basin entropy: a new tool to analyze uncertainty in dynamical systems, Sci. Rep., 6, 31416 (2016). @@ -426,5 +161,3 @@ You can find more examples in `src/examples` [4] A. Daza, A. Wagemakers and M. A. F. Sanjuán, Ascertaining when a basin is Wada: the merging method, Sci. Rep., 8, 9954 (2018). [5] A. Daza, A. Wagemakers, M. A. F. Sanjuán and J. A. Yorke, Testing for Basins of Wada, Sci. Rep., 5, 16579 (2015). - -[6] P. Menck, J. Heitzig, N. Marwan et al. How basin stability complements the linear-stability paradigm. Nature Phys 9, 89–92 (2013). diff --git a/src/Basins.jl b/src/Basins.jl index 8edd2f3..5c6172b 100644 --- a/src/Basins.jl +++ b/src/Basins.jl @@ -7,20 +7,8 @@ using ImageFiltering using NearestNeighbors using Combinatorics using Roots -using LsqFit -using RegionTrees -using DataStructures -export basin_entropy, detect_wada_grid_method, detect_wada_merge_method -export box_counting_dim -export draw_basin!, basin_map, basin_general_ds -export basin_stability, uncertainty_exponent, static_estimate - -export compute_saddle,compute_basin_precise -# Write your package code here. -include("basins_attraction.jl") -include("basin_entropy.jl") -include("uncertain_dim.jl") +export compute_saddle, detect_wada_grid_method, detect_wada_merge_method include("wada_detection.jl") include("straddle.jl") diff --git a/src/basin_entropy.jl b/src/basin_entropy.jl index dad6cb2..7cb4783 100644 --- a/src/basin_entropy.jl +++ b/src/basin_entropy.jl @@ -1,131 +1,3 @@ -""" - basin_entropy(basin; eps_x=20, eps_y=20) -This algorithm computes the basin entropy of a computed basin of attraction on a regular grid. -The function return the basin entropy and the boundary basin entropy. -[A. Daza, A. Wagemakers, B. Georgeot, D. Guéry-Odelin and M. A. F. Sanjuán, Basin entropy: -a new tool to analyze uncertainty in dynamical systems, Sci. Rep., 6, 31416, (2016).] -## Arguments -* `basin` : the matrix containing the information of the basin. - -## Keyword arguments -* `eps_x`, `eps_y` : define the size of the box that will be used to sample the basin -""" -function basin_entropy(basin; eps_x=20, eps_y=20) - r,c= size(basin) - vals = unique(basin) - S=Int16(length(vals)) - pn=zeros(Float64,1,S) - Sb=0 - Nb=0 - N=0 - for x in 1:eps_x:(r-eps_x+1), y in 1:eps_y:(c-eps_y+1) - I = CartesianIndices((x:x+eps_x-1,y:y+eps_y-1)) - box_values=[basin[k] for k in I] - N=N+1 - for (k,v) in enumerate(vals) - pn[k]=count(x->x==v,box_values)/length(box_values) - end - Nb = Nb + (length(unique(box_values))>1) - Sb = Sb + sum(entropy.(pn)) - end - return Sb/N, Sb/Nb -end - - -function entropy(p::Float64) - if p == 0. - h = 0. - else - h=p*log(1/p) - end - return h -end - -function basin_entropy(bsn::BasinInfo; eps_x=20, eps_y=20) - ind = findall(iseven.(bsn.basin) .== true) - # Not sure it is necessary to deepcopy - basin_test = deepcopy(bsn.basin) - # Remove attractors from the picture - for k in ind; basin_test[k] = basin_test[k]+1; end - basin_entropy(basin_test; eps_x=eps_x, eps_y=eps_y) -end - -""" - basin_stability(basin) -This algorithm computes the basin stability of a computed basin of attraction on a regular grid. -This function returns a vector with the relative volume of the basins. - -[P. Menck, J. Heitzig, N. Marwan et al. How basin stability complements the linear-stability paradigm. Nature Phys 9, 89–92 (2013).] - -## Arguments -* `basin` : the matrix containing the information of the basin. - -""" -function basin_stability(basin) - -N = length(unique(basin)) -v = zeros(1,N) - -for (k,b) in enumerate(unique(basin)) - v[k] = count(basin .== b)/length(basin) -end - -return v - -end - - -function basin_stability(bsn::BasinInfo) - ind = findall(iseven.(bsn.basin) .== true) - # Not sure it is necessary to deepcopy - basin_test = deepcopy(bsn.basin) - # Remove attractors from the picture - for k in ind; basin_test[k] = basin_test[k]+1; end - return basin_stability(basin_test) -end - - - -""" - fractal_test(bsn::BasinInfo) - -## Arguments -* `bsn` : - -## Keyword arguments -* - -""" -function fractal_test(bsn::BasinInfo) - -# disk radius -εₕ=5 - -# Sanity check. -if length(bsn.xg)/εₕ < 50 - @warn "Maybe the size of the grid is not fine enough." -end - -nx = length(bsn.xg) -ny = length(bsn.yg) -xi=bsn.xg[1]; xf=bsn.xg[end]; yi=bsn.yg[1]; yf=bsn.yg[end]; - - -# Now get the values in the boxes. -while Nb < 1e4 - -x = rand()*(xf-xi)+xi -y = rand()*(yf-yi)+yi -ind = find_in_range(x,y,εₕ) - -end - - -end - -function find_in_range(x,y,εₕ) - - -end +# All functions moved to Chaostools. diff --git a/src/basins_attraction.jl b/src/basins_attraction.jl index a63112c..4ca7e48 100644 --- a/src/basins_attraction.jl +++ b/src/basins_attraction.jl @@ -1,605 +1,4 @@ -PoincareMap = ChaosTools.PoincareMap -mutable struct BasinInfo{I,F,D,Q} - basin :: I - xg :: F - yg :: F - iter_f! :: Function - reinit_f! :: Function - get_u :: Function - current_color :: Int64 - next_avail_color :: Int64 - consecutive_match :: Int64 - consecutive_other_basins :: Int64 - prevConsecutives :: Int64 - prev_attr :: Int64 - prev_bas :: Int64 - prev_step :: Int64 - step :: Int64 - attractors :: D #Dict{Int16,Vector{Vector{Float64}}} - Na :: Int64 - visited :: Q - available :: Q -end -function Base.show(io::IO, bsn::BasinInfo) - println(io, "Basin of attraction structure") - println(io, rpad(" size : ", 14), size(bsn.basin)) - println(io, rpad(" Number of attractors found: ", 14), bsn.Na ) -end - -""" - basins_map2D(xg, yg, integ; kwargs...) -Compute an estimate of the basin of attraction on a two-dimensional plane using a map of the plane onto itself. -The dynamical system should be a discrete two dimensional system such as: - * Discrete 2D map. - * 2D poincaré map. - * A 2D stroboscopic map. - -[H. E. Nusse and J. A. Yorke, Dynamics: numerical explorations, Springer, New York, 1997] - -## Arguments -* `xg`, `yg` : 1-dim range vector that defines the grid of the initial conditions to test. -* `integ` : A integrator handle of the dynamical system. For a Poincaré map, the handle is a `pmap` -as defined in [`poincaremap`](@ref) - -## Keyword Arguments -* `T` : Period of the stroboscopic map. -* `Ncheck` : A parameter that sets the number of consecutives hits of an attractor before deciding the basin -of the initial condition. - - -## Example: - -```jl -using DynamicalSystems, ChaosTools -ds = Systems.rikitake(μ = 0.47, α = 1.0) -xg=range(-6.,6.,length=150); yg=range(-6.,6.,length=150) -pmap = poincaremap(ds, (3, 0.), Tmax=1e6; idxs = 1:2, rootkw = (xrtol = 1e-8, atol = 1e-8), reltol=1e-9) -bsn = basins_map2D(xg, yg, pmap) -``` - -""" -function basins_map2D(xg, yg, pmap::PoincareMap; Ncheck = 3) - reinit_f! = (pmap,y) -> _init_map(pmap, y, pmap.i) - get_u = (pmap) -> pmap.integ.u[pmap.i] - basin = draw_basin!(xg, yg, pmap, step!, reinit_f!, get_u, Ncheck) -end - - -function _init_map(pmap::PoincareMap, y, idxs) - u = zeros(1,length(pmap.integ.u)) - u[idxs] = y - # all other coordinates are zero - reinit!(pmap, u) -end -# nueva linea - - -function basins_map2D(xg, yg, integ; T=0., Ncheck = 2) - if T>0 - iter_f! = (integ) -> step!(integ, T, true) - else - iter_f! = (integ) -> step!(integ) - end - reinit_f! = (integ,y) -> reinit!(integ, y) - get_u = (integ) -> integ.u - - return draw_basin!(xg, yg, integ, iter_f!, reinit_f!, get_u, Ncheck) -end - - - - - - -""" - basins_general(xg, yg, integ; T=1., idxs=1:2) -Compute an estimate of the basin of attraction on a two-dimensional plane using a stroboscopic map. - -## Arguments -* `xg`, `yg` : 1-dim range vector that defines the grid of the initial conditions to test. -* `integ` : integrator handle of the dynamical system. - -## Keyword Arguments -* `dt` : Time step of the integrator. It recommended to use values above 1. -* `idxs = 1:D` : Optionally you can choose which variables to save. Defaults to the entire state. -* `Ncheck` : A parameter that sets the number of consecutives hits of an attractor before deciding the basin -of the initial condition. - -## Example: -```jl -using DynamicalSystems, ChaosTools -ds = Systems.magnetic_pendulum(γ=1, d=0.2, α=0.2, ω=0.8, N=3) -integ = integrator(ds, u0=[0,0,0,0], reltol=1e-9) -xg=range(-4,4,length=150) -yg=range(-4,4,length=150) -@time bsn = basins_general(xg, yg, integ; dt=1., idxs=1:2) -""" -function basins_general(xg, yg, integ; dt=1., idxs=1:2, Ncheck = 10) - i = typeof(idxs) <: Int ? i : SVector{length(idxs), Int}(idxs...) - iter_f! = (integ) -> step!(integ, dt, true) - reinit_f! = (integ,y) -> _init_ds(integ, y, i) - get_u = (integ) -> integ.u[i] - return draw_basin!(xg, yg, integ, iter_f!, reinit_f!,get_u, Ncheck) -end - - -function _init_ds(integ, y, idxs) - u = zeros(length(integ.u)) - u[idxs] = y - # all other coordinates are zero - reinit!(integ, u) -end - - -## Procedure described in H. E. Nusse and J. A. Yorke, Dynamics: numerical explorations, Springer, New York, 1997 -# The idea is to color the grid with the current color. When an attractor box is hit (even color), the initial condition is colored -# with the color of its basin (odd color). If the trajectory hits another basin 10 times in row the IC is colored with the same -# color as this basin. -function procedure!(bsn_nfo::BasinInfo, n, u, Ncheck::Int) - max_check = 60 - #next_c = bsn_nfo.basin[n] - next_c=get_data(bsn_nfo.basin, n) - bsn_nfo.step += 1 - - - if iseven(next_c) && bsn_nfo.consecutive_match < max_check - # check wether or not we hit an attractor (even color). Make sure we hit Ncheck consecutive times. - if bsn_nfo.prev_attr == next_c - bsn_nfo.prevConsecutives +=1 - else - # reset prevConsecutives - bsn_nfo.prev_attr = next_c - bsn_nfo.prevConsecutives =1 - return 0; - end - - if bsn_nfo.prevConsecutives >= Ncheck - # Wait if we hit the attractor a Ncheck times in a row just to check if it is not a nearby trajectory - c3 = next_c+1 - if Ncheck == 2 - # For maps we can color the previous steps as well. Every point of the trajectory lead - # to the attractor - find_and_replace!(bsn_nfo, bsn_nfo.current_color+1, c3) - else - # For higher dimensions we erase the past iterations. - find_and_replace!(bsn_nfo, bsn_nfo.current_color+1, 1) - end - reset_bsn_nfo!(bsn_nfo) - return c3 - end - end - - if next_c == 1 && bsn_nfo.consecutive_match < max_check - # uncolored box, color it with current odd color - set_data!(bsn_nfo.basin, n, bsn_nfo.current_color+1) - push!(bsn_nfo.visited,n) - bsn_nfo.consecutive_match = 0 - return 0 - elseif next_c == 1 && bsn_nfo.consecutive_match >= max_check - # Maybe chaotic attractor, perodic or long recursion. - # Color this box as part of an attractor - set_data!(bsn_nfo.basin, n, bsn_nfo.current_color) - store_attractor!(bsn_nfo, u) - bsn_nfo.consecutive_match = max_check - return 0 - elseif next_c == bsn_nfo.current_color + 1 - # hit a previously visited box with the current color, possible attractor? - if bsn_nfo.consecutive_match < max_check - bsn_nfo.consecutive_match += 1 - return 0 - else - set_data!(bsn_nfo.basin, n, bsn_nfo.current_color) - store_attractor!(bsn_nfo, u) - # We continue iterating until we hit again the same attractor. In which case we stop. - return 0; - end - elseif isodd(next_c) && 0 < next_c < bsn_nfo.current_color && bsn_nfo.consecutive_match < max_check && Ncheck == 2 - # hit a colored basin point of the wrong basin, happens all the time, we check if it happens - #10 times in a row or if it happens N times along the trajectory whether to decide if it is another basin. - bsn_nfo.consecutive_other_basins += 1 - - if bsn_nfo.prev_bas == next_c && bsn_nfo.prev_step == bsn_nfo.step-1 - bsn_nfo.prevConsecutives +=1 - bsn_nfo.prev_step += 1 - else - bsn_nfo.prev_bas = next_c - bsn_nfo.prev_step = bsn_nfo.step - bsn_nfo.prevConsecutives =1 - end - - if bsn_nfo.consecutive_other_basins > 60 || bsn_nfo.prevConsecutives > 10 - find_and_replace!(bsn_nfo, bsn_nfo.current_color+1, next_c) - reset_bsn_nfo!(bsn_nfo) - return next_c - end - return 0 - elseif iseven(next_c) && (max_check <= bsn_nfo.consecutive_match < 2*max_check) - # We make sure we hit the attractor 60 consecutive times - bsn_nfo.consecutive_match+=1 - return 0 - elseif iseven(next_c) && bsn_nfo.consecutive_match >= max_check*2 - # We have checked the presence of an attractor: tidy up everything and get a new box. - find_and_replace!(bsn_nfo, bsn_nfo.current_color+1, 1) - # store color - set_data!(bsn_nfo.basin, n, bsn_nfo.current_color) - store_attractor!(bsn_nfo, u) # store attractor - bsn_nfo.current_color = bsn_nfo.next_avail_color - bsn_nfo.next_avail_color += 2 - println("even and > max check ", next_c) - reset_bsn_nfo!(bsn_nfo) - return next_c+1; - else - return 0 - end -end - - -function store_attractor!(bsn_nfo::BasinInfo, u) - if haskey(bsn_nfo.attractors , bsn_nfo.current_color) - push!(bsn_nfo.attractors[bsn_nfo.current_color], u) # store attractor - else - bsn_nfo.attractors[bsn_nfo.current_color] = Dataset([SVector(u[1],u[2])]) # init dic - end -end - - -""" - draw_basin!(xg, yg, integ, iter_f!::Function, reinit_f!::Function) -Compute an estimate of the basin of attraction on a two-dimensional plane. This is a low level function, -for higher level functions see: `basin_poincare_map`, `basin_discrete_map`, `basin_stroboscopic_map` - -## Arguments -* `xg`, `yg` : 1-dim range vector that defines the grid of the initial conditions to test. -* `integ` : integrator handle of the dynamical system. -* `iter_f!` : function that iterates the map or the system, see step! from DifferentialEquations.jl and -examples for a Poincaré map of a continuous system. -* `reinit_f!` : function that sets the initial condition to test on a two dimensional projection of the phase space. -""" -function draw_basin!(xg, yg, integ, iter_f!::Function, reinit_f!::Function, get_u::Function, Ncheck) - - complete = 0; - - bsn_nfo = BasinInfo(ones(Int16, length(xg), length(yg)), xg, yg, iter_f!, reinit_f!, get_u, 2,4,0,0,0,1,1,0,0, - Dict{Int16,Dataset{2,Float64}}(),0,Vector{CartesianIndex}(),Vector{CartesianIndex}()) - - reset_bsn_nfo!(bsn_nfo) - - I = CartesianIndices(bsn_nfo.basin) - j = 1 - - while complete == 0 - # pick the first empty box - ind = 0 - for k in j:length(bsn_nfo.basin) - if bsn_nfo.basin[I[k]] == 1 - j = k - ind=I[k] - break - end - end - - if ind == 0 - # We are done - complete = 1 - break - end - - ni = ind[1]; mi = ind[2]; - x0 = xg[ni]; y0 = yg[mi] - - # Tentatively assign a color: odd is for basins, even for attractors. - # First color is one - bsn_nfo.basin[ni,mi] = bsn_nfo.current_color + 1 - u0=[x0, y0] - bsn_nfo.basin[ni,mi] = get_color_point!(bsn_nfo, integ, u0; Ncheck=Ncheck) - end - - bsn_nfo.Na = Int((bsn_nfo.current_color-2)/2) - - return bsn_nfo -end - - - -function get_color_point!(bsn_nfo::BasinInfo, integ, u0; Ncheck=2) - # This routine identifies the attractor using the previously defined basin. - # reinitialize integrator - bsn_nfo.reinit_f!(integ, u0) - reset_bsn_nfo!(bsn_nfo) - - done = 0; - inlimbo = 0 - - while done == 0 - old_u = bsn_nfo.get_u(integ) - bsn_nfo.iter_f!(integ) - new_u = bsn_nfo.get_u(integ) - - n = get_box(new_u, bsn_nfo) - - if !isnothing(n) # apply procedure only for boxes in the defined space - done = procedure!(bsn_nfo, n, new_u, Ncheck) - inlimbo = 0 - else - # We are outside the defined grid - inlimbo +=1 - end - - if inlimbo > 60 - done = check_outside_the_screen!(bsn_nfo, new_u, old_u, inlimbo) - end - end - - return done -end - - - - -function get_box(u, bsn_nfo::BasinInfo) - xg = bsn_nfo.xg; yg = bsn_nfo.yg; # aliases - xstep = (xg[2]-xg[1]) - ystep = (yg[2]-yg[1]) - xu=u[1] - yu=u[2] - if typeof(bsn_nfo.basin) <: Cell - if xu >= xg[1] && xu <= xg[end] && yu >= yg[1] && yu <= yg[end] - return findleaf(bsn_nfo.basin,u) - else - return nothing - end - end - n=0; m=0; - # check boundary - if xu >= xg[1] && xu <= xg[end] && yu >= yg[1] && yu <= yg[end] - n = Int(round((xu-xg[1])/xstep))+1 - m = Int(round((yu-yg[1])/ystep))+1 # +1 for 1 based indexing - else - return nothing - end - return CartesianIndex(n,m) -end - - -function check_outside_the_screen!(bsn_nfo::BasinInfo, new_u, old_u, inlimbo) - - if norm(new_u-old_u) < 1e-5 - find_and_replace!(bsn_nfo,bsn_nfo.current_color+1, 1) - reset_bsn_nfo!(bsn_nfo) - # this CI goes to a attractor outside the screen, set to -1 (even color) - return -1 # get next box - elseif inlimbo > 60*20 - find_and_replace!(bsn_nfo,bsn_nfo.current_color+1, 1) - reset_bsn_nfo!(bsn_nfo) - # this CI is problematic or diverges, set to -1 (even color) - return -1 # get next box - end - return 0 -end - - -function find_and_replace!(bsn_nfo::BasinInfo, old_c, new_c) - if typeof(bsn_nfo.basin) <: Matrix - while !isempty(bsn_nfo.visited) - ind = pop!(bsn_nfo.visited) - if bsn_nfo.basin[ind] == old_c - bsn_nfo.basin[ind] = new_c - end - end - # ind = (bsn_nfo.basin .== old_c) - # bsn_nfo.basin[ind] .= new_c - else - while !isempty(bsn_nfo.visited) - lf = pop!(bsn_nfo.visited) - if lf.data == old_c - lf.data = new_c - end - end - # for lf in allleaves(bsn_nfo.basin) - # if lf.data==old_c - # lf.data==new_c - # end - # end - - end -end - - - - -function reset_bsn_nfo!(bsn_nfo::BasinInfo) - #@show bsn_nfo.step - bsn_nfo.consecutive_match = 0 - bsn_nfo.consecutive_other_basins = 0 - bsn_nfo.prevConsecutives = 0 - bsn_nfo.prev_attr = 1 - bsn_nfo.prev_bas = 1 - bsn_nfo.prev_step = 0 - bsn_nfo.step = 0 -end - -function basins_map2D_tree(xg, yg, integ; T=0., Ncheck = 2, r_init=0.2, r_max=0.1) - if T>0 - iter_f! = (integ) -> step!(integ, T, true) - else - iter_f! = (integ) -> step!(integ) - end - reinit_f! = (integ,y) -> reinit!(integ, y) - get_u = (integ) -> integ.u - - return draw_basin_tree!(xg, yg, integ, iter_f!, reinit_f!, get_u, Ncheck, r_init, r_max) -end - - -import RegionTrees: AbstractRefinery, needs_refinement, refine_data - -struct MyRefinery <: AbstractRefinery - tolerance::Float64 -end - -# These two methods are all we need to implement -function needs_refinement(r::MyRefinery, cell) - return maximum(cell.boundary.widths) > r.tolerance -end - -function refine_data(r::MyRefinery, cell::Cell, indices) - return 1 -end - -function is_in_cell(cell::Cell, point::AbstractVector, indices) - or = cell.boundary.origin - wd = cell.boundary.widths - dv = cell.divisions - if point[1] > or[1] + wd[1] || point[1] < or[1] - return false - elseif point[2] > or[2] + wd[2] || point[2] < or[2] - return false - else - ind = [1,1] - if point[1] > dv[1] - ind[1]=2 - end - if point[2] > dv[2] - ind[2]=2 - end - if ind[1] == indices[1] && ind[2] == indices[2] - return true - else - return false - end - end -end - - -function refine_data(cell::Cell, indices, bsn_nfo::BasinInfo) - if iseven(cell.data) - att = cell.data - for p in bsn_nfo.attractors[att] - if is_in_cell(cell,p, indices) - #@show att,p,cell.boundary,indices - return att - end - end - end - return 1 -end - -function get_neighbor_cells(root::Cell,cell::Cell) - or = cell.boundary.origin - wd = cell.boundary.widths - dv = cell.divisions - W=findleaf(root,dv + SVector(wd[1],0)).data - E=findleaf(root,dv + SVector(-wd[1],0)).data - N=findleaf(root,dv + SVector(0,wd[2])).data - S=findleaf(root,dv + SVector(0,-wd[2])).data - #@show [W E N S cell.data] - return length(unique([W E N S cell.data])) > 1 -end - -function check_if_complete!(root::Cell, r_max, bsn_nfo::BasinInfo) - complete = true - refine_function = (cell, indices) -> refine_data(cell, indices, bsn_nfo) - leaf_stack = Vector{Cell}() - for lf in allleaves(root) - if maximum(lf.boundary.widths) > r_max && get_neighbor_cells(root,lf) - complete = false - push!(leaf_stack,lf) - end - end - - while !isempty(leaf_stack) - l = pop!(leaf_stack) - split!(l, refine_function) - for cl in children(l) - if cl.data == 1 - push!(bsn_nfo.available,cl) - end - end - end - - return complete -end - - - -function draw_basin_tree!(xg, yg, integ, iter_f!::Function, reinit_f!::Function, get_u::Function, Ncheck, r_init, r_max) - - complete = 0; - - xi=xg[1]; yi=yg[1]; xf=xg[end]; yf=yg[end] - - bsn_nfo = BasinInfo(Cell(SVector(xi, yi), SVector(xf-xi, yf-yi),1), xg, yg, iter_f!, reinit_f!, get_u, 2,4,0,0,0,1,1,0,0, - Dict{Int16,Dataset{2,Float64}}(),0,Vector{Cell}(),Vector{Cell}()) - - reset_bsn_nfo!(bsn_nfo) - - # Initial refinement of the grid. - r = MyRefinery(r_init) - @show bsn_nfo.basin - adaptivesampling!(bsn_nfo.basin, r) - - for lf in allleaves(bsn_nfo.basin) - push!(bsn_nfo.available,lf) - end - - while complete == 0 - - # Get next available candidate - lf = nothing - while !isempty(bsn_nfo.available) - l = pop!(bsn_nfo.available) - if l.data == 1 - lf = l - break - end - end - - - if isempty(bsn_nfo.available) - println("refinement!") - if check_if_complete!(bsn_nfo.basin, r_max, bsn_nfo) - complete=1 - break - end - # Get next available candidate - while !isempty(bsn_nfo.available) - l = pop!(bsn_nfo.available) - if l.data == 1 - lf = l - break - end - end - - println("end refinement! res=", lf.boundary.widths) - - end - - # Tentatively assign a color: odd is for basins, even for attractors. - # First color is one - lf.data = bsn_nfo.current_color + 1 - u0 = RegionTrees.center(lf) - lf.data = get_color_point!(bsn_nfo, integ, u0; Ncheck=Ncheck) - - end - bsn_nfo.Na = (bsn_nfo.current_color-2)÷2 - return bsn_nfo -end - -function get_data(bsn::Matrix,n) - return bsn[n] -end - -function get_data(bsn::Cell,n) - return n.data -end - -function set_data!(bsn::Matrix,n,c) - bsn[n]=c -end - -function set_data!(bsn::Cell,n,c) - n.data=c -end +## All functions in this file have been merged into ChaosTools!! diff --git a/src/examples/tst_ba_lorenz.jl b/src/examples/tst_ba_lorenz.jl index 0a4a436..c49bcaf 100644 --- a/src/examples/tst_ba_lorenz.jl +++ b/src/examples/tst_ba_lorenz.jl @@ -1,26 +1,15 @@ using Revise using Plots using DynamicalSystems -using Basins - +using OrdinaryDiffEq # Multistability, phase diagrams, and intransitivity in the Lorenz-84 low-order atmospheric circulation model # Chaos 18, 033121 (2008); https://doi.org/10.1063/1.2953589 -@inline @inbounds function lorenz84(u, p, t) - F = p[1]; G = p[2]; a = p[3]; b = p[4]; - x = u[1]; y = u[2]; z = u[3]; - dx = -y^2 -z^2 -a*x + a*F - dy = x*y - y - b*x*z +G - dz = b*x*y + x*z -z - return SVector{3}(dx, dy, dz) -end - - F=6.846; G=1.287; a=0.25; b=4.; -p= [F, G, a, b] -ds = ContinuousDynamicalSystem(lorenz84, rand(3), p) -xg=range(-1.,1.,length=200) -yg=range(-1.5,1.6,length=200) -pmap = poincaremap(ds, (3, 0.), Tmax=1e6; idxs = 1:2, rootkw = (xrtol = 1e-8, atol = 1e-8), reltol=1e-9) -@time bsn = basins_map2D(xg, yg, pmap) -plot(xg,yg,bsn.basin',seriestype=:heatmap) +ds = Systems.lorenz84(; F,G,a,b) +xg=range(-5.,5.,length = 500) +yg=range(-5.,5.,length = 500) +zg=range(-2.,2.,length = 10) +#pmap = poincaremap(ds, (3, 0.)) +@time bsn, att = basins_of_attraction((xg, yg, zg), ds) +plot(xg,yg,bsn[:,:,5]',seriestype=:heatmap) diff --git a/src/examples/tst_ba_mag_pend.jl b/src/examples/tst_ba_mag_pend.jl index 2f13c8a..e6fa5c0 100644 --- a/src/examples/tst_ba_mag_pend.jl +++ b/src/examples/tst_ba_mag_pend.jl @@ -1,46 +1,19 @@ -using Revise using Plots -using Basins using DynamicalSystems -# struct MagPendulum{T<:AbstractFloat} -# magnets::Vector{SVector{2, T}} -# end -# -# function (m::MagPendulum)(u, p, t) -# x, y, vx, vy = u -# γ, d, α, ω = p -# dx, dy = vx, vy -# dvx, dvy = @. -ω^2*(x, y) - α*(vx, vy) -# for ma in m.magnets -# δx, δy = (x - ma[1]), (y - ma[2]) -# D = sqrt(δx^2 + δy^2 + d^2) -# dvx -= γ*(x - ma[1])/D^5 -# dvy -= γ*(y - ma[2])/D^5 -# end -# return SVector(dx, dy, dvx, dvy) -# end -# -# function mag_pendulum(u = [sincos(0.12553*2π)..., 0, 0]; -# γ = 1.0, d = 0.3, α = 0.2, ω = 0.5, N = 3) -# m = MagPendulum([SVector(cos(2π*i/N), sin(2π*i/N)) for i in 1:N]) -# p = [γ, d, α, ω] -# ds = ContinuousDynamicalSystem(m, u, p) -# end - # https://cdn.ima.org.uk/wp/wp-content/uploads/2020/03/Chaos-in-the-Magnetic-Pendulum-from-MT-April-2020.pdf #ds = mag_pendulum(γ=1, d=0.5, α=0.175, ω=1., N=4) ds = Systems.magnetic_pendulum(γ=1, d=0.3, α=0.2, ω=0.5, N=4) -integ = integrator(ds, u0=[0,0,0,0], reltol=1e-9) -xg=range(-2,2,length=500) -yg=range(-2,2,length=500) -@time bsn = Basins.basins_general(xg, yg, integ; dt=1., idxs=1:2) +xg=range(-2, 2,length = 400) +yg=range(-2, 2,length = 400) +default_diffeq = (reltol = 1e-9, alg = Vern9()) +@time bsn,att = basins_of_attraction((xg, yg), ds; diffeq = default_diffeq) # Uncertainty exponent for these parameter and grid -bd = box_counting_dim(xg, yg, bsn) +v, e, bd = basins_fractal_dimension(bsn) α = 2 - bd -plot(xg,yg,bsn.basin',seriestype=:heatmap) +plot(xg, yg, bsn',seriestype = :heatmap) diff --git a/src/examples/tst_ba_more_chaos.jl b/src/examples/tst_ba_more_chaos.jl index 5c692f3..e273501 100644 --- a/src/examples/tst_ba_more_chaos.jl +++ b/src/examples/tst_ba_more_chaos.jl @@ -11,5 +11,5 @@ integ=integrator(ds) xg=range(-10.,10.,length=100) yg=range(-10.,10.,length=100) pmap = poincaremap(ds, (3, 0.), Tmax=1e6; idxs = 1:2, rootkw = (xrtol = 1e-8, atol = 1e-8), reltol=1e-9) -@time bsn = Basins.basins_map2D(xg, yg, pmap) -plot(xg,yg,bsn.basin',seriestype=:heatmap) +@time bsn,att = basins_of_attraction((xg, yg), pmap) +plot(xg, yg, bsn',seriestype = :heatmap) diff --git a/src/examples/tst_ba_newton.jl b/src/examples/tst_ba_newton.jl index 120ef0f..27771c3 100644 --- a/src/examples/tst_ba_newton.jl +++ b/src/examples/tst_ba_newton.jl @@ -1,10 +1,7 @@ using Revise using Plots -using Basins using DynamicalSystems - - function newton_map(dz,z, p, n) f(x) = x^p[1]-1 df(x)= p[1]*x^(p[1]-1) @@ -21,52 +18,35 @@ function newton_map_J(J,z0, p, n) return end ds = DiscreteDynamicalSystem(newton_map,[0.1, 0.2], [3] , newton_map_J) -integ = integrator(ds) -a=zeros(1,8) -b=zeros(1,8) -c=zeros(1,8) -d=zeros(1,8) +xg=range(-1.,1.,length=1000) +yg=range(-1.,1.,length=1000) -# for (k,res) in enumerate(range(100,800,step=100)) -# -# xg=range(-1.,1.,length=res) -# yg=range(-1.,1.,length=res) -# -# @time bsn=Basins.basins_map2D(xg, yg, integ) -# -# a[k],De,e,f=uncertainty_exponent(bsn,integ) -# b[k],e,f=Basins.static_estimate(xg,yg,bsn) -# c[k],e,f=Basins.ball_estimate(xg,yg,bsn) -# d[k]=2-box_counting_dim(xg,yg,bsn) -# end +@time bsn, att = basins_of_attraction((xg, yg), ds) -xg=range(-1.,1.,length=300) -yg=range(-1.,1.,length=300) - -@time bsn=Basins.basins_map2D(xg, yg, integ) - - -b=zeros(1,40) -c=zeros(1,40) - -precision = 1e-5 -for k in 1:40 - #a[k],De,e,f=uncertainty_exponent(bsn,integ) - b[k],e,f=Basins.static_estimate(xg,yg,bsn; precision=precision) - c[k],e,f=Basins.ball_estimate(xg,yg,bsn; precision=precision) - #d[k]=2-box_counting_dim(xg,yg,bsn) -end +plot(xg, yg, bsn',seriestype = :heatmap) -using Statistics -println(mean(b), " ", mean(c)) -#plot(xg,yg,bsn.basin',seriestype=:heatmap) -plot(b') -plot!(c') - - -#sa,sb = compute_saddle(integ, bsn, [1], [2,3]; N=100) - -#s=Dataset(sa) -#plot!(s[:,1],s[:,2],seriestype=:scatter) +# b=zeros(1,40) +# c=zeros(1,40) +# +# precision = 1e-5 +# for k in 1:40 +# #a[k],De,e,f=uncertainty_exponent(bsn,integ) +# b[k],e,f=Basins.static_estimate(xg,yg,bsn; precision=precision) +# c[k],e,f=Basins.ball_estimate(xg,yg,bsn; precision=precision) +# #d[k]=2-box_counting_dim(xg,yg,bsn) +# end +# +# +# using Statistics +# println(mean(b), " ", mean(c)) +# #plot(xg,yg,bsn.basin',seriestype=:heatmap) +# plot(b') +# plot!(c') +# +# +# #sa,sb = compute_saddle(integ, bsn, [1], [2,3]; N=100) +# +# #s=Dataset(sa) +# #plot!(s[:,1],s[:,2],seriestype=:scatter) diff --git a/src/examples/tst_ba_riddled.jl b/src/examples/tst_ba_riddled.jl index 5cfb34d..6cb2131 100644 --- a/src/examples/tst_ba_riddled.jl +++ b/src/examples/tst_ba_riddled.jl @@ -4,7 +4,7 @@ using Basins using DifferentialEquations using DynamicalSystems -# Blowout bifurcation with non-normal parameters in population dynamicsBernard Cazelle +# Blowout bifurcation with non-normal parameters in population dynamics Bernard Cazelle function franke_yabuku(dz,z, p, n) r=2.825; s=20.25; c1=1.2; c2=0.1; x = z[1]; y = z[2]; @@ -20,12 +20,11 @@ function franke_yabuku_J(J,z0, p, n) end ds = DiscreteDynamicalSystem(franke_yabuku,[0.1, 0.2], [] , franke_yabuku_J) -integ = integrator(ds) xg=range(8.,16.,length=100) yg=range(0.5,1.,length=100) # compute basin -@time basin=basins_map2D(xg, yg, integ) -plot(xg,yg,basin',seriestype=:heatmap) +@time basins, att = basins_of_attraction((xg, yg), ds) +plot(xg, yg, basins', seriestype = :heatmap) diff --git a/src/examples/tst_ba_rikitake.jl b/src/examples/tst_ba_rikitake.jl index 2443a14..f7c8f91 100644 --- a/src/examples/tst_ba_rikitake.jl +++ b/src/examples/tst_ba_rikitake.jl @@ -1,14 +1,22 @@ using Revise using Plots using DynamicalSystems -using DifferentialEquations -using Basins +using OrdinaryDiffEq μ=0.47 ds = Systems.rikitake(μ = μ, α = 1.0) -integ=integrator(ds) -xg=range(-6.,6.,length=200) -yg=range(-6.,6.,length=200) -pmap = poincaremap(ds, (3, 0.), Tmax=1e6; idxs = 1:2, rootkw = (xrtol = 1e-8, atol = 1e-8), reltol=1e-9) -@time bsn = Basins.basins_map2D(xg, yg, pmap) -plot(xg,yg,bsn.basin',seriestype=:heatmap) +xg = yg = range(-2.5,2.5,length = 100) +df_de = (;reltol = 1e-10, abstol = 1e-10, alg = Vern9(), maxiters = Int(1e12)) +pmap = poincaremap(ds, (3, 0.); rootkw = (xrtol = 1e-8, atol = 1e-8), diffeq = df_de) +bsn, att = basins_of_attraction((xg, yg), ds; mx_chk_loc_att = 100) +plot(xg, yg, bsn',seriestype = :heatmap) + +# This system is problematic: very very long transient (t > 2000 sometimes) that can be mistaken with attractors. + + + delete!(att, 3) # remove "wrong" attractor + delete!(att, 4) # remove "wrong" attractor + + xg = yg = range(-2.5,2.5,length = 400) + bsn, att = basins_of_attraction((xg, yg), ds; attractors = att) + plot(xg, yg, bsn',seriestype = :heatmap) diff --git a/src/examples/tst_ba_thomas.jl b/src/examples/tst_ba_thomas.jl index 555e826..27c02bb 100644 --- a/src/examples/tst_ba_thomas.jl +++ b/src/examples/tst_ba_thomas.jl @@ -1,18 +1,16 @@ using Revise using Plots using DynamicalSystems -using DifferentialEquations -using Basins +using OrdinaryDiffEq ### First version using high-level function. b=0.1665 #b=0.22 ds = Systems.thomas_cyclical(b = b) -#integ=integrator(ds, reltol=1e-9) xg=range(-6.,6.,length=200) yg=range(-6.,6.,length=200) -pmap = poincaremap(ds, (3, 0.), Tmax=1e6; idxs = 1:2, rootkw = (xrtol = 1e-8, atol = 1e-8), reltol=1e-9) -@time bsn = Basins.basins_map2D(xg, yg, pmap) +pmap = poincaremap(ds, (3, 0.), 1e6; diffeq = (;reltol = 1e-9, alg = Vern9())) +@time bsn, att = basins_of_attraction((xg, yg), pmap) -plot(xg,yg,bsn.basin',seriestype=:heatmap) +plot(xg, yg, bsn', seriestype = :heatmap) diff --git a/src/examples/tst_coupled_maps.jl b/src/examples/tst_coupled_maps.jl index f5ad8dd..60116dc 100644 --- a/src/examples/tst_coupled_maps.jl +++ b/src/examples/tst_coupled_maps.jl @@ -24,12 +24,12 @@ integ = integrator(ds) xg=range(0.,0.999999,length=100) yg=range(0.,0.999999,length=100) -@time bsn=Basins.basins_general(xg, yg, integ; dt=1) -plot(xg,yg,bsn.basin',seriestype=:heatmap) +@time bsn, att = basins_of_attraction((xg, yg), ds) +plot(xg,yg,bsn',seriestype=:heatmap) -# Estimation using box counting -ind = findall(iseven.(bsn.basin) .== true) -basin_test = deepcopy(bsn.basin) -[basin_test[k] =basin_test[k]+1 for k in ind ] -bd = box_counting_dim(xg,yg, basin_test) -α = 2 - bd +# # Estimation using box counting +# ind = findall(iseven.(bsn.basin) .== true) +# basin_test = deepcopy(bsn.basin) +# [basin_test[k] =basin_test[k]+1 for k in ind ] +# bd = box_counting_dim(xg,yg, basin_test) +# α = 2 - bd diff --git a/src/examples/tst_duffing_all.jl b/src/examples/tst_duffing_all.jl index 1be4029..9989fdd 100644 --- a/src/examples/tst_duffing_all.jl +++ b/src/examples/tst_duffing_all.jl @@ -1,8 +1,8 @@ using Revise -using Basins using Printf using Plots -using DynamicalSystems +using DynamicalSystems, DifferentialEquations +using Basins @inline @inbounds function duffing(u, p, t) d = p[1]; F = p[2]; omega = p[3] @@ -12,48 +12,26 @@ using DynamicalSystems end #F=0.3818791946308725; ω= 0.1966442953020134 -F=0.2771812080536913; ω=0.1; # smooth boundary -#ω=0.1617;F = 0.395 -#ω=0.3;F = 0.1 +# F=0.2771812080536913; ω=0.1; # smooth boundary +ω=0.1617;F = 0.395 +# ω=0.3;F = 0.1 ds = ContinuousDynamicalSystem(duffing, rand(2), [0.15, F, ω]) -integ_df = integrator(ds; reltol=1e-8, save_everystep=false) -xg = yg = range(-2.2,2.2,length=100) -@time bsn = Basins.basins_map2D(xg, yg, integ_df; T=2*pi/ω) +xg = yg = range(-2.2,2.2,length=400) +grid = (xg,yg) +@time bsn, att = basins_of_attraction(grid, ds; T=2*pi/ω, diffeq = (;reltol = 1e-9, alg = Vern9())) # Basin entropy -@show Sb,Sbb = basin_entropy(bsn; eps_x=20, eps_y=20) +@show Sb,Sbb = basin_entropy(bsn) + +# Basins fractal test +test_res, Sbb = basins_fractal_test(bsn; ε = 5, Ntotal = 10000) # Wada merge Haussdorff distances -@time max_dist,min_dist = detect_wada_merge_method(xg, yg, bsn) -epsilon = xg[2]-xg[1] -@show dmax = max_dist/epsilon -@show dmin = min_dist/epsilon +@show max_dist,min_dist = detect_wada_merge_method(bsn) # Wada grid -W = detect_wada_grid_method(integ_df, bsn; max_iter=8) +W = detect_wada_grid_method(grid, ds; attractors = att, basins = bsn, max_iter=8, T = 2*pi/ω, diffeq = (;reltol = 1e-9, alg = Vern9())) @show W[:,end] -# Uncertainty exponent for these parameter and grid -bd = box_counting_dim(xg, yg, bsn) -α = 2 - bd - -D = uncertainty_exponent(bsn, integ_df) -@show 2-D - -println("---------------") -println("---------------") -println("Basin Report: ") -println("---------------") -println("---------------") - -@printf("Basin entropy %.2f \n", Sb) -@printf("Boundary Basin Entropy: %.2f\n", Sbb) -@printf("Uncertainty exponent: α= %.2f\n", α ) -@printf("Box counting dim: bd= %.2f\n", bd) -@printf("Uncertainty dim estimator: d = %.2f\n", 2-D[1]) -@printf("Number of basins: %d\n", bsn.Na) -@printf("Merge Method: Max fattening parameter: %.2f\n", dmax) -@printf("Wada Grid Method: W_Na = %.2f\n ", W[end,end] ) - - -plot(xg, yg, bsn.basin', seriestype=:heatmap) +f,e,α = uncertainty_exponent(bsn) +@show α diff --git a/src/examples/tst_duffing_straddle.jl b/src/examples/tst_duffing_straddle.jl deleted file mode 100644 index 43ebd21..0000000 --- a/src/examples/tst_duffing_straddle.jl +++ /dev/null @@ -1,80 +0,0 @@ -using Revise -using Basins -using Plots -using DynamicalSystems -using DifferentialEquations - -# -# -# ω=1. -# F = 0.2 -# ds =Systems.duffing([0.1, 0.25]; ω = ω, f = F, d = 0.15, β = -1) -# integ_df = integrator(ds; alg=Tsit5(), reltol=1e-8, save_everystep=false) -# xg = range(-2.2,2.2,length=150) -# yg = range(-2.2,2.2,length=150) -# -# @time bsn = Basins.basins_map2D(xg, yg, integ_df; T=2*pi/ω) -# -# sa,sb = compute_saddle(integ_df, bsn, [1], [2]; N=1000) - -# plot(xg,yg,bsn.basin', seriestype=:heatmap) -# s = Dataset(sa) -# plot!(s[:,1],s[:,2],seriestype=:scatter) - -# io = open("myfile.txt", "w"); -# for v in s -# write(io, string(v[1], " ", v[2], "; \n")); -# end -# close(io); - - - - -# Forced pendulum - -# Equations of motion: -function forced_pendulum!(du, u, p, t) - d = p[1]; F = p[2]; omega = p[3] - du[1] = u[2] - du[2] = -d*u[2] - sin(u[1])+ F*cos(omega*t) -end - -# We have to define a callback to wrap the phase in [-π,π] -function affect!(integrator) - if integrator.u[1] < 0 - integrator.u[1] += 2*π - else - integrator.u[1] -= 2*π - end -end - -condition(u,t,integrator) = (integrator.u[1] < -π || integrator.u[1] > π) - -cb = DiscreteCallback(condition,affect!) - -#d, F ,w -F = 1.66 -ω = 1. -d=0.2 -p=[d, F, ω] -#p=[0.15, 0.2, 0.1] -df = ODEProblem(forced_pendulum!,rand(2),(0.0,20.0), p) -integ_df = init(df, alg=AutoTsit5(Rosenbrock23()); reltol=1e-9, abstol=1e-9, save_everystep=false, callback=cb) - -xres=200 -yres=200 - -# range for forced pend -xg = range(-pi,pi,length=xres) -yg = range(-2.,4.,length=yres) - -# compute basin -@time bsn = Basins.basins_map2D(xg, yg, integ_df; T=2*pi/ω) - -Na = length(unique(bsn.basin)) - -sa,sb = compute_saddle(integ_df, bsn, [1], [2,3]; N=1000) - -plot(xg,yg,bsn.basin', seriestype=:heatmap) -s = Dataset(sa) -plot!(s[:,1],s[:,2],seriestype=:scatter, markercolor=:blue) diff --git a/src/examples/tst_henon_all.jl b/src/examples/tst_henon_all.jl index 29e772d..ca8c469 100644 --- a/src/examples/tst_henon_all.jl +++ b/src/examples/tst_henon_all.jl @@ -7,20 +7,16 @@ using Printf ds = Systems.henon(zeros(2); a = 1.4, b = 0.3) -integ_df = integrator(ds) - -xres=200 -yres=200 - +res=200 # range for forced pend -xg = range(-2.,2.,length=xres) -yg = range(-2.,2.,length=yres) +xg = range(-2.,2.,length=res) +yg = range(-2.,2.,length=res) # compute basin -@time basin = Basins.basins_map2D(xg, yg, integ_df) +@time basins, att = basins_of_attraction((xg, yg), ds) # Basin entropy -@show Sb,Sbb = basin_entropy(bsn.basin; eps_x=20, eps_y=20) +@show Sb,Sbb = basin_entropy(basins; eps_x=20, eps_y=20) # Wada merge Haussdorff distances @time max_dist,min_dist = detect_wada_merge_method(xg, yg, bsn.basin) diff --git a/src/examples/tst_map_feudel.jl b/src/examples/tst_map_feudel.jl index 6550a5a..a999e8e 100644 --- a/src/examples/tst_map_feudel.jl +++ b/src/examples/tst_map_feudel.jl @@ -38,7 +38,7 @@ xg = range(0.0, 1.0, length = 250) integ.p[2] = 3.833 integ.p[1] = 0.0015 -@time bsn = Basins.basins_map2D(xg, θ, integ) +@time bsn, att = basins_of_attraction((xg, θ), ds) #@time bsn=Basins.basin_general_ds(xg, θ, integ; dt=1, Ncheck=10) @show bd = box_counting_dim(xg, θ, bsn) diff --git a/src/examples/tst_map_grbogi.jl b/src/examples/tst_map_grbogi.jl index 351cad7..c2e5836 100644 --- a/src/examples/tst_map_grbogi.jl +++ b/src/examples/tst_map_grbogi.jl @@ -20,13 +20,12 @@ function grebogi_map_J(J,z0, p, n) return end ds = DiscreteDynamicalSystem(grebogi_map,[1., -1.], [] , grebogi_map_J) -integ = integrator(ds) -θg=range(0,2π,length=600) -xg=range(-0.5,0.5,length=600) -@time bsn=Basins.basins_map2D(θg, xg, integ) -#@time bns2=ChaosTools.basin_map(θg, xg, integ) +θg = range(0, 2π, length = 600) +xg = range(-0.5, 0.5, length = 600) + +@time bsn, att = basins_of_attration((θg, xg), ds) #plot(θg,xg,bsn.basin',seriestype=:heatmap) diff --git a/src/examples/tst_map_mcdonald.jl b/src/examples/tst_map_mcdonald.jl index 94c1cc4..3870c3e 100644 --- a/src/examples/tst_map_mcdonald.jl +++ b/src/examples/tst_map_mcdonald.jl @@ -9,7 +9,6 @@ using BenchmarkTools # Steven W. MCDONALD, Celso GREBOGY Edward OTT and James A. YORKE # Physica D 17 1985 function chaotic_map(dz,z, p, n) - xn = z[1]; yn = z[2] λ₁=p[1];λ₂=p[2]; dz[1]=mod(λ₁*xn,1) diff --git a/src/examples/tst_pendulum_all.jl b/src/examples/tst_pendulum_all.jl index ad1d432..e824abd 100644 --- a/src/examples/tst_pendulum_all.jl +++ b/src/examples/tst_pendulum_all.jl @@ -1,26 +1,19 @@ using Revise -using DifferentialEquations +using OrdinaryDiffEq using DynamicalSystems using Basins using Plots using Printf -# Equations of motion: -function forced_pendulum(u, p, t) - @inbounds begin - d = p[1]; F = p[2]; omega = p[3] - du1 = u[2] - du2 = -d*u[2] - sin(u[1])+ F*cos(omega*t) - return SVector{2}(du1, du2) - end -end - # We have to define a callback to wrap the phase in [-π,π] function affect!(integrator) + uu = integrator.u if integrator.u[1] < 0 - integrator.u[1] += 2*π + set_state!(integrator, SVector(uu[1] + 2π, uu[2])) + u_modified!(integrator, true) else - integrator.u[1] -= 2*π + set_state!(integrator, SVector(uu[1] - 2π, uu[2])) + u_modified!(integrator, true) end end @@ -32,40 +25,31 @@ cb = DiscreteCallback(condition,affect!) F = 1.66 ω = 1. d=0.2 -p=[d, F, ω] -#p=[0.15, 0.2, 0.1] -df = ODEProblem(forced_pendulum,rand(2),(0.0,20.0), p) -integ_df = init(df, alg=AutoTsit5(Rosenbrock23()); reltol=1e-9, save_everystep=false, callback=cb) - -xres=200 -yres=200 - +ds = Systems.forced_pendulum([0.,0.]; ω = ω, f = F, d = d); # range for forced pend -xg = range(-pi,pi,length=xres) -yg = range(-2.,4.,length=yres) - -# compute basin -@time bsn = Basins.basins_map2D(xg, yg, integ_df; T=2*pi/ω) +res = 300 +xg = range(-pi, pi,length = res) +yg = range(-2., 4.,length = res) +grid = (xg,yg) +default_diffeq = (reltol = 1e-9, alg = Vern9(), callback = cb) +bsn, att = basins_of_attraction(grid, ds; T = 2*pi/ω, diffeq = default_diffeq) # Basin entropy -@show Sb,Sbb = basin_entropy(bsn; eps_x=20, eps_y=20) +@show Sb,Sbb = basin_entropy(bsn) + +# Basins fractal test +test_res, Sbb = basins_fractal_test(bsn; ε = 5, Ntotal = 10000) # Wada merge Haussdorff distances -@time max_dist,min_dist = detect_wada_merge_method(xg, yg, bsn) -epsilon = xg[2]-xg[1] -@show dmax = max_dist/epsilon -@show dmin = min_dist/epsilon +@show max_dist,min_dist = detect_wada_merge_method(bsn) # Wada grid -W = detect_wada_grid_method(integ_df, bsn; max_iter=8) +W = detect_wada_grid_method(grid, ds; attractors = att, basins = bsn, max_iter=8, T = 2*pi/ω, diffeq = default_diffeq) @show W[:,end] # Uncertainty exponent for these parameter and grid -bd = box_counting_dim(xg, yg, bsn) -α = 2 - bd - -D = uncertainty_exponent(bsn, integ_df) -@show 2-D +ε, N_ε ,α = uncertainty_exponent(bsn) +@show α println("---------------") @@ -73,14 +57,14 @@ println("---------------") println("Basin Report: ") println("---------------") println("---------------") - -@printf("Basin entropy %.2f \n", Sb) -@printf("Boundary Basin Entropy: %.2f\n", Sbb) -@printf("Uncertainty exponent: α= %.2f\n", α ) -@printf("Box counting dim: bd= %.2f\n", bd) -@printf("Uncertainty dim estimator: d = %.2f\n", 2-D[1]) -@printf("Number of basins: %d\n", bsn.Na) -@printf("Merge Method: Max fattening parameter: %.2f\n", dmax) -@printf("Wada Grid Method: W_Na = %.2f\n ", W[end,end] ) - -plot(xg, yg, bsn.basin', seriestype=:heatmap) +# +# @printf("Basin entropy %.2f \n", Sb) +# @printf("Boundary Basin Entropy: %.2f\n", Sbb) +# @printf("Uncertainty exponent: α= %.2f\n", α ) +# @printf("Box counting dim: bd= %.2f\n", bd) +# @printf("Uncertainty dim estimator: d = %.2f\n", 2-D[1]) +# @printf("Number of basins: %d\n", bsn.Na) +# @printf("Merge Method: Max fattening parameter: %.2f\n", dmax) +# @printf("Wada Grid Method: W_Na = %.2f\n ", W[end,end] ) + +plot(xg, yg, basins', seriestype=:heatmap) diff --git a/src/examples/tst_riddles_ott.jl b/src/examples/tst_riddles_ott.jl index ea234f3..6c8efdc 100644 --- a/src/examples/tst_riddles_ott.jl +++ b/src/examples/tst_riddles_ott.jl @@ -1,62 +1,46 @@ using Revise -using Basins -using DifferentialEquations +using OrdinaryDiffEq +using DynamicalSystems using Plots -# Equations of motion: E. Ott, et al. I Physica D 76 (1994) 384-410 -function forced_particle!(du, u, p, t) -γ=0.05 ; x̄ = 1.9 ; f₀=2.3 ; ω =3.5 -x₀=1. ; y₀=0.; - - x=u[1]; dx = u[2]; y=u[3]; dy = u[4]; - - du[1] = dx - du[2] = -γ*dx -(-4*x*(1-x^2) + y^2) + f₀*sin(ω*t)*x₀ - du[3] = dy - du[4] = -γ*dy -(2*y*(x+x̄)) + f₀*sin(ω*t)*y₀ -end - - -df = ODEProblem(forced_particle!,rand(4),(0.0,20.0), []) -integ = init(df, alg=Tsit5(); reltol=1e-9, save_everystep=false) - -xres=300 -yres=300 - -# range for forced pend -xg = range(0.,1.2,length=xres) -yg = range(0.,1.2,length=yres) - +ds = System.riddled_basins(zeros(4); γ=0.05, x̄ = 1.9, f₀=2.3, ω =3.5, x₀=1, y₀=0) +integ = init(ds; reltol=1e-9, save_everystep=false) +res=300 ω =3.5 - -# compute basin -#@time bsn = Basins.basin_general_ds(xg, yg, integ; dt=2*pi/ω, idxs=[1,3]) - -#plot(xg,yg,bsn.basin', seriestype=:heatmap) - -function escape_function(u) - - if abs(u[3]) > 20 && abs(u[4])>100 && u[3]*u[4]>0 - @show integ.t - return 2 - elseif abs(u[3]) < 1e-8 && abs(u[4])<1e-9 && u[3]*u[4]<0 - return 1 - end - return 0 -end - -basin_t = zeros(length(xg),length(yg)) -for (k,x) in enumerate(xg), (n,y) in enumerate(yg) - reinit!(integ,[x,0,y,0]) - step!(integ, 2π/ω*10,true) - while escape_function(integ.u) == 0 - step!(integ, 2π/ω,true) - integ.t > 200 && break - end - - basin_t[k,n]=escape_function(integ.u); -end - -using JLD - -@save "riddle.jld" basin_t +xg = range(-2,2,length=res) +yg = range(0.,2,length=res) +default_diffeq = (reltol = 1e-9, alg = Vern9()) +basins, att = basins_of_attraction((xg, yg), df; T = 2π/ω, diffeq = default_diffeq, horizon_limit = 10) + +# +# # compute basin +# #@time bsn = Basins.basin_general_ds(xg, yg, integ; dt=2*pi/ω, idxs=[1,3]) +# +# #plot(xg,yg,bsn.basin', seriestype=:heatmap) +# +# function escape_function(u) +# +# if abs(u[3]) > 20 && abs(u[4])>100 && u[3]*u[4]>0 +# @show integ.t +# return 2 +# elseif abs(u[3]) < 1e-8 && abs(u[4])<1e-9 && u[3]*u[4]<0 +# return 1 +# end +# return 0 +# end +# +# basin_t = zeros(length(xg),length(yg)) +# for (k,x) in enumerate(xg), (n,y) in enumerate(yg) +# reinit!(integ,[x,0,y,0]) +# step!(integ, 2π/ω*10,true) +# while escape_function(integ.u) == 0 +# step!(integ, 2π/ω,true) +# integ.t > 200 && break +# end +# +# basin_t[k,n]=escape_function(integ.u); +# end +# +# using JLD +# +# @save "riddle.jld" basin_t diff --git a/src/straddle.jl b/src/straddle.jl index 933009f..383334b 100755 --- a/src/straddle.jl +++ b/src/straddle.jl @@ -2,18 +2,20 @@ -function bisection_refine!(u_A, u_B, bsn_nfo, integ, basin_A, basin_B, tol) +function bisection_refine!(u_A, u_B, mapper, basin_A, basin_B, tol) + + bsn_nfo = deepcopy(mapper.bsn_nfo) + integ = deepcopy(mapper.integ) # shortcut functions function get_col(u0) - a = get_color_point!(bsn_nfo, integ, u0) - return iseven(a) ? Int(a/2) : Int((a-1)/2) + a = mapper(u0) end function get_col_next(u0) - bsn_nfo.reinit_f!(integ,u0) + bsn_nfo.complete_and_reinit!(integ,u0) bsn_nfo.iter_f!(integ) - a=get_col(bsn_nfo.get_u(integ)) + a = get_col(bsn_nfo.get_projected_state(integ)) return a end @@ -101,12 +103,12 @@ end """ - compute_saddle(bsn_nfo::BasinInfo, integ; max_iter=10) + compute_saddle(grid::Tuple, ds; bas_A = nothing, bas_B = nothing, N = 100, init_tol = 1e-6, kwargs...) The algorithm the saddle that lies in a boundary of the basin of attraction of a dynamical system. When the boundary is fractal, this set is known as the chaotic saddle and is responsible for the transient dynamics of the system. The saddle is computed with the Saddle-Straddle algorithm. It is necessary to define two `generalized basin`, that is we must separate the basin into two sets (see also keyword arguments). -[H. E. Nusse and J. A. Yorke, Dynamics: numerical explorations, Springer, New York, 2012] +[H. E. Nusse and J. A. Yorke, Dynamics: numerical explorations, Springer, New York, 1997] ## Arguments * `bsn_nfo` : structure that holds the information of the basin as well as the map function. This structure is set when the basin is first computed with `basin_stroboscopic_map` or `basin_poincare_map`. @@ -126,15 +128,21 @@ sa,sb = compute_saddle(bsn, integ_df, [1], [2,3],1000) ``` """ -function compute_saddle(integ, bsn_nfo::BasinInfo, bas_A, bas_B; N=100, init_tol = 1e-6) +function compute_saddle(grid::Tuple, ds; bas_A = nothing, bas_B = nothing, attractors = nothing, N = 100, init_tol = 1e-6, kwargs...) + + att = attractors; + if isnothing(attractors) + basins, att = basins_of_attraction(grid, ds; kwargs...) + end + mapper = AttractorMapper(ds; attractors = att, kwargs...) # shortcut functions function get_col(u0) - a = get_color_point!(bsn_nfo, integ, u0) - return iseven(a) ? Int(a/2) : Int((a-1)/2) + a = mapper(u0) end - Na = bsn_nfo.Na + Na = length(att) + # basic check if (Set(bas_A) ∪ Set(bas_B)) != Set(collect(1:Na)) @error "Generalized basins are not well defined" @@ -148,15 +156,16 @@ function compute_saddle(integ, bsn_nfo::BasinInfo, bas_A, bas_B; N=100, init_tol tol = init_tol - xg = bsn_nfo.xg; yg = bsn_nfo.yg; # aliases - + @show att # Set initial condition near attractors of the selected basins. - u_A = u_B = [0., 0.] - v = bsn_nfo.attractors[bas_A[1]*2] + u_A = u_B = zeros(length(att[1])) + v = att[bas_A[1]] u_A = v[1] - v = bsn_nfo.attractors[bas_B[1]*2] + v = att[bas_B[1]] u_B = v[1] + bsn_nfo = deepcopy(mapper.bsn_nfo) + integ = deepcopy(mapper.integ) saddle_series_A = Array{typeof(u_A),1}(undef, N) saddle_series_B = Array{typeof(u_A),1}(undef, N) @@ -167,17 +176,17 @@ function compute_saddle(integ, bsn_nfo::BasinInfo, bas_A, bas_B; N=100, init_tol while k<= N # Initial bisection - u_A_r,u_B_r = bisection_refine!(u_A, u_B, bsn_nfo, integ, bas_A, bas_B, tol) + u_A_r,u_B_r = bisection_refine!(u_A, u_B, mapper, bas_A, bas_B, tol) dist = norm(u_A_r-u_B_r) - bsn_nfo.reinit_f!(integ,u_A_r) + bsn_nfo.complete_and_reinit!(integ, u_A_r) bsn_nfo.iter_f!(integ) - u_A_it = deepcopy(bsn_nfo.get_u(integ)) # for some reason I have to deepcopy this vector + u_A_it = bsn_nfo.get_projected_state(integ) # for some reason I have to deepcopy this vector ca=get_col(u_A_it) - bsn_nfo.reinit_f!(integ,u_B_r) + bsn_nfo.complete_and_reinit!(integ,u_B_r) bsn_nfo.iter_f!(integ) - u_B_it = deepcopy(bsn_nfo.get_u(integ)) + u_B_it = bsn_nfo.get_projected_state(integ) cb=get_col(u_B_it) diff --git a/src/uncertain_dim.jl b/src/uncertain_dim.jl index d15efc2..4b08b0a 100644 --- a/src/uncertain_dim.jl +++ b/src/uncertain_dim.jl @@ -1,331 +1,2 @@ -""" - box_counting_dim(xg, yg, basin; kwargs...) -This function compute the box-counting dimension of the boundary. It is related to the uncertainty dimension. -[C. Grebogi, S. W. McDonald, E. Ott, J. A. Yorke, Final state sensitivity: An obstruction to predictability, Physics Letters A, 99, 9, 1983] - -## Arguments -* `basin` : the matrix containing the information of the basin. -* `xg`, `yg` : 1-dim range vector that defines the grid of the initial conditions to test. - -## Keyword arguments -* `kwargs...` these arguments are passed to `generalized_dim` see `ChaosTools.jl` for the docs. - -""" -function box_counting_dim(xg, yg, basin; kwargs...) - # get points coordinates in the boudary - bnd = get_boundary_filt(basin) - I1 = findall(bnd .== 1); - v = Dataset{2,eltype(xg)}() - for ind in I1; push!(v, [xg[ind[1]] , yg[ind[2]]]); end; - return generalized_dim(v; q=0, kwargs...) -end - -function box_counting_dim(xg, yg, bsn::BasinInfo; kwargs...) - ind = findall(iseven.(bsn.basin) .== true) - # Not sure it is necessary to deepcopy - basin_test = deepcopy(bsn.basin) - # Remove attractors from the picture - for k in ind; basin_test[k] = basin_test[k]+1; end - return box_counting_dim(xg, yg, basin_test; kwargs...) -end - - - - -""" - uncertainty_exponent(bsn::BasinInfo, integ; precision=1e-3) -> α,e,ε,f_ε -This function estimates the uncertainty exponent of the boundary. It is related to the uncertainty dimension. - -[C. Grebogi, S. W. McDonald, E. Ott, J. A. Yorke, Final state sensitivity: An obstruction to predictability, Physics Letters A, 99, 9, 1983] - -## Arguments -* `bsn` : structure that holds the information of the basin. -* `integ` : handle to the iterator of the dynamical system. - -## Keyword arguments -* `precision` variance of the estimator of the uncertainty function. - -""" -function uncertainty_exponent(bsn::BasinInfo, integ; precision=1e-3) - xg = bsn.xg; yg = bsn.yg; - nx=length(xg) - ny=length(yg) - xi = xg[1]; xf = xg[end]; - grid_res_x=xg[2]-xg[1] - - num_step=10 - N_u = zeros(1,num_step) # number of uncertain box - N = zeros(1,num_step) # number of boxes - ε = zeros(1,num_step) # resolution - - # First we compute a crude estimate to have an idea at which resolution - D,e,f = static_estimate(xg,yg,bsn); - println("First estimate: α=",D) - - # resolution in log scale - # estimate the minimum resolution to estimate a f_ε ∼ 10^-3 - min_ε = -3/D[1] - @show min_ε = max(min_ε,-7) # limit resolution to 10^-7 to avoid numerical inconsistencies. - max_ε = log10(grid_res_x*5) # max radius of the ball is roughly 5 pixels. - - r_ε = 10 .^ range(min_ε,max_ε,length=num_step) - - for (k,eps) in enumerate(r_ε) - Nb=0; Nu=0; μ=0; σ²=0; M₂=0; - completed = 0; - # Find uncertain boxes - while completed == 0 - k1 = rand(1:nx) - k2 = rand(1:ny) - - c1 = bsn.basin[k1,k2] - c2 = get_color_point!(bsn, integ, [xg[k1]+eps,yg[k2]]) - c3 = get_color_point!(bsn, integ, [xg[k1]-eps,yg[k2]]) - - if length(unique([c1,Int(c2),Int(c3)]))>1 - Nu = Nu + 1 - end - Nb += 1 - - # Welford's online average estimation and variance of the estimator - M₂ = wel_var(M₂, μ, Nu/Nb, Nb) - μ = wel_mean(μ, Nu/Nb, Nb) - σ² = M₂/Nb - # If the process is binomial, the variance of the estimator is Nu/Nb·(1-Nu/Nb)/Nb (p·q/N) - # simulations matches - #push!(sig,Nu/Nb*(1-Nu/Nb)/Nb) - - # Stopping criterion: variance of the estimator of the mean bellow 1e-3 - if Nu > 50 && σ² < precision - completed = 1 - @show Nu,Nb,σ² - end - - if Nu < 10 && Nb>10000 - # skip this value, we don't find the boundary - # corresponds to roughly f_ε ∼ 10^-4 - @warn "Resolution may be too small for this basin, skip value" - completed = 1 - Nu = 0 - end - - end - N_u[k]=Nu - N[k]=Nb - ε[k]=eps - end - - # uncertain function - f_ε = N_u./N - - # remove zeros in case there are: - ind = f_ε .> 0 - f_ε = f_ε[ind] - ε = ε[ind] - - # @show i1,D=linear_region(vec(log10.(ε)), vec(log10.(f_ε))) - # Dϵ=0. - # @show i1 - # get exponent from a linear region - @. model(x, p) = p[1]*x+p[2] - fit = curve_fit(model, vec(log10.(ε)), vec(log10.(f_ε)), [2., 2.]) - D = coef(fit) - Dϵ = estimate_errors(fit) - return D[1], Dϵ[1], vec(log10.(ε)),vec(log10.(f_ε)) -end - - - -function wel_var(M₂, μ, xₙ, n) - - μ₂ = μ + (xₙ - μ)/n - M₂ = M₂ + (xₙ - μ)*(xₙ - μ₂) - return M₂ -end - -function wel_mean(μ, xₙ, n) - return μ + (xₙ - μ)/n -end - - -function static_estimate(xg,yg,bsn; precision=1e-4) - - xg = bsn.xg; yg = bsn.yg; - y_grid_res=yg[2]-yg[1] - nx=length(xg) - ny=length(yg) - - num_step=10 - N_u = zeros(1,num_step) # number of uncertain box - N = zeros(1,num_step) # number of boxes - ε = zeros(1,num_step) # resolution - - # resolution in pixels - min_ε = 1; - max_ε = 10; - - r_ε = min_ε:max_ε - #@show r_ε = 10 .^ range(log10(min_ε),log10(max_ε),length=num_step) - for (k,eps) in enumerate(r_ε) - Nb=0; Nu=0; μ=0; σ²=0; M₂=0; - completed = 0; - # Find uncertain boxes - while completed == 0 - kx = rand(1:nx) - ky = rand(ceil(Int64,eps+1):floor(Int64,ny-eps)) - - indy = range(ky-eps,ky+eps,step=1) - c = [bsn.basin[kx,ky] for ky in indy] - - if length(unique(c))>1 - Nu = Nu + 1 - end - Nb += 1 - - # Welford's online average estimation and variance of the estimator - M₂ = wel_var(M₂, μ, Nu/Nb, Nb) - μ = wel_mean(μ, Nu/Nb, Nb) - σ² = M₂/Nb - - # Stopping criterion: variance of the estimator of the mean bellow 1e-3 - if Nu > 50 && σ² < precision - completed = 1 - #@show Nu,Nb,σ² - end - - end - N_u[k]=Nu - N[k]=Nb - ε[k]=eps*y_grid_res - end - - # uncertain function - f_ε = N_u./N - - # remove zeros in case there are: - ind = f_ε .> 0. - f_ε = f_ε[ind] - ε = ε[ind] - # get exponent - @. model(x, p) = p[1]*x+p[2] - fit = curve_fit(model, vec(log10.(ε)), vec(log10.(f_ε)), [2., 2.]) - D = coef(fit) - @show estimate_errors(fit) - #D = linear_region(vec(log10.(ε)), vec(log10.(f_ε))) - return D[1],vec(log10.(ε)), vec(log10.(f_ε)) - -end - - - -function ball_estimate(xg,yg,bsn; precision = 1e-4) - - xg = bsn.xg; yg = bsn.yg; - xf=xg[end];xi=xg[1];yf=yg[end];yi=yg[1]; - y_grid_res=yg[2]-yg[1] - nx=length(xg) - ny=length(yg) - - num_step=10 - N_u = zeros(1,num_step) # number of uncertain box - N = zeros(1,num_step) # number of boxes - ε = zeros(1,num_step) # resolution - - # resolution in pixels - min_ε = 1; - max_ε = 10; - @show r_ε = 10 .^ range(log10(min_ε),log10(max_ε),length=num_step) - #r_ε = (min_ε:max_ε)*y_grid_res - - for (k,eps) in enumerate(r_ε) - Nb=0; Nu=0; μ=0; σ²=0; M₂=0; - completed = 0; - c_ind = get_ind_in_circle(eps) - # Find uncertain boxes - while completed == 0 - - kx = rand()*((nx-eps)-(1+eps))+1+eps - ky = rand()*((ny-eps)-(1+eps))+1+eps - - # I = get_indices_in_range(kx, ky, eps, c_ind) - # c = [ bsn.basin[n[1],n[2]] for n in eachrow(I)] - - Ix,Iy = get_indices_in_range(kx, ky, eps, c_ind) - c = [ bsn.basin[n,m] for n in Ix, m in Iy] - - if length(unique(c))>1 - Nu = Nu + 1 - end - Nb += 1 - - # Welford's online average estimation and variance of the estimator - M₂ = wel_var(M₂, μ, Nu/Nb, Nb) - μ = wel_mean(μ, Nu/Nb, Nb) - σ² = M₂/Nb - - # Stopping criterion: variance of the estimator of the mean bellow 1e-3 - if Nu > 50 && σ² < precision - completed = 1 - @show Nu,Nb,σ² - end - - end - N_u[k]=Nu - N[k]=Nb - ε[k]=eps*y_grid_res - end - - # uncertain function - f_ε = N_u./N - - # remove zeros in case there are: - ind = f_ε .> 0. - f_ε = f_ε[ind] - ε = ε[ind] - # get exponent - @. model(x, p) = p[1]*x+p[2] - fit = curve_fit(model, vec(log10.(ε)), vec(log10.(f_ε)), [2., 2.]) - D = coef(fit) - @show estimate_errors(fit) - return D[1],vec(log10.(ε)), vec(log10.(f_ε)) - -end - - - -function get_indices_in_range(kx, ky, eps,c_ind) - - # center and scale - - #indx = ceil.(Int64,(kx-eps):1.:(kx+eps)) - #indy = ceil.(Int64,(ky-eps):1.:(ky+eps)) - indx = range(ceil.(Int64,kx-eps),ceil(Int64,kx+eps)-1,step=1) - indy = range(ceil.(Int64,ky-eps),ceil(Int64,ky+eps)-1,step=1) - - - #I = deepcopy(c_ind) - - - #Find indices in a circle of radius eps - #I[:,1] .= I[:,1] .+ ceil(Int64,kx) - #I[:,2] .= I[:,2] .+ ceil(Int64,ky) - #@show [[n[1],n[2]] for n in I] - - return indx,indy - #return I -end - - -function get_ind_in_circle(r) - #Find indices in a circle of radius eps - I = Array{Int64,1}[] - r_l = ceil(Int16,r) - for n in -r_l:r_l, m in -r_l:r_l - if n^2+m^2 ≤ r^2 - push!(I, [n; m]) - end - end - v= hcat(I...)' - return v -end +# Functions moved to Chaostools. diff --git a/src/wada_detection.jl b/src/wada_detection.jl index 6e82fae..be7bd30 100644 --- a/src/wada_detection.jl +++ b/src/wada_detection.jl @@ -17,7 +17,7 @@ function get_boundary_filt(basin) return res end -function get_list(basin,xg ,yg) +function get_list(basin) num_att = length(unique(basin)) @@ -29,7 +29,7 @@ function get_list(basin,xg ,yg) # original boundary bnd = get_boundary_filt(basin) I1=findall(bnd .== 1); - for ind in I1; push!(v_list , [xg[ind[1]]; yg[ind[2]]]); end + for ind in I1; push!(v_list , [ind[1]; ind[2]]); end #@show v_list push!(v_set,v_list) @@ -45,7 +45,7 @@ function get_list(basin,xg ,yg) I1=findall(bnd .== 1); # Get coordinates lists from matrix #push!(v_list , hcat([[xg[ind[1]]; yg[ind[2]]] for ind in I1 ]...)) - for ind in I1; push!(v_list , [xg[ind[1]]; yg[ind[2]]]); end + for ind in I1; push!(v_list , [ind[1]; ind[2]]); end push!(v_set,v_list) end @@ -56,39 +56,30 @@ end """ - detect_wada_merge_method(xg,yg, basin) -The algorithm gives the maximum and minimum Haussdorff distances between merged basins. These two distances can help to decide if the basin has the Wada property. + detect_wada_merge_method(basins) -> mx_dist, mn_dist +The algorithm gives the maximum and minimum Haussdorff distances between combinations +of merged basins in unit of pixels. These two distances can help to decide if +the basin has the Wada property. + +The maximum Haussdorff distance can be interpreted as the minimum Fattening parameter +of the boundaries you need to match all basins. See Ref. [A. Daza, A. Wagemakers and M. A. F. Sanjuán, Ascertaining when a basin is Wada: the merging method, Sci. Rep., 8, 9954 (2018)] ## Arguments -* `basin` : the matrix containing the information of the basin. -* `xg`, `yg` : 1-dim range vector that defines the grid of the initial conditions to test. +* `basins` : the matrix containing the information of the basin. ## Example ``` -max_dist,min_dist = detect_wada_merge_method(basin,xg,yg) -# grid resolution -epsilon = xg[2]-xg[1] -# if dmax is large then this is not Wada -@show dmax = max_dist/epsilon -@show dmin = min_dist/epsilon +max_dist,min_dist = detect_wada_merge_method(basins) ``` """ -function detect_wada_merge_method(xg,yg,bsn::BasinInfo) - ind = findall(iseven.(bsn.basin) .== true) - basin_test = deepcopy(bsn.basin) - for k in ind; basin_test[k] =basin_test[k]+1; end - return detect_wada_merge_method(xg,yg,basin_test) -end - -function detect_wada_merge_method(xg,yg,basin) +function detect_wada_merge_method(basins) - num_att = length(unique(basin)) - - v_list=get_list(basin, xg, yg) + num_att = length(unique(basins)) + v_list=get_list(basins) # compute distances using combbinations of 2 elements from a collection ind = combinations(1:num_att,2) @@ -131,7 +122,7 @@ end mutable struct ds_info{I} - bsn_nfo::BasinInfo # basin info for BA routine + bsn_nfo # basin info for BA routine integ::I # integrator end @@ -144,7 +135,7 @@ end """ detect_wada_grid_method(integ, bsn_nfo::BasinInfo; max_iter=10) The algorithm test for Wada basin in a dynamical system. It uses the dynamical system to look if all the atractors are represented in the boundary. - +Warning: only works for 2D (for now) [A. Daza, A. Wagemakers, M. A. F. Sanjuán and J. A. Yorke, Testing for Basins of Wada, Sci. Rep., 5, 16579 (2015)] ## Arguments @@ -155,28 +146,31 @@ The algorithm test for Wada basin in a dynamical system. It uses the dynamical s * `max_iter` : set the maximum depth of subdivisions to look for an atractor. The number of points doubles at each step. """ -function detect_wada_grid_method(integ, bsn_nfo::BasinInfo; max_iter=10) +function detect_wada_grid_method(grid::Tuple, ds; max_iter=10, attractors = nothing, basins = nothing, kwargs...) - ds_nfo = ds_info(bsn_nfo, integ) - num_att = bsn_nfo.Na - if findfirst(x->x==-1, bsn_nfo.basin) != nothing + if isnothing(attractors) || isnothing(basins) + basins, attractors = basins_of_attraction(grid, ds; kwargs...) + end + att = attractors; + mapper = AttractorMapper(ds; attractors = att, kwargs...) + + #ds_nfo = ds_info(bsn_nfo, integ) + num_att = length(att) + + if findfirst(x->x==-1, basins) != nothing @error "The basin contains escapes or undefined attractors, cannot test for Wada" return nothing end - + xg = grid[1] + yg = grid[2] # helper function to obtain coordinates - index_to_coord(p) = [bsn_nfo.xg[p[1]], bsn_nfo.yg[p[2]]] - - # We remove the atractors for this computation - tmp_bsn = deepcopy(bsn_nfo.basin) - ind = findall(iseven.(tmp_bsn) .== true) - [tmp_bsn[k] = tmp_bsn[k]+1 for k in ind ] + index_to_coord(p) = [xg[p[1]], yg[p[2]]] # obtain points in the boundary - bnd = get_boundary_filt(tmp_bsn) + bnd = get_boundary_filt(basins) p1_ind = findall(bnd .> 0) # initialize empty array of indices and collection of empty sets of colors @@ -186,7 +180,7 @@ function detect_wada_grid_method(integ, bsn_nfo::BasinInfo; max_iter=10) # Initialize matrices (step 1) for (k,p1) in enumerate(p1_ind) - p2, nbgs = get_neighbor_and_colors(tmp_bsn, [p1[1], p1[2]]) + p2, nbgs = get_neighbor_and_colors(basins, [p1[1], p1[2]]) if length(nbgs) > 1 # keep track of different colors and neighbor point push!(clr_mat[k],nbgs...) @@ -203,7 +197,7 @@ function detect_wada_grid_method(integ, bsn_nfo::BasinInfo; max_iter=10) pc1 = index_to_coord(p1_ind[k]) pc2 = index_to_coord(p2_ind[k]) # update number of colors - clr_mat[k]=divide_and_test_W(ds_nfo, pc1, pc2, n, clr_mat[k], num_att) + clr_mat[k]=divide_and_test_W(mapper, pc1, pc2, n, clr_mat[k], num_att) # update W matrix W[length(clr_mat[k]),n] += 1 end @@ -220,7 +214,7 @@ function detect_wada_grid_method(integ, bsn_nfo::BasinInfo; max_iter=10) end -function get_neighbor_and_colors(basin, p) +function get_neighbor_and_colors(basins, p) radius = 1 n=p[1] m=p[2] @@ -229,9 +223,9 @@ function get_neighbor_and_colors(basin, p) # check neihbors and collect basin colors for k=n-radius:n+radius, l=m-radius:m+radius try - push!(v,basin[k,l]) + push!(v,basins[k,l]) if k != n || l != m - if basin[n,m] != basin[k,l] + if basins[n,m] != basins[k,l] p2=CartesianIndex(k,l) end end @@ -244,7 +238,7 @@ end -function divide_and_test_W(ds_nfo, p1, p2, nstep, clrs, Na) +function divide_and_test_W(mapper, p1, p2, nstep, clrs, Na) if length(clrs) == Na return clrs @@ -261,7 +255,7 @@ function divide_and_test_W(ds_nfo, p1, p2, nstep, clrs, Na) # get colors and update color set for this box! for pnt in pnt_to_test - clr = get_color_point!(ds_nfo.bsn_nfo, ds_nfo.integ, pnt) + clr = mapper(pnt) push!(clrs, clr) if length(clrs) == Na break diff --git a/test/runtests.jl b/test/runtests.jl index 84ff5f0..ddbe4f0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,430 +1,7 @@ using Basins using DynamicalSystems using Test -using DifferentialEquations @testset "All tests" begin -@testset "Test basin_discrete_map" begin - v_bnd=[1.2727272727272727 1.1111111111111112; -1.6767676767676767 1.7575757575757576; --0.2222222222222222 0.30303030303030304; --0.6666666666666666 -1.595959595959596; --0.26262626262626265 0.42424242424242425; -1.1515151515151516 1.1919191919191918; -0.5454545454545454 -1.797979797979798; --1.4747474747474747 0.8282828282828283; -1.2727272727272727 -0.2222222222222222; -0.9494949494949495 1.5151515151515151; --1.2323232323232323 -0.020202020202020204; -1.1111111111111112 -0.5454545454545454; -0.5858585858585859 -1.878787878787879; -0.7878787878787878 -1.3131313131313131; --1.4747474747474747 0.9494949494949495; -1.7171717171717171 1.797979797979798; -0.5858585858585859 0.9090909090909091; -1.7171717171717171 1.7575757575757576; --1.1111111111111112 2.0; --0.7070707070707071 -1.5555555555555556; -1.0707070707070707 -0.6666666666666666; -1.0707070707070707 0.3838383838383838; --1.2727272727272727 0.020202020202020204; -1.1515151515151516 1.878787878787879; --1.393939393939394 0.7474747474747475; --0.9090909090909091 1.4747474747474747; -0.8282828282828283 -1.3131313131313131; -1.3535353535353536 1.3131313131313131; --0.8686868686868687 1.3131313131313131; -1.5151515151515151 1.0303030303030303; -0.9494949494949495 -0.9494949494949495; -1.5555555555555556 1.0303030303030303; -0.3434343434343434 0.46464646464646464; --0.5454545454545454 0.6666666666666666; -1.6767676767676767 1.7575757575757576; -1.3131313131313131 1.6363636363636365; -1.4343434343434343 1.4343434343434343; --1.2323232323232323 0.020202020202020204; -0.9090909090909091 -1.2323232323232323; -0.7474747474747475 1.1515151515151516; --1.6363636363636365 1.5555555555555556; -1.4747474747474747 1.5555555555555556; -1.2727272727272727 0.9090909090909091; --0.3434343434343434 0.42424242424242425; --0.7070707070707071 -1.4343434343434343; --1.2727272727272727 0.26262626262626265; -1.0707070707070707 1.7575757575757576; -0.46464646464646464 -2.0; -1.4343434343434343 1.5555555555555556; -0.30303030303030304 0.42424242424242425; --0.8282828282828283 -1.2727272727272727; --1.5151515151515151 1.0303030303030303; -1.3535353535353536 0.3838383838383838; --0.8686868686868687 1.3535353535353536; -1.4343434343434343 0.7070707070707071; --0.98989898989899 1.5555555555555556; -1.4747474747474747 0.6262626262626263; -1.4343434343434343 1.9595959595959596; --0.6666666666666666 0.7878787878787878; -1.1111111111111112 -0.7070707070707071; --1.2323232323232323 0.10101010101010101; --1.393939393939394 0.5050505050505051; -1.3535353535353536 1.6767676767676767; -0.8686868686868687 1.3535353535353536; -0.42424242424242425 0.5050505050505051; -0.8686868686868687 1.393939393939394; --0.5050505050505051 0.6262626262626263; -1.0707070707070707 1.6363636363636365; -1.3131313131313131 1.0303030303030303; -1.0303030303030303 -0.9494949494949495; -0.5454545454545454 -1.9191919191919191; --1.4747474747474747 0.98989898989899; -1.1919191919191918 0.6262626262626263; -1.0707070707070707 0.7474747474747475; -0.6262626262626263 -1.7171717171717171; -0.42424242424242425 0.5454545454545454; --1.2727272727272727 0.18181818181818182; --0.9494949494949495 1.595959595959596; --0.7878787878787878 -1.2323232323232323; -1.1515151515151516 0.5050505050505051; -0.9494949494949495 1.5151515151515151; -1.7575757575757576 1.7575757575757576; -0.46464646464646464 0.7070707070707071; -1.0303030303030303 2.0; --0.5858585858585859 0.7070707070707071; -1.1111111111111112 1.0303030303030303; -1.1919191919191918 -0.5050505050505051; -1.5151515151515151 0.7878787878787878; --0.5050505050505051 -1.8383838383838385; -1.1919191919191918 -0.2222222222222222; --0.9090909090909091 1.3535353535353536; -1.1111111111111112 0.9090909090909091; --1.595959595959596 1.595959595959596; -0.7070707070707071 -1.5151515151515151; -1.1111111111111112 1.878787878787879; --1.0707070707070707 -0.5050505050505051; --1.2727272727272727 0.020202020202020204; -1.0303030303030303 1.7575757575757576; --0.3838383838383838 -1.9191919191919191; -0.7474747474747475 -1.5151515151515151; -1.1919191919191918 -0.3838383838383838; -1.3131313131313131 0.9494949494949495; --1.1919191919191918 -0.1414141414141414; --0.9494949494949495 1.5555555555555556; -1.1111111111111112 0.7474747474747475; --0.10101010101010101 0.3434343434343434; -1.0707070707070707 1.9595959595959596; --1.0707070707070707 1.6767676767676767; -0.42424242424242425 0.5050505050505051; -1.1919191919191918 -0.3838383838383838; -0.8686868686868687 1.5555555555555556; -1.4343434343434343 0.5050505050505051; --0.6666666666666666 0.9090909090909091; -1.3131313131313131 1.6767676767676767; -1.0303030303030303 1.9595959595959596; -1.3535353535353536 1.9191919191919191; -1.4343434343434343 2.0; -1.1111111111111112 1.0303030303030303; -0.6262626262626263 -1.6767676767676767; -1.0707070707070707 1.5151515151515151; --0.9494949494949495 -0.9494949494949495; -0.18181818181818182 0.3434343434343434; --0.5858585858585859 0.7878787878787878; -1.0707070707070707 0.5858585858585859; -1.1919191919191918 0.6262626262626263; -0.6666666666666666 -1.7575757575757576; -1.3535353535353536 1.9595959595959596; --0.42424242424242425 0.5858585858585859; -1.1919191919191918 -0.3434343434343434; --1.6363636363636365 1.878787878787879; --1.1919191919191918 -0.1414141414141414; --0.9494949494949495 -0.9494949494949495; --1.1515151515151516 -0.3838383838383838; -1.1919191919191918 1.0707070707070707; -0.98989898989899 1.595959595959596; -1.0707070707070707 0.42424242424242425; --0.7070707070707071 -1.4747474747474747; --1.1919191919191918 0.020202020202020204; -0.5454545454545454 0.6666666666666666; -1.797979797979798 1.9191919191919191; -0.3434343434343434 0.5050505050505051; -1.7575757575757576 1.9595959595959596; -1.7575757575757576 1.878787878787879; --0.8686868686868687 1.2727272727272727; -1.3131313131313131 -0.06060606060606061; -1.1111111111111112 0.5050505050505051; -1.595959595959596 1.1919191919191918; --0.7878787878787878 1.1111111111111112; -1.1515151515151516 1.6767676767676767; --1.0707070707070707 1.9595959595959596; --0.9090909090909091 -0.98989898989899; --1.1919191919191918 -0.06060606060606061; -0.5050505050505051 0.6666666666666666; --0.46464646464646464 -1.878787878787879; --0.5050505050505051 0.6262626262626263; --0.5858585858585859 -1.6767676767676767; -0.9494949494949495 1.5555555555555556; -1.3535353535353536 1.7575757575757576; --0.7474747474747475 -1.4343434343434343; -1.1515151515151516 1.9191919191919191; --0.6262626262626263 -1.5151515151515151; -0.020202020202020204 0.30303030303030304; -0.98989898989899 -0.8686868686868687; -1.0303030303030303 0.30303030303030304; -1.393939393939394 1.2323232323232323; --1.6767676767676767 1.7171717171717171; -1.3131313131313131 1.2323232323232323; -0.98989898989899 -0.98989898989899; --1.6363636363636365 1.5151515151515151; --0.7070707070707071 -1.393939393939394; -1.3535353535353536 0.2222222222222222; -0.7474747474747475 1.0303030303030303; -1.1919191919191918 1.1919191919191918; -0.98989898989899 -1.0707070707070707; --0.6666666666666666 -1.5555555555555556; -0.5050505050505051 0.7070707070707071; --0.9090909090909091 -1.0707070707070707; -0.7070707070707071 -1.6767676767676767; -0.9494949494949495 0.30303030303030304; -1.0303030303030303 0.7474747474747475; -1.2727272727272727 0.10101010101010101; --0.9494949494949495 -1.0303030303030303; -0.8686868686868687 1.5151515151515151; -1.393939393939394 0.3434343434343434; -1.393939393939394 1.9595959595959596; -1.5555555555555556 1.9595959595959596; -1.4747474747474747 1.5151515151515151; --0.7878787878787878 -1.3535353535353536; -0.6666666666666666 -1.595959595959596; --1.4343434343434343 0.8282828282828283; -1.0303030303030303 1.7171717171717171; -1.0303030303030303 1.8383838383838385; -0.1414141414141414 0.30303030303030304; --1.0303030303030303 -0.5454545454545454; --0.42424242424242425 -1.9595959595959596; -0.9494949494949495 0.3838383838383838; -0.98989898989899 0.46464646464646464; -1.1919191919191918 1.3131313131313131; --1.0707070707070707 -0.6262626262626263; -1.6363636363636365 1.1919191919191918]; - ds = Systems.henon(zeros(2); a = 1.4, b = 0.3) - integ_df = integrator(ds) - xg = range(-2.,2.,length=100) - yg = range(-2.,2.,length=100) - bsn = Basins.basins_map2D(xg, yg, integ_df) - - # Compute boundary - ind = findall(iseven.(bsn.basin) .== true) - basin_test = bsn.basin - for k in ind; basin_test[k] =basin_test[k]+1; end - bnd = Basins.get_boundary_filt(basin_test) - I1 = findall(bnd .== 1); - v = hcat([[xg[ind[1]]; yg[ind[2]]] for ind in I1 ]...) - - hd = Basins.haussdorff_dist(v,v_bnd'); - @test (hd < 0.3) - - #@test abs(count(bsn_nfo.basin .== 3) - 4127) < 10 - #@test abs(count(bsn_nfo.basin .== -1) - 5730) < 10 -end - - -@testset "Test basin_entropy" begin - ω=1.; F = 0.2 - ds =Systems.duffing([0.1, 0.25]; ω = ω, f = F, d = 0.15, β = -1) - integ_df = integrator(ds; alg=Tsit5(), reltol=1e-8, save_everystep=false) - xg = range(-2.2,2.2,length=150) - yg = range(-2.2,2.2,length=150) - bsn = Basins.basins_map2D(xg, yg, integ_df; T=2*pi/ω) - Sb,Sbb = basin_entropy(bsn.basin; eps_x=20, eps_y=20) - @test (abs(trunc(Sb;digits=3) - 0.663) < 0.001) - @test (abs(trunc(Sbb;digits=3) - 0.663) < 0.001) -end - -@testset "Test Wada detection" begin - # Equations of motion: - function forced_pendulum!(du, u, p, t) - d = p[1]; F = p[2]; omega = p[3] - du[1] = u[2] - du[2] = -d*u[2] - sin(u[1])+ F*cos(omega*t) - end - # We have to define a callback to wrap the phase in [-π,π] - function affect!(integrator) - if integrator.u[1] < 0 - integrator.u[1] += 2*π - else - integrator.u[1] -= 2*π - end - end - condition(u,t,integrator) = (integrator.u[1] < -π || integrator.u[1] > π) - cb = DiscreteCallback(condition,affect!) - - F = 1.66; ω = 1.; d=0.2 - df = ODEProblem(forced_pendulum!,rand(2),(0.0,20.0), [d, F, ω]) - integ = init(df, alg=AutoTsit5(Rosenbrock23()); reltol=1e-9, save_everystep=false, callback=cb) - xg = range(-pi,pi,length=100); yg = range(-2.,4.,length=100) - bsn = Basins.basins_map2D(xg, yg, integ; T=2*pi/ω) - - # Wada merge Haussdorff distances - # First remove attractors - max_dist,min_dist = detect_wada_merge_method(xg,yg,bsn) - epsilon = xg[2]-xg[1] - dmax = max_dist/epsilon - dmin = min_dist/epsilon - @test (round(dmax) == 6.) - @test (round(dmin) == 2.) - - W = detect_wada_grid_method(integ, bsn; max_iter=5) - @test abs(trunc(W[3];digits=2)) ≥ 0.9 - -end - - -@testset "Test basin_stability" begin - ω=1.; F = 0.2 - ds =Systems.duffing([0.1, 0.25]; ω = ω, f = F, d = 0.15, β = -1) - integ_df = integrator(ds; alg=Tsit5(), reltol=1e-8, save_everystep=false) - xg = range(-2.2,2.2,length=150) - yg = range(-2.2,2.2,length=150) - bsn = Basins.basins_map2D(xg, yg, integ_df; T=2*pi/ω) - bs = basin_stability(bsn.basin) - @test abs((trunc(bs[1];digits=3) - 0.492)) < 0.005 - @test abs((trunc(bs[2];digits=3) - 0.507)) < 0.005 -end - - -@testset "Test box_counting_dimension" begin - ω=1.; F = 0.2 - ds =Systems.duffing([0.1, 0.25]; ω = ω, f = F, d = 0.15, β = -1) - integ_df = integrator(ds; alg=Tsit5(), reltol=1e-8, save_everystep=false) - xg = range(-2.2,2.2,length=150) - yg = range(-2.2,2.2,length=150) - bsn = Basins.basins_map2D(xg, yg, integ_df; T=2*pi/ω) - bd = box_counting_dim(xg, yg, bsn) - @test abs((trunc(bd;digits=2) - 1.9))<0.1 -end - -@testset "Test compute saddle straddle" begin - v = [-0.7024272244361562 0.7358943669783878; - 0.6884108286512226 0.2883122496731336; - 0.12248777086662276 0.17266855751576238; - 0.21416001598855067 0.1918853810113003; - -0.13345893435577877 -0.10181462673868892; - -0.24612685358870395 0.40609732169904256; - 0.5443525330629784 0.3637084477432412; - -0.8543699001914774 -0.5504191054878991; - -0.81424974452034 0.3876573128398122; - -0.3759850591885575 0.3029466252309591; - 1.2695931242637941 0.5649763379043485; - 0.12319548967885363 0.1988038193210813; - 0.08993268394324054 0.08763876246458675; - 0.5998194677317075 0.42267589676742806; - -0.7471048204006432 0.4173581434012174; - -0.3371249654066878 0.41666798624746265; - 0.7345532189737 0.3144499987176077; - 0.24865664622871408 0.28838520561505615; - -0.7047154958073086 -0.5213541656840109; - -0.9056118879414881 0.33965597519541424; - -0.4547496618520956 0.21744361580215943; - -0.2921283747093624 0.4214569617780614; - 0.6291128784123545 0.3716326220739068; - -0.6010738837992531 -0.43440268456863923; - -0.9652688621193279 0.24069456477396253; - -0.5501113728424021 0.15418201057091677; - -0.4754568961525241 0.10081275951313767; - -0.46112021536727277 0.09967906030053533; - -0.4391298368793252 0.11974265732897998; - -0.3662749654329448 0.20493692475101746; - -0.4603329215220823 0.6299851440826713; - 0.41313450130864204 0.28437379172110344; - -0.7299720705099347 -0.5197447840442747; - -0.890354901563265 0.3447067895447532; - -0.44687447282979653 0.22139280314316853; - -0.2966913365698902 0.45089092769960826; - 0.5501212667337514 0.3597852144857481; - -0.8200621499512287 -0.5412846282080438; - -0.8349858140921245 0.3758250755769921; - -0.39618264893265825 0.2735219498521972; - -0.9437331293087026 0.7574461204326174; - 0.9293611309012172 -0.15957232746111116; - 1.14599339117127 0.4412985174267645; - 1.2396517289691749 0.5188989178649657; - 0.8959919160009866 0.7090718967680737; - -1.3461532547680735 0.008447716620129747; - -0.6398916478847173 0.6224442163197923; - 0.8276869977268546 0.06600919221071469; - 1.0607570462732214 0.6160653104393968; - -0.20258432736829904 -0.10525350252623214; - -0.3349774815559638 0.19748108608150314; - -0.663004994554249 0.724099207141831; - 0.6512663688101122 0.31465908699556383; - -0.19176348994343262 -0.1154502893224026; - -0.33183913666489595 0.1997476599045356; - -0.7339332858876773 0.7450732722276884; - 0.7126484478736383 0.2607963013619859; - 0.3565945556966931 0.376716105526132; - -1.1424583292686683 -0.5181306982978956; - -0.6691798553091624 0.48524940958270907; - -0.7750640627426556 0.7447635013295555; - 0.7617537229889098 0.18145588387180991; - 0.7785640313451194 0.6364791644564665; - -1.3274874780888364 -0.13525661761177976; - -0.6135443463613546 0.580073323763926; - 0.8897837152389886 -0.06286426899098416; - 1.1440178842431123 0.4993844802023786; - 1.071146877796164 0.6852925880897688; - -1.0130674106682127 -0.5202505944779693; - -0.7194964747274598 0.40984327478863253; - -0.3337755688232476 0.43117917669428807; - 0.6995015529974392 0.3448638592935367; - -0.09956822741216353 -0.02900414846930688; - -0.8926894003760993 0.7683209764135135; - 0.8340505900639127 0.011479529745785512; - 1.0830174110506972 0.5892351694272199; - 0.21204311031642775 0.27421769347596214; - -0.5644430315888005 -0.4442700978076312; - -0.9833207542715674 0.22818523879153269; - -0.5545537845226697 0.1593390081784191; - -0.472231790321022 0.10490889580215614; - -0.44780564466696554 0.11274653830575736; - -0.39238113473288266 0.17168250535822036; - -0.2884018330542508 0.44862621970732874; - 0.5312147137503425 0.35474079716608226; - -0.8467395756009193 -0.5497951927484416; - -0.8188633312213228 0.38589323804591835; - -0.37959560576794266 0.29752006205944304; - 0.4696497379499575 0.5119045223664362; - -1.3699491532902666 0.0061234201083442685; - -0.7483772772691258 0.6806762370674923; - 0.8713539043026773 -0.04321198642007509; - 1.1294232129600181 0.527761159044075; - 0.8977550302213696 0.7036254552586797; - -1.3389562103162267 -0.04653014856135522; - -0.6203974969982126 0.6005519054452689; - 0.8447656866536019 0.030093293866878423; - 1.096993725067534 0.5832862071649934; - 0.3206215383741723 0.3693331736195564; - -1.1173678907969715 -0.5383392012987194; - -0.6806544022816443 0.47871124803080173; - -0.626554882949195 0.6976130759887653; - 0.6516324418974114 0.31993125496392566; - -0.2171788716252682 -0.13785785370964124; - -0.43298341824677933 0.08012326524258628; - -0.42802570899573394 0.12623349965496122; - -0.33932641866474456 0.2438970723609983; - -1.2365333502089975 -0.2700884633323442]; - - v = v'; - ω=1.; F = 0.2 - ds =Systems.duffing([0.1, 0.25]; ω = ω, f = F, d = 0.15, β = -1) - integ_df = integrator(ds; alg=AutoTsit5(Rosenbrock23()), reltol=1e-8, abstol=1e-8, save_everystep=false) - xg = range(-2.2,2.2,length=200) - yg = range(-2.2,2.2,length=200) - @time bsn = Basins.basins_map2D(xg, yg, integ_df; T=2*pi/ω) - sa,sb = compute_saddle(integ_df, bsn, [1], [2]; N=100) - s = hcat(sa...) - hd = Basins.haussdorff_dist(s,v) - @test (hd < 0.4) - -end - end