selen 0.15.5

//! Step 6.3: Subproblem Solving Strategies
//!
//! This module implements specialized solving strategies for different types of subproblems
//! created by the variable partitioning system. Each subproblem type gets an optimized
//! solving approach:
//!
//! - **Float subproblems**: Direct bounds optimization with interval arithmetic
//! - **Integer subproblems**: Enhanced constraint propagation with binary search
//! - **Hybrid coordination**: Managing solution combination and validation
//!
//! The goal is to achieve 10-100x performance improvements over monolithic solving
//! by applying the most appropriate algorithm to each subproblem type.

use crate::model::Model;
use crate::core::solution::Solution;
use crate::optimization::variable_partitioning::{VariablePartition, PartitionResult};
use crate::variables::{Var, VarId};
use crate::variables::domain::float_interval::precision_to_step_size;
use std::time::{Duration, Instant};
use std::collections::HashMap;

/// Specialized solver for float-only subproblems
#[derive(Debug, Clone)]
pub struct FloatSubproblemSolver {
    /// Precision for floating-point operations (decimal places)
    /// This controls the granularity of float solutions and ensures
    /// that subproblem solutions respect the model's precision settings
    precision_digits: i32,
    /// Timeout for solving operations
    timeout: Duration,
}

/// Specialized solver for integer-only subproblems  
#[derive(Debug, Clone)]
pub struct IntegerSubproblemSolver {
    /// Maximum search depth for integer solving
    max_depth: usize,
    /// Timeout for solving operations
    timeout: Duration,
}

/// Coordinator for managing multiple subproblem solutions
#[derive(Debug)]
pub struct SubproblemCoordinator {
    /// Float solver instance
    float_solver: FloatSubproblemSolver,
    /// Integer solver instance
    integer_solver: IntegerSubproblemSolver,
    /// Overall solving timeout
    global_timeout: Duration,
}

/// Result of solving a single subproblem
#[derive(Debug, Clone)]
pub struct SubproblemSolution {
    /// Variable assignments for this subproblem
    pub variable_assignments: HashMap<VarId, SubproblemValue>,
    /// Time taken to solve this subproblem
    pub solve_time: Duration,
    /// Whether the subproblem was solved successfully
    pub is_solved: bool,
    /// Number of variables in this subproblem
    pub variable_count: usize,
}

/// Value type for subproblem solutions
#[derive(Debug, Clone, PartialEq)]
pub enum SubproblemValue {
    /// Float value
    Float(f64),
    /// Integer value  
    Integer(i32),
}

/// Combined result from solving multiple subproblems
#[derive(Debug, Clone)]
pub struct CombinedSolution {
    /// All variable assignments from all subproblems
    pub all_assignments: HashMap<VarId, SubproblemValue>,
    /// Individual subproblem results
    pub subproblem_results: Vec<SubproblemSolution>,
    /// Total solving time across all subproblems
    pub total_time: Duration,
    /// Whether all subproblems were solved successfully
    pub is_complete: bool,
    /// Performance improvement over monolithic solving (estimated)
    pub speedup_factor: f64,
}

/// Errors that can occur during subproblem solving
#[derive(Debug, Clone, PartialEq)]
pub enum SubproblemSolvingError {
    /// Float subproblem solving failed
    FloatSolvingFailed(FloatSolvingError),
    /// Integer subproblem solving failed
    IntegerSolvingFailed(IntegerSolvingError),
    /// Timeout exceeded during solving
    TimeoutExceeded,
    /// Solution combination failed
    CombinationFailed(CombinationError),
    /// No subproblems to solve
    NoSubproblems,
}

/// Specific errors for float subproblem solving
#[derive(Debug, Clone, PartialEq)]
pub enum FloatSolvingError {
    /// No float variables in partition
    EmptyPartition,
    /// Variable is not a float type
    InvalidVariableType(VarId),
    /// Bounds are invalid (e.g., min > max)
    InvalidBounds(VarId),
    /// Numerical computation failed
    ComputationFailed(VarId),
}

/// Specific errors for integer subproblem solving
#[derive(Debug, Clone, PartialEq)]
pub enum IntegerSolvingError {
    /// No integer variables in partition
    EmptyPartition,
    /// Variable is not an integer type
    InvalidVariableType(VarId),
    /// Domain is empty
    EmptyDomain(VarId),
    /// Search depth exceeded
    DepthExceeded,
}

/// Specific errors for solution combination
#[derive(Debug, Clone, PartialEq)]
pub enum CombinationError {
    /// Conflicting variable assignments
    ConflictingAssignments(VarId),
    /// Missing required variable
    MissingVariable(VarId),
    /// Invalid solution structure
    InvalidStructure,
}

impl FloatSubproblemSolver {
    /// Create a new float subproblem solver
    pub fn new(precision_digits: i32) -> Self {
        Self {
            precision_digits,
            timeout: Duration::from_millis(1000), // 1 second default
        }
    }
    
    /// Set solving timeout
    pub fn with_timeout(mut self, timeout: Duration) -> Self {
        self.timeout = timeout;
        self
    }
    
    /// Solve a float-only subproblem using direct bounds optimization
    ///
    /// This method leverages the fact that float subproblems often have
    /// continuous domains and can be solved using interval arithmetic
    /// and bounds propagation rather than expensive search.
    pub fn solve_float_subproblem(
        &self,
        model: &Model,
        partition: &VariablePartition,
    ) -> Result<SubproblemSolution, SubproblemSolvingError> {
        let start_time = Instant::now();
        
        if partition.float_variables.is_empty() {
            return Err(SubproblemSolvingError::FloatSolvingFailed(
                FloatSolvingError::EmptyPartition
            ));
        }
        
        let mut assignments = HashMap::new();
        
        // For Step 6.3, we implement a simplified bounds-based approach
        // In a full implementation, this would use interval arithmetic and constraint propagation
        
        for &var_id in &partition.float_variables {
            if let Some(solution_value) = self.solve_single_float_variable(model, var_id)? {
                assignments.insert(var_id, SubproblemValue::Float(solution_value));
            }
        }
        
        let solve_time = start_time.elapsed();
        
        if solve_time > self.timeout {
            return Err(SubproblemSolvingError::TimeoutExceeded);
        }
        
        let is_solved = !assignments.is_empty();
        let variable_count = partition.float_variables.len();
        
        Ok(SubproblemSolution {
            variable_assignments: assignments,
            solve_time,
            is_solved,
            variable_count,
        })
    }
    
    /// Solve a single float variable using bounds analysis with precision awareness
    fn solve_single_float_variable(
        &self,
        model: &Model,
        var_id: VarId,
    ) -> Result<Option<f64>, SubproblemSolvingError> {
        let vars = model.get_vars();
        
        // Get the variable from the model
        let var = &vars[var_id];
        
        match var {
            Var::VarF(float_interval) => {
                // Use precision-aware calculations based on the solver's precision setting
                let step_size = precision_to_step_size(self.precision_digits);
                
                let min_val = float_interval.min;
                let max_val = float_interval.max;
                
                if min_val.is_finite() && max_val.is_finite() {
                    // Use midpoint as a reasonable solution, but round to solver's precision
                    let midpoint = (min_val + max_val) / 2.0;
                    // Round to the solver's step size, not the interval's
                    let solution = (midpoint / step_size).round() * step_size;
                    Ok(Some(solution))
                } else if min_val.is_finite() {
                    // Only lower bound, move one step from the minimum
                    let candidate = min_val + step_size;
                    let solution = (candidate / step_size).round() * step_size;
                    Ok(Some(solution))
                } else if max_val.is_finite() {
                    // Only upper bound, move one step from the maximum
                    let candidate = max_val - step_size;
                    let solution = (candidate / step_size).round() * step_size;
                    Ok(Some(solution))
                } else {
                    // Unbounded, use 0 as default but round to solver's precision
                    let solution = (0.0 / step_size).round() * step_size;
                    Ok(Some(solution))
                }
            },
            _ => Err(SubproblemSolvingError::FloatSolvingFailed(
                FloatSolvingError::InvalidVariableType(var_id)
            )),
        }
    }
}

impl Default for IntegerSubproblemSolver {
    fn default() -> Self {
        Self::new()
    }
}

impl IntegerSubproblemSolver {
    /// Create a new integer subproblem solver
    pub fn new() -> Self {
        Self {
            max_depth: 1000,
            timeout: Duration::from_millis(5000), // 5 seconds default
        }
    }
    
    /// Set maximum search depth
    pub fn with_max_depth(mut self, max_depth: usize) -> Self {
        self.max_depth = max_depth;
        self
    }
    
    /// Set solving timeout
    pub fn with_timeout(mut self, timeout: Duration) -> Self {
        self.timeout = timeout;
        self
    }
    
    /// Solve an integer-only subproblem using enhanced constraint propagation
    ///
    /// This method uses the existing CSP solving capabilities but optimized
    /// for integer-only domains where we can use more aggressive pruning.
    pub fn solve_integer_subproblem(
        &self,
        model: &Model,
        partition: &VariablePartition,
    ) -> Result<SubproblemSolution, SubproblemSolvingError> {
        let start_time = Instant::now();
        
        if partition.integer_variables.is_empty() {
            return Err(SubproblemSolvingError::IntegerSolvingFailed(
                IntegerSolvingError::EmptyPartition
            ));
        }
        
        let mut assignments = HashMap::new();
        
        // For Step 6.3, implement a simplified integer solving approach
        // In a full implementation, this would use the existing CSP solver with integer-specific optimizations
        
        for &var_id in &partition.integer_variables {
            if let Some(solution_value) = self.solve_single_integer_variable(model, var_id)? {
                assignments.insert(var_id, SubproblemValue::Integer(solution_value));
            }
        }
        
        let solve_time = start_time.elapsed();
        
        if solve_time > self.timeout {
            return Err(SubproblemSolvingError::TimeoutExceeded);
        }
        
        let is_solved = !assignments.is_empty();
        let variable_count = partition.integer_variables.len();
        
        Ok(SubproblemSolution {
            variable_assignments: assignments,
            solve_time,
            is_solved,
            variable_count,
        })
    }
    
    /// Solve a single integer variable using domain analysis
    fn solve_single_integer_variable(
        &self,
        model: &Model,
        var_id: VarId,
    ) -> Result<Option<i32>, SubproblemSolvingError> {
        let vars = model.get_vars();
        
        // Get the variable from the model
        let var = &vars[var_id];
        
        match var {
            Var::VarI(sparse_set) => {
                // For Step 6.3, use the middle value from the domain
                let min_val = sparse_set.min();
                let max_val = sparse_set.max();
                
                // Use midpoint as a reasonable solution
                let solution = (min_val + max_val) / 2;
                Ok(Some(solution))
            },
            _ => Err(SubproblemSolvingError::IntegerSolvingFailed(
                IntegerSolvingError::InvalidVariableType(var_id)
            )),
        }
    }
}

impl SubproblemCoordinator {
    /// Create a new subproblem coordinator
    pub fn new(precision_digits: i32) -> Self {
        Self {
            float_solver: FloatSubproblemSolver::new(precision_digits),
            integer_solver: IntegerSubproblemSolver::new(),
            global_timeout: Duration::from_millis(10000), // 10 seconds default
        }
    }
    
    /// Set global timeout for all solving operations
    pub fn with_global_timeout(mut self, timeout: Duration) -> Self {
        self.global_timeout = timeout;
        self
    }
    
    /// Solve all subproblems from a partition result
    ///
    /// This is the main entry point for Step 6.3. It coordinates solving
    /// of both float and integer subproblems and combines the results.
    pub fn solve_partitioned_problem(
        &self,
        model: &Model,
        partition_result: &PartitionResult,
    ) -> Result<CombinedSolution, SubproblemSolvingError> {
        let overall_start = Instant::now();
        let mut subproblem_results = Vec::new();
        let mut all_assignments = HashMap::new();
        
        // Check if we have any subproblems to solve
        if partition_result.float_partition.is_none() && partition_result.integer_partition.is_none() {
            return Err(SubproblemSolvingError::NoSubproblems);
        }
        
        // Solve float subproblem if it exists
        if let Some(float_partition) = &partition_result.float_partition {
            match self.float_solver.solve_float_subproblem(model, float_partition) {
                Ok(float_solution) => {
                    // Merge float assignments
                    for (var_id, value) in &float_solution.variable_assignments {
                        all_assignments.insert(*var_id, value.clone());
                    }
                    subproblem_results.push(float_solution);
                },
                Err(e) => return Err(e),
            }
        }
        
        // Solve integer subproblem if it exists
        if let Some(integer_partition) = &partition_result.integer_partition {
            match self.integer_solver.solve_integer_subproblem(model, integer_partition) {
                Ok(integer_solution) => {
                    // Merge integer assignments
                    for (var_id, value) in &integer_solution.variable_assignments {
                        all_assignments.insert(*var_id, value.clone());
                    }
                    subproblem_results.push(integer_solution);
                },
                Err(e) => return Err(e),
            }
        }
        
        let total_time = overall_start.elapsed();
        
        // Check global timeout
        if total_time > self.global_timeout {
            return Err(SubproblemSolvingError::TimeoutExceeded);
        }
        
        // Calculate performance improvement estimate
        let speedup_factor = self.estimate_speedup_factor(&subproblem_results, partition_result);
        
        let is_complete = subproblem_results.iter().all(|result| result.is_solved);
        
        Ok(CombinedSolution {
            all_assignments,
            subproblem_results,
            total_time,
            is_complete,
            speedup_factor,
        })
    }
    
    /// Estimate the speedup factor compared to monolithic solving
    fn estimate_speedup_factor(
        &self,
        subproblem_results: &[SubproblemSolution],
        partition_result: &PartitionResult,
    ) -> f64 {
        if subproblem_results.is_empty() {
            return 1.0;
        }
        
        // Calculate actual solving time
        let actual_time: Duration = subproblem_results.iter()
            .map(|result| result.solve_time)
            .sum();
        
        // Estimate monolithic solving time based on problem size
        // This is a heuristic: O(n^2) for mixed problems, where n is variable count
        let total_vars = partition_result.total_variables;
        let estimated_monolithic_time = Duration::from_micros(
            (total_vars * total_vars * 100) as u64 // 100 microseconds per variable^2
        );
        
        // Calculate speedup ratio
        if actual_time.as_nanos() > 0 {
            estimated_monolithic_time.as_nanos() as f64 / actual_time.as_nanos() as f64
        } else {
            10.0 // Default conservative estimate
        }
    }
    
    /// Convert combined solution to a standard Solution object
    pub fn to_solution(&self, combined: &CombinedSolution, _model: &Model) -> Option<Solution> {
        if !combined.is_complete {
            return None;
        }
        
        // For Step 6.3, create a simplified solution
        // In a full implementation, this would properly construct a Solution object
        
        // We'll return None for now since constructing a full Solution requires
        // more integration with the existing solution framework
        None
    }
}

impl std::fmt::Display for SubproblemSolvingError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            SubproblemSolvingError::FloatSolvingFailed(err) => {
                write!(f, "Float subproblem solving failed: {}", err)
            },
            SubproblemSolvingError::IntegerSolvingFailed(err) => {
                write!(f, "Integer subproblem solving failed: {}", err)
            },
            SubproblemSolvingError::TimeoutExceeded => {
                write!(f, "Solving timeout exceeded")
            },
            SubproblemSolvingError::CombinationFailed(err) => {
                write!(f, "Solution combination failed: {}", err)
            },
            SubproblemSolvingError::NoSubproblems => {
                write!(f, "No subproblems to solve")
            },
        }
    }
}

impl std::fmt::Display for FloatSolvingError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            FloatSolvingError::EmptyPartition => {
                write!(f, "No float variables in partition")
            },
            FloatSolvingError::InvalidVariableType(var_id) => {
                write!(f, "Variable {:?} is not a float type", var_id)
            },
            FloatSolvingError::InvalidBounds(var_id) => {
                write!(f, "Variable {:?} has invalid bounds", var_id)
            },
            FloatSolvingError::ComputationFailed(var_id) => {
                write!(f, "Computation failed for variable {:?}", var_id)
            },
        }
    }
}

impl std::fmt::Display for IntegerSolvingError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            IntegerSolvingError::EmptyPartition => {
                write!(f, "No integer variables in partition")
            },
            IntegerSolvingError::InvalidVariableType(var_id) => {
                write!(f, "Variable {:?} is not an integer type", var_id)
            },
            IntegerSolvingError::EmptyDomain(var_id) => {
                write!(f, "Variable {:?} has empty domain", var_id)
            },
            IntegerSolvingError::DepthExceeded => {
                write!(f, "Search depth exceeded")
            },
        }
    }
}

impl std::fmt::Display for CombinationError {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            CombinationError::ConflictingAssignments(var_id) => {
                write!(f, "Conflicting assignments for variable {:?}", var_id)
            },
            CombinationError::MissingVariable(var_id) => {
                write!(f, "Missing required variable {:?}", var_id)
            },
            CombinationError::InvalidStructure => {
                write!(f, "Invalid solution structure")
            },
        }
    }
}

impl std::error::Error for SubproblemSolvingError {}

/// Convenience function to solve a partitioned problem end-to-end
pub fn solve_with_partitioning(
    model: &Model,
    partition_result: &PartitionResult,
) -> Result<CombinedSolution, SubproblemSolvingError> {
    let coordinator = SubproblemCoordinator::new(model.float_precision_digits());
    coordinator.solve_partitioned_problem(model, partition_result)
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::model::Model;
    use crate::optimization::variable_partitioning::{VariablePartition, PartitionResult};

    #[test]
    fn test_float_solver_respects_precision() {
        // Test that FloatSubproblemSolver actually uses the precision setting
        
        // Create model with high precision (8 decimal places)
        let mut model = Model::with_float_precision(8);
        let var_id = model.float(0.0, 1.0).into();
        
        // Create partition with this float variable
        let partition = VariablePartition {
            float_variables: vec![var_id],
            integer_variables: vec![],
            constraint_count: 0,
        };
        
        // Create solver with the same precision
        let solver = FloatSubproblemSolver::new(8);
        
        // Solve the subproblem
        let result = solver.solve_float_subproblem(&model, &partition)
            .expect("Should solve float subproblem");
        
        assert!(result.is_solved);
        assert_eq!(result.variable_assignments.len(), 1);
        
        // Verify the solution value is within the precision bounds
        if let Some(SubproblemValue::Float(value)) = result.variable_assignments.get(&var_id) {
            // Value should be rounded to 8 decimal places
            let step_size = precision_to_step_size(8); // 0.00000001
            let rounded_value = (value / step_size).round() * step_size;
            let diff = (value - rounded_value).abs();
            assert!(diff < 1e-12, "Solution should be rounded to precision: {} vs {}", value, rounded_value);
        } else {
            panic!("Expected float value in solution");
        }
    }
    
    #[test]
    fn test_float_solver_different_precisions() {
        // Test that different precision settings produce differently rounded results
        
        let test_cases = vec![
            (1, 0.1),      // 1 decimal place
            (2, 0.01),     // 2 decimal places  
            (4, 0.0001),   // 4 decimal places
            (6, 0.000001), // 6 decimal places
        ];
        
        for (precision_digits, expected_step) in test_cases {
            let mut model = Model::with_float_precision(precision_digits);
            let var_id = model.float(0.0, 1.0).into();
            
            let partition = VariablePartition {
                float_variables: vec![var_id],
                integer_variables: vec![],
                constraint_count: 0,
            };
            
            let solver = FloatSubproblemSolver::new(precision_digits);
            let result = solver.solve_float_subproblem(&model, &partition)
                .expect("Should solve float subproblem");
            
            // Check that the step size matches expectations
            let actual_step = precision_to_step_size(precision_digits);
            let diff = (actual_step - expected_step).abs();
            assert!(diff < 1e-12, "Step size mismatch for precision {}: {} vs {}", 
                precision_digits, actual_step, expected_step);
            
            // Verify solution is properly rounded
            if let Some(SubproblemValue::Float(value)) = result.variable_assignments.get(&var_id) {
                let remainder = value % actual_step;
                assert!(remainder.abs() < 1e-12 || (actual_step - remainder).abs() < 1e-12,
                    "Value {} should be aligned to step size {} (remainder: {})", 
                    value, actual_step, remainder);
            }
        }
    }
    
    #[test]
    fn test_coordinator_uses_model_precision() {
        // Test that SubproblemCoordinator correctly passes model precision to float solver
        
        let mut model = Model::with_float_precision(3); // 3 decimal places
        let float_var = model.float(0.0, 10.0).into();
        let int_var = model.int(0, 100).into();
        
        // Create partition result with both types
        let partition_result = PartitionResult {
            float_partition: Some(VariablePartition {
                float_variables: vec![float_var],
                integer_variables: vec![],
                constraint_count: 0,
            }),
            integer_partition: Some(VariablePartition {
                float_variables: vec![],
                integer_variables: vec![int_var],
                constraint_count: 0,
            }),
            is_separable: true,
            total_variables: 2,
            total_constraints: 0,
        };
        
        // Coordinator should extract precision from model
        let coordinator = SubproblemCoordinator::new(model.float_precision_digits());
        let result = coordinator.solve_partitioned_problem(&model, &partition_result)
            .expect("Should solve partitioned problem");
        
        assert!(result.is_complete);
        assert_eq!(result.all_assignments.len(), 2);
        
        // Verify float solution respects 3-decimal precision
        if let Some(SubproblemValue::Float(value)) = result.all_assignments.get(&float_var) {
            let step_size = precision_to_step_size(3); // 0.001
            let remainder = value % step_size;
            assert!(remainder.abs() < 1e-12 || (step_size - remainder).abs() < 1e-12,
                "Float value {} should be aligned to 3-decimal precision (step {})", 
                value, step_size);
        }
    }
    
    #[test]
    fn test_precision_mismatch_handling() {
        // Test behavior when solver precision doesn't match model precision
        
        let mut model = Model::with_float_precision(6); // Model has 6 decimal places
        let var_id = model.float(0.0, 1.0).into();
        
        let partition = VariablePartition {
            float_variables: vec![var_id],
            integer_variables: vec![],
            constraint_count: 0,
        };
        
        // Create solver with different precision (2 decimal places)
        let solver = FloatSubproblemSolver::new(2);
        
        let result = solver.solve_float_subproblem(&model, &partition)
            .expect("Should still solve despite precision mismatch");
        
        // Should solve, but use solver's precision, not model's
        assert!(result.is_solved);
        
        if let Some(SubproblemValue::Float(value)) = result.variable_assignments.get(&var_id) {
            // Should be rounded to solver's 2-decimal precision, not model's 6-decimal
            let solver_step = precision_to_step_size(2); // 0.01
            let remainder = value % solver_step;
            assert!(remainder.abs() < 1e-12 || (solver_step - remainder).abs() < 1e-12,
                "Should use solver precision (2 decimals), not model precision (6 decimals)");
        }
    }
    
    #[test]
    fn test_convenience_function_precision_propagation() {
        // Test that solve_with_partitioning correctly propagates model precision
        
        let mut model = Model::with_float_precision(4);
        let var_id = model.float(-5.0, 5.0).into();
        
        let partition_result = PartitionResult {
            float_partition: Some(VariablePartition {
                float_variables: vec![var_id],
                integer_variables: vec![],
                constraint_count: 0,
            }),
            integer_partition: None,
            is_separable: true,
            total_variables: 1,
            total_constraints: 0,
        };
        
        let result = solve_with_partitioning(&model, &partition_result)
            .expect("Should solve with partitioning");
        
        assert!(result.is_complete);
        
        // Verify precision propagation
        if let Some(SubproblemValue::Float(value)) = result.all_assignments.get(&var_id) {
            let step_size = precision_to_step_size(4); // 0.0001
            let remainder = value % step_size;
            assert!(remainder.abs() < 1e-12 || (step_size - remainder).abs() < 1e-12,
                "Convenience function should propagate model precision correctly");
        }
    }
}