theory_core 0.1.0

//! SMT Solver Backend Integration (Stub Implementation)
//!
//! Provides integration with SMT solvers. Z3 integration is commented out
//! due to system dependencies - this is a stub implementation for Phase 2A demo.

use crate::ir::{ConstraintIR, TheoryTag};
use crate::surrogate::SolutionResult;
use anyhow::Result;
use serde::{Deserialize, Serialize};
use std::collections::HashMap;
use std::time::Instant;

/// Represents a bound value (inclusive or exclusive)
#[derive(Debug, Clone, Copy)]
enum Bound {
    Inclusive(i64),
    Exclusive(i64),
}

impl Bound {
    fn value(&self) -> i64 {
        match self {
            Bound::Inclusive(v) | Bound::Exclusive(v) => *v,
        }
    }
}

/// Represents bounds on a variable from inequality constraints
#[derive(Debug, Clone)]
struct VariableBounds {
    lower: Option<Bound>,
    upper: Option<Bound>,
    is_integer: bool,
}

impl VariableBounds {
    fn update_lower(&mut self, new_bound: Bound) {
        match (self.lower, new_bound) {
            (None, bound) => self.lower = Some(bound),
            (Some(Bound::Inclusive(old)), Bound::Inclusive(new)) if new > old => {
                self.lower = Some(Bound::Inclusive(new))
            }
            (Some(Bound::Exclusive(old)), Bound::Exclusive(new)) if new > old => {
                self.lower = Some(Bound::Exclusive(new))
            }
            (Some(Bound::Inclusive(old)), Bound::Exclusive(new)) if new >= old => {
                self.lower = Some(Bound::Exclusive(new))
            }
            (Some(Bound::Exclusive(old)), Bound::Inclusive(new)) if new > old => {
                self.lower = Some(Bound::Inclusive(new))
            }
            _ => {}
        }
    }

    fn update_upper(&mut self, new_bound: Bound) {
        match (self.upper, new_bound) {
            (None, bound) => self.upper = Some(bound),
            (Some(Bound::Inclusive(old)), Bound::Inclusive(new)) if new < old => {
                self.upper = Some(Bound::Inclusive(new))
            }
            (Some(Bound::Exclusive(old)), Bound::Exclusive(new)) if new < old => {
                self.upper = Some(Bound::Exclusive(new))
            }
            (Some(Bound::Inclusive(old)), Bound::Exclusive(new)) if new <= old => {
                self.upper = Some(Bound::Exclusive(new))
            }
            (Some(Bound::Exclusive(old)), Bound::Inclusive(new)) if new < old => {
                self.upper = Some(Bound::Inclusive(new))
            }
            _ => {}
        }
    }

    fn lower_value(&self) -> Option<i64> {
        self.lower.as_ref().map(|b| match b {
            Bound::Inclusive(v) => *v,
            Bound::Exclusive(v) => v + 1, // For integers, x > 5 means x >= 6
        })
    }

    fn upper_value(&self) -> Option<i64> {
        self.upper.as_ref().map(|b| match b {
            Bound::Inclusive(v) => *v,
            Bound::Exclusive(v) => v - 1, // For integers, x < 5 means x <= 4
        })
    }

    /// Format the bounds as a range string
    fn format_range(&self, var_name: &str) -> String {
        match (self.lower, self.upper) {
            (Some(lower), Some(upper)) => {
                let lower_op = match lower {
                    Bound::Inclusive(_) => "<=",
                    Bound::Exclusive(_) => "<",
                };
                let upper_op = match upper {
                    Bound::Inclusive(_) => "<=",
                    Bound::Exclusive(_) => "<",
                };
                format!(
                    "{} {} {} {} {}",
                    lower.value(),
                    lower_op,
                    var_name,
                    upper_op,
                    upper.value()
                )
            }
            (Some(lower), None) => {
                let op = match lower {
                    Bound::Inclusive(_) => "<=",
                    Bound::Exclusive(_) => "<",
                };
                format!("{} {} {}", lower.value(), op, var_name)
            }
            (None, Some(upper)) => {
                let op = match upper {
                    Bound::Inclusive(_) => "<=",
                    Bound::Exclusive(_) => "<",
                };
                format!("{} {} {}", var_name, op, upper.value())
            }
            (None, None) => format!("{} = any value", var_name),
        }
    }
}

/// SMT solver backend configuration
#[derive(Debug, Clone, Serialize, Deserialize)]
pub struct SmtBackendConfig {
    /// Timeout for solver calls in milliseconds
    pub timeout_ms: u64,
    /// Maximum number of solver instances to pool
    pub max_solver_pool: usize,
    /// Enable incremental solving
    pub incremental: bool,
    /// Z3 specific options
    pub z3_options: HashMap<String, String>,
}

impl Default for SmtBackendConfig {
    fn default() -> Self {
        let mut z3_options = HashMap::new();
        z3_options.insert("model".to_string(), "true".to_string());
        z3_options.insert("proof".to_string(), "false".to_string());

        Self {
            timeout_ms: 5000,
            max_solver_pool: 4,
            incremental: true,
            z3_options,
        }
    }
}

/// SMT solver backend (stub implementation)
pub struct SmtBackend {
    config: SmtBackendConfig,
}

impl SmtBackend {
    pub fn new(config: SmtBackendConfig) -> Result<Self> {
        Ok(Self { config })
    }

    /// Solve constraints using appropriate SMT solver (stub)
    pub async fn solve(&self, constraints: &ConstraintIR) -> Result<SolutionResult> {
        let start_time = Instant::now();

        // Route to appropriate solver based on theory tags
        let solver_type = self.select_solver(&constraints.theory_tags);

        let result = match solver_type {
            SolverType::Z3 => self.solve_with_stub(constraints).await?,
            SolverType::CVC5 => {
                // TODO: Implement CVC5 integration
                return Err(anyhow::anyhow!("CVC5 not yet implemented"));
            }
        };

        let solve_time = start_time.elapsed();

        // Update result with actual solving time
        let mut final_result = result;
        final_result.estimated_time_ms = solve_time.as_millis() as f64;

        Ok(final_result)
    }

    /// Select appropriate solver based on theory tags
    fn select_solver(&self, _theory_tags: &[TheoryTag]) -> SolverType {
        // For now, always use Z3 stub as it's the most general-purpose
        SolverType::Z3
    }

    /// Solve using stub implementation (placeholder for Z3)
    async fn solve_with_stub(&self, constraints: &ConstraintIR) -> Result<SolutionResult> {
        use crate::ir::{BinaryOp, Constraint, Expr};

        // Try to solve as inequality ranges first (if only inequalities present)
        if let Ok(solution) = self.solve_inequality_ranges(constraints) {
            return Ok(solution);
        }

        // Try to solve as a linear system first
        match self.solve_linear_system(constraints) {
            Ok(solution) => {
                // Verify the solution satisfies all constraints
                if self.verify_solution(constraints, &solution) {
                    return Ok(SolutionResult {
                        satisfiable: true,
                        assignment: solution,
                        confidence: 0.95,
                        estimated_time_ms: 0.0,
                        method: "LinearSystemSolver".to_string(),
                        winning_branch_id: Some("LinearSystemSolver".to_string()),
                        winning_strategy: Some("LinearAlgebra".to_string()),
                        winning_backend: Some("CPU".to_string()),
                        variable_ranges: None,
                    });
                }
            }
            Err(_) => {
                // Linear solver failed, will try non-linear
            }
        }

        // Try non-linear solver for small domains
        match self.solve_nonlinear_system(constraints) {
            Ok(solution) => {
                return Ok(SolutionResult {
                    satisfiable: true,
                    assignment: solution,
                    confidence: 0.95,
                    estimated_time_ms: 0.0,
                    method: "NonLinearSolver".to_string(),
                    winning_branch_id: Some("NonLinearSolver".to_string()),
                    winning_strategy: Some("BruteForceSearch".to_string()),
                    winning_backend: Some("CPU".to_string()),
                    variable_ranges: None,
                });
            }
            Err(_) => {
                // Non-linear solver failed, will try fallback
            }
        }

        // Fall back to simple constraint solver
        let mut assignment: HashMap<String, crate::surrogate::SolutionValue> = HashMap::new(); // Initialize variables with default values
        for (var_id, variable) in &constraints.variables {
            let var_name = format!("var_{}", var_id.0);
            let value = match &variable.domain {
                crate::ir::Domain::Boolean => crate::surrogate::SolutionValue::Bool(false),
                crate::ir::Domain::Integer { min, max: _ } => {
                    crate::surrogate::SolutionValue::Integer(min.unwrap_or(0))
                }
                crate::ir::Domain::Real { min, max: _ } => {
                    crate::surrogate::SolutionValue::Real(min.unwrap_or(0.0))
                }
                crate::ir::Domain::BitVector { width } => {
                    crate::surrogate::SolutionValue::BitVector {
                        value: 0,
                        width: *width,
                    }
                }
                crate::ir::Domain::String { max_length: _ } => {
                    crate::surrogate::SolutionValue::String("stub_value".to_string())
                }
                _ => crate::surrogate::SolutionValue::Integer(0),
            };
            assignment.insert(var_name, value);
        }

        // Process constraints to solve for variable values
        for constraint in &constraints.constraints {
            match constraint {
                Constraint::Equal { lhs, rhs } => {
                    // Try to solve: var = expr
                    if let Expr::Var(var_id) = lhs {
                        if let Ok(value) = self.evaluate_expr(rhs, &assignment) {
                            assignment.insert(format!("var_{}", var_id.0), value);
                        }
                    } else if let Expr::Var(var_id) = rhs {
                        if let Ok(value) = self.evaluate_expr(lhs, &assignment) {
                            assignment.insert(format!("var_{}", var_id.0), value);
                        }
                    }
                }
                _ => {
                    // Other constraint types not yet implemented in stub
                }
            }
        }

        // Verify the solution satisfies all constraints
        let satisfiable = self.verify_solution(constraints, &assignment);

        Ok(SolutionResult {
            satisfiable,
            assignment,
            confidence: if satisfiable { 0.9 } else { 0.0 },
            estimated_time_ms: 0.0, // Will be updated by caller
            method: "SMTStub".to_string(),
            winning_branch_id: Some("SMTStub".to_string()),
            winning_strategy: Some("SMT".to_string()),
            winning_backend: Some("CPU".to_string()),
            variable_ranges: None,
        })
    }

    /// Verify that an assignment satisfies all constraints
    fn verify_solution(
        &self,
        constraints: &ConstraintIR,
        assignment: &HashMap<String, crate::surrogate::SolutionValue>,
    ) -> bool {
        use crate::ir::Constraint;

        for constraint in &constraints.constraints {
            match constraint {
                Constraint::Equal { lhs, rhs } => {
                    let lhs_val = self.evaluate_expr(lhs, assignment);
                    let rhs_val = self.evaluate_expr(rhs, assignment);

                    match (lhs_val, rhs_val) {
                        (Ok(l), Ok(r)) => {
                            if !self.values_equal(&l, &r) {
                                return false;
                            }
                        }
                        _ => return false, // Can't evaluate, constraint not satisfied
                    }
                }
                Constraint::Less { lhs, rhs } => {
                    let lhs_val = self.evaluate_expr(lhs, assignment);
                    let rhs_val = self.evaluate_expr(rhs, assignment);

                    match (lhs_val, rhs_val) {
                        (Ok(l), Ok(r)) => {
                            if !self.value_less_than(&l, &r) {
                                return false;
                            }
                        }
                        _ => return false,
                    }
                }
                _ => {
                    // Other constraint types - assume satisfied for now
                }
            }
        }
        true
    }

    /// Check if two solution values are equal
    fn values_equal(
        &self,
        a: &crate::surrogate::SolutionValue,
        b: &crate::surrogate::SolutionValue,
    ) -> bool {
        use crate::surrogate::SolutionValue;
        match (a, b) {
            (SolutionValue::Integer(x), SolutionValue::Integer(y)) => x == y,
            (SolutionValue::Real(x), SolutionValue::Real(y)) => (x - y).abs() < 1e-9,
            (SolutionValue::Bool(x), SolutionValue::Bool(y)) => x == y,
            (SolutionValue::String(x), SolutionValue::String(y)) => x == y,
            _ => false,
        }
    }

    /// Check if value a < value b
    fn value_less_than(
        &self,
        a: &crate::surrogate::SolutionValue,
        b: &crate::surrogate::SolutionValue,
    ) -> bool {
        use crate::surrogate::SolutionValue;
        match (a, b) {
            (SolutionValue::Integer(x), SolutionValue::Integer(y)) => x < y,
            (SolutionValue::Real(x), SolutionValue::Real(y)) => x < y,
            _ => false,
        }
    }

    /// Solve constraints that are purely inequalities defining ranges
    /// Returns a SolutionResult with ranges instead of single values
    fn solve_inequality_ranges(&self, constraints: &ConstraintIR) -> Result<SolutionResult> {
        use crate::ir::Constraint;

        // Check if all constraints are inequalities
        let all_inequalities = constraints
            .constraints
            .iter()
            .all(|c| matches!(c, Constraint::Less { .. } | Constraint::LessEqual { .. }));

        if !all_inequalities || constraints.constraints.is_empty() {
            return Err(anyhow::anyhow!("Not all constraints are inequalities"));
        }

        // Extract bounds for each variable
        let mut var_bounds: HashMap<String, VariableBounds> = HashMap::new();

        // Initialize with variable domains
        for (var_id, variable) in &constraints.variables {
            let var_name = format!("var_{}", var_id.0);
            let bounds = match &variable.domain {
                crate::ir::Domain::Integer { min, max } => VariableBounds {
                    lower: min.map(|v| Bound::Inclusive(v)),
                    upper: max.map(|v| Bound::Inclusive(v)),
                    is_integer: true,
                },
                crate::ir::Domain::Real { min, max } => VariableBounds {
                    lower: min.map(|v| Bound::Inclusive(v as i64)),
                    upper: max.map(|v| Bound::Inclusive(v as i64)),
                    is_integer: false,
                },
                _ => return Err(anyhow::anyhow!("Unsupported domain type for inequalities")),
            };
            var_bounds.insert(var_name, bounds);
        }

        // Extract bounds from inequality constraints
        for constraint in &constraints.constraints {
            match constraint {
                Constraint::Less { lhs, rhs } => {
                    self.extract_bounds_from_inequality(lhs, rhs, false, &mut var_bounds)?;
                }
                Constraint::LessEqual { lhs, rhs } => {
                    self.extract_bounds_from_inequality(lhs, rhs, true, &mut var_bounds)?;
                }
                _ => {}
            }
        }

        // Create assignment with range information
        let mut assignment = HashMap::new();
        let mut ranges = HashMap::new();
        for (var_name, bounds) in &var_bounds {
            // For now, pick a value in the middle of the range
            let value = self.pick_value_from_bounds(bounds)?;
            assignment.insert(var_name.clone(), value);

            // Store range information for display
            ranges.insert(var_name.clone(), bounds.format_range(var_name));
        }

        // Verify the solution
        let satisfiable = self.verify_solution(constraints, &assignment);

        Ok(SolutionResult {
            satisfiable,
            assignment,
            confidence: if satisfiable { 0.95 } else { 0.0 },
            estimated_time_ms: 0.0,
            method: "InequalityRangeSolver".to_string(),
            winning_branch_id: Some("InequalityRangeSolver".to_string()),
            winning_strategy: Some("IntervalArithmetic".to_string()),
            winning_backend: Some("CPU".to_string()),
            variable_ranges: Some(ranges),
        })
    }

    /// Extract bounds from an inequality constraint
    fn extract_bounds_from_inequality(
        &self,
        lhs: &crate::ir::Expr,
        rhs: &crate::ir::Expr,
        inclusive: bool,
        var_bounds: &mut HashMap<String, VariableBounds>,
    ) -> Result<()> {
        use crate::ir::Expr;

        // Handle case: var < const or var <= const
        if let Expr::Var(var_id) = lhs {
            if let Expr::Const(val) = rhs {
                let var_name = format!("var_{}", var_id.0);
                if let Some(bounds) = var_bounds.get_mut(&var_name) {
                    let bound_val = self.const_to_i64(val)?;
                    let new_upper = if inclusive {
                        Bound::Inclusive(bound_val)
                    } else {
                        Bound::Exclusive(bound_val)
                    };
                    bounds.update_upper(new_upper);
                }
            }
        }

        // Handle case: const < var or const <= var
        if let Expr::Const(val) = lhs {
            if let Expr::Var(var_id) = rhs {
                let var_name = format!("var_{}", var_id.0);
                if let Some(bounds) = var_bounds.get_mut(&var_name) {
                    let bound_val = self.const_to_i64(val)?;
                    let new_lower = if inclusive {
                        Bound::Inclusive(bound_val)
                    } else {
                        Bound::Exclusive(bound_val)
                    };
                    bounds.update_lower(new_lower);
                }
            }
        }

        Ok(())
    }

    /// Convert a ConstValue to i64
    fn const_to_i64(&self, val: &crate::ir::ConstValue) -> Result<i64> {
        match val {
            crate::ir::ConstValue::Integer(i) => Ok(*i),
            crate::ir::ConstValue::Real(r) => Ok(*r as i64),
            _ => Err(anyhow::anyhow!("Cannot convert const to i64")),
        }
    }

    /// Pick a value from bounds (middle of range)
    fn pick_value_from_bounds(
        &self,
        bounds: &VariableBounds,
    ) -> Result<crate::surrogate::SolutionValue> {
        let lower = bounds.lower_value().unwrap_or(i64::MIN);
        let upper = bounds.upper_value().unwrap_or(i64::MAX);

        if lower > upper {
            return Err(anyhow::anyhow!("Invalid bounds: lower > upper"));
        }

        let middle = if lower == i64::MIN && upper == i64::MAX {
            0 // No bounds specified
        } else if lower == i64::MIN {
            upper - 1
        } else if upper == i64::MAX {
            lower + 1
        } else {
            (lower + upper) / 2
        };

        Ok(crate::surrogate::SolutionValue::Integer(middle))
    }

    /// Attempt to solve constraints as a linear system
    /// Returns Ok with solution if it's a solvable linear system, Err otherwise
    fn solve_linear_system(
        &self,
        constraints: &ConstraintIR,
    ) -> Result<HashMap<String, crate::surrogate::SolutionValue>> {
        use crate::ir::Constraint;

        // Extract linear equations: a*x + b*y + ... = c
        let mut equations: Vec<(Vec<(usize, i64)>, i64)> = Vec::new(); // (coefficients, constant)
        let num_vars = constraints.variables.len();

        for constraint in &constraints.constraints {
            if let Constraint::Equal { lhs, rhs } = constraint {
                // Try to extract linear equation
                if let Some(equation) = self.extract_linear_equation(lhs, rhs, num_vars) {
                    equations.push(equation);
                } else {
                    // Not a linear equation, can't solve as linear system
                    return Err(anyhow::anyhow!("Not a linear system"));
                }
            } else {
                return Err(anyhow::anyhow!("Non-equality constraint in system"));
            }
        }

        // Check if we have enough equations
        if equations.len() < num_vars {
            return Err(anyhow::anyhow!("Underdetermined system"));
        }

        // Solve using Gaussian elimination
        self.gaussian_elimination(&equations, num_vars)
    }

    /// Solve non-linear constraints using brute-force search for small integer domains
    /// This handles systems like x * y = 12, x^2 + y = 10, etc.
    fn solve_nonlinear_system(
        &self,
        constraints: &ConstraintIR,
    ) -> Result<HashMap<String, crate::surrogate::SolutionValue>> {
        // Only proceed if all variables have small integer domains
        let mut var_ranges: Vec<(String, std::ops::RangeInclusive<i64>)> = Vec::new();

        for (var_id, variable) in &constraints.variables {
            let var_name = format!("var_{}", var_id.0);
            if let crate::ir::Domain::Integer { min, max } = &variable.domain {
                let min_val = min.unwrap_or(-20);
                let max_val = max.unwrap_or(20);
                let range_size = (max_val - min_val + 1) as usize;

                // Limit search space to prevent excessive computation
                if range_size > 10000 {
                    return Err(anyhow::anyhow!("Domain too large for brute force"));
                }

                var_ranges.push((var_name, min_val..=max_val));
            } else {
                // Non-integer domains not yet supported
                return Err(anyhow::anyhow!(
                    "Non-integer domain not supported for non-linear solver"
                ));
            }
        }

        // Check if search space is reasonable
        let total_combinations: usize = var_ranges
            .iter()
            .map(|(_, range)| range.clone().count())
            .product();

        if total_combinations > 1000000 {
            return Err(anyhow::anyhow!(
                "Search space too large: {} combinations",
                total_combinations
            ));
        }

        // Try all combinations (brute force)
        self.brute_force_search(constraints, &var_ranges)
    }

    /// Brute force search through all combinations of variable values
    fn brute_force_search(
        &self,
        constraints: &ConstraintIR,
        var_ranges: &[(String, std::ops::RangeInclusive<i64>)],
    ) -> Result<HashMap<String, crate::surrogate::SolutionValue>> {
        // Recursive helper to generate all combinations
        fn generate_combinations(
            var_ranges: &[(String, std::ops::RangeInclusive<i64>)],
            current: &mut HashMap<String, crate::surrogate::SolutionValue>,
            index: usize,
        ) -> Vec<HashMap<String, crate::surrogate::SolutionValue>> {
            if index >= var_ranges.len() {
                return vec![current.clone()];
            }

            let (var_name, range) = &var_ranges[index];
            let mut results = Vec::new();

            for value in range.clone() {
                current.insert(
                    var_name.clone(),
                    crate::surrogate::SolutionValue::Integer(value),
                );
                results.extend(generate_combinations(var_ranges, current, index + 1));
            }

            results
        }

        let mut initial = HashMap::new();
        let combinations = generate_combinations(var_ranges, &mut initial, 0);

        // Test each combination
        for assignment in combinations {
            if self.verify_solution(constraints, &assignment) {
                return Ok(assignment);
            }
        }

        Err(anyhow::anyhow!("No solution found in search space"))
    }

    /// Extract coefficients from a linear equation like: x + y = 10 or x - y = 2
    /// Returns (coefficients, constant) where coefficients[i] is the coefficient for var_i
    fn extract_linear_equation(
        &self,
        lhs: &crate::ir::Expr,
        rhs: &crate::ir::Expr,
        num_vars: usize,
    ) -> Option<(Vec<(usize, i64)>, i64)> {
        let mut coefficients: Vec<i64> = vec![0; num_vars];
        let mut constant = 0i64;

        // Process LHS (left side gets positive coefficients)
        if !self.extract_coefficients(lhs, &mut coefficients, &mut constant, 1) {
            return None;
        }

        // Process RHS (right side gets negative coefficients, or moves to constant)
        if !self.extract_coefficients(rhs, &mut coefficients, &mut constant, -1) {
            return None;
        }

        // Convert to sparse representation
        let sparse_coeffs: Vec<(usize, i64)> = coefficients
            .iter()
            .enumerate()
            .filter(|(_, coef)| **coef != 0)
            .map(|(idx, &coef)| (idx, coef))
            .collect();

        Some((sparse_coeffs, -constant))
    }

    /// Extract coefficients from an expression (recursive)
    fn extract_coefficients(
        &self,
        expr: &crate::ir::Expr,
        coefficients: &mut [i64],
        constant: &mut i64,
        sign: i64,
    ) -> bool {
        use crate::ir::{BinaryOp, ConstValue, Expr};

        match expr {
            Expr::Const(ConstValue::Integer(val)) => {
                *constant += sign * val;
                true
            }
            Expr::Var(var_id) => {
                if var_id.0 < coefficients.len() {
                    coefficients[var_id.0] += sign;
                    true
                } else {
                    false
                }
            }
            Expr::Binary { op, lhs, rhs } => match op {
                BinaryOp::Add => {
                    self.extract_coefficients(lhs, coefficients, constant, sign)
                        && self.extract_coefficients(rhs, coefficients, constant, sign)
                }
                BinaryOp::Sub => {
                    self.extract_coefficients(lhs, coefficients, constant, sign)
                        && self.extract_coefficients(rhs, coefficients, constant, -sign)
                }
                BinaryOp::Mul => {
                    // Handle coefficient * variable (e.g., 2*x)
                    if let Expr::Const(ConstValue::Integer(coef)) = **lhs {
                        if let Expr::Var(var_id) = **rhs {
                            if var_id.0 < coefficients.len() {
                                coefficients[var_id.0] += sign * coef;
                                return true;
                            }
                        }
                    }
                    // Handle variable * coefficient (e.g., x*2)
                    if let Expr::Var(var_id) = **lhs {
                        if let Expr::Const(ConstValue::Integer(coef)) = **rhs {
                            if var_id.0 < coefficients.len() {
                                coefficients[var_id.0] += sign * coef;
                                return true;
                            }
                        }
                    }
                    false // Non-linear
                }
                _ => false, // Division, modulo, etc. are non-linear
            },
            _ => false,
        }
    }

    /// Solve a system of linear equations using Gaussian elimination
    fn gaussian_elimination(
        &self,
        equations: &[(Vec<(usize, i64)>, i64)],
        num_vars: usize,
    ) -> Result<HashMap<String, crate::surrogate::SolutionValue>> {
        // Convert to dense matrix [A|b]
        let n = equations.len().min(num_vars);
        let mut matrix: Vec<Vec<f64>> = vec![vec![0.0; num_vars + 1]; n];

        for (i, (coeffs, constant)) in equations.iter().take(n).enumerate() {
            for &(var_idx, coef) in coeffs {
                if var_idx < num_vars {
                    matrix[i][var_idx] = coef as f64;
                }
            }
            matrix[i][num_vars] = *constant as f64;
        }

        // Forward elimination
        for col in 0..n.min(num_vars) {
            // Find pivot
            let mut pivot_row = col;
            let mut max_val = matrix[col][col].abs();
            for row in (col + 1)..n {
                if matrix[row][col].abs() > max_val {
                    max_val = matrix[row][col].abs();
                    pivot_row = row;
                }
            }

            if max_val < 1e-10 {
                continue; // Singular or underdetermined
            }

            // Swap rows
            if pivot_row != col {
                matrix.swap(col, pivot_row);
            }

            // Eliminate column
            for row in (col + 1)..n {
                let factor = matrix[row][col] / matrix[col][col];
                for c in col..=num_vars {
                    matrix[row][c] -= factor * matrix[col][c];
                }
            }
        }

        // Back substitution
        let mut solution = vec![0.0; num_vars];
        for row in (0..n).rev() {
            // Find the pivot column for this row
            let mut pivot_col = None;
            for col in 0..num_vars {
                if matrix[row][col].abs() > 1e-10 {
                    pivot_col = Some(col);
                    break;
                }
            }

            if let Some(col) = pivot_col {
                let mut sum = matrix[row][num_vars];
                for c in (col + 1)..num_vars {
                    sum -= matrix[row][c] * solution[c];
                }
                solution[col] = sum / matrix[row][col];
            }
        }

        // Convert to result format
        let mut assignment = HashMap::new();
        for (idx, &val) in solution.iter().enumerate() {
            let var_name = format!("var_{}", idx);
            // Round if close to integer
            let int_val = val.round();
            if (val - int_val).abs() < 1e-9 {
                assignment.insert(
                    var_name,
                    crate::surrogate::SolutionValue::Integer(int_val as i64),
                );
            } else {
                assignment.insert(var_name, crate::surrogate::SolutionValue::Real(val));
            }
        }

        Ok(assignment)
    }

    /// Evaluate an expression to a concrete value
    fn evaluate_expr(
        &self,
        expr: &crate::ir::Expr,
        assignment: &HashMap<String, crate::surrogate::SolutionValue>,
    ) -> Result<crate::surrogate::SolutionValue> {
        use crate::ir::{BinaryOp, ConstValue, Expr};

        match expr {
            Expr::Const(val) => match val {
                ConstValue::Bool(b) => Ok(crate::surrogate::SolutionValue::Bool(*b)),
                ConstValue::Integer(i) => Ok(crate::surrogate::SolutionValue::Integer(*i)),
                ConstValue::Real(r) => Ok(crate::surrogate::SolutionValue::Real(*r)),
                ConstValue::String(s) => Ok(crate::surrogate::SolutionValue::String(s.clone())),
                ConstValue::BitVector { value, width } => {
                    Ok(crate::surrogate::SolutionValue::BitVector {
                        value: *value,
                        width: *width,
                    })
                }
            },
            Expr::Var(var_id) => {
                let var_name = format!("var_{}", var_id.0);
                assignment
                    .get(&var_name)
                    .cloned()
                    .ok_or_else(|| anyhow::anyhow!("Variable {} not found in assignment", var_name))
            }
            Expr::Binary { op, lhs, rhs } => {
                let lhs_val = self.evaluate_expr(lhs, assignment)?;
                let rhs_val = self.evaluate_expr(rhs, assignment)?;

                match (lhs_val, rhs_val) {
                    (
                        crate::surrogate::SolutionValue::Integer(l),
                        crate::surrogate::SolutionValue::Integer(r),
                    ) => {
                        let result = match op {
                            BinaryOp::Add => l + r,
                            BinaryOp::Sub => l - r,
                            BinaryOp::Mul => l * r,
                            BinaryOp::Div => {
                                if r == 0 {
                                    return Err(anyhow::anyhow!("Division by zero"));
                                }
                                l / r
                            }
                            BinaryOp::Mod => {
                                if r == 0 {
                                    return Err(anyhow::anyhow!("Modulo by zero"));
                                }
                                l % r
                            }
                            BinaryOp::Power => {
                                if r < 0 {
                                    return Err(anyhow::anyhow!(
                                        "Negative exponent not supported for integers"
                                    ));
                                }
                                if r > 100 {
                                    return Err(anyhow::anyhow!("Exponent too large"));
                                }
                                l.pow(r as u32)
                            }
                            _ => return Err(anyhow::anyhow!("Unsupported operation {:?}", op)),
                        };
                        Ok(crate::surrogate::SolutionValue::Integer(result))
                    }
                    (
                        crate::surrogate::SolutionValue::Real(l),
                        crate::surrogate::SolutionValue::Real(r),
                    ) => {
                        let result = match op {
                            BinaryOp::Add => l + r,
                            BinaryOp::Sub => l - r,
                            BinaryOp::Mul => l * r,
                            BinaryOp::Div => {
                                if r == 0.0 {
                                    return Err(anyhow::anyhow!("Division by zero"));
                                }
                                l / r
                            }
                            BinaryOp::Power => l.powf(r),
                            _ => return Err(anyhow::anyhow!("Unsupported operation {:?}", op)),
                        };
                        Ok(crate::surrogate::SolutionValue::Real(result))
                    }
                    _ => Err(anyhow::anyhow!("Type mismatch in binary operation")),
                }
            }
            Expr::Unary { op, operand } => {
                use crate::ir::UnaryOp;
                let val = self.evaluate_expr(operand, assignment)?;
                match (op, val) {
                    (UnaryOp::Negate, crate::surrogate::SolutionValue::Integer(i)) => {
                        Ok(crate::surrogate::SolutionValue::Integer(-i))
                    }
                    (UnaryOp::Negate, crate::surrogate::SolutionValue::Real(r)) => {
                        Ok(crate::surrogate::SolutionValue::Real(-r))
                    }
                    (UnaryOp::Sqrt, crate::surrogate::SolutionValue::Integer(i)) => {
                        if i < 0 {
                            return Err(anyhow::anyhow!("Square root of negative number"));
                        }
                        let result = (i as f64).sqrt();
                        // Return integer if it's a perfect square
                        let rounded = result.round();
                        if (result - rounded).abs() < 1e-9 {
                            Ok(crate::surrogate::SolutionValue::Integer(rounded as i64))
                        } else {
                            Ok(crate::surrogate::SolutionValue::Real(result))
                        }
                    }
                    (UnaryOp::Sqrt, crate::surrogate::SolutionValue::Real(r)) => {
                        if r < 0.0 {
                            return Err(anyhow::anyhow!("Square root of negative number"));
                        }
                        Ok(crate::surrogate::SolutionValue::Real(r.sqrt()))
                    }
                    (UnaryOp::Abs, crate::surrogate::SolutionValue::Integer(i)) => {
                        Ok(crate::surrogate::SolutionValue::Integer(i.abs()))
                    }
                    (UnaryOp::Abs, crate::surrogate::SolutionValue::Real(r)) => {
                        Ok(crate::surrogate::SolutionValue::Real(r.abs()))
                    }
                    _ => Err(anyhow::anyhow!("Unsupported unary operation {:?}", op)),
                }
            }
            _ => Err(anyhow::anyhow!("Unsupported expression type")),
        }
    }
}

/// Supported SMT solver types
#[derive(Debug, Clone, Copy)]
enum SolverType {
    Z3,
    CVC5,
}