alkahest-cas 3.12.0

//! Tests for the verification gate itself.
//!
//! These call [`residual_is_zero`] directly with hand-built candidates rather
//! than going through `dsolve`, because the point is to pin what the *gate*
//! does with a candidate — including candidates no implemented solving class
//! would ever produce.  Both directions are pinned: a correct solution of a
//! singular ODE must still verify, and a wrong candidate whose blow-up used to
//! be skipped must now be caught.

use super::*;
use crate::kernel::Domain;

/// `y`, `y'`, … named as `OdeInput` names them, plus an integration constant.
struct Fixture {
    pool: ExprPool,
    x: ExprId,
    y: ExprId,
    c1: ExprId,
}

impl Fixture {
    fn new() -> Self {
        let pool = ExprPool::new();
        let x = pool.symbol("x", Domain::Real);
        let y = pool.symbol("y", Domain::Real);
        let c1 = pool.symbol("C1", Domain::Real);
        Fixture { pool, x, y, c1 }
    }

    fn neg(&self, e: ExprId) -> ExprId {
        self.pool.mul(vec![self.pool.integer(-1_i32), e])
    }

    /// `x − n/d`.
    fn x_minus(&self, n: i32, d: i32) -> ExprId {
        let r = self.pool.rational(n, d);
        self.pool.add(vec![self.x, self.neg(r)])
    }

    /// `√(1/2 − x)` — real for `x < 1/2`, NaN for the `0.61` and `0.79`
    /// samples, which is what makes an expression containing it "singular
    /// there" for the sampler.
    fn sqrt_half_minus_x(&self) -> ExprId {
        let half = self.pool.rational(1_i32, 2_i32);
        let arg = self.pool.add(vec![half, self.neg(self.x)]);
        self.pool.func("sqrt", vec![arg])
    }

    /// `eˣ·e⁻ˣ`, which is `1` but which `simplify` does not collapse.
    ///
    /// Multiplying a term by this is how these tests force the gate past its
    /// symbolic-zero branch and onto the numeric sampler — the same artefact
    /// the real corpus residuals carry (`x⁻¹·eˣ·e⁻ˣ − x⁻¹` for `y'' − y = 1/x`).
    fn one_that_does_not_simplify(&self) -> ExprId {
        let ex = self.pool.func("exp", vec![self.x]);
        let emx = self.pool.func("exp", vec![self.neg(self.x)]);
        self.pool.mul(vec![ex, emx])
    }
}

// ---------------------------------------------------------------------------
// The soundness direction: a wrong candidate must not be certified
// ---------------------------------------------------------------------------

/// The reproducer for the bug this module was restructured to fix.
///
/// The ODE is `y' = 0`, which is regular at every `x`; its general solution is
/// `y = C1`.  The candidate adds `√(x − ½)·(x − 0.61)²·(x − 0.79)²`, which is
///
/// * NaN at the samples `0.11, 0.27, 0.43` (negative radicand), and
/// * a double root at `0.61` and `0.79`, so the residual `y'` is exactly zero
///   at the two samples that *do* evaluate.
///
/// The old gate skipped the three NaN samples, saw `2 × 3 = 6` finite samples
/// all `≈ 0`, cleared its `≥ 6` bar and returned `Ok(())` — a wrong `y(x)`
/// certified.  The candidate blows up where the ODE is perfectly regular, so
/// the classifier now books those samples as blow-ups and refuses.
#[test]
fn wrong_candidate_blowing_up_at_a_regular_point_is_rejected() {
    let f = Fixture::new();
    let (input, yp) = OdeInput::first_order(f.x, f.y, &f.pool);
    let input = input.with_equation(yp);

    let root = f.pool.func("sqrt", vec![f.x_minus(-1, -2)]);
    let sq = |n: i32, d: i32| f.pool.pow(f.x_minus(n, d), f.pool.integer(2_i32));
    let wrong = simp(
        f.pool
            .add(vec![f.c1, f.pool.mul(vec![root, sq(61, 100), sq(79, 100)])]),
        &f.pool,
    );

    let (residual, derivs) = build_residual(&input, wrong, &f.pool).expect("residual builds");
    let report = numeric_report(&input, residual, &derivs, &[f.c1], &[], &f.pool);

    // The shape of the trap: exactly six agreeing samples, no disagreement, and
    // nine samples that the old code discarded.
    assert_eq!(report.agree, 6, "report: {report}");
    assert_eq!(report.disagree, 0, "report: {report}");
    assert_eq!(report.blowup_at_regular_point, 9, "report: {report}");
    assert_eq!(report.skipped_singular_ode, 0, "report: {report}");
    assert!(!report.certifies(), "report: {report}");

    residual_is_zero(&input, wrong, &[f.c1], &f.pool)
        .expect_err("a candidate that blows up where the ODE is regular must not certify");
}

/// A candidate that simply disagrees at finite samples is still rejected — the
/// restructure must not have loosened the ordinary path.  `y = C1·eˣ` is not a
/// solution of `y' = 2y`.
#[test]
fn plain_finite_disagreement_is_still_rejected() {
    let f = Fixture::new();
    let (input, yp) = OdeInput::first_order(f.x, f.y, &f.pool);
    // y' − 2·y·(eˣ·e⁻ˣ) = 0, i.e. y' = 2y, spelled so the symbolic branch misses.
    let two_y = f.pool.mul(vec![
        f.pool.integer(-2_i32),
        f.y,
        f.one_that_does_not_simplify(),
    ]);
    let input = input.with_equation(f.pool.add(vec![yp, two_y]));

    let cand = simp(
        f.pool.mul(vec![f.c1, f.pool.func("exp", vec![f.x])]),
        &f.pool,
    );
    let (residual, derivs) = build_residual(&input, cand, &f.pool).expect("residual builds");
    let report = numeric_report(&input, residual, &derivs, &[f.c1], &[], &f.pool);
    assert_eq!(report.agree, 0, "report: {report}");
    assert_eq!(report.disagree, 15, "report: {report}");
    assert!(!report.certifies(), "report: {report}");
}

// ---------------------------------------------------------------------------
// The capability direction: a correct solution must still be certified
// ---------------------------------------------------------------------------

/// A correct solution of an ODE that is *singular over part of the sample
/// range* must still verify.
///
/// The equation is `y' − y − √(½ − x)·eˣ·e⁻ˣ + √(½ − x) = 0`.  The last two
/// terms cancel as functions (`eˣ·e⁻ˣ = 1`), so the equation is `y' = y` and
/// `y = C1·eˣ` is exactly right — but as *written* the equation does not
/// evaluate at all for `x > ½`, so the samples `0.61` and `0.79` carry no
/// information about the candidate and must be skipped, not counted against it.
#[test]
fn correct_solution_of_a_singular_ode_still_verifies() {
    let f = Fixture::new();
    let (input, yp) = OdeInput::first_order(f.x, f.y, &f.pool);
    let root = f.sqrt_half_minus_x();
    let eq = f.pool.add(vec![
        yp,
        f.neg(f.y),
        f.neg(f.pool.mul(vec![root, f.one_that_does_not_simplify()])),
        root,
    ]);
    let input = input.with_equation(eq);

    let cand = simp(
        f.pool.mul(vec![f.c1, f.pool.func("exp", vec![f.x])]),
        &f.pool,
    );

    // The premise: the ODE really is unevaluable past the branch point, and
    // fine before it.
    assert!(ode_is_regular_at(&input, 0.43, &HashMap::new(), &f.pool));
    assert!(!ode_is_regular_at(&input, 0.61, &HashMap::new(), &f.pool));

    let (residual, derivs) = build_residual(&input, cand, &f.pool).expect("residual builds");
    let report = numeric_report(&input, residual, &derivs, &[f.c1], &[], &f.pool);
    assert_eq!(report.agree, 9, "report: {report}");
    assert_eq!(report.disagree, 0, "report: {report}");
    assert_eq!(report.blowup_at_regular_point, 0, "report: {report}");
    assert_eq!(report.skipped_singular_ode, 6, "report: {report}");
    assert!(report.certifies(), "report: {report}");

    residual_is_zero(&input, cand, &[f.c1], &f.pool)
        .expect("a correct solution of a singular ODE must still verify");
}

/// When the candidate blows up *and* the ODE is singular at the same sample,
/// the sample is still information-free and must stay a skip.
///
/// `y' = y·eˣ·e⁻ˣ/√(½ − x)` has the correct solution `y = C1·e^{−2√(½ − x)}`.
/// Both the equation and the solution are NaN for `x > ½`, so the ordering in
/// [`classify_nonfinite`] — singular ODE first — is what keeps this verifying.
#[test]
fn candidate_blowup_where_the_ode_is_also_singular_stays_a_skip() {
    let f = Fixture::new();
    let (input, yp) = OdeInput::first_order(f.x, f.y, &f.pool);
    let root = f.sqrt_half_minus_x();
    let inv_root = f.pool.pow(root, f.pool.integer(-1_i32));
    let rhs = f
        .pool
        .mul(vec![f.y, f.one_that_does_not_simplify(), inv_root]);
    let input = input.with_equation(f.pool.add(vec![yp, f.neg(rhs)]));

    // y = C1·exp(−2·√(½ − x))
    let expo = f.pool.mul(vec![f.pool.integer(-2_i32), root]);
    let cand = simp(
        f.pool.mul(vec![f.c1, f.pool.func("exp", vec![expo])]),
        &f.pool,
    );

    assert!(!ode_is_regular_at(&input, 0.61, &HashMap::new(), &f.pool));

    let (residual, derivs) = build_residual(&input, cand, &f.pool).expect("residual builds");
    let report = numeric_report(&input, residual, &derivs, &[f.c1], &[], &f.pool);
    assert_eq!(report.disagree, 0, "report: {report}");
    assert_eq!(report.blowup_at_regular_point, 0, "report: {report}");
    assert_eq!(report.skipped_singular_ode, 6, "report: {report}");
    assert_eq!(report.agree, 9, "report: {report}");

    residual_is_zero(&input, cand, &[f.c1], &f.pool)
        .expect("correct solution, ODE singular past the branch point");
}

// ---------------------------------------------------------------------------
// The regularity probe
// ---------------------------------------------------------------------------

/// `ode_is_regular_at` must answer about the *equation*, not about one state:
/// a single unlucky probe value may make a perfectly regular equation
/// unevaluable, so one finite probe is enough to call it regular.
#[test]
fn regularity_probe_ignores_an_unlucky_state() {
    let f = Fixture::new();
    let (input, yp) = OdeInput::first_order(f.x, f.y, &f.pool);
    // y' − log(y) = 0: NaN for the negative probe, finite for the positive ones.
    let eq = f.pool.add(vec![yp, f.neg(f.pool.func("log", vec![f.y]))]);
    let input = input.with_equation(eq);
    assert!(ode_is_regular_at(&input, 0.43, &HashMap::new(), &f.pool));

    // y' − log(½ − x) = 0 is genuinely unevaluable past the branch point, for
    // every state.
    let (input2, yp2) = OdeInput::first_order(f.x, f.y, &f.pool);
    let half = f.pool.rational(1_i32, 2_i32);
    let arg = f.pool.add(vec![half, f.neg(f.x)]);
    let eq2 = f.pool.add(vec![yp2, f.neg(f.pool.func("log", vec![arg]))]);
    let input2 = input2.with_equation(eq2);
    assert!(ode_is_regular_at(&input2, 0.43, &HashMap::new(), &f.pool));
    assert!(!ode_is_regular_at(&input2, 0.61, &HashMap::new(), &f.pool));
}

/// An equation the sampler cannot evaluate at all (an unknown special function)
/// is never called regular, so a candidate is never condemned on the strength
/// of a probe that did not run.
#[test]
fn regularity_probe_declines_on_an_unknown_construct() {
    let f = Fixture::new();
    let (input, yp) = OdeInput::first_order(f.x, f.y, &f.pool);
    let ei = f.pool.func("Ei", vec![f.x]);
    let input = input.with_equation(f.pool.add(vec![yp, f.neg(ei)]));
    assert!(!ode_is_regular_at(&input, 0.43, &HashMap::new(), &f.pool));
}

// ---------------------------------------------------------------------------
// Named constants are not free parameters
// ---------------------------------------------------------------------------

/// `pi` is an ordinary symbol here, so the collector used to hand it to the
/// sampler and the gate evaluated candidates at `π = 1.7`.  Anything whose
/// correctness depends on π being π then disagrees at every sample — a correct
/// answer refused.
#[test]
fn pi_is_bound_to_its_value_rather_than_sampled() {
    let f = Fixture::new();
    let (input, yp) = OdeInput::first_order(f.x, f.y, &f.pool);
    let pi = f.pool.symbol("pi", Domain::Real);
    // y' − cos(π + π)·y·(eˣ·e⁻ˣ) = 0.  `cos(2π) = 1`, so this is `y' = y` and
    // `y = C1·eˣ` is exactly right — but only at the real π.  At `π = 1.7` the
    // coefficient is `cos(3.4) = −0.967` and the candidate is wrong.
    let coeff = f.pool.func("cos", vec![f.pool.add(vec![pi, pi])]);
    let eq = f.pool.add(vec![
        yp,
        f.neg(f.pool.mul(vec![coeff, f.y, f.one_that_does_not_simplify()])),
    ]);
    let input = input.with_equation(eq);
    let cand = simp(
        f.pool.mul(vec![f.c1, f.pool.func("exp", vec![f.x])]),
        &f.pool,
    );

    let (residual, _) = build_residual(&input, cand, &f.pool).expect("residual builds");
    // Non-vacuity: `simp` must not have folded `cos(2π)` away, or the case
    // would never reach the sampler at all.
    assert!(
        contains(residual, pi, &f.pool),
        "the residual must still mention pi: {}",
        f.pool.display(residual)
    );
    assert!(
        free_parameters(&input, &[residual], &[f.c1], &f.pool).is_empty(),
        "pi must not be collected as a free parameter"
    );
    // `pi` resolves without ever being in the environment.
    let sin_pi = f.pool.func("sin", vec![pi]);
    assert!(
        eval(sin_pi, &HashMap::new(), &f.pool)
            .expect("pi has a value")
            .abs()
            < 1e-12,
        "sin(pi) must evaluate to 0"
    );

    residual_is_zero(&input, cand, &[f.c1], &f.pool)
        .expect("a solution that is correct at the real pi must verify");
}

/// The imaginary unit is the other symbol that already denotes a number.  Bound
/// to a real sample it turns `e^{iπ}` into an ordinary real exponential and the
/// gate reports a disagreement that is entirely an artefact of the binding.
#[test]
fn the_imaginary_unit_is_not_sampled_either() {
    let f = Fixture::new();
    let (input, yp) = OdeInput::first_order(f.x, f.y, &f.pool);
    let pi = f.pool.symbol("pi", Domain::Real);
    let i = f.pool.imaginary_unit();
    // y' + e^{iπ}·y·(eˣ·e⁻ˣ) = 0.  `e^{iπ} = −1`, so this is `y' = y`.
    let coeff = f.pool.func("exp", vec![f.pool.mul(vec![i, pi])]);
    let eq = f.pool.add(vec![
        yp,
        f.pool.mul(vec![coeff, f.y, f.one_that_does_not_simplify()]),
    ]);
    let input = input.with_equation(eq);
    let cand = simp(
        f.pool.mul(vec![f.c1, f.pool.func("exp", vec![f.x])]),
        &f.pool,
    );

    let (residual, _) = build_residual(&input, cand, &f.pool).expect("residual builds");
    assert!(
        contains(residual, i, &f.pool) && contains(residual, pi, &f.pool),
        "the residual must still mention both constants: {}",
        f.pool.display(residual)
    );
    assert!(
        free_parameters(&input, &[residual], &[f.c1], &f.pool).is_empty(),
        "neither pi nor the imaginary unit is a free parameter"
    );
    residual_is_zero(&input, cand, &[f.c1], &f.pool)
        .expect("e^{i*pi} = -1 must be read as -1, not as a sampled product");
}

// ---------------------------------------------------------------------------
// Parameters are sampled at both signs
// ---------------------------------------------------------------------------

/// The soundness hole a positive-only [`PARAM_SETS`] left open.
///
/// The equation is `y' = √(k²)·y`, whose general solution is `y = C1·e^{|k|x}`.
/// The candidate `y = C1·e^{kx}` is **right for every `k > 0` and wrong for
/// every `k < 0`** — the single commonest way a symbolic-coefficient class goes
/// wrong, since `√(k²) → k` is the step it is tempted to take.  Sampled only at
/// positive `k` the two are indistinguishable and the wrong candidate was
/// certified; the negative rows separate them.
#[test]
fn a_candidate_right_only_for_positive_parameters_is_refused() {
    let f = Fixture::new();
    let (input, yp) = OdeInput::first_order(f.x, f.y, &f.pool);
    let k = f.pool.symbol("k", Domain::Real);
    let abs_k = f
        .pool
        .func("sqrt", vec![f.pool.pow(k, f.pool.integer(2_i32))]);
    let input = input.with_equation(f.pool.add(vec![yp, f.neg(f.pool.mul(vec![abs_k, f.y]))]));

    let wrong = simp(
        f.pool.mul(vec![
            f.c1,
            f.pool.func("exp", vec![f.pool.mul(vec![k, f.x])]),
        ]),
        &f.pool,
    );
    let (residual, derivs) = build_residual(&input, wrong, &f.pool).expect("residual builds");
    // Non-vacuity: `simp` must not have rewritten `√(k²)` to `k`, which would
    // make the candidate correct as written and the case meaningless.
    assert!(
        !is_symbolic_zero(residual, &f.pool),
        "the residual must not already be zero: {}",
        f.pool.display(residual)
    );
    let params = free_parameters(&input, &[residual], &[f.c1], &f.pool);
    assert_eq!(params, vec![k], "k is the one free parameter");

    // On the positive rows alone the candidate is indistinguishable from the
    // truth — this is what the old table could see, and why it certified.
    let positive_only = parametric_report_over(
        &PARAM_SETS[..3],
        &input,
        residual,
        &derivs,
        &[f.c1],
        &params,
        &f.pool,
    );
    assert_eq!(positive_only.inner.disagree, 0, "{positive_only}");
    assert!(
        positive_only.certifies(),
        "the positive rows on their own certify a wrong candidate: {positive_only}"
    );

    // With both signs the disagreement is visible and the gate refuses.
    let full = parametric_report(&input, residual, &derivs, &[f.c1], &params, &f.pool);
    assert!(full.has_counterevidence(), "{full}");
    residual_is_zero(&input, wrong, &[f.c1], &f.pool)
        .expect_err("a candidate wrong for every k < 0 must not be certified");

    // …and the correct answer still is.
    let right = simp(
        f.pool.mul(vec![
            f.c1,
            f.pool.func("exp", vec![f.pool.mul(vec![abs_k, f.x])]),
        ]),
        &f.pool,
    );
    residual_is_zero(&input, right, &[f.c1], &f.pool)
        .expect("y = C1*e^{|k|x} is the general solution and must still verify");
}