ripsolve 0.0.1

Rust Integer Programming Solver
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147

#![feature(nll)]
#[macro_use(s,stack)]
extern crate ndarray;
extern crate rand;

use ndarray::prelude::*;
use ndarray::Zip;
use ndarray::Array;
use rand::{Rng,thread_rng};
use std::time::Instant;

type Matrix = Array2<f32>;
type Vector = Array1<f32>;

#[derive(Debug, PartialEq, Default)]
struct Pivot {
    row: usize,
    col: usize
}

#[allow(non_snake_case)]
fn initial_tableau(c:&Vector, A:&Matrix, b:&Vector) -> Matrix {
    let mut tableau = Array::from_elem((A.rows()+1, A.cols()+1), 0.);
    
    tableau.slice_mut(s![..-1,..-1]).assign(A);
    tableau.slice_mut(s![..-1,  -1]).assign(b);
    tableau.slice_mut(s![-1,  ..-1]).assign(c);

    tableau
}

fn can_improve(tableau:&Matrix) -> bool {
    tableau.slice(s![-1, ..-1]).iter().any(|&x| x < 0.0)
}

fn find_pivot(tableau:&Matrix) -> Result<Pivot, &'static str>{
    // pick minimum index of most negative value in the last row
    let last_row = tableau.slice(s![-1, ..-1]);
    let mut vmax =  0.0;
    let mut col  =  std::usize::MAX;
    for (i, &value) in last_row.iter().enumerate() {
        if value < vmax {
            col = i;
            vmax = value;
        }
    }

    // pick row index minimizing the quotient
    let lc = tableau.cols() - 1;
    let mut qmax =  std::f32::INFINITY;
    let mut row = std::usize::MAX;
    for (i, r) in tableau.slice(s![..-1, ..]).outer_iter().enumerate() {
        if r[col] > 0.0 {
            let q = r[lc] / r[col];
            if qmax > q {
                qmax = q;
                row = i;
            }
        }
    }
    // if row is still std::usize::MAX, then solution is unbounded
    if row == std::usize::MAX {
        Err("Solution is unbounded")
    } else {
        Ok(Pivot{row:row, col:col})
    }
}

fn do_pivot(tableau:&mut Matrix, pivot:&Pivot) {
    let mut pivoting_row = tableau.subview_mut(Axis(0), pivot.row);
    pivoting_row /= pivoting_row[pivot.col];

    let pivoting_row = pivoting_row.to_owned(); 
    
    // subtract the pivot row from all other rows
    for (i, mut r) in tableau.outer_iter_mut().enumerate() {
        if i != pivot.row {
            let multiplier = r[pivot.col];
            Zip::from(&mut r)
                 .and(&pivoting_row)
                 .apply(|w, &x| {*w -= multiplier*x});
        }
    }
}

#[allow(non_snake_case)]
fn simplex(c:&Vector, A:&Matrix, b:&Vector, basis:&mut [usize]) -> Result<f32, &'static str>  {
    let mut tableau = initial_tableau(c, A, b);
    let mut count = 0;
    while can_improve(&tableau) {
        count += 1;
        // println!("Basis = {:?}\n", basis);
        let pivot = find_pivot(&tableau)?;
        // println!("{:?}", pivot);
        do_pivot(&mut tableau, &pivot);
        // println!("{}", tableau);
        basis[pivot.row] = pivot.col;
    }
    println!{"Iterations: {}", count};
    let lc = tableau.cols() - 1;
    let lr = tableau.rows() - 1;
    let last_col = tableau.subview_mut(Axis(1), lc);
    let cost = -last_col[lr];
    Ok(cost)
}

#[allow(non_snake_case)]
fn main() {
    let vars = 100; 
    let constraints = 50;
    let mut rng = thread_rng();
    // let c = array![-10., -12., -12., 0., 0., 0.];
    // let A = array![[1.,    2.,   2., 1., 0., 0.],
    //                [2.,    1.,   2., 0., 1., 0.],
    //                [2.,    2.,   1., 0., 0., 1.]];
    //let b = array![20., 20., 20.];
    

    // let c = array![-16.0, -12.0, -18.0, 0., 0., 0.];
    // let A = array![[2.,    2.,   1., 1., 0., 0.],
    //                [2.,    1.,   2., 0., 1., 0.],
    //                [1.,    2.,   1., 0., 0., 1.]];

    let b = Vector::zeros(constraints)+10.0;
    let mut a_p = Matrix::zeros((constraints, vars));
    for x in a_p.iter_mut() {
        *x = if rng.gen::<f32>() < 0.5 {1.0} else {2.0};
    }
    let A = stack![Axis(1), a_p, Matrix::eye(constraints)];

    let mut c_p = Vector::zeros(vars);
    for x in c_p.iter_mut() {
        *x = rng.gen_range::<f32>(-20., -10.);
    }
    let c = stack![Axis(0), c_p, Vector::zeros(constraints)];


    //println!("{}\n{}", A, c);
    let mut basis:Vec<usize> = (vars..vars+constraints).collect();

    let now = Instant::now();
    match simplex(&c, &A, &b, &mut basis) {
        Ok(cost) => println!("Cost is {}", cost),
        Err(e) => println!("{}", e)
    }
    println!("Elapsed time: {:?}", now.elapsed());
}