use crate::errors::{OptimizationError, Result};
use crate::problem::{ObjectiveType, Problem, Variable};
use crate::solvers::{
    Solution, SolutionStatus, Solver, SolverConfig, SolverResult, SolverStatistics,
};
use approx::AbsDiffEq;
use std::time::{Duration, Instant};

/// simplex tableau for solving LP problems
#[derive(Debug, Clone)]
pub struct SimplexTableau {
    /// tableau matrix (inc RHS column)
    pub matrix: Vec<Vec<f64>>,
    /// num original decision vars
    pub num_vars: usize,
    /// num constraints
    pub num_constraints: usize,
    /// indices of basic vars
    pub basic_variables: Vec<usize>,
    /// simplex phase (1 or 2)
    pub phase: u8,
}

/// simplex method solver
pub struct SimplexSolver {
    config: SolverConfig,
    /// pivot rule to use
    pivot_rule: PivotRule,
    /// use bland's rule to prevent cycling
    use_blands_rule: bool,
}

/// diff pivot rules for simplex method
#[derive(Debug, Clone, Copy)]
pub enum PivotRule {
    /// choose var w/ most negative reduced cost (dantzig's rule)
    MostNegative,
    /// choose first var w/ negative reduced cost
    FirstNegative,
    /// choose based on largest improvement per step
    SteepestEdge,
}

impl SimplexSolver {
    /// creates new simplex solver w/ default config
    pub fn new() -> Self {
        SimplexSolver {
            config: SolverConfig::default(),
            pivot_rule: PivotRule::MostNegative,
            use_blands_rule: false,
        }
    }

    /// Creates a new simplex solver with custom configuration
    pub fn with_config(config: SolverConfig) -> Self {
        SimplexSolver {
            config,
            pivot_rule: PivotRule::MostNegative,
            use_blands_rule: false,
        }
    }

    /// Sets the pivot rule
    pub fn with_pivot_rule(mut self, rule: PivotRule) -> Self {
        self.pivot_rule = rule;
        self
    }

    /// Enables Bland's rule to prevent cycling
    pub fn with_blands_rule(mut self) -> Self {
        self.use_blands_rule = true;
        self
    }

    /// Solves using the two-phase simplex method
    fn solve_two_phase(&self, problem: &Problem) -> SolverResult {
        let start_time = Instant::now();
        let mut iterations = 0;
        let mut basis_factorizations = 0;
        let mut degenerate_pivots = 0;

        // Convert to standard form
        let (standard_problem, initial_basic_vars, _) =
            problem
                .to_standard_form()
                .map_err(|e| OptimizationError::InvalidProblem {
                    reason: e.to_string(),
                })?;

        // Create initial tableau
        let mut tableau = SimplexTableau::from_problem(&standard_problem, initial_basic_vars)?;

        // Phase 1: Find initial basic feasible solution
        if tableau.has_artificial_variables(&standard_problem.variables) {
            tableau.phase = 1;
            let phase1_result = self.run_phase_1(
                &mut tableau,
                &standard_problem,
                &mut iterations,
                &mut degenerate_pivots,
            )?;
            basis_factorizations += 1;

            if !phase1_result {
                return Ok(Solution::infeasible(
                    SolverStatistics::new()
                        .with_iterations(iterations)
                        .with_solve_time(start_time.elapsed())
                        .with_simplex_info(basis_factorizations, degenerate_pivots)
                        .with_info("phase", "1"),
                ));
            }
        }

        // Phase 2: Optimize original objective
        tableau.phase = 2;
        self.setup_phase_2(&mut tableau, &problem.objective)?;
        basis_factorizations += 1;

        // Main simplex iterations
        loop {
            if start_time.elapsed() > self.config.max_time {
                return Ok(Solution::new(SolutionStatus::TimeLimit));
            }

            if iterations >= self.config.max_iterations {
                return Ok(Solution::new(SolutionStatus::IterationLimit));
            }

            if tableau.is_optimal(self.config.optimality_tolerance) {
                break;
            }

            let entering = self.find_entering_variable(&tableau)?;
            let leaving =
                tableau.find_leaving_variable(entering, self.config.feasibility_tolerance)?;

            if leaving.is_none() {
                return Ok(Solution::unbounded(
                    SolverStatistics::new()
                        .with_iterations(iterations)
                        .with_solve_time(start_time.elapsed())
                        .with_simplex_info(basis_factorizations, degenerate_pivots)
                        .with_info("maximize", &problem.objective.objective_type.to_string()),
                ));
            }

            let leaving = leaving.unwrap();

            // Check for degeneracy
            if tableau.matrix[leaving + 1][tableau.matrix[0].len() - 1].abs()
                < self.config.feasibility_tolerance
            {
                degenerate_pivots += 1;
            }

            tableau.pivot(leaving, entering);
            iterations += 1;

            if self.config.verbose && iterations % 100 == 0 {
                println!(
                    "Iteration {}: Objective = {:.6}",
                    iterations,
                    tableau.get_objective_value()
                );
            }
        }

        // Extract solution
        let solution = self.extract_solution(
            &tableau,
            problem,
            iterations,
            start_time.elapsed(),
            basis_factorizations,
            degenerate_pivots,
        )?;
        Ok(solution)
    }

    /// Runs Phase 1 of the simplex method
    fn run_phase_1(
        &self,
        tableau: &mut SimplexTableau,
        standard_problem: &Problem,
        iterations: &mut usize,
        degenerate_pivots: &mut usize,
    ) -> Result<bool> {
        // Set up auxiliary objective to minimize sum of artificial variables
        tableau.setup_phase_1_objective(&standard_problem.variables);

        loop {
            if *iterations >= self.config.max_iterations {
                return Err(OptimizationError::MaxIterationsReached {
                    max_iter: self.config.max_iterations,
                });
            }

            if tableau.is_optimal(self.config.optimality_tolerance) {
                break;
            }

            let entering = self.find_entering_variable(tableau)?;
            let leaving =
                tableau.find_leaving_variable(entering, self.config.feasibility_tolerance)?;

            if leaving.is_none() {
                // Auxiliary problem is unbounded, which means original is infeasible
                return Ok(false);
            }

            let leaving = leaving.unwrap();

            if tableau.matrix[leaving + 1][tableau.matrix[0].len() - 1].abs()
                < self.config.feasibility_tolerance
            {
                *degenerate_pivots += 1;
            }

            tableau.pivot(leaving, entering);
            *iterations += 1;
        }

        // Check if Phase 1 found a feasible solution
        let auxiliary_value = tableau.get_objective_value();
        Ok(auxiliary_value.abs() < self.config.feasibility_tolerance)
    }

    /// Sets up Phase 2 objective
    fn setup_phase_2(
        &self,
        tableau: &mut SimplexTableau,
        original_objective: &crate::problem::Objective,
    ) -> Result<()> {
        let objective = original_objective.to_minimization();

        // Replace objective row with original objective
        for j in 0..tableau.num_vars {
            tableau.matrix[0][j] = if j < objective.coefficients.len() {
                objective.coefficients[j]
            } else {
                0.0
            };
        }
        let num_cols = tableau.matrix[0].len();
        tableau.matrix[0][num_cols - 1] = -objective.constant;

        // Zero out basic variables in objective row
        for i in 0..tableau.num_constraints {
            let basic_var = tableau.basic_variables[i];
            let basic_coeff = tableau.matrix[0][basic_var];

            if basic_coeff.abs() > f64::EPSILON {
                for j in 0..tableau.matrix[0].len() {
                    tableau.matrix[0][j] -= basic_coeff * tableau.matrix[i + 1][j];
                }
            }
        }

        Ok(())
    }

    /// Finds the entering variable based on the pivot rule
    fn find_entering_variable(&self, tableau: &SimplexTableau) -> Result<usize> {
        let epsilon = self.config.optimality_tolerance;

        match self.pivot_rule {
            PivotRule::MostNegative => {
                let mut min_coeff = -epsilon;
                let mut entering = None;

                for j in 0..tableau.num_vars {
                    if tableau.matrix[0][j] < min_coeff {
                        min_coeff = tableau.matrix[0][j];
                        entering = Some(j);
                    }
                }

                entering.ok_or_else(|| OptimizationError::SolverError {
                    reason: "No entering variable found".to_string(),
                })
            }
            PivotRule::FirstNegative => {
                for j in 0..tableau.num_vars {
                    if tableau.matrix[0][j] < -epsilon {
                        return Ok(j);
                    }
                }
                Err(OptimizationError::SolverError {
                    reason: "No entering variable found".to_string(),
                })
            }
            PivotRule::SteepestEdge => {
                // For now, use most negative as approximation
                self.find_entering_variable(&tableau)
            }
        }
    }

    /// Extracts the final solution from the tableau
    fn extract_solution(
        &self,
        tableau: &SimplexTableau,
        original_problem: &Problem,
        iterations: usize,
        solve_time: Duration,
        basis_factorizations: usize,
        degenerate_pivots: usize,
    ) -> Result<Solution> {
        let mut solution_values = vec![0.0; original_problem.num_variables()];

        // Extract basic variable values
        for i in 0..tableau.num_constraints {
            let basic_var = tableau.basic_variables[i];
            if basic_var < original_problem.num_variables() {
                solution_values[basic_var] = tableau.matrix[i + 1][tableau.matrix[0].len() - 1];
            }
        }

        // Compute objective value for original problem
        let objective_value = match original_problem.objective.objective_type {
            ObjectiveType::Minimize => {
                -tableau.get_objective_value() + original_problem.objective.constant
            }
            ObjectiveType::Maximize => {
                tableau.get_objective_value() - original_problem.objective.constant
            }
        };

        // Compute dual variables (shadow prices)
        let dual_variables = self.compute_dual_variables(tableau, original_problem)?;
        let reduced_costs = self.compute_reduced_costs(tableau, original_problem)?;

        let statistics = SolverStatistics::new()
            .with_iterations(iterations)
            .with_solve_time(solve_time)
            .with_simplex_info(basis_factorizations, degenerate_pivots)
            .with_info("method", "two-phase simplex")
            .with_info("pivot_rule", &format!("{:?}", self.pivot_rule));

        Ok(
            Solution::optimal(solution_values, objective_value, statistics)
                .with_dual_info(dual_variables, reduced_costs),
        )
    }

    /// Computes dual variable values
    fn compute_dual_variables(
        &self,
        tableau: &SimplexTableau,
        problem: &Problem,
    ) -> Result<Vec<f64>> {
        let mut duals = vec![0.0; problem.num_constraints()];

        // For basic implementation, use tableau's dual information
        // This is a simplified version - full implementation would be more complex
        for i in 0..tableau.num_constraints.min(duals.len()) {
            // This is a placeholder - proper dual computation requires more work
            duals[i] = 0.0;
        }

        Ok(duals)
    }

    /// Computes reduced costs
    fn compute_reduced_costs(
        &self,
        tableau: &SimplexTableau,
        problem: &Problem,
    ) -> Result<Vec<f64>> {
        let mut costs = vec![0.0; problem.num_variables()];

        for j in 0..problem.num_variables().min(tableau.num_vars) {
            costs[j] = tableau.matrix[0][j];
        }

        Ok(costs)
    }
}

impl SimplexTableau {
    /// Creates a tableau from a problem in standard form
    pub fn from_problem(problem: &Problem, initial_basic_vars: Vec<usize>) -> Result<Self> {
        let num_constraints = problem.num_constraints();
        let total_vars = problem.num_variables();
        let num_cols = total_vars + 1; // +1 for RHS
        let num_rows = num_constraints + 1; // +1 for objective

        let mut matrix = vec![vec![0.0; num_cols]; num_rows];

        // Fill constraint rows
        for (i, constraint) in problem.constraints.iter().enumerate() {
            for (j, &coeff) in constraint.coefficients.iter().enumerate() {
                matrix[i + 1][j] = coeff;
            }
            matrix[i + 1][num_cols - 1] = constraint.rhs;
        }

        // Fill objective row (will be modified in phase setup)
        for (j, &coeff) in problem.objective.coefficients.iter().enumerate() {
            matrix[0][j] = coeff;
        }

        Ok(SimplexTableau {
            matrix,
            num_vars: total_vars,
            num_constraints,
            basic_variables: initial_basic_vars,
            phase: 1,
        })
    }

    /// Checks if the tableau represents an optimal solution
    pub fn is_optimal(&self, tolerance: f64) -> bool {
        self.matrix[0][0..self.num_vars]
            .iter()
            .all(|&coeff| coeff >= -tolerance)
    }

    /// Finds the leaving variable using the minimum ratio test
    pub fn find_leaving_variable(&self, entering: usize, tolerance: f64) -> Result<Option<usize>> {
        let mut min_ratio = f64::INFINITY;
        let mut leaving = None;

        for i in 1..=self.num_constraints {
            let pivot_element = self.matrix[i][entering];
            let rhs = self.matrix[i][self.matrix[0].len() - 1];

            if pivot_element > tolerance {
                let ratio = rhs / pivot_element;
                if ratio < min_ratio
                    || (ratio.abs_diff_eq(&min_ratio, tolerance) && leaving.is_none())
                {
                    min_ratio = ratio;
                    leaving = Some(i - 1);
                }
            }
        }

        Ok(leaving)
    }

    /// Performs the pivot operation
    pub fn pivot(&mut self, leaving_row: usize, entering_col: usize) {
        let pivot_row = leaving_row + 1; // Convert to matrix indices
        let pivot_element = self.matrix[pivot_row][entering_col];

        if pivot_element.abs() < f64::EPSILON {
            return; // Avoid division by near-zero
        }

        // Normalize pivot row
        for j in 0..self.matrix[0].len() {
            self.matrix[pivot_row][j] /= pivot_element;
        }

        // Eliminate column
        for i in 0..self.matrix.len() {
            if i != pivot_row {
                let multiplier = self.matrix[i][entering_col];
                for j in 0..self.matrix[0].len() {
                    self.matrix[i][j] -= multiplier * self.matrix[pivot_row][j];
                }
            }
        }

        // Update basic variables
        self.basic_variables[leaving_row] = entering_col;
    }

    /// Gets the current objective function value
    pub fn get_objective_value(&self) -> f64 {
        -self.matrix[0][self.matrix[0].len() - 1]
    }

    /// Checks if the tableau has artificial variables
    pub fn has_artificial_variables(&self, variables: &[Variable]) -> bool {
        variables.iter().any(|var| var.name.starts_with('a'))
    }

    /// Sets up the Phase 1 auxiliary objective
    pub fn setup_phase_1_objective(&mut self, variables: &[Variable]) {
        // Zero out objective row
        for j in 0..self.matrix[0].len() {
            self.matrix[0][j] = 0.0;
        }

        // Set coefficients for artificial variables to 1
        for (j, var) in variables.iter().enumerate() {
            if var.name.starts_with('a') && j < self.num_vars {
                self.matrix[0][j] = 1.0;
            }
        }

        // Eliminate basic artificial variables from objective
        for i in 0..self.num_constraints {
            let basic_var = self.basic_variables[i];
            if basic_var < variables.len() && variables[basic_var].name.starts_with('a') {
                let coeff = self.matrix[0][basic_var];
                if coeff.abs() > f64::EPSILON {
                    for j in 0..self.matrix[0].len() {
                        self.matrix[0][j] -= coeff * self.matrix[i + 1][j];
                    }
                }
            }
        }
    }
}

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

impl Solver for SimplexSolver {
    fn name(&self) -> &str {
        "Simplex"
    }

    fn solve(&mut self, problem: &Problem) -> SolverResult {
        self.solve_two_phase(problem)
    }

    fn solve_with_config(&mut self, problem: &Problem, config: &SolverConfig) -> SolverResult {
        self.config = config.clone();
        
        let mut working_problem = problem.clone();
        let mut presolve_result = None;
        let mut scaling_result = None;
        
        // apply presolve if enabled
        if config.presolve {
            let presolve = crate::utils::Presolve::new();
            match presolve.apply(&working_problem) {
                Ok(result) => {
                    if config.verbose {
                        println!("Presolve applied: {} vars -> {}, {} constraints -> {}",
                                result.statistics.original_variables,
                                result.statistics.presolved_variables,
                                result.statistics.original_constraints,
                                result.statistics.presolved_constraints);
                    }
                    working_problem = result.problem.clone();
                    presolve_result = Some(result);
                }
                Err(e) => {
                    if config.verbose {
                        eprintln!("Presolve failed: {}, continuing without presolve", e);
                    }
                }
            }
        }
        
        // apply scaling if enabled  
        if config.scaling {
            let scaling = crate::utils::Scaling::new(crate::utils::ScalingMethod::Equilibration);
            match scaling.apply(&working_problem) {
                Ok(result) => {
                    if config.verbose {
                        println!("Scaling applied, condition number: {:.2e} -> {:.2e}",
                                result.statistics.condition_number_before,
                                result.statistics.condition_number_after);
                    }
                    working_problem = result.problem.clone();
                    scaling_result = Some(result);
                }
                Err(e) => {
                    if config.verbose {
                        eprintln!("Scaling failed: {}, continuing without scaling", e);
                    }
                }
            }
        }
        
        // solve the (possibly preprocessed) problem
        let mut solution = self.solve_two_phase(&working_problem)?;
        
        // unscale solution if scaling was applied
        if let Some(scaling_res) = scaling_result {
            let scaling = crate::utils::Scaling::new(crate::utils::ScalingMethod::Equilibration);
            solution.variables = scaling.unscale_solution(&solution.variables, &scaling_res.column_scales);
        }
        
        // map solution back to original variables if presolve was applied
        if let Some(_presolve_res) = presolve_result {
            // TODO: implement solution mapping from presolved back to original variables
            // for now, we'll just use the solution as-is
        }
        
        Ok(solution)
    }

    fn can_solve(&self, problem: &Problem) -> bool {
        let stats = problem.statistics();
        !stats.is_integer // Can only solve linear (not integer) problems
    }

    fn default_config(&self) -> SolverConfig {
        SolverConfig::default()
            .with_option("method", "two-phase")
            .with_option("pivot_rule", "most_negative")
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::problem::*;
    use approx::assert_abs_diff_eq;

    fn create_test_problem() -> Problem {
        ProblemBuilder::new()
            .add_non_negative_variable("x")
            .unwrap()
            .add_non_negative_variable("y")
            .unwrap()
            .add_constraint_by_name(&[("x", 1.0), ("y", 1.0)], ConstraintType::LessEqual, 10.0)
            .unwrap()
            .add_constraint_by_name(
                &[("x", 2.0), ("y", -1.0)],
                ConstraintType::GreaterEqual,
                2.0,
            )
            .unwrap()
            .objective_by_name(&[("x", 3.0), ("y", 2.0)], ObjectiveType::Maximize)
            .unwrap()
            .build()
            .unwrap()
    }

    #[test]
    fn test_simplex_solver_creation() {
        let solver = SimplexSolver::new();
        assert_eq!(solver.name(), "Simplex");
        assert!(solver.can_solve(&create_test_problem()));
    }

    #[test]
    fn test_simplex_solve() {
        let problem = create_test_problem();
        let mut solver = SimplexSolver::new();

        let solution = solver.solve(&problem).unwrap();

        assert!(solution.is_optimal());
        assert_eq!(solution.variables.len(), 2);

        // The optimal solution should be approximately x=4, y=6 with objective value 24
        assert_abs_diff_eq!(solution.variables[0], 4.0, epsilon = 1e-6);
        assert_abs_diff_eq!(solution.variables[1], 6.0, epsilon = 1e-6);
        assert_abs_diff_eq!(solution.objective_value, 24.0, epsilon = 1e-6);
    }

    #[test]
    fn test_simplex_with_config() {
        let problem = create_test_problem();
        let mut solver = SimplexSolver::new();
        let config = SolverConfig::new().with_max_iterations(100).verbose();

        let solution = solver.solve_with_config(&problem, &config).unwrap();
        assert!(solution.is_optimal());
    }

    #[test]
    fn test_infeasible_problem() {
        let problem = ProblemBuilder::new()
            .add_non_negative_variable("x")
            .unwrap()
            .add_constraint_by_name(&[("x", 1.0)], ConstraintType::LessEqual, -1.0)
            .unwrap()
            .objective_by_name(&[("x", 1.0)], ObjectiveType::Maximize)
            .unwrap()
            .build()
            .unwrap();

        let mut solver = SimplexSolver::new();
        let solution = solver.solve(&problem).unwrap();

        assert!(solution.is_infeasible());
    }

    #[test]
    fn test_unbounded_problem() {
        let problem = ProblemBuilder::new()
            .add_non_negative_variable("x")
            .unwrap()
            .add_non_negative_variable("y")
            .unwrap()
            .add_constraint_by_name(&[("x", 1.0), ("y", -1.0)], ConstraintType::LessEqual, 0.0)
            .unwrap()
            .objective_by_name(&[("x", 1.0), ("y", 1.0)], ObjectiveType::Maximize)
            .unwrap()
            .build()
            .unwrap();

        let mut solver = SimplexSolver::new();
        let solution = solver.solve(&problem).unwrap();

        assert!(solution.is_unbounded());
    }
}
