Documentation

BFGS_periodic(FJ::Function, x::AbstractVector, q::Integer;
              maxiter::Integer=50)

Find a periodic orbit using BFGS from the Optim package.

Arguments:

• FJ: A function that returns the map and its derivative F, dFdx = FJ(x)
• x: An initial point to find the orbit from
• q: The period of the orbit
• maxiter=50: The maximum number of optimization steps allows

Output:

• xs: A periodic trajectory of length d
• res: An Optim return object, (can check convergence with Optim.converged(res))

newton_periodic(FJ::Function, x::AbstractVector, q::Integer;
                maxiter::Integer=50, rtol::Number=1e-8, verbose::Bool=false)

Find a periodic orbit using Newton's method with line search

Arguments:

• FJ: A function that returns the map and its derivative F, dFdx = FJ(x)
• x: An initial point to find the orbit from
• q: The period of the orbit
• maxiter=50: The maximum number of optimization steps allows
• rtol=1e-8: The residual tolerance required
• verbose=false: If true, outputs information about convergence

Output:

• xs: A periodic trajectory of length d
• converged: A flag that indicates whether the orbit converged in maxiter steps

abstract type InvariantCircle

An abstract type representing a chain of invariant circles zᵢ where zᵢ₊₁ = F(zᵢ) for i<=p and z₀ = F(zₚ₋₁). The main implementation of this type of FourierCircle, but other potential implementations of circles could exist.

Implementations should have the following functions for immediate portability:

Basic stuff

• Initializer InvariantCircleImplementation(stuff)
• Similar Base.similar(z::InvariantCircle)

Getters and Setters
- Get number of unknown variables per circle get_N(z::InvariantCircle)
- Get period of circle get_p(z::InvariantCircle)

• Get parameters (except τ) get_a(z::InvariantCircle)
• Set parameters (except τ) set_a!(z::InvariantCircle, a::AbstractArray)
• Get rotation number τ in [0,2π) get_τ(z::InvariantCircle)
• Set τ set_τ!(z::InvariantCircle, τ::Number)

Evaluation related routines

• Evaluate the circle (see (z::InvariantCircle)(θ, i_circle)) evaluate(z::InvariantCircle, θ::AbstractVector; i_circle::Integer=1)
• Evaluate the derivative of the circle w.r.t. θ deval(z::InvariantCircle, θ::AbstractVector; i_circle::Integer=1)
• Get a basis Φ for z evaluated at Nθ equispaced points (requires linearity) get_Φ(z::InvariantCircle, Nθ::Integer; τ::Number=0.0)
• Evaluate the invariant circle on Φ function grid_eval!(x::AbstractArray, z::InvariantCircle, Φ::AbstractArray; i_circle=0)
• Evaluate the derivative of the invariant circle on Φ grid_deval!(dx::AbstractArray, z::InvariantCircle, Φ::AbstractArray; i_circle=0)

Constraints (for Newton iteration)

• Constrain the value at 0 fixed_θ0_constraint(z::InvariantCircle)

Other routines (useful for continuation)

• Get average radius LinearAlgebra.average_radius(z::InvariantCircle)
• Get the area of the invariant circle area(z::InvariantCircle; i_circle=1, Ns = 100)

evaluate(z::InvariantCircle, θ::Number; i_circle::Integer=1)

Evaluate the i_circleth circle in z at the point θ

deval(z::InvariantCircle, θ::Number; i_circle::Integer=1)

Evaluate the derivative of the i_circleth circle in z at the point θ

shifted_eval(z::InvariantCircle, θ::AbstractVector; i_circle::Integer=1)

Evaluate the i_circleth circle in z at the point θ, shifted by τ

get_circle_residual(F::Function, z::InvariantCircle, Nθ::Integer)

Get the KAM-like residual of a chain of p invariant circles.
> Rᵢ₁ = z₁(θᵢ+τ) - F(z_p(θᵢ))
> Rᵢⱼ = zⱼ(θᵢ) - F(zⱼ₋₁(θᵢ)) for 2<=j<=p

Arguments:

• F: Symplectic map
• z: Invariant circle
• Nθ: Number of θ points where the residual is taken

gn_circle(FJ::Function, z::InvariantCircle, Nθ::Integer;
          maxiter::Integer=10, rtol::Number=1e-8, verbose::Bool=false,
          monitor::Function=(z, rnorm) -> nothing,
          constraint::Function=fixed_θ0_constraint, λ::Number=0)

Find an invariant circle using Gauss-Newton with linesearch. Implemented with dense linear algebra (no FFTs yet). Returns the number of iterations taken to converge. If an evaluation fails, returns -1. If the routine fails to converge, returns -2

Arguments:

• z::InvariantCircle: An initial connecting orbit guess, see linear_initial_connecting_orbit
• FJ::Function: Function defined on R² with signature F(x), J(x) = FJ(x) where F is the symplectic map and J = dF/dx
• Nθ: Number of quadrature nodes for the least squares
• maxiter = 10: Maximum number of Gauss Newton steps
• rtol = 1e-8: Stopping tolerance

FourierCircle <: InvariantCircle

Constructors:

• FourierCircle(Na::Integer; a::AbstractArray=[], τ::Number=0., p::Integer=1)

Fourier invariant circle data structure (see InvariantCircle), used to compute a circle z that is invariant of a map Fᵖ(z), where p is the period. This is done by storing all of the circles Fⁱ(z) for 0<=i<p. The circles are stored as Fourier series, i.e. z(θ) = a₀ + ∑₁ᴺᵃ Aⱼ [cos(jθ), sin(jθ)]ᵀ The coefficients of the invariant circle can be accessed and set in chunks via array indexing. That is, You can get and set a₀ of the ith circle using z[0, i] You can get and set Aⱼ of the ith circle using z[j, i]

get_Na(z::FourierCircle)

Get number of Fourier modes

get_p(z::FourierCircle)

Get number of islands

Base.length(z::FourierCircle)

Get total number of parameters, excluding τ

get_a0(z::FourierCircle; i_circle::Integer=1)

Get constant term of Fourier series of the i_circleth circle in island chain

set_a0!(z::FourierCircle, a0::AbstractArray; i_circle::Integer=1)

Set constant term of Fourier series of the i_circleth circle in island chain

get_Am(z::FourierCircle, m::Integer; i_circle::Integer=1)

Get mth Fourier coefficients of the i_circleth circle in island chain

set_Am!(z::FourierCircle, m::Integer, a::AbstractArray; i_circle::Integer=1)

Set mth Fourier coefficients of the i_circleth circle in island chain

Base.getindex(z::FourierCircle, i_A::Integer, i_circle::Integer)

Get the coefficients of z. z[0, i_circle] gets the constant term. z[i_A, i_circle] gets higher coefficients

Base.setindex!(z::FourierCircle, X::AbstractArray, i_A::Integer,
               i_circle::Integer)

Set the coefficients of z. z[0, i_circle] sets the constant term. z[i_A, i_circle] sets higher coefficients

get_τ(z::FourierCircle)

Get the rotation number (the circle is 2π periodic)

set_τ!(z::FourierCircle, τ::Number)

Set the rotation number (the circle is 2π periodic)

average_radius(z::FourierCircle)

Norm defined as the L2 norm of the radius, i.e. ‖z‖² = 1/2πp ∫ ||z(θ) - a₀||^2 dθ = 1/2p ∑ₖ Tr(AₖᵀAₖ) If i_circle = 0, takes the norm over all circles. If i_circle != 0, finds the average radius of a single circle.

evaluate(z::FourierCircle, θ::AbstractVector{T}; i_circle::Integer=1) where {T}

Evaluate the i_circleth circle in the chain at a vector of points θ

deval(z::FourierCircle, θ::AbstractVector; i_circle::Integer=1)

Evaluate the derivative of the i_circleth circle in the chain at a vector of points θ

deriv(z::FourierCircle)

Return the derivative of the FourierCircle, wrapped in a circle object

area(z::FourierCircle; i_circle=1)

Get the area of the circle.

area(c::ConnectingOrbit; origin::AbstractArray=[0.0,0.0],
     i_p::Integer=1, rtol=1e-8)

Get the "area" of a connecting orbit ∫₋₁¹ (x dy/ds - y dy/ds) ds

circle_linear!(z::FourierCircle, FJ::Function, x::AbstractArray, h::Number)

Fill z's Fourier coefficients using a linearization about a stable periodic orbit. Useful for seeding continuation.

Arguments:

• z: An invariant circle
• FJ: A function that returns the map and its derivative F, dFdx = FJ(x)
• x: A periodic orbit of F
• h: The average radius of the first linearized invariant circle

FourierTorus <: InvariantTorus

Constructors:

• FourierTorus(d::Integer, Na::Vector{<:Integer}; a::AbstractArray=[], τ::AbstractVector=[0.], p::Integer=1)

Fourier invariant torus data structure (see InvariantTorus), used to compute a torus z that is invariant of a map Fᵖ(z), where p is the period. This is done by storing all of the tori Fⁱ(z) for 0<=i<p. The tori are stored as Fourier series, i.e. z(θ) = a₀ + ∑₁ᴺᵃ Aⱼ [cos(jθ), sin(jθ)]ᵀ The coefficients of the invariant torus can be accessed and set in chunks via array indexing. That is, You can get and set a₀ of the ith torus using z[0, i] You can get and set Aⱼ of the ith torus using z[j, i]

evaluate_on_grid(tor::FourierTorus, thetavecs::AbstractVector)

Evaluate the torus on a grid defined by the vector-of-vectors thetavecs.

kam_residual(tor::FourierTorus, F::Function, thetavecs::AbstractVector)

Compute the function F(tor(θ))-tor(θ+τ) on the grid defined by the vector-of-vectors thetavecs. If the invariant torus was computed with a non-identity observable h, then F needs to be given by h ∘ F ∘ h^-1.

kam_residual_rnorm(tor::FourierTorus, F::Function, thetavecs::AbstractVector)

Compute a relative error of the KAM residual |F(tor(θ))-tor(θ+τ)|/|tor|. If the invariant torus was computed with a non-identity observable h, then F needs to be given by h ∘ F ∘ h^-1.

ConnectingOrbit(Na::Integer, p::Integer)

Initialize a connecting orbit consisting of p segments cⱼ : [0,1] -> R². The endpoints cⱼ(0) and cⱼ(1) are both hyperbolic perodic orbits. The connecting orbits connect these, acting as the outer boundary of an island. The connecting orbits are parameterized by Chebyshev polynomials of degree Na.

See linear_initial_connecting_orbit and gn_connecting! for functions to initialize and compute a connecting orbit.

get_Na(c::ConnectingOrbit)

Get the polynomial degree of the connecting orbit.

get_p(c::ConnectingOrbit)

Get the period of the connecting orbit

get_am(c::ConnectingOrbit, m::Integer; i_p::Integer=1)

Get the mth Chebyshev coefficient of the i_pth connecting orbit

set_am!(c::ConnectingOrbit, m::Integer, am::AbstractArray; i_p::Integer=1)

Set the mth Chebyshev coefficient of the i_pth connecting orbit

Base.getindex(c::ConnectingOrbit, i_A::Integer, i_p::Integer)

Wrapper for get_am.

Base.setindex!(c::ConnectingOrbit, X::AbstractArray, i_A::Integer, i_p::Integer)

Wrapper for set_am!

evaluate(c::ConnectingOrbit, x::AbstractArray; i_p::Integer=1)

Evaluate the i_pth connecting orbit at a set of points x[j] in [0,1]

evaluate(c::ConnectingOrbit, x::Number; i_p::Integer=1)

Evaluate the i_pth connecting orbit at a point x in [0,1]

deval(c::ConnectingOrbit, x::Number; i_p::Integer=1)

Evaluate the derivative of the i_pth connecting orbit at a point x in [0,1]

area(c::ConnectingOrbit; origin::AbstractArray=[0.0,0.0],
     i_p::Integer=1, rtol=1e-8)

Get the "area" of a connecting orbit ∫₋₁¹ (x dy/ds - y dy/ds) ds

linear_initial_connecting_orbit(x0, x1, Na::Integer)

Initialize a ConnectingOrbit of degree Na using the two hyperbolic orbits x0 and x1 of the form xn = [xn0, F(xn0), ..., F^(p-1)(xn0)].

gn_connecting!(c::ConnectingOrbit, FJ::Function, ends::AbstractArray;
               Ns = 100, maxiters = 10, rtol = 1e-8, verbose=false)

Find a connecting orbit using Gauss Newton with linesearch.

Arguments:

• c::ConnectingOrbit: An initial connecting orbit guess, see linear_initial_connecting_orbit
• FJ::Function: Function defined on R² with signature F(x), J(x) = FJ(x) where F is the symplectic map and J = dF/dx
• ends::AbstractArray: A 2p × 2 matrix containing the end periodic orbits [x00, x10; F(x00), F(x10); ... ; F^(p-1)(x00), F^(p-1)(x10)]
• Ns = 100: Number of quadrature nodes for the least squares
• maxiters = 10: Maximum number of Gauss Newton steps
• rtol = 1e-8: Stopping tolerance

KernelLabel

Constructors:

• KernelLabel(x::AbstractArray, c::AbstractVector, σ::Number; kernel::Symbol=:SquaredExponential)
• KernelLabel(x::AbstractArray, c::AbstractVector, σ::AbstractArray; kernel::Symbol=:SquaredExponential)

An approximately invariant kernel label function.

Constructor elements:

• x: The d × 2N array of interpolating points, where x[:, n+N] = F(x[:, n]) for n<=N and some symplectic map F.
• c: The length 2N vector of coefficients of the kernel function
• σ: The length scale of the kernel (if a vector, length scale is different in different directions)
• kernel: The type of kernel used

The current supported kernels are:

• :SquaredExponential: The squared exponential kernel K(x, y) = exp(-|x-y|²)
• :InverseMultiquadric: The β=1 Inverse Multiquadric kernel K(x, y) = (1+|x-y|²)^(-1/2)
• :FourierSE: The Fourier × Cartesian squared exponential kernel K(x, y) = exp(-sin²(x₁-y₁) - (x₂-y₂)²)
• :Fourier: The Fourier × Fourier squared exponential kernel K(x, y) = exp(-sin²(x₁-y₁) - sin²(x₂-y₂))

get_matrix(k::KernelLabel, x::AbstractArray)

Get the kernel matrix K[i,j] = K(k.x[:,i], x[:,j])

evaluate(k::KernelLabel, x::AbstractArray)

Evaluate the kernel matrix at the columns of x

window_weight(xs::AbstractArray, lims::AbstractVector, α::Number;
                   ind::Integer=2)

A simple boundary weighting function for kernel_eigs and kernel_bvp. Returns a weight that is approximately 1 outside of lims and 0 inside via the sum of two logistic functions.

Arguments:

• xs: A d × N array of points, where d is the size of the phase space and N is the number of points.
• lims: A 2-vector giving the interval where the weight is approximately 0
• α: The length scale of the logistic functions (typically small relative to the size of the domain)
• ind=2: The index over which the window is applied. Defaults to 2 for maps on T×R (such as the standard map)

Output:

• w: A N-vector with the window weights

rectangular_window_weight(xs::AbstractArray, xlims::AbstractVector,
                          ylims::AbstractVector, α::Number)

A simple boundary weighting function for kernel_eigs and kernel_bvp. Returns a weight that is approximately 1 outside of xlims × ylims and 0 inside via the a function of logistic functions.

Arguments:

• xs: A d × N array of points, where d is the size of the phase space and N is the number of points.
• xlims and ylims: 2-vectors giving the rectangle where the weight is approximately 0
• α: The length scale of the logistic functions (typically small relative to the size of the domain)

Output:

• w: A N-vector with the window weights

kernel_sample_F(F::Function, N::Integer, xb::AbstractVector,
                yb::AbstractVector)

Sobol sample N points in the rectangle xb × yb. Then, evaluate F at each point. Input can be used for kernel_eigs and kernel_bvp

Arguments:

• F: The symplectic map from the state space to itself
• N: The number of points to be sampled
• xb × yb: The 2-vector bounds of the rectangle to be sampled

Output:

• xs: A 2 × 2N array of samples compatable with kernel_eigs and kernel_bvp, of the form [x_1, F(x_1), x_2, F(x_2), ...]

kernel_eigs(xs::AbstractArray, ϵ::Number, nev::Integer, σ::Number,
            boundary_weights::Vector; kernel::Symbol=:SquaredExponential,
            zero_mean = false, check = 1)

Solve the the invariant eigenvalue problem, given by the Rayleigh quotient
> min_c (‖GKc‖² + ‖Kc‖²_bd + ϵ‖c‖²_k)/‖Kc‖²
where

• ‖GKc‖² is a norm penalizing invariance (and possible a non-zero mean)
• ‖Kc‖²_bd is a norm penalizing boundary violation
• ‖c‖²_k is the smoothing kernel norm
• ‖Kc‖² is the ℓ² norm of the points

The eigenvalue problem is solved via Arpack.jl

Arguments:

• xs: interpolation points of size d × 2N, where xs[:, N+1:2N] = F.(xs[:, 1:N])
• ϵ: Amount of regularization
• nev: Number of eigenvalues to find
• σ: Kernel width
• boundary_weights: Boundary weighting vector, should be positive and O(1) at points x where one wants |k(x)| << 1
• kernel: Type of kernel to interpolate (see KernelLabel)
• zero_mean = false: Set to true to add a constraint that encourages k to have a zero mean. Useful when xs are sampled on an invariant set and boundary_weights=0
• check = 1: See Arpack.eigs. If 1, return all of the converged eigenvectors and eigenvalues. If 0, throw an error instead.

Output:

• λs: The eigenvalues
• vs: The eigenvectors
• k: The kernel label. Use set_c!(k, vs[:,n]) to load the nth eigenvector into k. By default, k stores the lowest order eigenvector.

kernel_bvp(xs::AbstractArray, ϵ::Number, σ::Number,
           boundary_weights::AbstractVector,
           boundary_values::AbstractVector;
           kernel::Symbol=:SquaredExponential, residuals::Bool=true)

Solve the the invariant boundary value least-squares problem
> min_c ‖GKc‖² + ‖Kc - h_{bd}‖²_bd + ϵ‖c‖²_k = R_inv + R_bd + R_eps
where

• ‖GKc‖² is a norm penalizing invariance (and possible a non-zero mean)
• ‖Kc - h_{bd}‖²_bd is a norm penalizing the function from violating the boundary condition
• ‖c‖²_k is the smoothing kernel norm

Arguments:

• xs: interpolation points of size d × 2N, where xs[:, N+1:2N] = F.(xs[:, 1:N])
• ϵ: Amount of regularization
• σ: Kernel width
• boundary_weights: Length 2N boundary weighting vector, should be positive and O(1) at points x where one wants |k(x)| << 1
• boundary_values: Length 2N boundary value vector, indicating the value the function should take at each point
• kernel=:SquaredExponential: Type of kernel to interpolate (see KernelLabel)
• residuals=true: True if you want the problem to return residuals.

Output:

• k: The kernel function

(if residuals=true)

• R: The total residual
• R_bd: The boundary condition residual
• R_inv: The invariance residual
• R_eps: The smoothness residual

get_energies(k::KernelLabel; W = 0. * I)

Get the relevant energies of a kernel label k with boundary weighting matrix W.

Output energies:

• EK = ‖c‖²_k is the smoothing kernel norm
• EInv = ‖GKc‖² is a norm penalizing invariance
• Ebd = ‖Kc‖²_W is a norm penalizing boundary violation
• EL2 = ‖Kc‖² is the ℓ² norm of the points

kernel_birkhoff(xs::AbstractArray, fs::AbstractVector, ϵ::Number, μ::Number;
                σ::Number=0., kernel::Symbol=:SquaredExponential,
                boundary_points::AbstractVector = [])

This function is mostly deprecated. The solutions to the infinitely discretized limit of this problem (the Birkhoff average) don't live in the native space. So, the results can be odd, hard to interpret, and wiggly. Use at your own risk.

Find the "Birkhoff average" of a function using the kernel approach. Solves the least-squares problem
> min_c μ⁻¹(‖GKc‖² + ‖Kc‖²_bd) + ϵ‖c‖²_k + ‖Kc - f‖²
where

• ‖GKc‖² is a norm penalizing invariance
• ‖c‖²_k is the smoothing kernel norm
• ‖Kc - f‖² is a least-squares fit norm
• ‖Kc‖²_bd is a norm penalizing boundary violation

Arguments:

• xs: interpolation points of size d × 2N, where xs[:, N+1:2N] = F.(xs[:, 1:N])
• fs: function values at points of size N
• ϵ: Amount of regularization (can be set to zero)
• μ: Weighting of invariance penalization to fit (should be small, e.g. 1e-6)
• σ: Scale of the problem
• kernel: Type of kernel to interpolate (see KernelLabel)
• boundary_points: A list of indices of points on the boundary

Output:

• k: A KernelLabel object

ContFrac(a::AbstractArray)

A continued fraction in (0,1]. The member a is a vector of the continued fraction coefficients so that ω = 1/(a[1] + 1/(a[2] + ...))

big_cont_frac_eval(c::ContFrac)

Find a high precision representation of the continued fraction as a BigFloat. To set the precision, use precision.

evaluate(c::ContFrac)

Find a floating point representation of the continued fraction.

partial_frac(c::ContFrac, n::Integer)

Find the nth convergent of the continued fraction c as a rational number.

denoms(c::ContFrac)

Return a list of the denominators of the continued fractions of c

big_partial_frac(c::ContFrac, n::Integer)

Find the nth convergent of the continued fraction c as a rational number as BigInts. Set the precision with precision.

big_denoms(c::ContFrac)

Return a list of the denominators of the continued fractions of c as BigInts. Set the precision with precision.

wba_weight(d, N)

Weights used in the weighted Birkhoff average, returned as a diagonal matrix. Returns as Diagonal([w_1, ..., w_1, ..., w_N, ..., w_N]), where each coefficient is repeated d times (this is used for vector MPE and RRE as a least-squares weighting)

weighted_birkhoff_average(hs::AbstractMatrix)

Finds a weighted Birkhoff average of a sequence of vector observations. The array input is assumed to be of size d × N, where the average is performed over the second index.

weighted_birkhoff_average(hs::AbstractVector)

Finds a weighted Birkhoff average of a sequence of scalar observations.

doubling_birkhoff_average(h::Function, F::Function, x0::AbstractVector;
                          tol::Number = 1e-10, T_init::Integer=10,
                          T_max::Integer=320)

Find the weighted Birkhoff ergodic average of an observable function h over a trajectory of the map F adaptively.

Arguments:

• h: A function from Rᵈ to Rⁿ, where d is the dimension of the state space and n is the dimension of the observation
• F: The (symplectic) map: a function from Rᵈ to Rᵈ
• x0: The initial point of the trajectory in Rᵈ
• tol: The tolerance by which convergence is judged. If the average does not change by more than tol over a doubling, the function returns
• T_init: The initial length of the trajectory considered
• T_max: The maximum trajectory length considered

Output:

• ave: The average of h over the trajectory
• xs: The trajectory
• hs: The value of h on the trajectory
• conv_flag: true iff the averages converged

mutable struct BRREsolution

See also: adaptive_birkhoff_extrapolation), birkhoff_extrapolation, get_w0!, adaptive_get_torus!, save_rre, load_rre

The solution object for the Birkhoff RRE Framework. Contains:

adaptive_birkhoff_extrapolation) and birkhoff_extrapolation output: - h: Observable function (not included in save and load) - F: Symplectic map, assumed to act from R^n -> R^n for some n (not included in save and load) - D: Output dimension of observable - xs: Trajectory used of points in phase space - hs: Trajectory of observables - K: Number of Birkhoff RRE filter unknowns - c: Filter coefficients - h_ave: Birkhoff average - resid_rre: Reduced Rank Extrapolation residual

get_w0! output: - d: Torus dimension - w0: Rotation vector (of form [1//p, w1, ..., wd], where p is the island period) - resid_w0: Log posterior residual for rotation vector

adaptive_get_torus! output: - tor: The output FourierTorus object - resid_tor: The validation residual for the torus

save_rre(file::AbstractString, sol::BRREsolution)

Save an RRE solution object to a file (not including the evaluation and observation function F and h). Currently only saves to JLD2 files.

load_rre(file::AbstractString)

Load an RRE solution object from a file (saved by save_rre).

vector_rre_backslash(x::AbstractArray, K::Integer; ϵ=0.0)

Applies Birkhoff vector RRE to a sequence x_n = x[:, n]

Arguments:

• x: The sequence
• K: The number of unknowns in the filter
• ϵ: A regularization parameter (typically zero gives best results)

Output:

• c: The learned filter
• sums: The partial sums (Birkhoff averages) obtained by applying c to the sequence
• resid: The residuals of the least-squares problem. Taking the norm gives an overall measure of the fit.

birkhoff_extrapolation(h::Function, F::Function, x0::AbstractVector,
                       N::Integer, K::Integer; 
                       x_prev::Union{AbstractArray,Nothing}=nothing,
                       weighted::Bool=false)

As input, takes an initial point x0, a symplectic map F, and an observable h (can choose the identity as a default h = (x)->x). Then, the method

1. Computes a time series xs[:,n+1] = Fⁿ(x0)
2. Computes observable series hs[:,n] = h(xs[:,n])
3. Performs sequence extrapolation (RRE or MPE) on hs to obtain a model c of length 2K+1 (with K unknowns), the extrapolated value applied to each window, and a residual
4. Returns a BRREsolution object

Use x_prev if you already know part of the sequence, but do not know the whole thing.

adaptive_birkhoff_extrapolation(h::Function, F::Function,
                x0::AbstractVector; rtol::Number=1e-10, Kinit = 20,
                Kmax = 100, Kstride=20, Nfactor::Integer=5)

Adaptively applies birkhoff_extrapolation to find a good enough filter length K, where "good enough" is defined by the rtol optional argument.

Arguments:

• h: The observable function (often the identity function h = (x)->x works)
• F: The symplectic map
• x0: The initial point of the trajectory
• Kinit, Kmax, Kstride: Increases filter size K as Kinit:Kstride:Kmax. For more complicated and higher dimensional tori, typically it is better to increase these values.
• rtol: Required tolerance for convergence (lower is better; inexact maps often require a looser tolerance)
• Nfactor: How rectangular the extrapolation least-squares problem is. Must be >=1.

Outputs:

• sol: A BRREsolution object

get_w0!(sol::BRREsolution, Nw0::Integer; Nsearch::Integer=30, Nsearchcutoff::Number=1e-6, 
    gridratio::Number=0., Nkz::Integer=50, sigma_w::Number = 1e-10, rattol::Number = 1e-8, 
    maxNisland::Number = 10)

Find the rotation vector from an output of Birkhoff RRE. Stores the output in the input sol object.

Input:

• sol: The output of a Birkhoff RRE extrapolation
• Nw0: The dimension of the torus
• (Nsearch,gridratio): Nsearch is the number of independent roots of c (i.e. ±ω are the same root) that we consider for choosing ω0. gridratio gives the size of grid we use for the discrete optimization to find ω0. If the torus is dominated by higher harmonics, then one or both of these parameters may need to be increased. The default values (chosen by gridratio==0.) are 2 for the 1D case and 5 for greater than 1D.
• sigma_w: Frequency accuracy used in the Bayesian inference
• Nsearchcutoff: The cutoff prominence of the Fourier modes used to put a maximum on Nsearch. Lowering this value may result in a rotation vector when otherwise it was not found, but it may reduce accuracy.
• Nkz=50: Number of frequencies to use for finding the KZ basis. This value likely does not need to be tuned by the user.

adaptive_get_torus!(sol::BRREsolution; max_size::Number=2000, validation_fraction::Number=0.05, 
                    ridge_factor::Number=10)

Fit an invariant torus from a trajectory. The output is saved to the sol input. See kam_residual for an a posteriori validation of the circle.

Inputs:

• sol: An input BRREsolution object. It is assumed that the rotation vector has already been computed using get_w0!.
• max_size: Maximum number of Fourier modes to use for the torus parameterization
• validation_fraction: Fraction of the trajectory which is used for computing the validation error of the trajectory.
• ridge_factor: Height of the ridges that the discrete gradient descent optimization can traverse.

sum_stats(sums)

Given a sequence of sums applied to a filter (an output of Birkhoff RRE extrapolation), find the average and standard deviation of the sums. Can be used as a measure of how "good" the filter is.

get_sum_ave(hs, c)

Given a signal hs (of size either d×N or N×1) and a filter c of size M×1, compute the average of the M-N windows of the filter applied to the signal.

lyapunov_exponent(x0::AbstractVector, FJ::Function, N::Integer)

Compute the lyapunov exponent of the map FJ after N iterates of the map.

Input:

• x0: Initial point
• FJ: Returns the symplectic map evaluation F and Jacobian J = dF/dx, i.e. has the function signature F, J = FJ(x)
• N: The number of iterations of the map to be used

Output:

• λ: The Lyapunov exponent λ = 1/N * log(σ_max(dF^N(x)/dx))

lyapunov_exponents(x0::AbstractVector, FJ::Function, N::Integer)

Compute the lyapunov exponents of the map FJ over N iterates of the map.

Input:

• x0: Initial point
• FJ: Returns the symplectic map evaluation F and Jacobian J = dF/dx, i.e. has the function signature F, J = FJ(x)
• N: The number of iterations of the map to be used

Output:

• λs: The Lyapunov exponent λ = 1/n * log(σ_max(dF^n(x)/dx)) for 1<=n<=N

standard_map_F(k)

Returns the Chirikov standard map with parameter k.

standard_map_F(K::AbstractMatrix, delta::AbstractVector)

Returns the a higher dimensional standard map with matrix K and offset delta,i.e. y+ = y - (K/(2π))*sin(2π(x+δ)) x+ = x + y+ where (x,y) take up the (1:d,d+1:2d) entries of the input.

standard_map_FJ(k)

Returns a function FJ that returns the Chirikov standard map and its derivative (i.e. F, dFdx = FJ(x)) with parameter k.

standard_map_FJ(K::AbstractMatrix, delta::AbstractVector)

Returns a function FJ that returns the higher dimensional standard map and its derivative (i.e. F, dFdx = FJ(x)) with parameter k (see standard_map_F(::AbstractMatrix,::AbstractVector)).

polar_map(d, z0)

Returns the d-dimensional polar map
> h:(θ,z)->((z-z0)cos(2πθ), (z-z0)sin(2πθ))
as well as its derivative HJ, its inverse hinv, and its inverse derivative HJinv. Useful for applying extrapolation methods to maps on Tᵈ×Rᵈ.

polar_map([d];z0 = -0.5)

Returns the polar map
> h:(θ,z)->((z-z0)cos(2πθ), (z-z0)sin(2πθ))
as well as its derivative HJ, its inverse hinv, and its inverse derivative HJinv. Useful for applying extrapolation methods to maps on T×R. Default value of z0 is useful for the standard map on T×[0,1] with k=0.7.

Plots.plot(z::InvariantCircle; kwargs...)

Creates a plot of the invariant circle z. See Plots.plot!(z::InvariantCircle) for a list of keyword arguments.

Plots.plot(p::Plots.Plot, z::InvariantCircle; kwargs...)

Plots the invariant circle z on p.

kwargs:

• N::Integer: Number of θ points used to plot each circle
• i_circle::Integer: Which period of an island chain to plot. Default value of 0 plots all circles of an island chain.
• label, color, linewidth, linestyle: see Plots.jl

parametric_plot(z::InvariantCircle; N::Integer=200, i_circle::Integer=1,
                linewidth=1, linestyle=:solid, label1="x", label2="y",
                xlabel="θ", plot_min_dθs=true, markersize=5)

Parametric plot of an invariant circle.

kwargs:

• N::Integer: Number of θ points used to plot each circle
• i_circle::Integer: Which period of an island chain to plot

Plots.plot(c::ConnectingOrbit; kwargs...)

Creates a plot of the connecting orbit c. See Plots.plot!(c::ConnectingOrbit) for a list of keyword arguments.

Plots.plot!(p::Plots.Plot, c::ConnectingOrbit; kwargs...)

Creates a plot of the connecting orbit c.

kwargs:

• N::Integer: Number of points used to plot each orbit
• i_circle::Integer: Which period of a connecting orbit chain to plot. Default value of 0 plots all connecting orbits of a island chain.
• label, color, linewidth: see Plots.jl

CairoMakie.lines!(ax, z::InvariantCircle; N::Integer=100, color=nothing, i_circle::Integer=0, linewidth=1)

Plot the invariant circle z on the CairoMakie axis ax.

Arguments:

• ax: CairoMakie Axis object
• z: The circle in R² to be plotted
• N: Number of points to plot
• i_circle: Which invariant circle of an island to plot. If 0, plot all
• color, linewidth: see CairoMakie.lines!

lines_periodic!(ax, z::InvariantCircle, hinv::Function; N::Integer=100,
                color=nothing, i_circle::Integer=0, linewidth=1)

Useful for plotting with invariant circles on the torus. I.e., if F : T×R→T×R, and one finds an invariant circle of z(θ+τ) = (h∘F)(z(θ)) where h : T×R→R², this plots h⁻¹∘z, the original invariant circle.

Arguments:

• ax: CairoMakie Axis object
• z: The circle in R² to be plotted
• hinv: The map to the torus by the real numbers h⁻¹ : R²→T×R
• N: Number of points to plot
• i_circle: Which invariant circle of an island to plot. If 0, plot all
• color, linewidth, label: see CairoMakie.lines!

plot_on_grid(x::AbstractVector, y::AbstractVector, k::KernelLabel;
             kwargs...)

Create a filled contour plot of the kernel label k on the x × y grid.

kwargs:

• balance=true: If true, makes the maximum value of the color scale equal to the minimum
• size, fontsize: See CairoMakie.Figure
• xlabel, ylabel, title: See CairoMakie.Axis
• levels, linewidth: See CairoMakie.contour!
• clabel: See CairoMakie.Colorbar

Output:

• f: The CairoMakie Figure
• f_grid: The data used to make the figure

poincare_plot(xb::AbstractVector, yb::AbstractVector, F::Function,
              Ninit::Integer, Niter::Integer; size=(800, 800),
              fontsize=25, xlabel="x", ylabel="y", xlims = nothing,
              ylims=nothing, markersize=3, title="Poincare Plot")

Create a Poincare plot of a 2D map F in a rectangular region xb×yb. The Poincare plot uses Ninit trajectories of length Niter.

Output:

• f: The figure
• xs: The trajectories used to make the figure