API reference

analyticity_factor(residual::AbstractArray{<:Real}; marg = true, radius = 1)
analyticity_factor(x, q, v0; marg = true, radius = 1, o = nothing)

Return the blockwise residual correction 1 / (1 + sum(norm(residual_S))/radius), where S ranges over nonempty subsets of parties. Marginals occupy the last index of each axis; exclude the fixed all-identity coordinate. Norms are Euclidean norms of vectorised blocks. Without marginals, return 1 / (1 + norm(residual)/radius). radius must be finite and strictly positive, and certify a local ball about the noise o in this same norm, within the normalised correlation affine space. The default unit radius applies to white noise; it is not generally valid for other noise. Locality of the ball is the caller's responsibility.

The second form uses the residual v0*q + (1-v0)*o - x, with white-noise centre by default. The input x must be a normalised local point and, with marginals, the all-identity coordinates of x, q, and o must equal one. The function does not check locality or normalisation. The corrected finite visibility is the returned factor times v0, along the same noise line: the corrected point is o + factor*(v0*q + (1-v0)*o - o). For white noise and q = shrinking_target(p, eta), the global white-noise visibility is additionally multiplied by prod(eta) (or eta^N for a common factor).

Pass full correlation tensors, not probabilities or symmetry-reduced coordinates. Floating-point evaluation is a numerical estimate; rigorous certificates require upper bounds on residual norms and lower bounds on shrinking factors.

Calls the lazy pairwise blended conditional gradient algorithm from Frank-Wolfe package.

The supplied p is always the finite target: measurement-shrinking compensation is explicit through shrinking_target. Returned inequalities and β refer to that finite target; shr2 scales only the displayed lower bounds.

Arguments:

• p: a correlation/probability tensor of order N.

Returns:

• x: a correlation/probability tensor of order N, the output of the Frank-Wolfe algorithm,
• ds: a deterministic strategy, the atom returned by the last LMO,
• primal: ½|x-(v₀*p+(1-v₀)*o)|²,
• dual_gap: ⟨x-(v₀*p+(1-v₀)*o), x-ds⟩,
• active_set: all deterministic strategies used for the decomposition of the last iterate x, contains fields weights, atoms, and x,
• M: a Bell inequality, meaningful only if the dual gap is small enough,
• β: the visibility bound along the line from o to p, reliable only if the last LMO is exact,
• status: the termination status returned by FrankWolfe.

Optional arguments:

• o: same type as p, corresponds to the noise to be added, by default the center of the polytope,
• prob: a boolean, indicates if p is a correlation or probability array,
• marg: a boolean, indicates if p contains marginals (by convention in the last index of each dimension),
• v0: the visibility used to make a nonlocal p closer to the local polytope,
• epsilon: the tolerance, used as a stopping criterion (when the primal value or the dual gap go below its value), by default 10Base.rtoldefault(T),
• shortcut: if positive, the ratio between primal and dual gap for early termination,
• verbose: an integer, indicates the level of verbosity from 0 to 4,
• shr2: a squared measurement shrinking factor (or one per party), used to scale displayed corrected correlation bounds. The usual product of shrinking factors applies to white noise with shrinking_target; for inverse-shrunk target and noise use shr2 = 1,
• radius: finite positive local radius about o in the residual block norm, default 1 (certified for white noise). For other noise the caller must supply a valid radius; see analyticity_factor,
• mode: an integer, 0 is for the heuristic LMO, 1 for the enumeration LMO,
• nb: an integer, number of random tries in the LMO, if heuristic, by default 10^2,
• TL: type of the last call of the LMO,
• mode_last: an integer, mode of the last call of the LMO, -1 for no last call,
• nb_last: an integer, number of random tries in the last LMO, if heuristic, by default 10nb,
• sym: a boolean, indicates if the symmetry of the input should be used, by default automatic choice,
• callback_interval: an integer, print interval if verbose = 3,
• seed: an integer, the initial random seed.

dyadic_certificate(result, target; v0, q = nothing, o = nothing, radius = 1,
                   weight_bits = 52, residual_bits = 40, backend = :blas,
                   row_block = 512, column_block = 512, atom_block = 4096)

Certify a finite dichotomic correlation target after a bell_frank_wolfe run, for any number of parties, with or without marginals. result may be the full solver result, its fifth element (the active set), or an ActiveSetStorage. The number of parties and marginal convention are inferred from the storage. Probability tensors and nondeterministic atoms are not supported.

Supply the intended EXACT target tensor and visibility, e.g. v0=69//100. Float inputs for the target, visibility, noise, or radius are rejected rather than silently treating an approximate physical tensor as exact. The solver's visibility and target are not stored in its result and cannot be inferred. With marginals, the identity setting is last on every axis.

Round the saved weights with dyadic_weights, reconstruct the local mixture with dyadic_orbit_sums, and bound the full residual with dyadic_residual_bound. No cached iterate or approximate residual is trusted. An optional q must describe a locality-preserving uniform orbit average; this is the CALLER'S premise and is not verified. Without q, the raw stored strategies are used, even if the solver ran in a symmetry subspace.

White noise and its unit block-norm local ball are the defaults. For other exact noise o, supply a certified positive rational radius in the same norm. The returned rational nu = radius/(radius + residual.norm_upper) guarantees locality on the same noise line at visibility_lower = nu*v0. No measurement-shrinking factor is applied by this form.

Return a named tuple containing visibility_lower, finite_visibility, nu, v0, radius, marg, the rounded weights, exact orbit sums, residual bounds, and the original stored ax and supplied q for replay. Inputs are not modified; ax and q are referenced, not copied.

res = bell_frank_wolfe(Float64.(p); v0=0.69, marg=true)
certificate = dyadic_certificate(res, p; v0=69//100)  # p is exact

See the keyword-only form for rounding measurements and applying shrinking.

dyadic_certificate(result; measurements, target = nothing, q = nothing,
                   v0, shr2 = nothing, coordinate_bits = 40,
                   shrinking_bits = 60, kwargs...)

Round measurement directions and certify a white-noise visibility for all measurements simulated by their convex hulls. measurements is a common m × d matrix, or a tuple/vector of one matrix per party. Marginal identity settings must NOT appear in these matrices. Party counts and the marginal convention come from the solver's correlation active set.

Each matrix is rounded by dyadic_measurements. target is a function receiving a TUPLE of the resulting exact rational matrices and returning the FULL exact finite correlation tensor for the intended state. This explicit builder is required for multipartite targets and for targets with marginals: the state cannot be recovered from the measurement directions or solver result. Its returned all-identity entry must be 1. When omitted for a bipartite target without marginals, use the fast lazy Gram target vertices[1]*vertices[2]'. The physical interpretation of the tensor and its measurement coordinates is the caller's responsibility; no state is inferred or numerically certified.

shr2 is a certified EXACT squared hull-radius lower bound for the NEW rounded points, either common to all parties or one per party. If omitted, compute each with shrinking_squared_exact. To reuse an old certified geometry, round it separately and use shrinking_squared_transfer. Floating-point shrinking estimates are rejected. Hulls include antipodes (outcome relabelling).

For marginals, choose exact rational lower bounds eta[n] on the square roots of shr2[n], compensate the target with shrinking_target(target, eta), and multiply its certified finite visibility by prod(eta). Without marginals, use a rational lower bound on sqrt(prod(shr2)) directly. Square roots already rational are preserved exactly; otherwise round DOWN to shrinking_bits binary places. This avoids applying a full-correlation shrinking exponent to uncompensated marginal blocks. shrinking_bits must be positive.

Return the finite certificate's fields, with visibility_lower now including shrinking, plus per-party measurements, eta2, eta, and shrinking_product. finite_visibility still refers to the supplied/compensated finite target. Other keywords are passed to the exact-target form. For nonwhite noise or inverse-shrunk targets, use that form and supply the exact noise line explicitly. The q premise is unchanged: the caller guarantees locality preservation.

res = bell_frank_wolfe(v*v'; v0=0.69)
cert = dyadic_certificate(res; measurements=v, q, v0=69//100)

For a general state, use target = vertices -> exact_correlations(vertices); that function must evaluate the state on the rounded vertices, not return the old floating-point tensor. All reported certificate bounds are exact rationals.

dyadic_measurements(vertices; bits = 40)

Approximate the normalised rows of a real m × d matrix by dyadic points in the closed unit ball. Return (numerators, denominator), with Int64 entries and common denominator 2^bits. bits must lie in 1:50.

Float64 normalisation proposes candidates only. Their squared norms are checked in integer arithmetic. An outside candidate is contracted by a rounded up integer square root, then truncated towards zero, guaranteeing sum(abs2, numerators[i, :]) ≤ denominator^2 exactly. Rows need not initially be unit length, but must be nonzero, finite, and representable in Float64. The input is not modified. Shrinking factors must refer to these NEW points; use shrinking_squared_exact or shrinking_squared_transfer.

dyadic_orbit_sums(ax, counts, q = nothing; marg = false, backend = :blas,
                 row_block = 512, column_block = 512, atom_block = 4096)

Reconstruct and optionally orbit-average an exact local correlation mixture. ax[n] is the party's BitMatrix of deterministic signs (atoms in rows, settings in columns; true means +1). counts are nonnegative integer weights with positive sum at most 2^52. With marg=true, append the identity setting of value +1 to every party; it is NOT stored in ax.

q is a full tensor of consecutive positive orbit labels, or nothing for no projection. The caller must ensure that uniform averaging over these labels preserves locality. Only dimensions, labels, and the singleton all-identity orbit are checked; no symmetry discovery or orbit verification is performed.

Return (numerators, multiplicities, denominator, dims, marg). The local entry at a full index i is numerators[q[i]] / (denominator*multiplicities[q[i]]); without q, use the linear index instead. Orbit numerators are Int128.

For any number of parties, group the first half as matrix rows and the rest as columns. Products of their signs remain ±1. Blocked ordinary IEEE binary64 GEMM is exact since every product and partial sum is an integer of magnitude at most sum(counts) ≤ 2^52. This assumes standard CPU BLAS arithmetic, not reduced-precision or Strassen multiplication. backend=:integer uses direct integer summation instead. The block sizes control workspace, not accuracy.

dyadic_residual_bound(target, orbit; q = nothing, marg = orbit.marg,
                      v0 = 1, o = nothing, bits = 40)

Bound the FULL residual v0*target + (1-v0)*o - local using integer floors, where orbit is returned by dyadic_orbit_sums, with the SAME q. target and optional noise o must be full arrays of exact integers/rationals (lazy AbstractArrays are allowed). v0 must be exact and in [0,1]. White noise is the default. With marginals, the last coordinate of each axis is the identity; the all-identity entry of target, noise, and local must be 1.

For F=2^bits, floor each target/noise-line entry and each local orbit value on the 1/F grid. If their integer difference is e, the true absolute error is at most (|e|+1)/F. Square these bounds and sum within each nonempty correlation block, then round each square root UP with integer arithmetic. bits must lie in 1:52; arbitrary-size integers avoid accumulator overflow.

Return exact rational squared_upper (the full Frobenius norm squared bound), norm_upper (the sum of block norm bounds used by analyticity_factor), and block_squared_upper, block_norm_upper, squared_numerators, and scale. Without marginals there is one block. Otherwise blocks are ordered by the nonzero bit masks of measured parties, party 1 being the least significant bit. Every full entry is evaluated: the target need NOT respect the symmetry of q. The bound includes grid-rounding slack even when the residual is exactly zero.

dyadic_weights(weights; bits = 52)

Normalise nonnegative finite weights and round them to integer counts summing to 2^bits, using the largest-remainder rule. bits must lie in 1:52. Floating inputs are interpreted as their exact binary values before rounding. The input is not modified; zero weights and a nonunit initial sum are allowed, but the sum must be positive.

Return (numerators, denominator, l1_change). Dividing the Int64 numerators by the common denominator gives an exact probability distribution. l1_change is the exact rational distance from the normalised input weights. The residual certificate uses the rounded mixture directly, so this change must not be added to its error bound a second time.

local_bound_correlation(FC; marg = false, mode = 0)

Compute the local bound of a Bell inequality in correlator notation parametrised by FC. No symmetry detection is implemented yet, used mostly for pedagogy and tests. The default mode=0 is heuristic; use mode=1 for exact enumeration.

local_bound_probability(FP; mode = 0)

Compute the local bound of a Bell inequality in probability notation parametrised by FP. No symmetry detection is implemented yet, used mostly for pedagogy and tests. The default mode=0 is heuristic; use mode=1 for exact enumeration.

move_marg(FC::AbstractArray, sense::Int = -1)

Change convention for the placement of marginals. By default, converts from first to last index. If sense=1, convert back from last to first index.

nonlocality_threshold(p::Array, lower_bound = 0, upper_bound = 1)

Compute the nonlocality threshold of the probability/correlation tensor p.

Returns:

• lower_bound: a lower bound including the blockwise analyticity correction for white-noise correlation searches; an approximate finite-scenario bound otherwise,
• upper_bound: a (heuristic) upper bound on the nonlocality threshold of p
• local_model: a decomposition of the tensor p with visibility lower_bound (up to a distance 2√epsilon),
• bell_inequality: a (heuristic) Bell inequality corresponding to upper_bound.

The local model or Bell inequality is nothing if the search terminates before finding the corresponding certificate. With analyticity correction, the returned local model is mixed with white noise to match the corrected visibility. This remains an approximate decomposition; the residual correction certifies locality without expanding its extra atoms. Floating-point evaluation does not provide interval-certified bounds.

The search and returned bounds concern the supplied finite target. To compensate for measurement shrinking, pass shrinking_target(p, eta) and multiply the returned lower bound by prod(eta) (or eta^ndims(p) for a common factor). shr2 only affects displayed bounds.

Optional arguments:

• upper: whether to start from the upper bound or the lower bound, false by default
• digits: number of digits of lower_bound, 3 by default,
• analyticity: apply the residual correction, by default for correlation tensors with white noise. Set to false to return the uncorrected numerical bracket. digits controls that bracket; the corrected lower bound can be smaller.
• radius: local radius used by analyticity_factor and displayed bounds. With nonwhite noise, only the displayed bounds are corrected; returned bounds and models retain the uncorrected finite-scenario convention.
• for the other optional arguments, see bell_frank_wolfe.

pythagorean_approximation(vec; epsilon = 1.0e-16)

Approximate the unit rows of an m × d real matrix by Rational{BigInt} rows of exactly unit squared norm, for any positive dimension d. The input is not modified. Entries smaller than the nonnegative finite epsilon are first set to zero. Rows must remain normalised up to floating-point precision after this cleanup. Exact rational unit rows are preserved.

Uses the recursive rational half-angle parametrisation from _normalize!. For BigFloat input, run at the desired precision (e.g. inside setprecision).

shrinking_squared(vec; verbose = true)
shrinking_squared(vecs; verbose = true)

Estimate the squared shrinking factor of an m × d Bloch matrix, including antipodal rows. For a vector of matrices, return the minimum squared factor. This uses numerical norms even for rational vertices; use shrinking_squared_exact for an exact rational result, or shrinking_squared_transfer to transfer a certified lower bound.

shrinking_squared_exact(vec; antipodal = true, verbose = true)
shrinking_squared_exact(numerators, denominator; kwargs...)
shrinking_squared_exact(vecs; kwargs...)

Compute the exact squared inradius about the origin of the convex hull of the rows of an m × d matrix. Include antipodal vertices by default, as in shrinking_squared. Entries must be integers or rationals, and all rows must lie in the unit ball. The hull must be full dimensional with the origin strictly inside. The result is a Rational{BigInt}.

Enumerate facets using Polyhedra with rational coordinates and minimise beta^2 / sum(abs2, normal); no square root or floating-point geometry enters the result. Floating-point inputs are rejected: convert or approximate them explicitly first. The two-argument form takes integer numerators and a positive common integer denominator. For a vector of matrices, return the minimum of their squared shrinking factors. verbose prints a numerical inradius only.

shrinking_squared_transfer(old_vertices, old_eta2, new_vertices; bits = 60)
shrinking_squared_transfer(old_vertices, old_eta2, numerators, denominator; bits = 60)

Transfer a previously certified squared shrinking lower bound old_eta2 to paired new vertices, without enumerating their hull. Both matrices must have identical sizes, exact integer/rational entries, and rows inside the unit ball. old_eta2 must be an exact integer or rational in (0, 1]; its validity for the old hull is the caller's responsibility. Use the same antipodal convention for both hulls, and pair corresponding rows. A common positive integer denominator may be supplied for the new integer numerators.

If delta is the largest distance between paired vertices, support functions change by at most delta, hence eta_new >= sqrt(old_eta2) - delta. Round the first square root down and the second up to multiples of 2^-bits, using integer arithmetic only. Return their positive squared difference as a Rational{BigInt}. Throw if this bound is nonpositive (increasing bits may help when rounding is responsible). Identical vertices preserve old_eta2. This is a certified lower bound, not in general the exact new inradius.

shrinking_target!(q, p, eta; marg = true)

Write the compensated target into q; see shrinking_target. The destination must have the same axes as p. Exact aliasing (q === p) is supported; other overlapping views are not.

shrinking_target(p::AbstractArray{<:Real}, eta; marg = true)
shrinking_target!(q, p, eta; marg = true)

Compensate a full correlation tensor for local measurement shrinking. eta is a common shrinking factor, or one factor per party, in (0, 1]. With marginals in the last index of each axis, multiply each nonempty correlator block p_S by prod(eta[n] for n outside S). Preserve the all-identity coordinate. Without marginals, the target is unchanged.

If the compensated target is local at visibility v, measurement shrinking gives white-noise visibility prod(eta) * v for the original tensor, or eta^N * v for a common factor and N parties. Apply this function before symmetry reduction. The finite target need not itself be physical.

The in-place form supports q === p; other overlapping views are not supported. The destination element type must represent the scaled values.