simplex シンプレックス法

struct Simplex {
    tableau: Vec<Vec<f64>>, // シンプレックス表
    m: usize,               // 制約式の数
    n: usize,               // 変数の数
}

impl Simplex {
    // A: 制約行列, b: 右辺定数ベクトル, c: 目的関数の係数ベクトル
    // offset: 目的関数の定数項
    fn new(a: Vec<Vec<f64>>, b: Vec<f64>, c: Vec<f64>, offset: f64) -> Self {
        let m = a.len();
        let n = c.len();
        // 表のサイズ: (行: m + 1) x (列: n + m + 1)
        // 列構成: [変数x | スラック変数s | 右辺RHS]
        let mut tableau = vec![vec![0.0; n + m + 1]; m + 1];

        for i in 0..m {
            for j in 0..n {
                tableau[i][j] = a[i][j];
            }
            tableau[i][n + i] = 1.0; // スラック変数の係数
            tableau[i][n + m] = b[i]; // RHS
        }

        for j in 0..n {
            tableau[m][j] = -c[j]; // 目的関数 (最大化のため符号反転)
        }
        tableau[m][n + m] = offset; // 定数項をRHSの初期値にセット

        Simplex { tableau, m, n }
    }

    fn solve(&mut self) -> Result<f64, String> {
        loop {
            // 1. ピボット列の選択 (目的関数行で最も負の要素)
            let mut pivot_col = 0;
            let mut min_val = 0.0;
            for j in 0..(self.n + self.m) {
                if self.tableau[self.m][j] < min_val {
                    min_val = self.tableau[self.m][j];
                    pivot_col = j;
                }
            }

            // 最適性の判定: 負の係数がなければ終了
            if min_val >= -1e-9 {
                break;
            }

            // 2. ピボット行の選択 (最小比テスト)
            let mut pivot_row = None;
            let mut min_ratio = f64::INFINITY;

            for i in 0..self.m {
                if self.tableau[i][pivot_col] > 1e-9 {
                    let ratio = self.tableau[i][self.n + self.m] / self.tableau[i][pivot_col];
                    if ratio < min_ratio {
                        min_ratio = ratio;
                        pivot_row = Some(i);
                    }
                }
            }

            let p_row = match pivot_row {
                Some(r) => r,
                None => return Err("問題は非有界（Unbounded）です。".to_string()),
            };

            // 3. ピボット操作 (ガウスの消去法)
            self.pivot(p_row, pivot_col);
        }
        Ok(self.tableau[self.m][self.n + self.m])
    }

    fn pivot(&mut self, row: usize, col: usize) {
        let divisor = self.tableau[row][col];
        for j in 0..=(self.n + self.m) {
            self.tableau[row][j] /= divisor;
        }

        for i in 0..=self.m {
            if i != row {
                let multiplier = self.tableau[i][col];
                for j in 0..=(self.n + self.m) {
                    self.tableau[i][j] -= multiplier * self.tableau[row][j];
                }
            }
        }
    }

    fn print_tableau(&self) {
        println!("--- Tableau ---");
        for row in &self.tableau {
            for val in row {
                print!("{:>8.2} ", val);
            }
            println!();
        }
    }
}

fn main() {
    println!("Solving 0-1 Integer Programming Problem:");
    println!("Maximize: Z = 3x1 + 4x2 + x3 + 2x4");
    println!("Subject to: 2x1 + 3x2 + x3 + 3x4 <= 4");
    println!("x_i in {{0, 1}} for i in {{1, 2, 3, 4}}\n");

    let c = vec![3.0, 4.0, 1.0, 2.0]; // 係数 (x1, x2, x3, x4)
    let a = vec![2.0, 3.0, 1.0, 3.0]; // 制約の係数
    let b = 4.0;                     // 制約の右辺

    // 1. 0-1 整数計画問題として解く (全探索)
    let mut max_z_int = -1.0;
    let mut best_x_int = vec![0; 4];

    for i in 0..(1 << 4) {
        let mut current_z = 0.0;
        let mut current_weight = 0.0;
        let mut x = vec![0; 4];
        for j in 0..4 {
            if (i >> j) & 1 == 1 {
                current_z += c[j];
                current_weight += a[j];
                x[j] = 1;
            }
        }

        if current_weight <= b {
            if current_z > max_z_int {
                max_z_int = current_z;
                best_x_int = x;
            }
        }
    }

    println!("--- 0-1 Integer Solution ---");
    println!("Optimal Z: {:.1}", max_z_int);
    for (i, val) in best_x_int.iter().enumerate() {
        println!("x{} = {}", i + 1, val);
    }
    println!();

    // 2. 線形緩和問題 (Continuous Relaxation) をシンプレックス法で解く
    println!("--- Continuous Relaxation (Simplex Method with xi <= 1) ---");
    let a_matrix = vec![
        vec![2.0, 3.0, 1.0, 3.0], // 主制約
        vec![1.0, 0.0, 0.0, 0.0], // x1 <= 1
        vec![0.0, 1.0, 0.0, 0.0], // x2 <= 1
        vec![0.0, 0.0, 1.0, 0.0], // x3 <= 1
        vec![0.0, 0.0, 0.0, 1.0], // x4 <= 1
    ];
    let b_vec = vec![4.0, 1.0, 1.0, 1.0, 1.0];
    let c_vec = vec![3.0, 4.0, 1.0, 2.0];
    let offset = 0.0;

    let mut simplex = Simplex::new(a_matrix, b_vec, c_vec, offset);
    match simplex.solve() {
        Ok(res) => {
            println!("Optimal Z (Upper Bound): {:.4}", res);
            println!("\nFinal Tableau for Relaxation:");
            simplex.print_tableau();
        }
        Err(e) => println!("Error: {}", e),
    }
}