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

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.

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.

vector_mpe_backslash(x::AbstractArray, K::Integer)

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

vector_mpe_backslash(x::AbstractArray, K::Integer)

Applies Birkhoff vector MPE to a sequence x_n = x[:, n] using the LSQR algorithm. This currently does not have preconditioning, and therefore is less accurate than vector_mpe_backslash.

Arguments:

• x: The sequence
• K: The number of unknowns in the filter
• c0: The initial guess of
• atol, btol: Tolerances. See IterativeSolvers.lsqr!

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

birkhoff_extrapolation(h::Function, F::Function, x0::AbstractVector,
                       N::Integer, K::Integer; iterative::Bool=false,
                       x_prev::Union{AbstractArray,Nothing}=nothing,
                       rre::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 as c, sums, resid, xs, hs[, history] where history is a diagnostic from the iterative solver that only returns when iterative=true

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, iterative::Bool=false,
                Nfactor::Integer=1, rre::Bool=true)

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 (can choose the identity as a default h = (x)->x)
• F: The symplectic map
• x0: The initial point of the trajectory
• rtol: Required tolerance for convergence (inexact maps often require a looser tolerance)
• Kinit: The length of the initial filter
• Kmax: The maximum allowed filter size
• Kstride: The amount K increases between applications of birkhoff_extrapolation
• iterative: Whether to use an iterative method to solve the Hankel system in the extrapolation step
• Nfactor: How rectangular the extrapolation algorithm is. Must be >=1.
• rre: Turn to true to use reduced rank extrapolation instead of minimal polynomial extrapolation.

Outputs:

• c: Linear model / filter
• sums: The extrapolated value applied to each window
• resid: The least squares residual
• xs: A time series x[:,n] = Fⁿ(x[:,0])
• hs: The observations h[:,n] = h(x[:,n])
• rnorm: The norm of resid
• K: The final degree of the filter
• history: Returned if iterative=true. The history of the final LSQR iteration

sum_stats(sums)

Given a sequence of sums applied to a filter (an output of invariant circle 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.

get_circle_info(hs::AbstractArray, c::AbstractArray; rattol::Number=1e-8,
                ratcutoff::Number=1e-4, max_island_d::Integer=30)

Get a Fourier representation of an invariant circle from the observations hs and the learned filter c (see adaptive_invariant_circle_model and invariant_circle_model to find the filter). See get_circle_residual for an a posteriori validation of the circle.

Optional Arguments:

• rattol: Roots are judged to be rational if |ω-m/n|<rattol
• ratcutoff: Relative prominence needed by a linear mode to qualify as "important" for deciding whether the sequence is an island
• max_island_d: Maximum denominator considered for islands.
• modetol: Tolerance (typically less than 1) used for determining the number of Fourier modes. Higher numbers result in more Fourier modes at the expense of robustness of the least-squares system. Decrease it (e.g. to 0.5) for noisy data

Output:

• z: An invariant circle of type FourierCircle

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_FJ(k)

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

polar_map(;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