Struct BytecodeTape 

pub struct BytecodeTape<F: Float> { /* private fields */ }

Available on crate feature bytecode only.

A bytecode tape that can be re-evaluated at different inputs.

Created via crate::api::record. After recording, call forward to re-evaluate and reverse to compute adjoints.

Implementations

impl<F: Float> BytecodeTape<F>

pub fn forward(&mut self, inputs: &[F])

Re-evaluate the tape at new inputs (forward sweep).

Overwrites values in-place — no allocation.

pub fn forward_nonsmooth(&mut self, inputs: &[F]) -> NonsmoothInfo<F>

Forward sweep with nonsmooth branch tracking.

Calls forward to evaluate the tape, then scans for nonsmooth operations and records which branch was taken at each one.

Tracked operations:

• Abs, Min, Max — kinks with nontrivial subdifferentials
• Signum, Floor, Ceil, Round, Trunc — step-function discontinuities (zero derivative on both sides, tracked for proximity detection only)

Returns crate::NonsmoothInfo containing all kink entries in tape order.

pub fn forward_into(&self, inputs: &[F], values_buf: &mut Vec<F>)

Forward evaluation into an external buffer.

Reads opcodes, constants, and argument indices from self, but writes computed values into values_buf instead of self.values. This allows parallel evaluation of the same tape at different inputs without cloning.

impl<F: Float> BytecodeTape<F>

pub fn jacobian(&mut self, inputs: &[F]) -> Vec<Vec<F>>

Compute the full Jacobian of a multi-output tape via reverse mode.

Performs m reverse sweeps (one per output). Returns J[i][j] = ∂f_i/∂x_j.

pub fn vjp_multi(&mut self, inputs: &[F], weights: &[F]) -> Vec<F>

Vector-Jacobian product for a multi-output tape.

Computes wᵀ · J where J is the Jacobian. More efficient than computing the full Jacobian when only the weighted combination is needed.

pub fn jacobian_limiting( &mut self, inputs: &[F], forced_signs: &[(u32, i8)], ) -> Vec<Vec<F>>

Compute a Jacobian with forced branch choices at specified tape indices.

For each (tape_index, sign) in forced_signs, the reverse sweep uses forced_reverse_partials instead of the standard partials at that index.

This is the building block for Clarke subdifferential enumeration.

pub fn clarke_jacobian( &mut self, inputs: &[F], tol: F, max_active_kinks: Option<usize>, ) -> Result<(NonsmoothInfo<F>, Vec<Vec<Vec<F>>>), ClarkeError>

Compute the Clarke generalized Jacobian via limiting Jacobian enumeration.

1. Runs forward_nonsmooth to detect all kink operations and their branches.
2. Identifies “active” kinks (|switching_value| < tol).
3. Enumerates all 2^k sign combinations for the k active kinks.
4. For each combination, computes a limiting Jacobian via forced reverse sweeps.

Returns the nonsmooth info and a vector of limiting Jacobians.

Errors

Returns crate::ClarkeError::TooManyKinks if the number of active kinks exceeds the limit (default 20, overridden by max_active_kinks).

pub fn jacobian_forward(&self, x: &[F]) -> Vec<Vec<F>>

Dense Jacobian via forward mode (one forward-tangent pass per input).

More efficient than reverse mode when num_inputs < num_outputs.

Custom ops are handled exactly (first order) via CustomOp::eval_dual, which propagates the tangent at the evaluation point x — so the Jacobian is correct at any x, not a linearization frozen at the recording inputs. (Second-order paths — hessian_vec, sparse_hessian_vec — still reject custom ops, since their linearizing sweep is only first-order-accurate.)

pub fn jacobian_cross_country(&mut self, inputs: &[F]) -> Vec<Vec<F>>

Dense Jacobian via cross-country (vertex) elimination.

Builds a linearized DAG from the tape, then eliminates intermediate vertices in Markowitz order. For functions where m ≈ n and the graph has moderate connectivity, this can require fewer operations than either pure forward mode (n passes) or reverse mode (m passes).

impl<F: Float> BytecodeTape<F>

pub fn dead_code_elimination(&mut self)

Remove unreachable nodes from the tape.

pub fn dead_code_elimination_for_outputs(&mut self, active_outputs: &[u32])

Eliminate dead code, keeping only the specified outputs alive.

Like dead_code_elimination but seeds reachability only from active_outputs. After compaction, output_indices contains only the active outputs (remapped), and output_index is set to the first active output.

Panics

Panics if active_outputs is empty or contains an out-of-range index.

pub fn cse(&mut self)

Common subexpression elimination.

Deduplicates identical (OpCode, arg0, arg1) triples, normalising argument order for commutative ops. Finishes with a DCE pass to remove the now-dead duplicates.

pub fn optimize(&mut self)

Run all tape optimizations: CSE followed by DCE.

In debug builds, validates internal consistency after optimization.

impl<F: Float> BytecodeTape<F>

pub fn reverse_seeded(&self, seeds: &[F]) -> Vec<F>

Reverse sweep with weighted seeds for multiple outputs.

Computes ∑_i weights[i] * ∂output_i/∂x — a vector-Jacobian product.

Returns the gradient with respect to all inputs (length num_inputs).

pub fn reverse(&self, seed_index: u32) -> Vec<F>

Reverse sweep: compute adjoints seeded at the output.

Returns the full adjoint vector (length = num_variables).

pub fn reverse_from(&self, values: &[F], seed_index: u32) -> Vec<F>

Reverse sweep reading from an external values buffer.

Like reverse but reads primal values from values instead of self.values. Pair with forward_into for parallel evaluation.

pub fn gradient(&mut self, inputs: &[F]) -> Vec<F>

Forward + reverse: compute the gradient at new inputs.

Returns only the input adjoints (indices 0..num_inputs).

pub fn gradient_with_buf( &mut self, inputs: &[F], adjoint_buf: &mut Vec<F>, ) -> Vec<F>

Like gradient but reuses a caller-provided buffer for the adjoint vector, avoiding allocation on repeated calls.

pub fn gradient_batch(&mut self, inputs: &[&[F]]) -> Vec<Vec<F>>

Evaluate the gradient at multiple input points.

Returns one gradient vector per input point.

impl<F: Float> BytecodeTape<F>

pub fn detect_sparsity(&self) -> SparsityPattern

Detect the structural sparsity pattern of the Hessian.

Walks the tape forward propagating input-dependency bitsets. At nonlinear operations, marks cross-pairs as potential Hessian interactions.

pub fn detect_jacobian_sparsity(&self) -> JacobianSparsityPattern

Detect the structural sparsity pattern of the Jacobian.

Walks the tape forward propagating input-dependency bitsets (first-order). For each output, determines which inputs it depends on.

pub fn sparse_hessian(&self, x: &[F]) -> (F, Vec<F>, SparsityPattern, Vec<F>)

Compute a sparse Hessian using structural sparsity detection and graph coloring.

Returns (value, gradient, pattern, hessian_values) where hessian_values[k] corresponds to (pattern.rows[k], pattern.cols[k]).

For problems with sparse Hessians, this requires only chromatic_number HVP calls instead of n, which can be dramatically fewer for banded or sparse interaction structures.

pub fn sparse_hessian_vec<const N: usize>( &self, x: &[F], ) -> (F, Vec<F>, SparsityPattern, Vec<F>)

Batched sparse Hessian: packs N colors per sweep using DualVec.

Reduces the number of forward+reverse sweeps from num_colors to ceil(num_colors / N). Each sweep processes N colors simultaneously.

Custom ops limitation: For tapes containing custom ops, this method uses first-order chain rule (linearized partials). For exact second-order derivatives through custom ops, use sparse_hessian instead, which calls CustomOp::eval_dual / CustomOp::partials_dual.

Returns (value, gradient, pattern, hessian_values).

pub fn sparse_jacobian( &mut self, x: &[F], ) -> (Vec<F>, JacobianSparsityPattern, Vec<F>)

Compute a sparse Jacobian using structural sparsity detection and graph coloring.

Auto-selects forward-mode (column compression) or reverse-mode (row compression) based on which requires fewer sweeps.

Returns (output_values, pattern, jacobian_values).

pub fn sparse_jacobian_forward( &mut self, x: &[F], ) -> (Vec<F>, JacobianSparsityPattern, Vec<F>)

Sparse Jacobian via forward-mode (column compression).

pub fn sparse_jacobian_reverse( &mut self, x: &[F], ) -> (Vec<F>, JacobianSparsityPattern, Vec<F>)

Sparse Jacobian via reverse-mode (row compression).

pub fn sparse_jacobian_with_pattern( &mut self, x: &[F], pattern: &JacobianSparsityPattern, colors: &[u32], num_colors: u32, forward_mode: bool, ) -> (Vec<F>, Vec<F>)

Sparse Jacobian with a precomputed sparsity pattern and coloring.

Skips re-detection of sparsity on repeated calls. Use column_coloring colors for forward mode or row_coloring colors for reverse mode. The forward flag selects the mode.

pub fn sparse_jacobian_vec<const N: usize>( &mut self, x: &[F], ) -> (Vec<F>, JacobianSparsityPattern, Vec<F>)

Batched sparse Jacobian: packs N colors per forward sweep using DualVec.

Reduces the number of forward sweeps from num_colors to ceil(num_colors / N).

pub fn sparse_hessian_with_pattern( &self, x: &[F], pattern: &SparsityPattern, colors: &[u32], num_colors: u32, ) -> (F, Vec<F>, Vec<F>)

Sparse Hessian with a precomputed sparsity pattern and coloring.

Skips re-detection on repeated calls (e.g. in solver loops).

impl<F: Float> BytecodeTape<F>

pub fn forward_tangent<T: NumFloat>(&self, inputs: &[T], buf: &mut Vec<T>)

Forward sweep with tangent-carrying numbers. Reads opcodes and constants from self, writing results into buf. Does not mutate the tape.

Generic over T: NumFloat so it works with both Dual<F> and DualVec<F, N>.

Custom-op accuracy

Custom operations use recording-time primals (self.values) for their first-order linearization. If the tape has been re-evaluated at different inputs via forward() but self.values was not updated to match the tangent inputs, the custom-op linearization point will be stale, producing O(||x - x_record||) errors in the tangent output. For exact derivatives through custom ops, use the specializations forward_tangent_dual (Dual<F>) or forward_tangent_dual2 (Dual<Dual<F>>), which route through CustomOp::eval_dual / CustomOp::partials_dual.

pub fn forward_tangent_dual2( &self, inputs: &[Dual<Dual<F>>], buf: &mut Vec<Dual<Dual<F>>>, )

Forward sweep specialized for Dual<Dual<F>> (forward-over-forward), propagating exact first- and second-order information through custom ops.

The inner dual level is evaluated with CustomOp::eval_dual at the current inputs, and the outer tangent applies the chain rule with partials from CustomOp::partials_dual — duals whose tangent components carry the partials’ own derivatives, so the custom op’s curvature lands in the eps.eps component exactly as the built-in arms produce it. Accuracy therefore matches what the op implements: ops overriding eval_dual/partials_dual are exact to second order; ops relying on the trait defaults degrade to constant partials (zero curvature) but are still evaluated at the current point. The generic forward_tangent would instead linearize around recording-time primals, which both drops curvature and is stale away from the recording point.

pub fn forward_tangent_dual(&self, inputs: &[Dual<F>], buf: &mut Vec<Dual<F>>)

Forward sweep specialized for Dual<F>, calling CustomOp::eval_dual so custom ops propagate exact first-order tangents at the current inputs. The generic forward_tangent instead linearizes custom ops around recording-time primals, which is stale away from the recording point.

pub fn hvp(&self, x: &[F], v: &[F]) -> (Vec<F>, Vec<F>)

Hessian-vector product via forward-over-reverse.

Returns (gradient, H·v) where both are Vec<F> of length num_inputs. The tape is not mutated.

pub fn hvp_with_buf( &self, x: &[F], v: &[F], dual_vals_buf: &mut Vec<Dual<F>>, adjoint_buf: &mut Vec<Dual<F>>, ) -> (Vec<F>, Vec<F>)

Like hvp but reuses caller-provided buffers to avoid allocation on repeated calls (e.g. inside hessian).

pub fn hessian(&self, x: &[F]) -> (F, Vec<F>, Vec<Vec<F>>)

Full Hessian matrix via n Hessian-vector products.

Returns (value, gradient, hessian) where hessian[i][j] = ∂²f/∂x_i∂x_j. The tape is not mutated.

pub fn hessian_vec<const N: usize>(&self, x: &[F]) -> (F, Vec<F>, Vec<Vec<F>>)

Full Hessian matrix via batched forward-over-reverse.

Processes ceil(n/N) batches instead of n individual HVPs, computing N Hessian columns simultaneously.

Custom ops limitation: For tapes containing custom ops, this method uses first-order chain rule (linearized partials). For exact second-order derivatives through custom ops, use hessian instead, which calls CustomOp::eval_dual / CustomOp::partials_dual.

pub fn third_order_hvvp( &self, x: &[F], v1: &[F], v2: &[F], ) -> (Vec<F>, Vec<F>, Vec<F>)

Third-order directional derivative: ∑_{jk} (∂³f/∂x_i∂x_j∂x_k) v1_j v2_k.

Given directions v1 and v2, computes:

• gradient: ∇f(x)
• hvp: H(x) · v1 (Hessian-vector product)
• third: (∂/∂v2)(H · v1) (third-order tensor contracted with v1 and v2)

Uses Dual<Dual<F>> (nested dual numbers): inner tangent for v1, outer tangent for v2.

impl<F: Float> BytecodeTape<F>

pub fn gradient_par(&self, inputs: &[F]) -> Vec<F>

Parallel gradient: forward + reverse using external buffers.

Takes &self instead of &mut self, enabling shared access across threads.

pub fn jacobian_par(&self, inputs: &[F]) -> Vec<Vec<F>>

Parallel Jacobian: one reverse sweep per output, parallelized.

Returns J[i][j] = ∂f_i/∂x_j.

pub fn hessian_par(&self, x: &[F]) -> (F, Vec<F>, Vec<Vec<F>>)

Parallel Hessian: one HVP per column, parallelized over columns.

Returns (value, gradient, hessian).

pub fn sparse_hessian_par( &self, x: &[F], ) -> (F, Vec<F>, SparsityPattern, Vec<F>)

Parallel sparse Hessian: parallelized over colors.

Returns (value, gradient, pattern, hessian_values).

pub fn sparse_jacobian_par( &self, x: &[F], ) -> (Vec<F>, JacobianSparsityPattern, Vec<F>)

Parallel sparse Jacobian: parallelized over colors.

Auto-selects forward (column compression) or reverse (row compression) based on num_outputs vs num_inputs.

pub fn gradient_batch_par(&self, inputs: &[&[F]]) -> Vec<Vec<F>>

Evaluate the gradient at multiple input points in parallel.

Uses forward_into + reverse_from with per-thread buffers.

pub fn hessian_batch_par( &self, inputs: &[&[F]], ) -> Vec<(F, Vec<F>, Vec<Vec<F>>)>

Compute Hessian at multiple input points in parallel.

Returns (value, gradient, hessian) for each input point.

impl<F: Float> BytecodeTape<F>

pub fn taylor_grad<const K: usize>( &self, x: &[F], v: &[F], ) -> (Taylor<F, K>, Vec<Taylor<F, K>>)

Forward-reverse Taylor pass for gradient + higher-order directional adjoints.

Builds Taylor inputs x_i(t) = x_i + v_i * t (with zero higher coefficients), runs forward_tangent, then reverse_tangent to get Taylor-valued adjoints.

Returns (output, adjoints) where:

• output is the Taylor expansion of f along direction v
• adjoints[i].coeff(0) = ∂f/∂x_i (gradient)
• adjoints[i].coeff(1) = Σ_j (∂²f/∂x_i∂x_j) v_j (HVP)
• adjoints[i].derivative(k) = k-th order directional adjoint

For K=2, the HVP component is equivalent to hvp. For K≥3, yields additional higher-order information in the same pass.

Like hvp, takes &self and does not call forward(x) before the Taylor pass. Custom ops will use primal values from recording time.

pub fn taylor_grad_with_buf<const K: usize>( &self, x: &[F], v: &[F], fwd_buf: &mut Vec<Taylor<F, K>>, adj_buf: &mut Vec<Taylor<F, K>>, ) -> (Taylor<F, K>, Vec<Taylor<F, K>>)

Like taylor_grad but reuses caller-provided buffers to avoid allocation on repeated calls.

pub fn ode_taylor_step<const K: usize>(&self, y0: &[F]) -> Vec<Taylor<F, K>>

Compute the Taylor expansion of the ODE solution y(t) to order K.

Given a tape representing the right-hand side f: R^n → R^n of the ODE y' = f(y), and an initial condition y(0) = y0, computes the Taylor coefficients y_0, y_1, ..., y_{K-1} such that y(t) ≈ y_0 + y_1·t + y_2·t² + ... + y_{K-1}·t^{K-1}.

The tape must have num_outputs == num_inputs (autonomous ODE: f maps R^n → R^n).

Returns one Taylor<F, K> per state variable. Use Taylor::eval_at to evaluate at a step size h, or inspect coefficients for error estimation.

pub fn ode_taylor_step_with_buf<const K: usize>( &self, y0: &[F], buf: &mut Vec<Taylor<F, K>>, ) -> Vec<Taylor<F, K>>

Like ode_taylor_step but reuses a caller-provided buffer to avoid allocation on repeated calls.

impl<F: Float> BytecodeTape<F>

pub fn new() -> Self

Create an empty bytecode tape.

pub fn with_capacity(est_ops: usize) -> Self

Create a bytecode tape with pre-allocated capacity.

pub fn new_input(&mut self, value: F) -> u32

Register a new input variable. Returns its index.

pub fn push_const(&mut self, value: F) -> u32

Register a scalar constant. Returns its index.

pub fn push_op(&mut self, op: OpCode, arg0: u32, arg1: u32, value: F) -> u32

Record an operation. Returns the result index.

Constant folding: if all operands point to Const entries (not Input), the operation is replaced by a single Const with the already-computed value.

Algebraic simplification: identity patterns that are exact under IEEE 754 for every future re-evaluation point (x * 1 → x, x / 1 → x, x + (-0.0) → x, x - (+0.0) → x) alias the live operand. Patterns whose result depends on the re-evaluable operand’s future value (x * 0, x - x, x / x, x + (+0.0)) are not folded: the tape’s defining contract is re-evaluation at new inputs, and freezing a recording-time result would silently mask singularities (x / x replayed at x = 0 is NaN, not 1) or flip signed zeros.

pub fn push_powi(&mut self, arg0: u32, exp: i32, value: F) -> u32

Record a powi operation. The i32 exponent is stored in arg_indices[1].

Constant folding: if the operand is a Const, emit Const instead.

Algebraic simplification: x^0 → 1 (guarded), x^1 → x, x^(-1) → Recip(x) (cheaper unary dispatch).

pub fn register_custom(&mut self, op: Arc<dyn CustomOp<F>>) -> CustomOpHandle

Register a custom operation. Returns a handle for use with crate::BReverse::custom_unary and crate::BReverse::custom_binary.

pub fn push_custom_unary( &mut self, arg0: u32, handle: CustomOpHandle, value: F, ) -> u32

Record a unary custom op. arg_indices = [arg0, callback_idx].

pub fn push_custom_binary( &mut self, arg0: u32, arg1: u32, handle: CustomOpHandle, value: F, ) -> u32

Record a binary custom op. arg_indices = [arg0, callback_idx], second operand stored in custom_second_args.

pub fn set_output(&mut self, index: u32)

Mark the output variable.

pub fn output_value(&self) -> F

Get the output value (available after forward() or initial recording).

pub fn output_index(&self) -> usize

Index of the (single) output variable.

Use this with the buffer produced by forward_tangent to read the output: buf[tape.output_index()].

pub fn num_inputs(&self) -> usize

Number of input variables.

pub fn num_ops(&self) -> usize

Number of operations (including inputs and constants).

pub fn set_outputs(&mut self, indices: &[u32])

Mark multiple output variables.

When set, num_outputs, output_values, jacobian, and vjp_multi become available. Single-output methods (output_index, gradient, etc.) continue to work using the first output.

pub fn num_outputs(&self) -> usize

Number of output variables. Returns 1 in single-output mode.

pub fn output_values(&self) -> Vec<F>

Get all output values (available after forward() or initial recording).

In single-output mode, returns a single-element vector.

pub fn all_output_indices(&self) -> &[u32]

Indices of all output entries in the tape buffer.

For multi-output tapes, returns all registered output indices. For single-output tapes, returns a single-element slice.

pub fn validate(&self) -> Result<(), TapeValidationError>

Check the tape’s structural invariants.

A well-formed tape — anything produced by crate::api::record, the push_* builders, or optimize — always passes. Err means the tape data was corrupted (most likely deserialized from tampered or truncated bytes); evaluating such a tape would panic on an out-of-bounds slot or silently produce meaningless derivatives from uncomputed slots.

Checked invariants:

• opcodes, arg_indices, and values each have exactly num_variables entries;
• the Input opcodes are exactly the first num_inputs entries;
• output_index and every entry of output_indices name a real tape slot;
• every operand index references a strictly earlier slot (the tape is stored in topological order), unary ops carry the UNUSED sentinel in their second slot, and custom-op callback indices are registered;
• custom_second_args keys name Custom ops and their values reference strictly earlier slots.

pub fn opcodes_slice(&self) -> &[OpCode]

Slice view of all opcodes in the tape.

pub fn arg_indices_slice(&self) -> &[[u32; 2]]

Slice view of all argument index pairs [arg0, arg1].

pub fn values_slice(&self) -> &[F]

Slice view of all primal values in the tape.

pub fn num_variables_count(&self) -> usize

Total number of tape entries (inputs + constants + operations).

pub fn has_custom_ops(&self) -> bool

Returns true if the tape contains any custom operations.

Trait Implementations

impl<F: Float> Default for BytecodeTape<F>

fn default() -> Self

Returns the “default value” for a type.

impl<'de, F: Float + Deserialize<'de>> Deserialize<'de> for BytecodeTape<F>

fn deserialize<D: Deserializer<'de>>(deserializer: D) -> Result<Self, D::Error>

Deserialize this value from the given Serde deserializer.

impl<F: Float + Serialize> Serialize for BytecodeTape<F>

fn serialize<S: Serializer>(&self, serializer: S) -> Result<S::Ok, S::Error>

Serialize this value into the given Serde serializer.