Usage

The base of the numerical set oriented methods of this framework are BoxSet (the discretization of a set of boxes) and BoxMap (the discretization of a map). Thus, in the following, we will have a closer look at the two and other useful things to know when using GAIO.jl.

To create a Box given its center point c = (c_1, c_2, ..., c_d) as well as its "radius" in every axis direction r = (r_1, r_2, ..., c_d), simply type

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

This creates a set $Q = [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 by calling

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

BoxPartition

Most algorithms in GAIO.jl revolve around a partition of the domain $Q$ into small boxes. To create an $n_1 \times \ldots \times n_d$ - element equidistant grid of boxes, we can pass the tuple $n = (n_1, \ldots, n_d)$ into the function BoxPartition

julia> n = (4, 2)(4, 2)
julia> P = BoxPartition(Q, n)4 x 2 - element BoxPartition

BoxPartitions use a cartesian indexing structure to be memory-efficient. These indices are accessed and used through the API:

julia> x = (0.2, 0.1)(0.2, 0.1)
julia> key = point_to_key(P, x)    # x is some point in the domain Q(1, 1)
julia> box = key_to_box(P, 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(P, x)    # performs both above functionsBox{2, Float64}:
   center: [0.125, 0.25]
   radius: [0.125, 0.25]

TreePartition

For partitions of $Q$ into variably sized boxes, one can use TreePartition:

julia> P2 = TreePartition(Q)TreePartition of depth 1

A TreePartition uses a binary tree structure to store a partition of the domain. Every Box of a TreePartition can be split using the command

julia> subdivide!(P2)TreePartition of depth 2
julia> subdivide!(P2)TreePartition of depth 3
julia> subdivide!(P2)TreePartition of depth 4

The axis direction along which to subdivide cycles through with the depth, i.e. subdividing at depth 1 splits along dimension 1, subdividing at depth d+1 splits along dimension 1 again.

The TreePartition created above is equivalent to a 4x2 BoxPartition. One can retrieve this using

julia> P3 = BoxPartition(P2)4 x 2 - element BoxPartition

TreePartitions use indices of the type (depth, cartesian_index) where cartesian_index is the equivalent index of a BoxPartition with the same size as a TreePartition subdivided depth times. In other words,

julia> key_to_box( P, (1, 1) ) == key_to_box( P2, (4, (1, 1)) )true
julia> key_to_box( P, (4, 2) ) == key_to_box( P2, (4, (4, 2)) )true

BoxSet

The core idea behind GAIO.jl is to approximate a subset of the domain via a collection of small boxes. To construct BoxSets, there are two main options: getting all boxes in the partition, or locating a box surrounding a point $x \in Q$

julia> B = cover(P, x)    # one box surrounding the point x1 - element BoxSet in 4 x 2 - element BoxPartition
julia> B = cover(P, :)    # set of all boxes in P8 - element BoxSet in 4 x 2 - element BoxPartition

One can also create a Boxset from an iterable of Boxes. This will cover every element of the iterable with boxes from P:

julia> x1 = (0.2, 0.1)(0.2, 0.1)
julia> box1 = point_to_box(P, x1)Box{2, Float64}:
   center: [0.125, 0.25]
   radius: [0.125, 0.25]
julia> x2 = (0.3, 0.6)(0.3, 0.6)
julia> box2 = point_to_box(P, x2)Box{2, Float64}:
   center: [0.375, 0.75]
   radius: [0.125, 0.25]
julia> B = cover(P, [box1, box2])2 - element BoxSet in 4 x 2 - element BoxPartition

BoxSet is a highly memory-efficient way of storing boxes. However, should you want to access the boxes or their internal data, this can be done via iteration:

julia> for box in B
           center, radius = box
           # do something
       end
       
       # get an array of boxes
julia> arr_of_boxes = collect(B)
       
       # get an array of box centers2-element Vector{Box{2, Float64}}:
 [0.0, 0.25) × [0.0, 0.5)
 [0.25, 0.5) × [0.5, 1.0)
julia> arr_of_centers = collect(box.center for box in B)
       
       # get an array of box radii2-element Vector{SVector{2, Float64}}:
 [0.125, 0.25]
 [0.375, 0.75]
julia> arr_of_radii = collect(box.radius for box in B)
       
       # (memory-efficiently) create a matrix where each center is a column2-element Vector{SVector{2, Float64}}:
 [0.125, 0.25]
 [0.125, 0.25]
julia> mat_of_centers = reinterpret(reshape, eltype(arr_of_centers[1]), arr_of_centers)2×2 reinterpret(reshape, Float64, ::Vector{SVector{2, Float64}}) with eltype Float64:
 0.125  0.375
 0.25   0.75

BoxMap

A BoxMap is a function which maps boxes to boxes. Given a pointmap f : ℝᵈ → ℝᵈ which accepts an SVector from StaticArrays.jl (or just an NTuple) and returns the same,

julia> function f(u)   # the Baker transformation
           x, y = u
           if x < 0.5
               (2x, y/2)
           else
               (2x - 1, y/2 + 1/2)
           end
       endf (generic function with 1 method)

initialize the corresponding BoxMap F by

julia> F = BoxMap(f, Q)SampledBoxMap with 2 sample points

This will generate a BoxMap which attempts to calculate setwise images of f. One can also use a DynamicalSystem from DynamicalSystems.jl:

julia> using DynamicalSystems: DiscreteDynamicalSystem
julia> using StaticArrays
julia> dynamic_rule(u::Vec, p, t) where Vec = Vec( f(u) ) # DynamicalSystems.jl syntaxdynamic_rule (generic function with 1 method)
julia> u0 = SA_F64[0, 0]   # this is dummy data2-element SVector{2, Float64} with indices SOneTo(2):
 0.0
 0.0
julia> p0 = SA_F64[0, 0]   # parameters, in case they are needed2-element SVector{2, Float64} with indices SOneTo(2):
 0.0
 0.0
julia> system = DiscreteDynamicalSystem(dynamic_rule, u0, p0)2-dimensional DeterministicIteratedMap
 deterministic: true
 discrete time: true
 in-place:      false
 dynamic rule:  dynamic_rule
 parameters:    [0.0, 0.0]
 time:          0
 state:         [0.0, 0.0]
julia> F̃ = BoxMap(system, Q)SampledBoxMap with 2 sample points

The same works for a continuous dynamical system.

julia> using DynamicalSystems: ContinuousDynamicalSystem
julia> using OrdinaryDiffEq
julia> function lorenz_dudx(u::Vec, p, t, Δt) where Vec
           x,y,z = u
           σ,ρ,β = p
       
           dudx = Vec(( σ*(y-x), ρ*x-y-x*z, x*y-β*z ))
           return Δt * dudx
       endlorenz_dudx (generic function with 1 method)
julia> Δt = 0.20.2
julia> lorenz(u, p=p0, t=0) = lorenz_dudx(u, p, t, Δt)lorenz (generic function with 3 methods)
julia> u0 = SA_F64[0, 0, 0]3-element SVector{3, Float64} with indices SOneTo(3):
 0.0
 0.0
 0.0
julia> p0 = SA_F64[10, 28, 0.4]3-element SVector{3, Float64} with indices SOneTo(3):
 10.0
 28.0
  0.4
julia> diffeq = (alg = RK4(), dt = 1/20, adaptive = false)(alg = RK4(; stage_limiter! = trivial_limiter!, step_limiter! = trivial_limiter!, thread = static(false),), dt = 0.05, adaptive = false)
julia> lorenz_system = ContinuousDynamicalSystem(lorenz, u0, p0; diffeq)3-dimensional CoupledODEs
 deterministic: true
 discrete time: false
 in-place:      false
 dynamic rule:  lorenz
 ODE solver:    RK4
 ODE kwargs:    (dt = 0.05, adaptive = false)
 parameters:    [10.0, 28.0, 0.4]
 time:          0.0
 state:         [0.0, 0.0, 0.0]
julia> Q̄ = Box((0,0,25), (30,30,30))Box{3, Float64}:
   center: [0.0, 0.0, 25.0]
   radius: [30.0, 30.0, 30.0]
julia> F̄ = BoxMap(lorenz_system, Q̄)SampledBoxMap with 2 sample points

By default, GAIO.jl will try to adaptively choose test points to compute setwise images by approximating the (local) Lipschitz constant of the map. There are many other types of BoxMap discretizations available, see the section on BoxMaps for more information.

Using BoxMap

Now, one can map a BoxSet via the BoxMap F by simply calling F as a function

julia> C = F(B)7 - element BoxSet in 4 x 2 - element BoxPartition

where the output C is also a BoxSet.

For long running computations, GAIO.jl can also display a progress meter

julia> using ProgressMeter
julia> C = F(B; show_progress = true)7 - element BoxSet in 4 x 2 - element BoxPartition

(Adding a progress meter adds a little bit of overhead, so for super short computations like the above it isn't recommended)

TransferOperator

The Perron-Frobenius operator (or transfer operator) [5] is discretized in GAIO.jl using the TransferOperator type. To initialize a TransferOperator that acts on a subdomain of $Q$, type

julia> B = cover(P, :)8 - element BoxSet in 4 x 2 - element BoxPartition
julia> T = TransferOperator(F, B)   # T operates on the domain covered by the box set B8 x 8 TransferOperator over 4 x 2 - element BoxPartition with 33 stored entries:

⎡⢿⣇⠶⠆⎤
⎣⠘⠛⠘⡿⎦

In this case, the codomain is generated automatically. This is not always ideal (e.g. in eigenvalue calculations), so the codomain should often be specified as the third argument

julia> T = TransferOperator(F, B, B)8 x 8 TransferOperator over 4 x 2 - element BoxPartition with 33 stored entries:

⎡⢻⣧⠩⠥⎤
⎣⢤⣄⢤⣗⎦

Again, a progress meter can be displayed for long computations

julia> using ProgressMeter
julia> T = TransferOperator(F, B, B; show_progress = true)8 x 8 TransferOperator over 4 x 2 - element BoxPartition with 33 stored entries:

⎡⢻⣧⠩⠥⎤
⎣⢤⣄⢤⣗⎦

To convert this to the underlying transfer matrix described in [6], one can simply call the sparse function from SparseArrays

julia> using SparseArrays
julia> sparse(T)8×8 SparseMatrixCSC{Float64, Int64} with 33 stored entries:
 0.333333  0.166667  0.166667   ⋅        0.5  0.25  0.214286    ⋅
 0.333333  0.166667  0.166667   ⋅         ⋅    ⋅     ⋅          ⋅
  ⋅        0.166667  0.166667  0.333333   ⋅   0.25  0.214286   0.428571
  ⋅        0.166667  0.166667  0.333333   ⋅    ⋅     ⋅          ⋅
  ⋅         ⋅         ⋅         ⋅         ⋅    ⋅    0.214286    ⋅
  ⋅         ⋅         ⋅         ⋅         ⋅    ⋅    0.214286   0.428571
 0.333333  0.166667  0.166667   ⋅        0.5  0.25  0.0714286   ⋅
  ⋅        0.166667  0.166667  0.333333   ⋅   0.25  0.0714286  0.142857

To find an approximate invariant measure over B use the eigs function from Arpack.jl. All keyword arguments from Arpack.eigs are supported.

julia> # for the Baker trafo, the Lebesgue measure
       # - i.e. the constant-weight measure - is invariant
       λ, ev = eigs(T);
julia> λ3-element Vector{ComplexF64}:
   0.9999999999999991 + 0.0im
  0.49999999999999967 + 0.0im
 -0.37896131637649966 + 0.0im
julia> ev3-element Vector{BoxMeasure{Box{2, Float64}, Tuple{Int64, Int64}, ComplexF64, BoxPartition{2, Float64, Int64}, OrderedCollections.OrderedDict{Tuple{Int64, Int64}, ComplexF64}}}:
 BoxMeasure in 4 x 2 - element BoxPartition with 8 boxes in its suport
 BoxMeasure in 4 x 2 - element BoxPartition with 8 boxes in its suport
 BoxMeasure in 4 x 2 - element BoxPartition with 8 boxes in its suport
julia> μ = ev[1]   # ev is an array of measures, grab the first oneBoxMeasure in 4 x 2 - element BoxPartition with 8 boxes in its suport

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

BoxMeasure

The second output of eigs(T) is a discrete measure, a BoxMeasure. This measure 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 4 x 2 - element BoxPartition with 8 boxes in its suport

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

julia> ν = T'μBoxMeasure in 4 x 2 - element BoxPartition with 8 boxes in its suport

One can evaulate a BoxMeasure on a BoxSet

julia> B8 - element BoxSet in 4 x 2 - element BoxPartition
julia> μ(B)-2.6065137370063387 + 0.0im

Similarly, one can integrate a function with respect to a BoxMeasure by calling

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

Marginal distributions can be accessed using the marginal function

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

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

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

Since a measure $\mu$ 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 4 x 2 - element BoxPartition with 8 boxes in its suport

For multiple BoxMeasures, the concatenation operator ∘ can be applied to each one using julia's broadcasting functionality

julia> real_ev = real .∘ ev3-element Vector{BoxMeasure{Box{2, Float64}, Tuple{Int64, Int64}, Float64, BoxPartition{2, Float64, Int64}, OrderedCollections.OrderedDict{Tuple{Int64, Int64}, Float64}}}:
 BoxMeasure in 4 x 2 - element BoxPartition with 8 boxes in its suport
 BoxMeasure in 4 x 2 - element BoxPartition with 8 boxes in its suport
 BoxMeasure in 4 x 2 - element BoxPartition with 8 boxes in its suport

Similarly, finite signed measures can be given a vector space structure. This is also supported in GAIO.jl

julia> ν + μBoxMeasure in 4 x 2 - element BoxPartition with 8 boxes in its suport
julia> 2ν - μ/2BoxMeasure in 4 x 2 - element BoxPartition with 8 boxes in its suport

A BoxMeasure is implemented by a dictionary, mapping boxes to weights

Iterators.take(μ, 3)

To access this structure oneself one can call

P = μ.partition
key_val_pairs = pairs(μ)

Graphs of Boxes

One could equivalently view the transfer operator as a weighted directed graph. That is, a transfer matrix in GAIO.jl is the (transposed) weighted adjacency matrix for a graph. This graph can be constructed 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 4 x 2 - element BoxPartition, and default weight 1.0

See the MetaGraphsNext documentation for how to interface with this data type.

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(B)

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.