Extension interface

split_block(ratio, c11, c12, c22, s1::SplitScale, s2::SplitScale) -> (rho1, rho2)

The member of the block's cover family the split selects, from the two channels each slot carries.

Parent values fix the split covariantly except at blocks where a vanishing value meets a degenerate feasible segment, and there no covariant local choice exists at all: the equalizing point runs to an endpoint that is itself zero or infinite. So the values are read wherever they carry a usable scale and the declared input scales stand in wherever they do not, which keeps covariance where it is attainable and asks for units only where it is not. value_bounds_split is the test and ratio its threshold; see split_ratio.

The guard compares two differently scaled quantities and is therefore not itself invariant: a rescaling moves a block's value ratio, and one sitting near ratio can cross it. Nothing narrower is available — the one degeneracy test a rescaling leaves alone is an exactly vanishing pair, and a value exactly zero at one point is of order the displacement at a nearby one, so that test catches a block at one evaluation point and not at its neighbors.

Cost: the guard is three comparisons on values already recorded, and the declared scales are built by the first block that asks for them rather than in advance. A program whose blocks all keep a usable value pair therefore pays nothing for the fallback; one that reaches a guarded block pays one extra O(L) forward sweep, once. guards_any_block answers which case a tape is in without running a cover.

split_ratio(ratio, T::Type) -> r

The degeneracy threshold a cover over element type T runs at, from the ratio keyword its entry point was given, in M = magtype(T). nothing means sqrt(eps(M)).

ratio is how lopsided a value pair must be before value_bounds_split counts it as carrying no scale, and it must lie in [0, 1): at 1 every pair is condemned. It caps the split asymmetry the value channel can express, since a multiplication block aims at rho1 = sqrt(c12*z2/z1) and so at a ratio of sqrt(ratio) between the two references.

The default resolves against M rather than being stored as a constant because the quantity it bounds is representability. eps(M) is the floor below which one parent value is indistinguishable from zero at the scale of the other; sqrt(eps(M)) is the half-the-digits margin above that floor, capping the asymmetry at eps(M)^(-1/4). A constant chosen at one precision is silently wrong at another — 1e-8 sits below eps(Float32), which would collapse the guard to an exactly-vanishing test on a Float32 program, and that test does not survive an ordinary perturbation of the evaluation point.

guards_any_block(t::Tape; ratio=nothing) -> Bool

Whether any block of t could reach the declared-scale fallback, decided by the same value_bounds_split test the split itself applies and under the same conditions in which cover_slots! consults a scale at all.

This is how a caller learns whether the unit-scale default is load-bearing for their program. On a tape that answers false the cover reads parent values everywhere, is exactly covariant under an input rescaling, and pays no extra sweep. On one that answers true the cover imports the units its inputs are expressed in at the blocks where no covariant local choice exists, and stays covariant only when the declared scales are rescaled alongside.

The test reads recorded values and second partials, never an adjoint and never an intermediate gradient. The adjoint scales a block's three entries together, so it changes neither the value ratio nor which diagonal entries are nonzero, except where it underflows one of them to zero. Reading it that way overcounts — a block whose adjoint vanishes contributes no cover term for the fallback to reach — and, in that underflow corner alone, can undercount. Ignoring intermediate gradients overcounts too: a block with a dead_slot is covered diagonally and never consults a scale at all.

So this predicts which programs will ask for the declared scales without running a cover; the covers themselves build those scales on demand and never call it, and they are the authority on what a given tape does.

sensitivity_scales(w, t::Tape) -> sigma

One abs_forward sweep seeded with the declared input scales w, so sigma_j >= sum_a w[a] * abs(∂v_j/∂x_a) with equality when no path cancellation occurs. w = nothing means unit scales, which asserts that the inputs share units.

These are the references the split falls back on where a block's parent values carry no scale. They transform exactly as the values do under a rescaling of the intermediate variables, so the covariance argument of balanced_split applies to them unchanged; under an input rescaling x -> D*x the cover stays exactly covariant when the scales are rescaled alongside, w -> D*w.

Cost: one O(L) forward sweep, paid once per cover and only by a program that reaches the fallback.

split_scales(w, t::Tape) -> scales

The per-node scale references local_cover_slots reads, one SplitScale per node: a lazy view of magnitude.(t.vals) for the value channel, paired with the declared scales w propagates, which are built only if some block asks for them. Obtained once per cover sweep and shared by every node.

SplitScale(value, declared)

Both scale references split_block may read at one slot: the parent value, and the means of reaching the declared scale that stands in for it. They travel together because the choice between them belongs to a block, not to a node — a variable whose value carries scale in one product can be the vanishing operand of another — so a single scale per node cannot express the split.

declared is whatever declared_scale resolves: the number itself where a backend already has it, and a DeferredScale where the sweep that produces it should run only if some block asks. dead_slot reads the same field, since a backend that can name a slot's declared scale can say whether the slot's intermediate gradient vanishes.

dead_slot(s::SplitScale) -> Bool

Whether the slot's intermediate gradient ∇v is identically zero.

Such a slot contributes no outer product to the Hessian H = sum_k vbar_k * sum_{i,j} C^k_{ij} * ∇v_{p_i} * ∇v_{p_j}', so every entry of its block that involves it is multiplied away and needs no cover. cover_slots! drops those entries, which is what keeps the cover from importing a live slot's scale into a term that vanishes.

The test is exact: a tolerance would either declare a live slot dead, leaving a real Hessian entry uncovered, or compare two quantities with different units.

A slot reported by a bare number carries no liveness of its own, and the number is read as the whole answer — which a declared scale is only while it cannot underflow. The backends here supply a reference that reports liveness separately: DeferredScale over a tape, ForwardScale in a propagation.

forward_scale(x::FVar) -> SplitScale

The scale references split_block reads at x — the forward analogue of split_scales, which reads a whole tape at once.

Both channels are already on the variable: the value, and the ForwardScale behind it. So a propagation splits its local covers by exactly the rule a tape does, and the two backends are comparable coordinate by coordinate.

Nothing is deferred here, unlike on a tape. sigma obeys the same recursion as the tape's sensitivity_scales sweep, sigma_k = sum_i magnitude(d_i) * sigma_{p_i}, but it is a single scalar per variable rather than a pass over recorded nodes — so carrying it costs O(1) per elemental whether or not any block asks for it, against the O(L) sweep a tape would have to run.

Tape{T}

A Wengert tape recorded while a function runs on TVar inputs, holding its values and partials in the element type T. eltype(t) reports T: it is the working precision of the whole construction, and a value or weight a caller supplies is converted into it. Signed results such as gradients are returned in T; covers and other nonnegative bounds are returned in magtype(T), which is T itself for a real floating-point type, and computed in boundtype(T).

T may be a type whose values stand for sets of numbers, such as an interval type from IntervalArithmetic.jl. The values and partials recorded over a box then enclose those at every point of the box, magnitude takes the largest absolute value over each set, and the covers computed from the tape hold everywhere on the box: the sweeps round their bounds upward (see boundtype), so this holds exactly rather than up to rounding. Program control flow that the set does not decide, such as a comparison between overlapping intervals, raises an error while recording.

The first ninputs nodes represent the function inputs. Later nodes are appended in topological order, so every parent precedes its child. Node k stores its value vals[k], parent indices parents[k], first partials ds[k][i] = ∂φ_k/∂v_{p_i}, and the second partials Cs[k], a Curvature holding the symmetric ∂²φ_k/∂v_{p_i}∂v_{p_j} in whatever form the elemental supplies: UnaryCurvature or BinaryCurvature for a unary or binary elemental, PairedCurvature for an inner-product node (see dotnode), DeclaredCurvature for a declared sub-map output (see submapnode), and NoCurvature for a node with no second-order term.

Input nodes and constant output nodes have no parents and no second partials.

TVar(t::Tape{T}, idx::Integer)
TVar(c::Real)
TVar{T}(c::Real)

A scalar used while recording a function, carrying a value of type T. The CoverVar backend that records onto a Tape.

TVar(t, idx) refers to node idx on t and takes the tape's element type. TVar(c) creates a constant with a value but no tape; this supports typed literals such as zero(eltype(x)) in generic numerical code, which reach TVar{T}(c) and so stay in the program's precision. Operations involving only constants are folded immediately. Operations that mix a constant and a recorded variable store the constant as a coefficient.

Examples

julia> t = Tape([0.5, 2.0]);

julia> x = TVar(t, 1);

julia> DerivativeCovers.isconstant(x), DerivativeCovers.value(x)
(false, 0.5)

julia> c = TVar(3.0);   # a constant: a value and no tape

julia> DerivativeCovers.isconstant(c * TVar(4.0)), DerivativeCovers.value(c * TVar(4.0))
(true, 12.0)

tape_function(f, x0::AbstractVector) -> (tape, outidx)      # scalar output
tape_function(f, x0::AbstractVector) -> (tape, outidxs)     # vector output

Record f at x0 and return its Tape and output node index or indices.

f must accept a vector of TVars and use operations supported by TVar. The second return value is:

• a single TVar output: the Int index outidx of the output node;
• an AbstractVector of TVars: the Vector{Int} of output-node indices outidxs, one per entry. Input nodes and repeated outputs are preserved. A constant output receives its own node and therefore contributes a zero Jacobian row.

A function whose outputs are all constant raises an error, as does a recording that check_values rejects.

Examples

julia> f(x) = exp(x[1]) + x[1] * x[2]^2;

julia> t, outidx = tape_function(f, [0.5, 2.0]);

julia> t.ninputs, outidx, t.vals[outidx]
(2, 6, 3.648721270700128)

julia> r(x) = [x[1] * x[2], exp(x[1]) - x[3], x[2] + 2x[3]];

julia> tr, outidxs = tape_function(r, [0.5, 2.0, -1.0]);

julia> outidxs
3-element Vector{Int64}:
 4
 6
 8

FVar{T, S, M, SM}

A scalar of type T that propagates a Hessian cover forward. This is the tape-free CoverVar backend.

S and SM are the storage of the intermediate-gradient and cover channels, and M = boundtype(T) the type the cover and the declared scale are propagated in (see ForwardContext). M and SM follow from T and S, so FVar{T, S} names the variable type as well. A constant carries the storage too, because generic numerical code builds typed literals from eltype, which is all it has to reach the backend by.

For each variable, the propagation holds the value, the intermediate gradient ∇v, a squared cover of ∇²v (see ForwardChannels), the scalar declared scale sigma the split falls back on, and whether ∇v is identically zero (see ForwardScale). Working memory is O(n) per live variable rather than per recorded node, so a program whose intermediates go out of scope costs O(n × live width) and never materializes a tape.

Build one with forward_cover rather than directly.

CoverVar{T}

A scalar that carries a value of type T through a program so that a cover construction can observe the program's elementals as it runs. TVar records them onto a Tape.

The elementals of overloads.jl are written against this type rather than against any one backend. Each computes its value, first partials and second partials and hands them to one of four recorders, which is where the backends differ. A concrete subtype V defines:

• value and isconstant;
• constantvar(::Type{V}, c). Generic numerical code obtains typed literals through zero, one, oneunit and convert of V, all of which route here, so a constant must carry every type parameter V has;
• whichever of unarynode, binarynode, dotnode and linearnode it supports. An unsupported recorder makes the elementals that reach it unavailable; a subtype supporting none of them cannot run a program.

A constant has no node behind it in any backend: an elemental folds a constant operand into the node that consumes it.

value(x)

The numeric value carried by x: the value of a CoverVar, constant or not, or x itself for a plain Real. The Real fallback lets generic code read the value of either a program variable or a plain number without a branch, detaching a quantity from the backend.

isconstant(x) -> Bool

Whether x is a constant: a value with no node behind it, as returned by constantvar and by zero/one/oneunit of the element type. A constant contributes coefficients to the nodes it enters and has no node of its own.

constantvar(::Type{V}, c::Real) -> V

The constant of backend V carrying value c and nothing else. See CoverVar.

nnodes(t::Tape)

Number of nodes recorded on t, inputs included. This is the tape length L.

natoms(t::Tape)

Number of cover atoms recorded on t: the independent rank-1 second-order terms its nodes expose, summed over the tape (see atom_slots). A unary or binary elemental with nonzero second partials exposes one atom, an inner-product node one per parent pair, and a linear node none.

This is the group count at which cover_hybrid reaches the quadrature merge. It counts what the tape records, so a node whose adjoint or curvature vanishes at the recorded point is counted even though it seeds no channel.

noderef(x::TVar) -> (tape, idx)

The tape and node index behind x. A constant has neither: every elemental folds a constant operand into one that records, so reaching here with a constant means a caller bypassed that folding.

pushnode!(t::Tape, val, parents, d, C) -> TVar

Append a node with value val, parent indices parents, first partials d, and second partials C, and return a TVar referring to it.

val and d are converted to the tape's element type. C is stored as it is given, and every Curvature form the sweeps resolve statically holds that same element type, so an elemental should build its curvature in eltype(t) too.

unarynode(x::TVar, val, d, c) -> TVar

Record a node φ(x) with value val, φ'(x) = d, and φ''(x) = c.

A zero c (one whose magnitude is zero) is stored as NoCurvature rather than as a zero UnaryCurvature: the node contributes no second-order term either way, and that form says so without an entry to load and test.

binarynode(x::TVar, y::TVar, val, d1, d2, c11, c12, c22) -> TVar

Record a node φ(x, y) with value val, first partials (d1, d2), and second partials [c11 c12; c12 c22].

An all-zero block is stored as NoCurvature, for the reason unarynode stores a zero c that way. Addition and subtraction of two tape variables are the common case.

dotnode(ws::AbstractVector{<:TVar}, xs::AbstractVector{<:TVar}) -> TVar

Record an inner-product node y = sum_j w_j*x_j. Parent slots are ordered (w_1, ..., w_q, x_1, ..., x_q), first partials are (x_1, ..., x_q, w_1, ..., w_q), and the second-partial matrix is PairedCurvature(q). Parent slots are independent even when ws and xs share entries (e.g. dotnode(xs, xs)): the decomposition is over slots, not values.

linearnode(coeffs::AbstractVector{<:Real}, xs::AbstractVector{<:TVar}) -> TVar

Record a node y = sum_j coeffs[j]*xs[j] for constant coeffs. First partials are the coefficients converted to the tape's element type; the second-partial matrix is empty, since a linear combination of constants and TVars has no second-order term.

constantnode(t::Tape, c::Real) -> TVar

Record a node with value c, no parents, and empty partials. Its absolute sensitivity is zero, so its Jacobian row vanishes and it contributes no atom to any cover.

Elementals fold a constant operand into the node that consumes it rather than giving it a node (see TVar), so this is reached only at the boundary where a program's outputs become node indices and a constant needs an index to be one: a fixed offset among a program's outputs, such as the unit residual that carries an additive constant of a least-squares loss.

atom_slots(C::Curvature, nslots) -> atoms

The groups of a node's parent slots that the quadrature merge of cover_hessian combines separately. One atom spanning every slot is correct whenever the node's cover is a single rank-1 term. A PairedCurvature(q) instead forms q atoms, one per pair (j, q+j), matching the block-diagonal structure of its second partials: the atoms it contributes are independent rank-1 terms, not one combined term for the whole node. A NoCurvature forms none.

The atoms must partition the slots carrying cover weight: a slot left out of every atom drops its contribution from the quadrature merge and from cover_hybrid, and a slot named twice counts it twice.

check_finite_partials(site::CoverSite, av, c::Tuple)

Throw unless the adjoint weight av of the elemental at site and its covering entries c are all finite, reporting the values around it.

The covering invariant is a set of comparisons against c, so it cannot fire meaningfully when c is not finite: every comparison against NaN is false and no finite rho covers an infinite entry, so the test would reject a split that is not at fault. A non-finite c says instead that the program is not twice-differentiable where it was evaluated — a division by zero, or a power, logarithm or root at its singularity. A non-finite av says the same of a node downstream, whose partials the adjoint has already passed through.

check_values(t::Tape)
check_values(site::CoverSite, val, ds, cs)

Throw if the values recorded on t show that the program left the domain of one of its elementals. tape_function calls this once recording finishes. Over numbers there is nothing to check, since a numeric elemental outside its domain raises an error or returns a non-finite value the cover checks report. A set-valued element type, whose elementals can instead return a result flagged as undefined on part of the set, extends this to look for that flag.

The second form checks one elemental: its value val, its first partials ds and its second partials cs, with site naming it in the error. A propagation that keeps no tape, such as forward_cover, calls it at each elemental.

Curvature

The second partials ∂²φ_k/∂v_{p_i}∂v_{p_j} of one tape node, in whatever form that node's elemental can supply them. Tape stores one per node.

Cover constructions never need the entries themselves, only bounds on their magnitudes, so a curvature may declare a bound instead of recording a matrix (DeclaredCurvature). The sweeps reach every form through the methods below, which is what lets a new elemental introduce a new form without touching them.

A new subtype must define:

• cover_slots: the node's rank-1 cover, slot-aligned with its parents. It reads its scale references and its error context through a CoverSite rather than off a tape, so one definition serves every backend.
• bound_contribution!: the node's term of the elementwise bound B.
• atom_slots: the parent-slot groups the quadrature merge combines separately. The fallback treats every slot as one atom, which is correct for any curvature whose cover is a single rank-1 term.

A subtype whose entries are recorded exactly may also define:

• curvature_apply!: the second-order term of a Hessian-vector product.
• has_exact_partials: false suppresses hvp on any tape containing the node. The fallback is true, so a subtype that declares a bound must override it.

getindex and size are introspection, not part of the contract: a curvature that records no entries has none to report.

NoCurvature() <: Curvature

The curvature of a node with no second-order term: an input node, a linear n-ary node such as linearnode, or an elemental whose curvature at the recorded point is zero. It contributes to no cover and to no bound, so the sweeps skip it by dispatch rather than by inspecting a matrix of zeros.

UnaryCurvature(c) <: Curvature

The second partial φ''(x) of a unary elemental, held inline as the single number it is.

A tape stores its curvatures in a container with abstract element type, so every form is reached through a reference; holding the entry inline rather than in a 1×1 matrix saves the second reference on the way to it, once per node per sweep, and an allocation per node while recording.

BinaryCurvature(c11, c12, c22) <: Curvature

The second partials [c11 c12; c12 c22] of a binary elemental, slot-aligned with the node's two parents and held inline. Symmetry is structural: one entry stands for both off-diagonal slots.

Inline for the reason UnaryCurvature is.

DenseCurvature(C::AbstractMatrix) <: Curvature

Recorded second partials as an explicit symmetric matrix, slot-aligned with the node's parents.

This is the general form, for a curvature with no structure to exploit and no fixed arity. The elementals of overloads.jl do not use it: a unary or binary block is UnaryCurvature or BinaryCurvature, which hold their entries inline.

The rank-1 cover of a block larger than 2×2 is not a single split, so cover_slots rejects one; the elementwise bound and the Hessian-vector product read any arity. An elemental with more parents whose curvature does have structure supplies a form that carries it, as PairedCurvature and DeclaredCurvature do.

PairedCurvature(q::Int) <: Curvature

Curvature of an inner-product node with q pairs of parent slots ordered (w_1, ..., w_q, x_1, ..., x_q) (see dotnode). Entry (i, j) is 1 when i and j are paired (abs(i - j) == q), the sole nonzero second partials ∂²y/∂w_j∂x_j of y = sum_j w_j*x_j, and 0 otherwise.

The 2q × 2q matrix is never materialized: it has 2q nonzeros, and its block-diagonal structure lets every operation run over the q pairs directly.

DeclaredCurvature(s::AbstractVector) <: Curvature

Declared bound on the second partials of a sub-map output (see submapnode): abs(C[i, j]) <= s[i]*s[j], with the true C not recorded. s is nonnegative and slot-aligned with the node's parents. On a Tape{T}, s is a vector of magtype(T).

Cover constructions consume second partials only through bounds on abs.(C), so the declaration substitutes for a recorded matrix everywhere except hvp, which needs the signed entries and fails fast — hence has_exact_partials is false. Covers built over such a node are valid conditional on the declaration, which is the caller's assertion: nothing on the tape can verify it.

cover_slots(C::Curvature, ratio, site, av, sl; check=true) -> Union{Nothing,Vector{M}}
cover_slots!(rho::Vector, C::Curvature, ratio, site, av, sl; check=true) -> Bool

The rank-1 cover of one elemental, slot-aligned with its operands: nonnegative rho with rho[i]*rho[j] >= av*magnitude(C[i,j]), where av is the weight the elemental's second-order term carries (the absolute adjoint magnitude(vbar[k]) in a reverse sweep, 1 in a forward one). Returns nothing when the elemental contributes no second-order term.

Entries involving a slot whose intermediate gradient vanishes are exempt: they are multiplied by zero in the Hessian, so covering them would only import the other slot's scale into a term that is not there. dead_slot is the test, and it reads the slot's entry of sl.

rho, av and the entries of sl are magnitudes, in M = magtype(eltype(site)). rho must be a Vector{M}: local_cover_slots asserts that type so that the sweeps, which call it once per node, need not box the result. The entries av*magnitude(C[i,j]) and the products that cover them are formed in boundtype(eltype(site)), so where that rounds upward the invariant holds exactly rather than up to rounding.

The mutating form resizes rho to the elemental's slot count and fills it, returning whether there is a second-order term at all and leaving rho untouched when there is not. Every construction lives there, so that a sweep can hand it one buffer and let no node allocate for its slot weights; cover_slots allocates a fresh rho and forwards.

split_block breaks the tie at off-diagonal-active blocks; ratio is its degeneracy threshold and sl its scale references, one per slot and in slot order — not indexed by node, so that a backend holding no nodes can supply them. check verifies the covering invariant against the recorded entries; a curvature that declares a bound rather than recording entries has nothing to check against.

site is a CoverSite, used for the element type and for error reporting. This is the per-elemental kernel shared by local_covers, the CoverOnly sweep, and the forward propagation; see local_cover_slots.

cover_slots(C::Curvature, ratio, site, av, sl; check=true) -> Union{Nothing,Vector{M}}
cover_slots!(rho::Vector, C::Curvature, ratio, site, av, sl; check=true) -> Bool

The rank-1 cover of one elemental, slot-aligned with its operands: nonnegative rho with rho[i]*rho[j] >= av*magnitude(C[i,j]), where av is the weight the elemental's second-order term carries (the absolute adjoint magnitude(vbar[k]) in a reverse sweep, 1 in a forward one). Returns nothing when the elemental contributes no second-order term.

Entries involving a slot whose intermediate gradient vanishes are exempt: they are multiplied by zero in the Hessian, so covering them would only import the other slot's scale into a term that is not there. dead_slot is the test, and it reads the slot's entry of sl.

rho, av and the entries of sl are magnitudes, in M = magtype(eltype(site)). rho must be a Vector{M}: local_cover_slots asserts that type so that the sweeps, which call it once per node, need not box the result. The entries av*magnitude(C[i,j]) and the products that cover them are formed in boundtype(eltype(site)), so where that rounds upward the invariant holds exactly rather than up to rounding.

The mutating form resizes rho to the elemental's slot count and fills it, returning whether there is a second-order term at all and leaving rho untouched when there is not. Every construction lives there, so that a sweep can hand it one buffer and let no node allocate for its slot weights; cover_slots allocates a fresh rho and forwards.

split_block breaks the tie at off-diagonal-active blocks; ratio is its degeneracy threshold and sl its scale references, one per slot and in slot order — not indexed by node, so that a backend holding no nodes can supply them. check verifies the covering invariant against the recorded entries; a curvature that declares a bound rather than recording entries has nothing to check against.

site is a CoverSite, used for the element type and for error reporting. This is the per-elemental kernel shared by local_covers, the CoverOnly sweep, and the forward propagation; see local_cover_slots.

bound_contribution!(B, C::Curvature, av, ps, utilde)

Accumulate node k's term of the elementwise bound into B: av * sum_{i,j} magnitude(C[i,j]) * utilde[ps[i]] * utilde[ps[j]]', where av is the absolute adjoint and utilde the forward absolute sensitivities. See hessian_bound.

A DeclaredCurvature contributes its declared bound s[i]*s[j] >= abs(C[i,j]), which factorizes, so its term is the single rank-1 outer product w*w' with w = sum_i s[i]*utilde[ps[i]].

curvature_apply!(vbardot, C::Curvature, ps, bk, dv)

Accumulate the second-order term of a Hessian-vector product into vbardot: vbardot[ps[i]] += bk * sum_j C[i,j] * dv[ps[j]], where bk is the node's adjoint and dv the forward tangents. See hvp.

Defined only for a curvature that records its entries exactly; has_exact_partials reports which do.

has_exact_partials(C::Curvature) -> Bool

Whether C records its second partials exactly, so an exact second-order product can be formed through it. hvp refuses to run on a tape carrying any node that answers false.

CoverSite

Where a local cover is being built, carried by cover_slots for error reporting.

The split needs the element type to work in and enough context to say which elemental it failed on, and nothing else about the backend. Naming that little through a site is what lets one split serve a tape sweep and a tapeless forward pass alike.

A subtype defines Base.eltype, site_label and site_values.

TapeSite(t::Tape, k::Integer)

The CoverSite naming node k of t, so that a cover built there reports which node it failed on and what the values around it were.

site_label(site::CoverSite) -> AbstractString

How to name the elemental in an error message, such as "node 7".

site_values(site::CoverSite) -> AbstractString

The values surrounding the elemental, for the error that reports where a program stopped being twice-differentiable: "node value 1.5, parent values (0.3, 2.0)" and the like.

submap(F, xs::AbstractVector{<:CoverVar}; scales=nothing, ratio=nothing) -> Vector{<:CoverVar}

Record F as a single declared call on the tape of xs: evaluate declare_submap at the parent values and embed the result via submapnode. The recording tape gains one node per output of F — independent of how many operations F performs internally — at the price that exact Hessian-vector products (hvp) become unavailable on that tape. A scalar-valued F yields a one-element vector. scales and ratio are passed to declare_submap and govern only the declaration's curvature covers.

The declaration is recorded on the parent values themselves, so on a tape whose values stand for sets, such as intervals over a box, it holds over the same sets as the rest of the tape.

The declaration it embeds is dense in the sub-map's outputs and arguments, so this keeps the outer tape on one tape at diagnostic tier. compose_jacobian_cover and compose_hessian_cover cover the same composite from covers of the factors alone.

declare_submap(F, x0::AbstractVector; scales=nothing, ratio=nothing) -> (; y, J, S)

The declaration submapnode consumes, distilled from F's own scratch tape at x0: exact values y = F(x0), the exact signed Jacobian J from per-output reverse sweeps, and curvature covers S whose row k is the sum-merge cover of ∇²F_k, so abs.(∇²F_k) <= S[k, :]*S[k, :]' by cover_hessian's guarantee. F may return a single TVar (declared as one output) or a vector of them; a linear output gets an exactly zero row of S.

Costs O(q*L_F) time — one q-channel adjoint sweep over F's tape plus one absolute sweep per output — once per declaration, and O(q*n) memory for J and S, which places it at diagnostic tier alongside submapnode; the scratch tape is discarded afterwards. The declaration transforms covariantly under rescalings of F's arguments and outputs because each ingredient inherits the covariance of the underlying sweeps, including the zero-parent-value fallback exception of the local tie-break.

scales and ratio govern the split behind the rows of S, as in cover_hessian. The scratch tape's inputs are F's arguments, so scales is expressed in those units, not in the recording tape's. The two splits are otherwise independent: the recording tape consumes S as a fixed bound, so its covers stay valid whatever choice stands behind the declaration.

y and J are in the element type of x0 (as Tape assigns it) and S in its magtype. The declaration holds wherever x0 does: at a point for numeric x0, and at every point of the box for an interval x0, whose y and J then enclose the values and Jacobian over the box and whose S covers the curvature there. A tape whose values stand for sets accepts only a declaration of the second kind (see submapnode).

submapnode(xs::AbstractVector{<:TVar}, y::AbstractVector{<:Real}, J::AbstractMatrix{<:Real},
           S::Union{Nothing,AbstractMatrix{<:Real}}=nothing) -> Vector{TVar}

Record a declared sub-map call y = F(v) as one node per output, without retracing F's internals: node k carries the value y[k], the arguments xs as parents, the exact signed Jacobian row J[k, :] (J[k, i] = ∂F_k/∂v_i at the parent values) as its first partials, and the declared curvature cover S[k, :] as a DeclaredCurvature: abs(∂²F_k/∂v_i∂v_j) <= S[k, i]*S[k, j]. S = nothing, or an all-zero row of S, declares that output linear in the arguments (NoCurvature).

Gradients and every first-order sweep through the call are exact. Hessian covers are valid conditional on the declaration, which is the caller's assertion — nothing on the tape can verify it — and exact Hessian-vector products are unavailable, so hvp throws on any tape containing a declared call. For the covers to stay scale-covariant the declaration must transform as y -> C*y, J -> C*J*inv(D), S -> sqrt.(C)*S*inv(D) under a rescaling v -> D*v of the arguments and a rescaling y -> C*y of the sub-map outputs (a symmetric cover scales as the square root of its output's scale); declare_submap builds declarations with both properties. The covers are then invariant to the output gauge C, which cancels between each output's adjoint and its declared curvature.

Both matrices are dense and are held per node, so a declared call costs O(q*p) memory and declare_submap one adjoint sweep per output to build: diagnostic tier. Where neither is affordable, compose_jacobian_cover and compose_hessian_cover combine covers of the two factors with no joint tape and no dense J. That route bounds the composition's first-order contribution by a rank-1 pair rather than using J itself, which loosens the J'∇²g J term of the chain rule and leaves the declared curvature term unchanged.

On a tape whose element type stands for sets (its magtype differs from it), y and J must be of such a type too, enclosing the values and Jacobian over the sets the parents range over, and S must cover the curvature over them; a declaration with a pointwise y or J is an ArgumentError. submap meets this by declaring on the parent values themselves.

compose_jacobian_cover(a2, b2, a1, b1) -> (; a, b)

Cover of a product of covered factors: if abs.(J2) <= a2*b2' (m × q) and abs.(J1) <= a1*b1' (q × n), then abs.(J2*J1) <= (b2⋅a1) * a2*b1', so the returned pair a = sqrt(b2⋅a1)*a2, b = sqrt(b2⋅a1)*b1 covers J2*J1 — the Jacobian of a composition F2 ∘ F1 wherever each factor's cover holds at the matching point. For the composite cover to hold over a set of points x, such as a box, (a1, b1) must hold at every x in the set and (a2, b2) at every F1(x). The inner product b2⋅a1 is invariant under rescalings of the input, middle, and output spaces (b2 and a1 carry inverse middle-space units), and its square root is split onto both sides so the rule commutes with a rescaling of either factor.

Equality holds entrywise for sign-definite rank-1 factors covered exactly; in general the rule pays the same signed-cancellation looseness as any magnitude-level bound on a product.

The arithmetic is that of the inputs' element type. Pass them as Upward vectors for a composite that holds exactly rather than up to rounding, as covers over a box should.

compose_hessian_cover(hg, gamma, a, b, S=nothing) -> (; h_sum, h_quad)

Symmetric covers of the Hessian of a composite f = g ∘ F at x, with y = F(x), from covers of the pieces:

• hg: symmetric cover of ∇²g(y) (abs.(∇²g) <= hg*hg'), one entry per output of F;
• gamma: entrywise bound gamma[k] >= abs(∂g/∂y_k) at y;
• (a, b): Jacobian cover of F at x (abs.(J_F) <= a*b');
• S: curvature covers with abs.(∇²F_k) <= S[k, :]*S[k, :]' (as declare_submap returns); nothing declares F linear.

The result covers ∇²f at every x where all four inputs hold. At a single point, exact gradient magnitudes qualify as gamma. For a cover over a set of points, such as a box, each input must hold over the whole set: (a, b) and S at every x in it, and hg and gamma at every F(x); the magnitudes of an interval gradient enclosing ∂g/∂y over the set qualify as gamma. The arithmetic is that of the inputs' element type; pass them as Upward vectors for a composite that holds exactly rather than up to rounding, as covers over a box should.

The chain rule ∇²f = J_F'*∇²g*J_F + Σ_k (∂g/∂y_k)*∇²F_k bounds entrywise as abs.(∇²f) <= (a⋅hg)² * b*b' + Σ_k gamma[k] * S[k,:]*S[k,:]' — a sum of rank-1 magnitude atoms u_0 = (a⋅hg)*b and u_k = sqrt(gamma[k])*S[k,:]. Both merges are returned: the sum merge h_sum = Σ_α u_α (cross terms are nonnegative slack) and the quadrature merge h_quad[i] = sqrt(Σ_α u_α[i]²) (Cauchy–Schwarz), with h_quad <= h_sum coordinatewise.

For g = 0.5*norm(r)^2 this generic rule gives norm(a, 1)*b for the Gauss–Newton term: the sum-merge Gram cover with each row of the surrogate Ũ replaced by its rank-1 bound a[k]*b. It is therefore at least the tape's gramcover(t, outidxs) = Ũ'*ones(m), which is the default h_gn of cover_nls and in turn at least its per-residual quadrature (ngroups=nothing); NLS consumers should prefer cover_nls.

magnitude(x)

The largest absolute value x can take: the quantity the cover sweeps accumulate wherever they need the size of a value, partial or adjoint. For a number this is abs(x).

A type whose values stand for a set of numbers, such as an interval, extends this with the largest absolute value over the set. A program recorded on such values then yields covers that hold at every point of the set.

magtype(T)

The type of magnitude applied to a T, which is the type the cover sweeps compute in for a Tape{T}. It is T for a real floating-point type and the real type of a complex one.

boundtype(T)

The type the cover sweeps over a Tape{T} accumulate their bounds in. For a type whose values are points this is magtype(T), and the covers hold up to the rounding of the sweeps. A type whose values stand for sets, such as an interval, opts in to rigorous covers by returning an Upward type, so that every sum, product and square root the sweeps take rounds upward and the covers hold for every point of the set. The covers are returned in magtype(T) either way.

Upward{F}(x)

A nonnegative magnitude of floating-point type F whose arithmetic rounds upward: +, *, /, sqrt, hypot and integer powers return a value no smaller than the exact result. Each elementary operation rounds to nearest and then steps to nextfloat, except where the result is exact (an operand that is zero, or a factor or divisor that is one); x <= nextfloat(fl(x)) holds for any rounding to nearest, so this needs no control over the hardware rounding mode. hypot and the powers compose these operations, so they may step more than once.

A cover is a bound built from magnitudes by sums, products and square roots, each increasing in its operands, so rounding every step upward makes the computed cover at least the exact one. Subtraction is not increasing in its second operand and throws, and transcendental functions are not defined; a computation that reaches one fails rather than silently losing the bound. Division rounds upward too, and is an upper bound only where its divisor is fixed rather than computed, such as a chosen weight.

Mixing an Upward with a plain real promotes to Upward, so an operation rounds upward whenever either operand is one. Converting an Upward{F} to F is exact.

boundtype selects this type for the sweeps over a tape whose values stand for sets.

decided_sign(x)

1 if x is positive, -1 if it is negative, and 0 otherwise: at zero, and for a value whose sign cannot be determined, such as NaN. A type whose values stand for sets extends this to return 1 or -1 only when every member of the set has that sign, so that a set containing zero, or members of both signs, gives 0.

abs of a program variable uses this to decide whether the kink of abs at zero is reached.

decided_one(c)
decided_zero(c)

Whether the constant c is one (resp. zero): isone(c) and iszero(c) for a number. A type whose values stand for sets extends these to return true only for the set containing just that number, and false, never an error, when the comparison cannot be decided.

Operations with a constant use these to skip work an identity makes unnecessary, such as recording a node for x * c with a unit c, or the second partial of x^c with c zero or one. A false answer keeps the general form, which is always valid.

separation(x, y)

How far apart x and y are: abs(x - y) for numbers. A type whose values stand for sets extends this with the smallest distance between members of the two sets, which is zero whenever they overlap. The consistency checks of cover_hessian use it to compare two computations of the same gradient.