API Reference

struct TestInput

Container for test parameters and objective function.

Fields:

• dim::Int: Problem dimension
• center::Vector{Float64}: Center point of search region
• GN::Union{Int,Nothing}: Grid size (optional)
• prec::Union{Tuple{Float64,Float64},Nothing}: Precision parameters (α,δ)
• tolerance::Union{Float64,Nothing}: Convergence tolerance
• noise::Union{Tuple{Float64,Float64},Nothing}: Noise parameters
• sample_range::Union{Float64,Vector{Float64},Nothing}: Sampling radius around center (scalar or per dimension)
• reduce_samples::Union{Float64,Nothing}: Sample reduction factor
• degree_max::Union{Int, Nothing}: Maximum polynomial degree
• objective: Callable objective (Function or callable struct like TolerantObjective)

Constructor(T::TestInput, degree; kwargs...) -> ApproxPoly

Construct a polynomial approximation of the objective function using discrete least squares.

This is the main entry point for creating polynomial approximations in Globtim. The function samples the objective on a tensorized grid of Chebyshev or Legendre nodes and fits a polynomial of the specified degree.

Arguments

• T::TestInput: Test input specification containing the objective function and domain
• degree::Int: Maximum degree of the polynomial approximation

Keyword Arguments

• verbose::Int=0: Verbosity level (0=silent, 1=basic info, 2=detailed)
• basis::Symbol=:chebyshev: Basis type (:chebyshev or :legendre)
• precision::PrecisionType=Float64Precision: Precision type for coefficients (default; use RationalPrecision for exact/symbolic paths — to_exact_monomial_basis, msolve)
• normalized::Bool=false: Whether to normalize the polynomial
• power_of_two_denom::Bool=false: Use power-of-two denominators for rationals
• grid::Union{Nothing,Matrix{Float64}}=nothing: Pre-generated grid matrix (rows are points)

Returns

• ApproxPoly: Polynomial approximation object containing:
  • coeffs: Coefficient matrix
  • nrm: L2-norm approximation error over the domain
  • scale_factor: Scaling factors used for the domain
  • Additional metadata about the approximation

Notes

• The approximation error (pol.nrm) provides a measure of approximation quality
• Higher degrees generally reduce approximation error but increase computational cost
• Chebyshev basis is the default (well-understood approximation properties; avoids ill-conditioning of monomial basis)
• The function automatically handles both uniform and non-uniform domain scaling

Examples

# 1D function example with scalar input
f1 = x -> sin(x)
TR = TestInput(f1, dim=1, center=[0.0], sample_range=10.0)
pol = Constructor(TR, 8)
println("L2-norm error: ", pol.nrm)

# Basic usage with default Chebyshev basis
f = Deuflhard
TR = TestInput(f, dim=2, center=[0.0, 0.0], sample_range=1.0)
pol = Constructor(TR, 8)
println("L2-norm error: ", pol.nrm)

# Using Legendre basis with higher verbosity
pol = Constructor(TR, 10, basis=:legendre, verbose=1)

# High precision approximation
pol = Constructor(TR, 12, normalized=true)

# Using a pre-generated anisotropic grid
grid_aniso = generate_anisotropic_grid([10, 5], basis=:chebyshev)
grid_matrix = convert_to_matrix_grid(vec(grid_aniso))
pol_aniso = Constructor(TR, 0, grid=grid_matrix)  # degree ignored when grid provided

solve_polynomial_system(x, n, d, coeffs; solver=:hc, kwargs...) -> Vector{Vector{Float64}} or Tuple

Critical point finder using polynomial system solving.

Find all critical points of a polynomial approximation by solving the gradient system ∇p(x) = 0. Two solver backends are available:

• :hc (default) — HomotopyContinuation.jl: numerical algebraic geometry via homotopy continuation. Fast, handles large systems, but may lose paths (miss solutions).
• :msolve — msolve binary: exact Gröbner basis computation over ℚ with certified real root isolation. Guaranteed to find all real solutions (no path loss), but may be slower for large systems.

Arguments

• x: Polynomial variables (from DynamicPolynomials)
• n::Int: Number of variables (dimension)
• d::Int: Polynomial degree
• coeffs: Coefficient matrix from polynomial approximation

Keyword Arguments

• solver::Symbol=:hc: Solver backend (:hc or :msolve)
• basis::Symbol=:chebyshev: Basis type (:chebyshev or :legendre)
• precision::PrecisionType=RationalPrecision: Precision type for coefficients
• normalized::Bool=true: Whether to use normalized basis polynomials
• power_of_two_denom::Bool=false: For rational precision, ensures denominators are powers of 2
• return_system::Bool=false: If true, also return the polynomial system information (:hc only)
• msolve_threads::Int=1: Number of threads for msolve (:msolve only)

Returns

• If return_system=false: Vector{Vector{Float64}} — Real solutions within [-1,1]ⁿ
• If return_system=true: Tuple containing:
  • Solutions vector
  • Tuple of (polynomial system, HC system, total solution count)

Notes

• Only returns real solutions within the domain [-1,1]ⁿ
• Complex solutions and solutions outside the domain are filtered out
• The number of solutions can vary significantly based on the polynomial degree

Examples

using DynamicPolynomials
@polyvar x[1:2]

# Using HomotopyContinuation (default)
crit_pts = solve_polynomial_system(x, 2, 8, coeffs)

# Using msolve (exact, no path loss)
crit_pts = solve_polynomial_system(x, 2, 8, coeffs; solver=:msolve)

# msolve with 4 threads
crit_pts = solve_polynomial_system(x, 2, 8, coeffs; solver=:msolve, msolve_threads=4)

solve_polynomial_system(x, pol::ApproxPoly; kwargs...)

Convenience method that automatically extracts dimension and degree from an ApproxPoly object.

Arguments

• x: Polynomial variables (from DynamicPolynomials)
• pol::ApproxPoly: Polynomial approximation object
• kwargs...: Additional keyword arguments passed to the main method

Returns

Same as the main solve_polynomial_system method.

Example

f = x -> sin(x)
TR = TestInput(f, dim=1, center=[0.0], sample_range=10.)
pol = Constructor(TR, 8)
@polyvar x
solutions = solve_polynomial_system(x, pol)  # No need to specify dim and degree

process_crit_pts(
    real_pts::Vector{<:AbstractVector},
    f::Function,
    TR::TestInput;
    skip_filtering::Bool = false,
    kwargs...
)::DataFrame

Process critical points in n-dimensional space and return a DataFrame. Points are automatically filtered to the [-1,1]^n hypercube (unless skip_filtering is true) and transformed according to the TestInput parameters.

Arguments

• real_pts: Vector of points in n-dimensional space
• f: Function to evaluate at each point
• TR: TestInput struct containing dimension, center, and sample range information
• skip_filtering: If true, skips the [-1,1] bounds filtering (default: false)
• kwargs...: Additional arguments for future extensions

Returns

• DataFrame with columns x1, x2, ..., xn (for n dimensions) and z (function values)

adaptive_refine(f, bounds::Vector{Tuple{Float64, Float64}},
                degree; kwargs...)

Main adaptive refinement loop with parallel processing.

Arguments

• f: Callable to approximate (any callable works, including TolerantObjective)
• bounds: Domain bounds as vector of (min, max) tuples
• degree: Polynomial degree (or degree specification)

Keyword Arguments

• l2_tolerance::Float64=1e-6: Target L2 error tolerance
• max_depth::Int=10: Maximum subdivision depth
• max_leaves::Int=1000: Maximum number of leaf subdomains
• optimize_cuts::Bool=true: Whether to optimize cut positions
• parallel::Bool=true: Whether to use CPU parallel processing
• basis::Symbol=:chebyshev: Basis type
• verbose::Bool=false: Print progress information
• phase_callback::Union{Function,Nothing}=nothing: Called with (f, :refine, 0) at start. Use with TolerantObjective to set solver tolerances per phase.
• enable_p_refinement::Bool=false: Enable hp-refinement (try higher degree before splitting)
• max_degree::Int=40: Maximum polynomial degree for p-refinement
• degree_step::Int=6: Degree increment per p-refinement step
• cond_threshold::Float64=1e14: Maximum Vandermonde condition number for p-refinement
• anisotropic_degree::Union{Nothing,NamedTuple}=nothing: opt-in E2 ρk-driven anisotropic degree (jl9z.7). nothing ⇒ isotropic (behavior unchanged). When set, each non-converged leaf's base fit is used as the E2 probe: per-axis Bernstein radii ρk are estimated and the leaf is refit at a data-adaptive per-dim degree (small on smooth/sloppy axes floored at floor_degree, larger on rough axes). Fires at most once per leaf; further refinement preserves the anisotropic shape. The NamedTuple is splatted into choose_per_dim_degree_lsfit — e.g. (; c=4.0, floor_degree=2, max_degree=8). Three keys are reserved for the Stage-2 active-subspace fallback and stripped before the splat: active_fallback::Bool=false — when the LS ρk spectrum is blind on EVERY axis (Stage 1 would blast the leaf to isotropic `maxdegree), probe the gradient covariance instead; an unambiguous spectrum rotates the leaf to its active frame with spectrum-derived per-dim degrees, an ambiguous one falls through to the ordinary bump/split ladder.fallbackncells::Int=3,fallbackh::Float64=0.01— probe grid (costs(n+1)·ncellsⁿ` extra objective evaluations on the leaf, only when triggered).

Returns

• SubdivisionTree with refined subdomains

Example

using Globtim

f(x) = sum(x.^2)
bounds = [(-1.0, 1.0), (-1.0, 1.0)]

# Default: relative L2 tolerance (0.03 = reliable CP recovery threshold)
tree = adaptive_refine(f, bounds, 4)

# Explicit relative tolerance
tree = adaptive_refine(f, bounds, 4; l2_tolerance=0.01, tolerance_mode=:relative)

# Absolute tolerance (backward compatible)
tree = adaptive_refine(f, bounds, 4; l2_tolerance=1e-4, tolerance_mode=:absolute)

# hp-adaptive: start at degree 10, bump up to 40 before splitting
tree = adaptive_refine(f, bounds, 10; l2_tolerance=1e-4, tolerance_mode=:absolute,
                       enable_p_refinement=true, max_degree=40, degree_step=6)

two_phase_refine(f, bounds::Vector{Tuple{Float64, Float64}},
                 degree; kwargs...)

Two-phase adaptive refinement: coarse balancing pass, then accuracy refinement.

Phase 1 subdivides until errors are relatively balanced (no stragglers). Phase 2 refines to meet the final tolerance.

Arguments

• f: Callable to approximate (any callable works, including TolerantObjective)
• bounds: Domain bounds
• degree: Polynomial degree

Keyword Arguments

• coarse_tolerance::Float64: Phase 1 tolerance (should be looser than fine)
• fine_tolerance::Float64: Phase 2 final tolerance
• balance_threshold::Float64=3.0: Phase 1 stops when max/min error ratio < this
• max_depth::Int=10: Maximum subdivision depth
• max_leaves::Int=1000: Maximum leaves
• parallel::Bool=true: Use CPU parallel processing
• basis::Symbol=:chebyshev: Basis type
• verbose::Bool=false: Print progress
• phase_callback::Union{Function,Nothing}=nothing: Called with (f, :coarse, 0) at Phase 1 start and (f, :fine, 0) at Phase 2 start. Use with TolerantObjective to switch solver tolerances between phases.
• enable_p_refinement::Bool=false: Enable hp-refinement (try higher degree before splitting)
• max_degree::Int=40: Maximum polynomial degree for p-refinement
• degree_step::Int=6: Degree increment per p-refinement step
• cond_threshold::Float64=1e14: Maximum Vandermonde condition number for p-refinement
• logger::Union{Metrics.MetricsLogger,Nothing}=nothing: When non-nothing, emit JSONL phase rows (two_phase_refine_start, two_phase_phase1_done, two_phase_phase2_start, two_phase_refine_done) plus one leaf row per processed leaf with stage="tree-build". Mirror of adaptive_refine's logger kwarg.
• leaf_extra_fn::Union{Function,Nothing}=nothing: Optional (sd, leaf_id) -> Dict that supplies per-leaf data the library can't know (HC-side stats, predicate counters); the returned entries are folded into each leaf row's extra.

Returns

• SubdivisionTree with refined subdomains

enhanced_bfgs_refinement(initial_points::Vector{Vector{Float64}},
                       initial_values::Vector{Float64},
                       orthant_labels::Vector{String},
                       objective_function::Function,
                       config::BFGSConfig = BFGSConfig();
                       expected_minimum::Union{Vector{Float64}, Nothing} = nothing)

Perform enhanced BFGS refinement with comprehensive hyperparameter tracking.

Arguments

• initial_points: Vector of starting points for refinement
• initial_values: Function values at initial points
• orthant_labels: Labels identifying orthants/regions
• objective_function: The objective function to minimize
• config: BFGSConfig with hyperparameters
• expected_minimum: Expected global minimum for distance tracking (optional)

Returns

• Vector{BFGSResult}: Detailed results for each refinement

refine_with_enhanced_bfgs(df::DataFrame, objective_function::Function,
                         config::BFGSConfig = BFGSConfig();
                         expected_minima::Union{Vector{Vector{Float64}}, Nothing} = nothing)

Apply enhanced BFGS refinement to critical points in a DataFrame.

Arguments

• df: DataFrame with critical points (columns x1, x2, ..., z)
• objective_function: The objective function
• config: BFGSConfig with hyperparameters
• expected_minima: Known global minima for comparison (optional)

Returns

• DataFrame: Enhanced DataFrame with BFGS refinement results

determine_convergence_reason(result::Optim.OptimizationResults, tolerance_used::Float64, config::BFGSConfig)

Analyze Optim result to determine why optimization stopped.

Arguments

• result: Optimization result from Optim.jl
• tolerance_used: The tolerance that was used
• config: BFGSConfig structure

Returns

• Symbol: Convergence reason (:gradient, :ftol, :xtol, :iterations, etc.)

classify_critical_points(hessians::Vector{Matrix{Float64}}; 
                       tol_zero=1e-8, tol_pos=1e-8, tol_neg=1e-8)::Vector{Symbol}

Classify critical points based on Hessian eigenvalue structure.

Arguments

• hessians: Vector of Hessian matrices
• tol_zero: Tolerance for zero eigenvalues (degeneracy detection)
• tol_pos: Tolerance for positive eigenvalues
• tol_neg: Tolerance for negative eigenvalues

Returns

Vector{Symbol}: Classification for each point (:minimum, :maximum, :saddle, :degenerate, :error)

compute_hessians(f, points::Matrix{Float64})::Vector{Matrix{Float64}}

Compute Hessian matrices at specified points using ForwardDiff automatic differentiation.

Arguments

• f: Objective function to analyze (any callable)
• points: Matrix where each row is a point (npoints × ndims)

Returns

Vector{Matrix{Float64}}: Hessian matrix for each point

analyze_basins(df::DataFrame, df_min::DataFrame, n_dims::Int, tol_dist::Float64)

Analyze basin of attraction properties for each unique minimizer.

Returns

Tuple{Vector{Int}, Vector{Float64}, Vector{Int}}:

• Basin sizes (point count for each minimizer)
• Average convergence steps for each minimizer
• Region coverage count for each minimizer

evaluate(poly::ApproxPoly, x::AbstractVector{<:Real})::Float64

Evaluate polynomial approximation at point x.

The point x should be in the original (unscaled) domain. The function internally scales x by poly.scale_factor to map to the [-1,1]^n reference domain where the orthogonal basis polynomials are defined.

Arguments

• poly::ApproxPoly: The polynomial approximation object
• x::AbstractVector{<:Real}: Point at which to evaluate (in original domain)

Returns

• Float64: Value of the polynomial approximation at x

Example

poly = MainGenerate(f, 2, (:one_d_for_all, 8), 0.05, 0.95, 1.5, 1.0)
val = evaluate(poly, [0.5, 0.3])

evaluate(poly::ApproxPoly, X::AbstractMatrix{<:Real})::Vector{Float64}

Evaluate polynomial approximation at multiple points.

Arguments

• poly::ApproxPoly: The polynomial approximation object
• X::AbstractMatrix{<:Real}: Points as rows (npoints × ndims)

Returns

• Vector{Float64}: Values at each point

gradient(poly::ApproxPoly, x::AbstractVector{<:Real})::Vector{Float64}

Compute gradient of polynomial approximation at point x using automatic differentiation.

The gradient is computed in the original (unscaled) domain coordinates.

Arguments

• poly::ApproxPoly: The polynomial approximation object
• x::AbstractVector{<:Real}: Point at which to compute gradient (in original domain)

Returns

• Vector{Float64}: Gradient vector ∇p(x)

Example

poly = MainGenerate(f, 2, (:one_d_for_all, 8), 0.05, 0.95, 1.5, 1.0)
grad = gradient(poly, [0.5, 0.3])

relative_l2_error(pol::ApproxPoly) -> Float64

Compute the relative L2 approximation error: ||f - p||_L2 / ||f||_L2.

The absolute L2 error (pol.nrm) is a quadrature-weighted norm of the residual on [-1,1]^n. This function normalizes it by the same weighted norm of the function values, giving a dimensionless ratio in [0, 1] for a good approximation.

Returns NaN if the function norm is zero (constant zero function).

sparsify_polynomial(pol::ApproxPoly, threshold::Real; mode=:relative, preserve_indices=[])

Set small coefficients to zero while tracking L²-norm impact.

Arguments

• pol: ApproxPoly to sparsify
• threshold: Threshold for zeroing coefficients
• mode: :relative or :absolute thresholding
• preserve_indices: Indices of coefficients to preserve

Returns

• NamedTuple with fields:
  • polynomial: Sparsified ApproxPoly
  • sparsity: Fraction of non-zero coefficients
  • zeroed_indices: Indices of zeroed coefficients
  • l2_ratio: L²-norm ratio (sparsified/original)
  • original_nnz: Original number of non-zero coefficients
  • new_nnz: New number of non-zero coefficients

analyze_sparsification_tradeoff(pol::ApproxPoly; thresholds=[1e-6, 1e-8, 1e-10, 1e-12])

Analyze sparsity vs accuracy tradeoffs for different thresholds.

Arguments

• pol: Polynomial to analyze
• thresholds: Array of thresholds to test

Returns

• Array of results for each threshold

truncate_polynomial_adaptive(poly, threshold::Real; relative::Bool=false)

Truncate polynomial coefficients using extended precision for accurate threshold comparison. This function is designed to work well with AdaptivePrecision polynomials.

Arguments

• poly: DynamicPolynomials.Polynomial with extended precision coefficients
• threshold::Real: Truncation threshold
• relative::Bool: If true, threshold is relative to largest coefficient

Returns

• Truncated polynomial with small coefficients removed
• Statistics about the truncation (number of terms removed, etc.)

to_exact_monomial_basis(pol::ApproxPoly; variables=nothing)

Convert a polynomial from orthogonal basis (Chebyshev/Legendre) to monomial basis using exact arithmetic.

Arguments

• pol::ApproxPoly: Polynomial approximation from Globtim
• variables: Array of polynomial variables (created automatically if not provided)

Returns

• DynamicPolynomials.Polynomial: Polynomial in monomial basis with exact coefficients

Example

TR = TestInput(x -> sin(x[1]), dim=1, center=[0.0], sample_range=1.0)
pol = Constructor(TR, 10, basis=:chebyshev)
@polyvar x
mono_poly = to_exact_monomial_basis(pol, variables=[x])

exact_polynomial_coefficients(f::Function, dim::Int, degree::Int; kwargs...)

Convenience function to get exact monomial coefficients directly from a function.

Arguments

• f::Function: Function to approximate
• dim::Int: Dimension of the input
• degree::Int: Maximum polynomial degree
• basis::Symbol = :chebyshev: Basis to use (:chebyshev or :legendre)
• center::Vector = zeros(dim): Center of approximation domain
• sample_range::Real = 1.0: Radius of approximation domain
• tolerance::Real = 0.5: Tolerance for approximation
• precision = Float64Precision: Arithmetic precision

Returns

• DynamicPolynomials.Polynomial: Polynomial in monomial basis

Example

f = x -> x[1]^2 + x[2]^2
mono_poly = exact_polynomial_coefficients(f, 2, 4, basis=:chebyshev)

generate_grid(grid_spec::Union{Int,Vector{Int}}, n_dims::Union{Int,Nothing}=nothing; basis=:chebyshev)

Unified interface for generating isotropic or anisotropic grids.

Arguments

• grid_spec: Either:
  • Int: Number of points per dimension (isotropic grid)
  • Vector{Int}: Number of points for each dimension (anisotropic grid)
• n_dims: Number of dimensions (only needed for isotropic case)
• basis: Node type (:chebyshev, :legendre, or :uniform)

Examples

# Isotropic 3D grid with 10 points per dimension
grid = generate_grid(9, 3)

# Anisotropic 3D grid
grid = generate_grid([9, 5, 3])

generate_anisotropic_grid(grid_sizes::Vector{Int}; basis=:chebyshev)

Generate an anisotropic grid with different number of points in each dimension.

Arguments

• grid_sizes::Vector{Int}: Number of points in each dimension (will generate grid_sizes[i] + 1 points in dimension i)
• basis::Symbol=:chebyshev: Choice of basis for node generation (:chebyshev, :legendre, or :uniform)

Returns

• Array of SVectors containing the grid points

Examples

# 2D grid with 5 points in x and 10 points in y
grid = generate_anisotropic_grid([4, 9], basis=:chebyshev)

# 3D grid with different resolution per axis
grid = generate_anisotropic_grid([10, 5, 3], basis=:legendre)

compute_l2_norm(poly::AbstractPolynomial, domain::AbstractDomain; n_points=20)

Compute the L²-norm of a polynomial over a given domain using discrete approximation.

Arguments

• poly: Polynomial in monomial basis
• domain: Domain of integration (box domain [-a,a]ⁿ supported)
• n_points: Number of grid points per dimension (default: 20)

Returns

• L²-norm value (numerical approximation)

compute_l2_norm_quadrature(f::Function, n_points::Vector{Int}, basis::Symbol=:chebyshev)

Compute the L2 norm of a function using Gaussian quadrature.

Arguments

• f: Function to compute L2 norm for. Should accept a vector input.
• n_points: Number of quadrature points in each dimension
• basis: Type of polynomial basis (:chebyshev, :legendre, :uniform)

Returns

• L2 norm value

Example

f = x -> exp(-(x[1]^2 + x[2]^2))
l2_norm = compute_l2_norm_quadrature(f, [10, 10], :chebyshev)

integrate_monomial(exponents::Vector{Int}, domain::BoxDomain)

Analytically integrate a monomial over a box domain.

Arguments

• exponents: Vector of exponents for each variable
• domain: Box domain [-a,a]ⁿ

Returns

• Integral value

Example

integrate_monomial([2, 0], BoxDomain(2, 1.0))  # ∫∫ x² dy dx over [-1,1]²