sublinear_solver/solver/

neumann.rs

//! Neumann series solver for asymmetric diagonally dominant systems.
//!
//! This module implements a sublinear-time solver based on the Neumann series
//! expansion (I - M)^(-1) = Σ M^k, optimized for diagonally dominant matrices.

use crate::matrix::Matrix;
use crate::types::{Precision, ErrorBounds, ErrorBoundMethod, MemoryInfo};
use crate::error::{SolverError, Result};
use crate::solver::{
    SolverAlgorithm, SolverState, SolverOptions, StepResult, utils
};
use alloc::{vec::Vec, string::String};

/// Neumann series solver implementation.
///
/// Solves systems of the form Ax = b by reformulating as (I - M)x = D^(-1)b
/// where M = I - D^(-1)A and D is the diagonal of A.
///
/// The solution is computed using the Neumann series:
/// x = (I - M)^(-1) D^(-1) b = Σ_{k=0}^∞ M^k D^(-1) b
///
/// For diagonally dominant matrices, ||M|| < 1, ensuring convergence.
#[derive(Debug, Clone)]
pub struct NeumannSolver {
    /// Maximum number of series terms to compute
    max_terms: usize,
    /// Tolerance for series truncation
    series_tolerance: Precision,
    /// Enable adaptive term selection
    adaptive_truncation: bool,
    /// Precompute and cache matrix powers
    cache_powers: bool,
}

impl NeumannSolver {
    /// Create a new Neumann series solver.
    ///
    /// # Arguments
    /// * `max_terms` - Maximum number of series terms (default: 50)
    /// * `series_tolerance` - Tolerance for series truncation (default: 1e-8)
    ///
    /// # Example
    /// ```
    /// use sublinear_solver::NeumannSolver;
    ///
    /// let solver = NeumannSolver::new(16, 1e-8);
    /// ```
    pub fn new(max_terms: usize, series_tolerance: Precision) -> Self {
        Self {
            max_terms,
            series_tolerance,
            adaptive_truncation: true,
            cache_powers: true,
        }
    }

    /// Create a solver with default parameters.
    pub fn default() -> Self {
        Self::new(50, 1e-8)
    }

    /// Create a solver optimized for high precision.
    pub fn high_precision() -> Self {
        Self {
            max_terms: 100,
            series_tolerance: 1e-12,
            adaptive_truncation: true,
            cache_powers: true,
        }
    }

    /// Create a solver optimized for speed.
    pub fn fast() -> Self {
        Self {
            max_terms: 20,
            series_tolerance: 1e-6,
            adaptive_truncation: false,
            cache_powers: false,
        }
    }

    /// Configure adaptive truncation.
    pub fn with_adaptive_truncation(mut self, enable: bool) -> Self {
        self.adaptive_truncation = enable;
        self
    }

    /// Configure matrix power caching.
    pub fn with_power_caching(mut self, enable: bool) -> Self {
        self.cache_powers = enable;
        self
    }
}

/// State for the Neumann series solver.
#[derive(Debug, Clone)]
pub struct NeumannState {
    /// Problem dimension
    dimension: usize,
    /// Current solution estimate
    solution: Vec<Precision>,
    /// Right-hand side vector (D^(-1)b)
    rhs: Vec<Precision>,
    /// Original (un-scaled) right-hand side `b`. Kept so `update_residual`
    /// can compute the residual of `A x = b` correctly — the previous
    /// implementation used `rhs = D⁻¹b` and the resulting "residual norm"
    /// was the residual of `A x = D⁻¹b`, a different equation, so
    /// convergence checks fired too late at larger n (visible as the
    /// bench divergence at n=64).
    original_rhs: Vec<Precision>,
    /// Current residual
    residual: Vec<Precision>,
    /// Residual norm
    residual_norm: Precision,
    /// Diagonal scaling matrix (D^(-1))
    diagonal_inv: Vec<Precision>,
    /// Iteration matrix M = I - D^(-1)A (not cloneable)
    #[allow(dead_code)]
    iteration_matrix: Option<Vec<Vec<Precision>>>,
    /// Cached matrix powers M^k
    matrix_powers: Vec<Vec<Precision>>,
    /// Current series term
    current_term: Vec<Precision>,
    /// Number of series terms computed
    terms_computed: usize,
    /// Number of matrix-vector operations
    matvec_count: usize,
    /// Previous solution for convergence checking
    previous_solution: Option<Vec<Precision>>,
    /// Series convergence indicator
    series_converged: bool,
    /// Error bounds estimation
    error_bounds: Option<ErrorBounds>,
    /// Memory usage tracking
    memory_usage: MemoryInfo,
    /// Target tolerance
    tolerance: Precision,
    /// Maximum allowed terms
    max_terms: usize,
    /// Series truncation tolerance
    series_tolerance: Precision,
}

impl NeumannState {
    /// Create a new Neumann solver state.
    fn new(
        matrix: &dyn Matrix,
        b: &[Precision],
        options: &SolverOptions,
        solver_config: &NeumannSolver,
    ) -> Result<Self> {
        let dimension = matrix.rows();

        if !matrix.is_square() {
            return Err(SolverError::InvalidInput {
                message: "Matrix must be square for Neumann series".to_string(),
                parameter: Some("matrix_dimensions".to_string()),
            });
        }

        if b.len() != dimension {
            return Err(SolverError::DimensionMismatch {
                expected: dimension,
                actual: b.len(),
                operation: "neumann_initialization".to_string(),
            });
        }

        // Check diagonal dominance
        if !matrix.is_diagonally_dominant() {
            return Err(SolverError::MatrixNotDiagonallyDominant {
                row: 0, // Would need to compute actual row
                diagonal: 0.0,
                off_diagonal_sum: 0.0,
            });
        }

        // Extract diagonal and compute D^(-1)
        let mut diagonal_inv = vec![0.0; dimension];
        for i in 0..dimension {
            if let Some(diag_val) = matrix.get(i, i) {
                if diag_val.abs() < 1e-14 {
                    return Err(SolverError::InvalidSparseMatrix {
                        reason: format!("Zero or near-zero diagonal element at position {}", i),
                        position: Some((i, i)),
                    });
                }
                diagonal_inv[i] = 1.0 / diag_val;
            } else {
                return Err(SolverError::InvalidSparseMatrix {
                    reason: format!("Missing diagonal element at position {}", i),
                    position: Some((i, i)),
                });
            }
        }

        // Compute scaled RHS: D^(-1)b
        let rhs: Vec<Precision> = b.iter()
            .zip(&diagonal_inv)
            .map(|(&b_val, &d_inv)| b_val * d_inv)
            .collect();

        // Initialize solution. `compute_next_term` adds term k = M^k · D⁻¹b
        // starting at k=0, so the solution must start at zero — otherwise the
        // k=0 term `D⁻¹b` is double-counted and a 2×2 diagonally dominant
        // toy system that should converge to ~[1, 1] ends up at ~[2, 2].
        // (Caller-supplied initial guesses are still honoured.)
        let solution = if let Some(ref initial) = options.initial_guess {
            if initial.len() != dimension {
                return Err(SolverError::DimensionMismatch {
                    expected: dimension,
                    actual: initial.len(),
                    operation: "initial_guess".to_string(),
                });
            }
            initial.clone()
        } else {
            vec![0.0; dimension]
        };

        let residual = vec![0.0; dimension];
        let current_term = rhs.clone();

        let matrix_powers = if solver_config.cache_powers {
            Vec::with_capacity(solver_config.max_terms)
        } else {
            Vec::new()
        };

        let memory_usage = MemoryInfo {
            current_usage_bytes: dimension * 8 * 5, // Rough estimate
            peak_usage_bytes: dimension * 8 * 5,
            matrix_memory_bytes: 0, // TODO: compute actual matrix memory
            vector_memory_bytes: dimension * 8 * 5,
            workspace_memory_bytes: 0,
            allocation_count: 5,
            deallocation_count: 0,
        };

        Ok(Self {
            dimension,
            solution,
            rhs,
            original_rhs: b.to_vec(),
            residual,
            residual_norm: Precision::INFINITY,
            diagonal_inv,
            iteration_matrix: None,
            matrix_powers,
            current_term,
            terms_computed: 0,
            matvec_count: 0,
            previous_solution: None,
            series_converged: false,
            error_bounds: None,
            memory_usage,
            tolerance: options.tolerance,
            max_terms: solver_config.max_terms,
            series_tolerance: solver_config.series_tolerance,
        })
    }

    /// Compute the next term in the Neumann series.
    fn compute_next_term(&mut self, matrix: &dyn Matrix) -> Result<()> {
        if self.terms_computed >= self.max_terms {
            return Ok(());
        }

        // For k=0: term = D^(-1)b (already in current_term)
        // For k>0: term = M * previous_term
        if self.terms_computed > 0 {
            self.apply_iteration_matrix(matrix)?;
        }

        // Add current term to solution: x += M^k * D^(-1)b
        for (sol, &term) in self.solution.iter_mut().zip(&self.current_term) {
            *sol += term;
        }

        self.terms_computed += 1;

        // Check series convergence
        let term_norm = utils::l2_norm(&self.current_term);
        if term_norm < self.series_tolerance {
            self.series_converged = true;
        }

        Ok(())
    }

    /// Apply the iteration matrix M = I - D^(-1)A to current term.
    fn apply_iteration_matrix(&mut self, matrix: &dyn Matrix) -> Result<()> {
        // Compute M * current_term = current_term - D^(-1) * A * current_term
        let mut temp_vec = vec![0.0; self.dimension];

        // Step 1: temp_vec = A * current_term
        matrix.multiply_vector(&self.current_term, &mut temp_vec)?;
        self.matvec_count += 1;

        // Step 2: temp_vec = D^(-1) * temp_vec
        for (temp, &d_inv) in temp_vec.iter_mut().zip(&self.diagonal_inv) {
            *temp *= d_inv;
        }

        // Step 3: current_term = current_term - temp_vec
        for (curr, &temp) in self.current_term.iter_mut().zip(&temp_vec) {
            *curr -= temp;
        }

        Ok(())
    }

    /// Update the residual and its norm.
    ///
    /// Computes `r = A·x − b` against the ORIGINAL RHS, not the scaled
    /// `D⁻¹b` we keep around for Neumann iteration. The previous
    /// implementation compared against `self.rhs = D⁻¹b`, which means
    /// the "residual" being driven to zero was for the equation
    /// `A·x = D⁻¹b`, not the system we're actually solving. That made
    /// convergence checks late and caused the solver to diverge at
    /// larger n where the scaled-residual heuristic stopped tracking
    /// the true residual (visible in the bench harness at n ≥ 64).
    fn update_residual(&mut self, matrix: &dyn Matrix) -> Result<()> {
        // r = A·x
        matrix.multiply_vector(&self.solution, &mut self.residual)?;
        self.matvec_count += 1;

        // r ← A·x − b   (against the original, un-scaled RHS)
        for (r, &b_val) in self.residual.iter_mut().zip(self.original_rhs.iter()) {
            *r -= b_val;
        }

        self.residual_norm = utils::l2_norm(&self.residual);
        Ok(())
    }

    /// Estimate error bounds based on series truncation.
    fn estimate_error_bounds(&mut self) -> Result<()> {
        if !self.series_converged || self.terms_computed == 0 {
            return Ok(());
        }

        // Estimate ||M||_2 from the computed terms
        let mut matrix_norm_estimate = 0.0;
        if self.terms_computed > 1 {
            let term_ratio = utils::l2_norm(&self.current_term) /
                           utils::l2_norm(&self.rhs);
            matrix_norm_estimate = term_ratio.powf(1.0 / (self.terms_computed - 1) as Precision);
        }

        if matrix_norm_estimate < 1.0 {
            // Error bound for geometric series truncation
            let remaining_sum_bound = matrix_norm_estimate.powi(self.terms_computed as i32) /
                                    (1.0 - matrix_norm_estimate);
            let error_bound = remaining_sum_bound * utils::l2_norm(&self.rhs);

            self.error_bounds = Some(ErrorBounds::upper_bound_only(
                error_bound,
                ErrorBoundMethod::NeumannTruncation,
            ));
        }

        Ok(())
    }
}

impl SolverState for NeumannState {
    fn residual_norm(&self) -> Precision {
        self.residual_norm
    }

    fn matvec_count(&self) -> usize {
        self.matvec_count
    }

    fn error_bounds(&self) -> Option<ErrorBounds> {
        self.error_bounds.clone()
    }

    fn memory_usage(&self) -> MemoryInfo {
        self.memory_usage.clone()
    }

    fn reset(&mut self) {
        self.solution.fill(0.0);
        self.residual.fill(0.0);
        self.residual_norm = Precision::INFINITY;
        self.current_term = self.rhs.clone();
        self.terms_computed = 0;
        self.matvec_count = 0;
        self.previous_solution = None;
        self.series_converged = false;
        self.error_bounds = None;
        self.matrix_powers.clear();
    }
}

impl SolverAlgorithm for NeumannSolver {
    type State = NeumannState;

    fn initialize(
        &self,
        matrix: &dyn Matrix,
        b: &[Precision],
        options: &SolverOptions,
    ) -> Result<Self::State> {
        NeumannState::new(matrix, b, options, self)
    }

    fn step(&self, state: &mut Self::State) -> Result<StepResult> {
        // Save previous solution for convergence checking
        state.previous_solution = Some(state.solution.clone());

        // Compute next term in Neumann series
        // We need access to the original matrix, but we don't have it in state
        // This is a design issue - we need to store matrix reference or pass it
        // For now, return an error indicating we need matrix access
        return Err(SolverError::AlgorithmError {
            algorithm: "neumann".to_string(),
            message: "Matrix reference needed for iteration - design limitation".to_string(),
            context: vec![],
        });

        // TODO: Fix this by either storing matrix ref in state or changing interface
        // state.compute_next_term(matrix)?;
        // state.update_residual(matrix)?;

        // if self.adaptive_truncation {
        //     state.estimate_error_bounds()?;
        // }

        // if state.series_converged || state.terms_computed >= state.max_terms {
        //     Ok(StepResult::Converged)
        // } else {
        //     Ok(StepResult::Continue)
        // }
    }

    fn is_converged(&self, state: &Self::State) -> bool {
        // Check multiple convergence criteria
        let residual_converged = state.residual_norm <= state.tolerance;
        let series_converged = state.series_converged;
        let max_terms_reached = state.terms_computed >= state.max_terms;

        // Converged if residual is small enough or series has converged
        residual_converged || (series_converged && !max_terms_reached)
    }

    fn extract_solution(&self, state: &Self::State) -> Vec<Precision> {
        state.solution.clone()
    }

    fn update_rhs(&self, state: &mut Self::State, delta_b: &[(usize, Precision)]) -> Result<()> {
        // Update the scaled RHS: D^(-1)(b + Δb)
        for &(index, delta) in delta_b {
            if index >= state.dimension {
                return Err(SolverError::IndexOutOfBounds {
                    index,
                    max_index: state.dimension - 1,
                    context: "rhs_update".to_string(),
                });
            }

            // Apply diagonal scaling to the update
            let scaled_delta = delta * state.diagonal_inv[index];
            state.rhs[index] += scaled_delta;

            // For incremental solving, we'd need to adjust the solution
            // This is a simplified implementation
            state.solution[index] += scaled_delta;
        }

        // Reset series computation state
        state.current_term = state.rhs.clone();
        state.terms_computed = 0;
        state.series_converged = false;

        Ok(())
    }

    fn algorithm_name(&self) -> &'static str {
        "neumann"
    }

    /// Custom solve implementation that provides matrix access to steps.
    fn solve(
        &self,
        matrix: &dyn Matrix,
        b: &[Precision],
        options: &SolverOptions,
    ) -> Result<crate::solver::SolverResult> {
        let mut state = self.initialize(matrix, b, options)?;
        let mut iterations = 0;

        #[cfg(feature = "std")]
        let start_time = std::time::Instant::now();

        while !self.is_converged(&state) && iterations < options.max_iterations {
            // Save previous solution
            state.previous_solution = Some(state.solution.clone());

            // Compute next Neumann series term
            state.compute_next_term(matrix)?;

            // Update residual every few iterations (expensive)
            if iterations % 5 == 0 {
                state.update_residual(matrix)?;
            }

            // Estimate error bounds if requested
            if options.compute_error_bounds && self.adaptive_truncation {
                state.estimate_error_bounds()?;
            }

            iterations += 1;

            // Check for numerical issues
            if !state.residual_norm.is_finite() {
                return Err(SolverError::NumericalInstability {
                    reason: "Non-finite residual norm".to_string(),
                    iteration: iterations,
                    residual_norm: state.residual_norm,
                });
            }

            // Early termination if series converged
            if state.series_converged {
                break;
            }
        }

        // Final residual computation
        state.update_residual(matrix)?;

        let converged = self.is_converged(&state);
        let solution = self.extract_solution(&state);
        let residual_norm = state.residual_norm();

        // Check for convergence failure
        if !converged && iterations >= options.max_iterations {
            return Err(SolverError::ConvergenceFailure {
                iterations,
                residual_norm,
                tolerance: options.tolerance,
                algorithm: self.algorithm_name().to_string(),
            });
        }

        let mut result = if converged {
            crate::solver::SolverResult::success(solution, residual_norm, iterations)
        } else {
            crate::solver::SolverResult::failure(solution, residual_norm, iterations)
        };

        // Add optional data if requested
        if options.collect_stats {
            #[cfg(feature = "std")]
            {
                let total_time = start_time.elapsed().as_millis() as f64;
                let mut stats = crate::types::SolverStats::new();
                stats.total_time_ms = total_time;
                stats.matvec_count = state.matvec_count();
                result.stats = Some(stats);
            }
        }

        if options.compute_error_bounds {
            result.error_bounds = state.error_bounds();
        }

        Ok(result)
    }
}

#[cfg(all(test, feature = "std"))]
mod tests {
    use super::*;
    use crate::matrix::SparseMatrix;

    #[test]
    fn test_neumann_solver_creation() {
        let solver = NeumannSolver::new(16, 1e-8);
        assert_eq!(solver.max_terms, 16);
        assert_eq!(solver.series_tolerance, 1e-8);
        assert!(solver.adaptive_truncation);
        assert!(solver.cache_powers);

        let fast_solver = NeumannSolver::fast();
        assert_eq!(fast_solver.max_terms, 20);
        assert!(!fast_solver.cache_powers);
    }

    #[test]
    fn test_neumann_solver_simple_system() {
        // Create a simple 2x2 diagonally dominant system
        let triplets = vec![
            (0, 0, 4.0), (0, 1, 1.0),
            (1, 0, 1.0), (1, 1, 3.0),
        ];
        let matrix = SparseMatrix::from_triplets(triplets, 2, 2).unwrap();
        let b = vec![5.0, 4.0];

        let solver = NeumannSolver::new(20, 1e-8);
        let options = SolverOptions::default();

        let result = solver.solve(&matrix, &b, &options);

        // The system should solve successfully
        match result {
            Ok(solution) => {
                assert!(solution.converged);
                // Expected solution: x = [1, 1] (approximately)
                // 4*1 + 1*1 = 5 ✓
                // 1*1 + 3*1 = 4 ✓
                let x = solution.solution;
                assert!((x[0] - 1.0).abs() < 0.1);
                assert!((x[1] - 1.0).abs() < 0.1);
            },
            Err(e) => {
                // Currently expected due to design limitation
                println!("Expected error: {:?}", e);
            }
        }
    }

    #[test]
    fn test_neumann_not_diagonally_dominant() {
        // Create a non-diagonally dominant matrix
        let triplets = vec![
            (0, 0, 1.0), (0, 1, 3.0),
            (1, 0, 2.0), (1, 1, 1.0),
        ];
        let matrix = SparseMatrix::from_triplets(triplets, 2, 2).unwrap();
        let b = vec![4.0, 3.0];

        let solver = NeumannSolver::new(20, 1e-8);
        let options = SolverOptions::default();

        let result = solver.solve(&matrix, &b, &options);

        // Should fail due to lack of diagonal dominance
        assert!(result.is_err());
        if let Err(SolverError::MatrixNotDiagonallyDominant { .. }) = result {
            // Expected error
        } else {
            panic!("Expected MatrixNotDiagonallyDominant error");
        }
    }

    #[test]
    fn test_neumann_state_initialization() {
        let triplets = vec![(0, 0, 2.0), (1, 1, 3.0)];
        let matrix = SparseMatrix::from_triplets(triplets, 2, 2).unwrap();
        let b = vec![4.0, 6.0];
        let solver = NeumannSolver::default();
        let options = SolverOptions::default();

        let state = solver.initialize(&matrix, &b, &options).unwrap();

        assert_eq!(state.dimension, 2);
        assert_eq!(state.diagonal_inv, vec![0.5, 1.0/3.0]);
        assert_eq!(state.rhs, vec![2.0, 2.0]); // D^(-1)b
        assert_eq!(state.terms_computed, 0);
        assert!(!state.series_converged);
    }
}