Usage

The base of the set oriented framework are BoxSets (as a discretization of state space) and BoxMaps (as a discretization of a map). In the following, we will have a closer look at the two concepts and other useful things to know when using GAIO.jl.

Box

To create a Box with center c$\in\mathbb{R}^d$ and radius r$\in\mathbb{R}^d$, type

julia> using GAIO
julia> c, r = (0.5, 0.5), (0.5, 0.5)((0.5, 0.5), (0.5, 0.5))
julia> box = Box(c, r)Box{2, Float64}:
   center: [0.5, 0.5]
   radius: [0.5, 0.5]

This creates the set $box = [c_1 - r_1, c_1 + r_1) \times \ldots \times [c_d - r_d, c_d + r_d)$. Conversely, one can get back the vectors c and r from a box $box$ by

julia> c, r = boxBox{2, Float64}:
   center: [0.5, 0.5]
   radius: [0.5, 0.5]

BoxGrid

Most algorithms in GAIO.jl revolve around a partition $\scr P$ of some domain $X\subset\mathbb{R}^d$ into small boxes. Often, the domain X is a Box and the partition 𝒫 is a grid of boxes on X. To create an $n_1 \times \ldots \times n_d$-element grid of boxes on a box X, pass X and the tuple n = (n_1,...,n_d) to the function BoxGrid

julia> X = box     # domainBox{2, Float64}:
   center: [0.5, 0.5]
   radius: [0.5, 0.5]
julia> n = (4, 2)(4, 2)
julia> 𝒫 = BoxGrid(X, n)4 x 2 - element BoxGrid

julia> x = (0.2, 0.1)(0.2, 0.1)
julia> key = point_to_key(𝒫, x)    # x is some point in the domain Q(1, 1)
julia> box = key_to_box(𝒫, key)    # cover the point x with a box from PBox{2, Float64}:
   center: [0.125, 0.25]
   radius: [0.125, 0.25]
julia> box = point_to_box(𝒫, x)    # performs both above functionsBox{2, Float64}:
   center: [0.125, 0.25]
   radius: [0.125, 0.25]

BoxTree

For partitions of X into variably sized boxes, one can use BoxTree:

julia> 𝒫_tree = BoxTree(X)BoxTree of depth 1

A BoxTree uses a binary tree structure to store a box partition of the domain. One can refine a partition by bisecting all boxes in 𝒫 along the $i$-th coordinate direction through

julia> i = 1   # for example, the first coord. direction1
julia> subdivide!(𝒫_tree, i)BoxTree of depth 2

BoxSet

A core idea of GAIO.jl is to approximate a subset of the domain $X$ via a collection of boxes. To construct a BoxSet, there are two typical options: retrieving all boxes from some partition, or locating a box surrounding a (or some) point(s) in $X$:

julia> ℬ = cover(𝒫, :)    # set of all boxes from 𝒫8 - element BoxSet in 4 x 2 - element BoxGrid
julia> x = (0.2, 0.4)(0.2, 0.4)
julia> ℬ = cover(𝒫, x)    # set of one box containing the point x1 - element BoxSet in 4 x 2 - element BoxGrid
julia> xs = [ rand(2) for _=1:10 ]10-element Vector{Vector{Float64}}:
 [0.006073523657513791, 0.6730544888011307]
 [0.8490039514387125, 0.26339689105327313]
 [0.45289051604537045, 0.8892310358963075]
 [0.36108367686430576, 0.13013644681561976]
 [0.7152533176189678, 0.11347461252374103]
 [0.3806613785030396, 0.8509468785463707]
 [0.21675983605382088, 0.7998822501202824]
 [0.0035151572025705624, 0.8900386881142215]
 [0.045454173697344835, 0.8751375744162035]
 [0.9365771517608589, 0.3091801720691635]
julia> ℬ = cover(𝒫, xs)   # set of boxes containing the points in xs5 - element BoxSet in 4 x 2 - element BoxGrid

You can access the boxes or their internal data via iteration

julia> for box in ℬ
           center, radius = box
           # do something
       end
       
       # get an array of boxes
julia> arr_of_boxes = collect(ℬ)
       
       # get an array of box centers5-element Vector{Box{2, Float64}}:
 [0.0, 0.25) × [0.5, 1.0)
 [0.5, 0.75) × [0.0, 0.5)
 [0.75, 1.0) × [0.0, 0.5)
 [0.25, 0.5) × [0.5, 1.0)
 [0.25, 0.5) × [0.0, 0.5)
julia> arr_of_centers = collect(box.center for box in ℬ)
       
       # get an array of box radii5-element Vector{SVector{2, Float64}}:
 [0.125, 0.75]
 [0.625, 0.25]
 [0.875, 0.25]
 [0.375, 0.75]
 [0.375, 0.25]
julia> arr_of_radii = collect(box.radius for box in ℬ)5-element Vector{SVector{2, Float64}}:
 [0.125, 0.25]
 [0.125, 0.25]
 [0.125, 0.25]
 [0.125, 0.25]
 [0.125, 0.25]

BoxMap

A box map is a function which maps a BoxSet to a BoxSet. Given a map $f : \mathbb{R}^d\to \mathbb{R}^d$

julia> f((x,y)) = (1-1.4*x^2+y, 0.3*x)   # the Hénon mapf (generic function with 1 method)

and some box X as domain, you construct a BoxMap by

julia> F = BoxMap(f, X)BoxMap over [0.0, 1.0) × [0.0, 1.0)

By default, GAIO.jl will try to adaptively choose sample points in a set to compute the image of the set by approximating the (local) Lipschitz constant of the map $f$. There are many other types of BoxMap discretizations available, see the section on BoxMaps.

We can now map a BoxSet ℬ via the BoxMap F by

julia> 𝒞 = F(ℬ)6 - element BoxSet in 4 x 2 - element BoxGrid

where the output 𝒞 is also a BoxSet.

For long running computations, you can also display a progress meter

julia> using ProgressMeter
julia> 𝒞 = F(ℬ; show_progress = true)6 - element BoxSet in 4 x 2 - element BoxGrid

TransferOperator

The Perron-Frobenius operator (or transfer operator) [[5]] associated to a map $f$ can be approximated using the TransferOperator type. To construct a (discretized) TransferOperator from a BoxMap $F$ on the domain BoxSet ℬ, you type

julia> 𝒫 = BoxGrid(X, (20,20))     # 20 x 20 so there's more going on20 x 20 - element BoxGrid
julia> ℬ = cover(𝒫, :)             # cover the whole partition400 - element BoxSet in 20 x 20 - element BoxGrid
julia> T = TransferOperator(F, ℬ, ℬ)400 x 400 TransferOperator over 20 x 20 - element BoxGrid with 1392 stored entries:

⎡⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎤
⎢⠠⠀⠄⠠⠀⠀⠄⠢⠄⠀⠐⠂⡄⠠⠂⠀⠠⠠⠢⠀⠤⠤⠐⠀⠀⠤⠀⠀⠀⠤⠤⠀⠂⠀⠠⠀⠀⠴⠆⡠⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⎥
⎢⠀⠀⠂⠀⠂⠀⡀⠂⡐⢐⠀⠀⠂⠒⠂⠀⠐⠒⢀⠐⢀⢂⠀⠐⠂⢒⠒⠀⠂⡀⣀⢀⠀⠀⢀⠀⠀⢀⡀⢀⎥
⎢⠀⠀⠁⠀⠀⠀⠀⠀⠀⠐⠀⢀⢁⢂⠀⠀⠐⠀⠀⠀⠀⡂⣘⠂⢀⡀⡐⡀⠀⢀⡀⡀⢀⣀⢀⢀⠀⠀⠀⠁⎥
⎢⠀⠠⡄⠡⡄⠀⠀⡉⢨⠀⠬⠠⠀⠠⡅⠄⠠⠄⠁⢨⠄⠡⠥⠨⡠⠥⠅⠄⡥⠠⠁⠥⡩⠤⠥⠈⠀⠈⠁⠀⎥
⎢⢐⣓⠃⠐⠀⠀⠀⠐⠀⡃⠀⡠⠂⡓⠀⠀⠐⠐⠄⠀⠂⠂⢼⠀⠠⠆⠀⠄⠀⠰⡂⠀⠀⢘⡀⢠⠀⠂⠀⠑⎥
⎢⠠⢤⠄⢁⢄⠀⢈⢀⢈⠄⡀⠌⠈⣄⡀⢄⠀⢠⡁⢨⡀⣁⡩⢈⠌⣁⠅⠁⣄⠈⡍⢁⢈⠁⡁⠈⡀⠠⡀⠀⎥
⎢⣬⢨⡄⡄⣤⠀⠀⠀⠀⡄⡌⡅⣄⢁⢀⠀⢠⣀⡀⡀⢠⢠⡎⢠⢀⡂⢠⢢⠁⢠⠁⡤⡄⣠⢂⢰⢠⢈⢡⣀⎥
⎢⠐⢐⡄⢀⡂⠀⠀⡀⢀⠔⣐⠂⠄⢂⡀⠂⢐⡐⡂⠐⡀⢂⠲⢂⡂⣀⠒⠀⣀⠀⠀⠐⡒⠠⠀⠂⠀⠐⠂⠔⎥
⎢⣉⢨⡀⡁⠌⠀⡀⢀⢀⡅⢁⢄⣥⢠⠀⡈⠨⣁⡁⣀⡠⣌⡥⢨⢀⡀⡬⣬⢁⣈⣄⣭⣡⣭⡠⢈⠡⢈⠈⠁⎥
⎢⠀⠰⠄⠠⠂⠀⠠⠠⠠⠂⠔⠂⠀⠰⠄⠂⠐⠠⠄⠀⠔⠠⠀⠠⠀⠦⠒⠒⠴⠀⠀⠒⠐⠒⠐⠀⠐⠀⠀⠐⎥
⎢⣒⢲⡂⡖⢒⠁⠚⠒⠚⠄⡒⠊⡤⠒⠂⠀⠰⢗⠗⢪⠓⢖⡴⢰⠘⡖⠶⠤⠐⢚⠞⢤⣎⣤⠑⢀⠀⢐⣂⣐⎥
⎢⠉⠀⠅⠫⠨⠀⠢⠸⠁⠨⠁⠀⠗⠜⠌⠠⠨⠧⠊⠒⠭⠯⠪⠰⠐⠨⠨⠨⠀⠊⠯⠐⠒⠈⠪⠈⠄⠕⠎⠈⎥
⎢⠚⠐⡃⠒⡙⠀⠔⡰⢄⠠⠃⠁⠂⡢⠈⠑⠘⠒⠣⢄⠨⠐⠄⠀⡠⠑⡘⠐⡑⠜⠴⠡⠢⠋⠨⠘⠂⠲⠖⠪⎥
⎢⠛⠺⠅⠳⣐⠀⠲⠴⠣⠨⢦⡀⠃⢺⠆⡰⢀⠺⡞⠂⢮⠸⠈⠠⠐⡿⢸⡸⢶⠚⠽⣀⠇⡛⢨⢀⠆⠸⠄⢙⎥
⎢⠩⣥⡅⠩⠡⠄⠠⠩⠈⡅⠁⠄⠄⡭⠁⠰⠨⠍⠆⠈⠅⠕⢨⠈⠁⠹⠙⠐⠁⠈⡵⠀⠡⢨⠀⠀⠂⠌⠀⠩⎥
⎢⢰⢒⡆⢒⢲⠀⢒⢐⠂⠒⢄⡂⡂⡒⠢⠊⢐⣒⡓⠒⠂⢂⢒⠐⠒⢐⠚⠈⢃⠢⡒⠐⠒⠠⠀⠠⠁⣒⡂⠐⎥
⎢⣐⡐⡀⡐⡐⠀⠐⢐⠀⡀⠒⠂⡂⣂⠂⡀⠀⢒⠂⢀⡂⢐⢐⠀⡀⡐⡀⠀⠐⢀⠀⠀⡂⣀⠀⡀⠀⢐⠂⣀⎥
⎣⠀⠀⠄⠀⠀⠀⠀⠄⠠⠨⠠⠀⠀⠠⠄⠀⠠⠤⠄⠀⠀⠄⠀⠠⠄⠤⠤⠀⠄⠀⠀⠈⠄⠀⠀⠀⠀⠀⠀⠀⎦

Again, a progress meter can be displayed with the additional keyword argument show_progress = true.

Formally, T is a linear (in fact, a Markov) operator on the space of steps functions on the box set ℬ. Of particular interest are often certain eigenvectors. For example, an eigenvector at the eigenvalue 1 is an (approximate) invariant measure, which characterizes the long term behaviour of $f$ according to Birkhoff's ergodic theorem.

julia> λ, ev = eigs(T)(ComplexF64[0.6024149797975805 + 0.0im, -0.293057626130105 + 0.0im, 0.10206207261597958 + 0.0im], BoxMeasure{Box{2, Float64}, Tuple{Int64, Int64}, ComplexF64, BoxGrid{2, Float64, Int64}, OrderedCollections.OrderedDict{Tuple{Int64, Int64}, ComplexF64}}[BoxMeasure in 20 x 20 - element BoxGrid with 400 boxes in its suport, BoxMeasure in 20 x 20 - element BoxGrid with 400 boxes in its suport, BoxMeasure in 20 x 20 - element BoxGrid with 400 boxes in its suport], 3, 1, 20, [1.3480855942419456e-7, -2.5156626722598564e-8, -9.08398430559264e-8, 4.397254943258681e-8, 1.723191910927363e-7, 1.7245427819635085e-8, -1.1679581193354769e-7, 8.706224017991166e-8, -7.292880156621694e-9, 7.971422625606301e-10  …  -9.408707674202613e-8, 2.186307394178859e-8, -8.495220030844801e-7, -1.992036070301162e-5, 1.3587462373948192e-7, 5.527422893691085e-8, -4.9106579924733965e-8, 7.430098308410094e-8, 1.6611550650774443e-8, 1.8986434129292708e-7])
julia> μ = ev[1]   # ev is an array of measures, grab the first oneBoxMeasure in 20 x 20 - element BoxGrid with 400 boxes in its suport

This can also be done with the adjoint Koopman operator T'.

BoxMeasure

The second output of eigs(T) is a vector of (discrete) measures, BoxMeasures. A BoxMeasure is absolutely continuous w.r.t. the volume (i.e. Lebesgue) measure and its density is piecewise constant on the boxes of the domain ℬ. One can let T act on a BoxMeasure simply through multiplication

julia> ν = T*μBoxMeasure in 20 x 20 - element BoxGrid with 90 boxes in its suport

Of course, the same holds for the the Koopman operator as well.

julia> ν = T'*μBoxMeasure in 20 x 20 - element BoxGrid with 182 boxes in its suport

One can evaulate a BoxMeasure on an arbitrary BoxSet ℬ

julia> μ(ℬ)6.335723408296507 + 0.0im

Similarly, one can integrate a function with respect to a BoxMeasure via

julia> sum(x -> sin(x[1] + 2x[2]), μ)5.036354373963653 + 0.0im

Marginal distributions can be accessed using the marginal function

julia> marginal(μ; dim=1)BoxMeasure in 20 - element BoxGrid with 20 boxes in its suport

The measures can also be associated with a (Lebesgue) density

julia> p = density(μ)(::GAIO.var"#eval_density#162"{BoxMeasure{Box{2, Float64}, Tuple{Int64, Int64}, ComplexF64, BoxGrid{2, Float64, Int64}, OrderedCollections.OrderedDict{Tuple{Int64, Int64}, ComplexF64}}, BoxGrid{2, Float64, Int64}}) (generic function with 1 method)

Since a measure μ is a function defined over measurable sets, composite measures $g \circ \mu$ are well-defined for functions $g : \mathbb{R} \to \mathbb{R}$ (or $g : \mathbb{C} \to \mathbb{C}$). This is supported in GAIO.jl for BoxMeasures

julia> η = exp ∘ μBoxMeasure in 20 x 20 - element BoxGrid with 400 boxes in its suport

Finite signed measures can be given a vector space structure. This is also supported:

julia> ν + μBoxMeasure in 20 x 20 - element BoxGrid with 400 boxes in its suport
julia> 2ν - μ/2BoxMeasure in 20 x 20 - element BoxGrid with 400 boxes in its suport

Graphs of Boxes

One can equivalently view the transfer operator as a weighted directed graph. That is, the matrix of a TrensferOperator is the (transposed) weighted adjacency matrix for a graph. This graph can be constructed explicitely using the MetaGraphsNext.jl package

julia> using Graphs, MetaGraphsNextWARNING: using Graphs.density in module Main conflicts with an existing identifier.
julia> G = MetaGraph(T)Meta graph based on a SimpleDiGraph{Int64} with vertex labels of type Tuple{Int64, Int64}, vertex metadata of type Nothing, edge metadata of type Float64, graph metadata given by 20 x 20 - element BoxGrid, and default weight 1.0

See also the Graphs and MetaGraphsNext documentation.

Plotting

GAIO.jl offers both Plots or Makie for plotting. To plot a BoxSet or a BoxMeasure, simply choose either Plots or a Makie backend, eg. GLMakie, and call plot on a BoxSet or BoxMeasure

using GLMakie: plot

plot(ℬ)
plot(μ)

Plotting works with all the functionality of either package. This means you can set box plots as subplots, add colorbars, etc., using the Plots or Makie interface. For an example, see examples/invariant_measure_2d.jl.