//! `slp` is a Linear Programming Solver.
//!
//! To see the usage docs, visit [here](https://docs.rs/crate/slp/).
//!
//! ## An example
//!
//! ```rust
//! fn main() {
//!     let input = "
//! 3 2 # 3 is number of constraints and 2 is number of variables
//! 2 3 # coefficients of objective function to be maximized: 2x1 + 3x2
//! 2 1 18 # Constraint 1: 2x1 +  x2 <= 18
//! 6 5 60 # Constraint 2: 6x1 + 5x2 <= 60
//! 2 5 40 # Constraint 3: 2x1 + 5x2 <= 40
//! 1 0 -1
//! # x1 >= 0 and x2 >= 0 are always assumed
//!     ";
//!     let mut lp = slp::LP::new_from_buf_reader(&mut input.as_bytes());
//!     let solution = lp.solve();
//!     assert_eq!(solution, slp::Solution::Optimal(28.0, vec![5.0, 6.0]));
//!     match solution {
//!         slp::Solution::Infeasible => println!("INFEASIBLE"),
//!         slp::Solution::Unbounded => println!("UNBOUNDED"),
//!         slp::Solution::Optimal(obj, model) => {
//!             println!("OPTIMAL {}", obj);
//!             print!("SOLUTION");
//!             for v in model {
//!                 print!(" {}", v);
//!             }
//!             println!();
//!         }
//!     }
//! }
//! ```

/// Represents an LP instance
pub struct LP {
    n_constraints: usize,
    n_vars: usize,
    basic_indicies: Vec<usize>,
    tableau: Vec<Vec<f32>>, // Row major format
}

/// Solution to an LP instance as returned by
/// [the solve method](struct.LP.html#method.solve) of
/// [an LP instance](struct.LP.html).
#[derive(Debug, PartialEq)]
pub enum Solution {
    Infeasible,
    Unbounded,
    /// The first value is the optimal value of the objective and
    /// the second value is the assignment.
    Optimal(f32, Vec<f32>),
}

impl LP {
    /// Input is taken from stdin.
    ///
    /// Example input format for
    /// ```txt
    /// max        2x1 + 3x2
    /// subject to 2x1 +  x2 <= 18
    ///            6x1 + 5x2 <= 60
    ///            2x1 + 5x2 <= 40
    ///             x1       >= 0
    ///             x2       >= 0
    /// ```
    /// is
    /// ```txt
    /// 3 2 # 3 is number of constraints and 2 is number of variables
    /// 2 3 # coefficients of objective function to be maximized: 2x1 + 3x2
    /// 2 1 18 # Constraint 1: 2x1 +  x2 <= 18
    /// 6 5 60 # Constraint 2: 6x1 + 5x2 <= 60
    /// 2 5 40 # Constraint 3: 2x1 + 5x2 <= 40
    /// # x1 >= 0 and x2 >= 0 are always assumed
    /// ```
    /// `#` is used for comments.
    pub fn new_from_stdin() -> Self {
        LP::new_from_buf_reader(&mut std::io::stdin().lock())
    }

    /// Input is taken from a file.
    ///
    /// See [new_from_stdin](struct.LP.html#method.new_from_std) function for input format.
    pub fn new_from_file(input_file: &str) -> Self {
        let file = std::fs::File::open(input_file).expect("File not found");
        LP::new_from_buf_reader(&mut std::io::BufReader::new(file))
    }

    /// Input is taken from a reader that implements
    /// [`std::io::BufRead`](https://doc.rust-lang.org/std/io/trait.BufRead.html).
    ///
    /// See [new_from_stdin](struct.LP.html#method.new_from_std) function for input format.
    pub fn new_from_buf_reader<F>(mut reader: &mut F) -> LP
    where
        F: std::io::BufRead,
    {
        let mut tableau: Vec<Vec<f32>> = Vec::new();
        let mut basic_indicies = Vec::new();
        let parsed: Vec<usize> = get_numbers_from_buf_reader(&mut reader, 2);
        let n_constraints = parsed[0];
        let n_vars = parsed[1];
        let mut obj: Vec<f32> = get_numbers_from_buf_reader(&mut reader, n_vars);
        for i in obj.iter_mut() {
            *i = -*i;
        }
        for _ in 0..=n_constraints {
            obj.push(0.0);
        }
        for i in 0..n_constraints {
            let parsed = get_numbers_from_buf_reader(&mut reader, n_vars + 1);
            let mut row = vec![];
            for &coeff in parsed[0..n_vars].iter() {
                row.push(coeff);
            }
            for j in 0..n_constraints {
                row.push(if i == j { 1.0 } else { 0.0 });
            }
            row.push(parsed[n_vars]);
            tableau.push(row);
            basic_indicies.push(n_vars + i);
        }
        tableau.push(obj);
        LP {
            n_constraints,
            n_vars,
            basic_indicies,
            tableau,
        }
    }

    /// Solves the LP.
    ///
    /// Uses naive version of simplex method.
    ///
    /// Returns [a solution](enum.Solution.html).
    pub fn solve(&mut self) -> Solution {
        let b_col = self.n_vars + self.n_constraints;
        let mut phase_one_req = false;
        for i in 0..self.n_constraints {
            if self.tableau[i][b_col] < 0.0 {
                phase_one_req = true;
                break;
            }
        }
        if phase_one_req {
            let mut tableau = vec![];
            for i in 0..self.n_constraints {
                let mut row = vec![];
                for j in 0..b_col {
                    row.push(self.tableau[i][j]);
                }
                row.push(-1.0);
                row.push(self.tableau[i][b_col]);
                tableau.push(row);
            }
            let mut auxi_obj = vec![0.0; b_col + 2];
            auxi_obj[b_col] = 1.0;
            tableau.push(auxi_obj);
            let mut most_neg = 0;
            for i in 1..self.n_constraints {
                if tableau[i][b_col + 1] < tableau[most_neg][b_col + 1] {
                    most_neg = i;
                }
            }
            LP::pivot(&mut tableau, b_col, most_neg);
            let mut auxi_basic_indicies = self.basic_indicies.clone();
            auxi_basic_indicies[most_neg] = b_col;
            let mut auxi_lp = LP {
                n_constraints: self.n_constraints,
                n_vars: self.n_vars + 1,
                basic_indicies: auxi_basic_indicies,
                tableau,
            };
            let sol = auxi_lp.simplex();
            match sol {
                Solution::Infeasible => return Solution::Infeasible,
                Solution::Unbounded => return Solution::Unbounded,
                Solution::Optimal(obj, _) => {
                    if obj != 0.0 {
                        return Solution::Infeasible;
                    }
                    for i in 0..self.n_constraints {
                        for j in 0..b_col - 1 {
                            self.tableau[i][j] = auxi_lp.tableau[i][j];
                        }
                        self.tableau[i][b_col] = auxi_lp.tableau[i][b_col + 1];
                    }
                    for (i, b) in self.basic_indicies.iter_mut().enumerate() {
                        *b = auxi_lp.basic_indicies[i];
                    }
                    for i in 0..self.n_constraints {
                        if self.basic_indicies[i] == b_col + 1 {
                            continue;
                        }
                        let multipler = self.tableau[self.n_constraints][self.basic_indicies[i]];
                        for j in 0..self.tableau[self.n_constraints].len() {
                            self.tableau[self.n_constraints][j] -= multipler * self.tableau[i][j];
                        }
                    }
                }
            }
        }
        self.simplex()
    }

    fn simplex(&mut self) -> Solution {
        let b_col = self.n_vars + self.n_constraints;
        loop {
            let mut entering_var = 0;
            for (i, &v) in self.tableau[self.n_constraints].iter().enumerate() {
                if v < 0.0 && i < b_col {
                    if v < self.tableau[self.n_constraints][entering_var] {
                        entering_var = i;
                    }
                }
            }
            if self.tableau[self.n_constraints][entering_var] >= 0.0 {
                let mut model = Vec::new();
                for i in 0..self.n_vars {
                    let mut found = self.n_constraints;
                    for (j, &v) in self.basic_indicies.iter().enumerate() {
                        if i == v {
                            found = j;
                            break;
                        }
                    }
                    if found == self.n_constraints {
                        model.push(0.0);
                    } else {
                        model.push(self.tableau[found][b_col]);
                    }
                }
                break Solution::Optimal(self.tableau[self.n_constraints][b_col], model);
            }
            let mut leaving_var = 0;
            for i in 0..self.n_constraints {
                if self.tableau[i][entering_var] > 0.0 {
                    if self.tableau[leaving_var][entering_var] <= 0.0
                        || self.tableau[i][b_col] / self.tableau[i][entering_var]
                            < self.tableau[leaving_var][b_col]
                                / self.tableau[leaving_var][entering_var]
                    {
                        leaving_var = i;
                    }
                }
            }
            if self.tableau[leaving_var][entering_var] <= 0.0 {
                break Solution::Unbounded;
            }
            LP::pivot(&mut self.tableau, entering_var, leaving_var);
            self.basic_indicies[leaving_var] = entering_var;
        }
    }

    fn pivot(tableau: &mut Vec<Vec<f32>>, entering_var: usize, leaving_var: usize) {
        let pivot_coeff = tableau[leaving_var][entering_var];
        for v in tableau[leaving_var].iter_mut() {
            *v /= pivot_coeff;
        }
        for k in 0..tableau.len() {
            if k != leaving_var {
                let multiplier = tableau[k][entering_var];
                for i in 0..tableau[k].len() {
                    tableau[k][i] -= multiplier * tableau[leaving_var][i];
                }
            }
        }
    }
}

use std::fmt::Debug;
use std::str::FromStr;

fn get_numbers_from_buf_reader<T, F>(reader: &mut F, n: usize) -> Vec<T>
where
    T: FromStr,
    T::Err: Debug,
    F: std::io::BufRead,
{
    let mut line = String::new();
    let parsed = loop {
        reader.read_line(&mut line).expect("Invalid input");
        let parsed: Vec<T> = parse_line(&line);
        if parsed.len() != 0 {
            break parsed;
        } else {
            line.clear();
        }
    };
    if parsed.len() != n {
        panic!("Invalid input: {:?}", line);
    }
    parsed
}

fn parse_line<T>(line: &str) -> Vec<T>
where
    T: FromStr,
    T::Err: Debug,
{
    let mut line_without_comments = String::new();
    if let Some(index) = line.find('#') {
        line_without_comments = String::from(line.split_at(index).0);
    }
    let parsed: Vec<T> = line_without_comments
        .trim()
        .split_whitespace()
        .map(|s| s.parse::<T>().unwrap())
        .collect();
    parsed
}