oxiproj-core 0.1.2

//! SMT-backed transform verification via OxiZ.
//!
//! Uses OxiZ's Linear Real Arithmetic (QF_LRA) solver to formally verify
//! properties of coordinate-transform parameters.
//!
//! # Feature
//!
//! This module is compiled only when the `oxiz` feature is enabled.

#![cfg(feature = "oxiz")]

use num_rational::Rational64;
use oxiz_core::ast::TermManager;
use oxiz_solver::{Solver, SolverResult};

/// Result of an SMT verification query.
#[derive(Debug, Clone, PartialEq)]
pub enum VerificationResult {
    /// The property holds for all inputs in the specified domain.
    Verified,
    /// The solver found a counterexample violating the property.
    Counterexample {
        /// The value of the first witness variable.
        x: f64,
        /// The value of the second witness variable (or 0.0 if unused).
        y: f64,
    },
    /// The solver timed out or could not determine satisfiability.
    Unknown,
}

/// Verify that a Helmert scale factor (in PPM) is within geodetically
/// reasonable bounds (|scale_ppm| ≤ 100 PPM).
///
/// Uses QF_LRA to assert ¬(−100 ≤ s ≤ 100) and checks for unsat (property
/// holds) vs sat (counterexample exists).
///
/// # Arguments
///
/// * `scale_ppm` — scale factor in parts per million; e.g. 1.0 means 1 PPM.
pub fn verify_helmert_scale_bounds(scale_ppm: f64) -> VerificationResult {
    let mut solver = Solver::new();
    let mut tm = TermManager::new();
    solver.set_logic("QF_LRA");

    let c100 = tm.mk_real(Rational64::new(100, 1));
    let cn100 = tm.mk_real(Rational64::new(-100, 1));

    // Approximate scale_ppm as rational with 6 decimal places
    let scale_num = (scale_ppm * 1_000_000.0).round() as i64;
    let s_val = tm.mk_real(Rational64::new(scale_num, 1_000_000));

    // Assert: s < -100 OR s > 100  (the negation of the property)
    let lt = tm.mk_lt(s_val, cn100);
    let gt = tm.mk_gt(s_val, c100);
    let violation = tm.mk_or([lt, gt]);
    solver.assert(violation, &mut tm);

    match solver.check(&mut tm) {
        SolverResult::Unsat => VerificationResult::Verified,
        SolverResult::Sat => VerificationResult::Counterexample {
            x: scale_ppm,
            y: 0.0,
        },
        SolverResult::Unknown => VerificationResult::Unknown,
    }
}

/// Verify that a Jacobian determinant `det` is non-zero (above `threshold`)
/// at a given evaluation point, confirming local invertibility.
///
/// Uses QF_LRA: asserts |det| ≤ threshold and checks unsat (invertible) vs
/// sat (singular or near-singular).
///
/// # Arguments
///
/// * `det` — Jacobian determinant at the evaluation point.
/// * `threshold` — minimum absolute value considered non-zero (e.g. 1e-10).
pub fn verify_invertible_at(det: f64, threshold: f64) -> VerificationResult {
    let mut solver = Solver::new();
    let mut tm = TermManager::new();
    solver.set_logic("QF_LRA");

    // Encode det with 6-digit precision
    let det_num = (det * 1_000_000.0).round() as i64;
    let det_val = tm.mk_real(Rational64::new(det_num, 1_000_000));

    // Encode threshold with 12-digit precision to handle small values like 1e-10
    let thr_num_raw = (threshold * 1_000_000_000_000.0).round() as i64;
    let thr_num = if thr_num_raw == 0 { 1 } else { thr_num_raw };
    let thr_val = tm.mk_real(Rational64::new(thr_num, 1_000_000_000_000));
    let neg_thr = tm.mk_neg(thr_val);

    // Assert violation: -threshold ≤ det ≤ threshold  (i.e. |det| ≤ threshold)
    let ge = tm.mk_ge(det_val, neg_thr);
    let le = tm.mk_le(det_val, thr_val);
    let singular = tm.mk_and([ge, le]);
    solver.assert(singular, &mut tm);

    match solver.check(&mut tm) {
        SolverResult::Unsat => VerificationResult::Verified,
        SolverResult::Sat => VerificationResult::Counterexample {
            x: det,
            y: threshold,
        },
        SolverResult::Unknown => VerificationResult::Unknown,
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn helmert_scale_within_bounds() {
        assert_eq!(
            verify_helmert_scale_bounds(0.0),
            VerificationResult::Verified
        );
        assert_eq!(
            verify_helmert_scale_bounds(50.0),
            VerificationResult::Verified
        );
        assert_eq!(
            verify_helmert_scale_bounds(-99.9),
            VerificationResult::Verified
        );
    }

    #[test]
    fn helmert_scale_out_of_bounds() {
        let result = verify_helmert_scale_bounds(200.0);
        assert!(matches!(result, VerificationResult::Counterexample { .. }));
        let result = verify_helmert_scale_bounds(-150.0);
        assert!(matches!(result, VerificationResult::Counterexample { .. }));
    }

    #[test]
    fn invertible_jacobian() {
        assert_eq!(
            verify_invertible_at(1.0, 1e-10),
            VerificationResult::Verified
        );
        assert_eq!(
            verify_invertible_at(-0.5, 1e-10),
            VerificationResult::Verified
        );
    }

    #[test]
    fn singular_jacobian() {
        let result = verify_invertible_at(0.0, 1e-10);
        assert!(matches!(result, VerificationResult::Counterexample { .. }));
    }
}