column_range/
column_range.rs

1//! # Reducing column ranges of a sparse integer matrix
2//!
3//! We would like to reduce the sum of all column ranges of a sparse integer matrix. More
4//! specifically, this concerns a sparse matrix in the sense that each `(row, column)` combination
5//! may or may not contain a value, but if it contains a value, it can be zero.
6//!
7//! This problem formulation introduces a non-negative variable for each row, which will be
8//! subtracted from each value in the corresponding row. Then, a variable for the minimum of each
9//! column, as well as a variable for the maximum of each column.
10//!
11//! Two constraints are introduced for each `(row, column)` combination that contains a value (that
12//! value might be zero):
13//!
14//! * `column minimum + row subtraction amount <= matrix value`
15//! * `column maximum + row subtraction amount >= matrix value`
16//!
17//! None of the "subtraction amount" variables are in the initial basis, but all of the column
18//! minimums and column maximums are, as well as the slack variable for each of the constraints that
19//! don't already contain a pivot from the column minimum or maximum.
20use num_traits::Zero;
21use relp_num::RationalBig;
22use relp_num::RB;
23
24use relp::algorithm::OptimizationResult;
25use relp::algorithm::two_phase::matrix_provider::matrix_data::MatrixData;
26use relp::algorithm::two_phase::matrix_provider::MatrixProvider;
27use relp::algorithm::two_phase::phase_two;
28use relp::algorithm::two_phase::strategy::pivot_rule::FirstProfitable;
29use relp::algorithm::two_phase::tableau::inverse_maintenance::carry::basis_inverse_rows::BasisInverseRows;
30use relp::algorithm::two_phase::tableau::inverse_maintenance::carry::Carry;
31use relp::algorithm::two_phase::tableau::inverse_maintenance::InverseMaintainer;
32use relp::algorithm::two_phase::tableau::kind::non_artificial::NonArtificial;
33use relp::algorithm::two_phase::tableau::Tableau;
34use relp::data::linear_algebra::matrix::{ColumnMajor, MatrixOrder};
35use relp::data::linear_algebra::vector::{DenseVector, Vector};
36use relp::data::linear_program::elements::VariableType;
37use relp::data::linear_program::general_form::Variable;
38
39fn main() {
40    // Needs at least one `Some` in each row.
41    // Values should be large enough, such that the minimal value in the matrix will always be
42    // non-negative.
43    let input_matrix = [
44        [Some(3), Some(3), Some(3)],
45        [None,    Some(3), Some(3)],
46        [Some(1), Some(2), Some(3)],
47    ];
48
49    let m = input_matrix.len();
50    let n = input_matrix[0].len();
51
52    let column_extreme = |f: fn(_) -> Option<i32>| (0..n)
53        .map(|j| f((0..m).filter_map(move |i| input_matrix[i][j])).unwrap())
54        .collect::<Vec<_>>();
55    let column_min = column_extreme(Iterator::min);
56    let column_max = column_extreme(Iterator::max);
57
58    // Variables
59    let subtraction_amount = (0..m).map(|_| Variable {
60            variable_type: VariableType::Continuous,
61            cost: RB!(0),
62            lower_bound: Some(RB!(0)),
63            upper_bound: None,
64            shift: RB!(0),
65            flipped: false,
66        })
67        .collect::<Vec<_>>();
68    let column_minimum = (0..n).map(|_| Variable {
69            variable_type: VariableType::Continuous,
70            // We subtract this variable from a corresponding variable to compute a column range.
71            cost: RB!(-1),
72            lower_bound: Some(RB!(0)),
73            upper_bound: None,
74            shift: RB!(0),
75            flipped: false,
76        })
77        .collect::<Vec<_>>();
78    let column_maximum = (0..n).map(|_| Variable {
79            variable_type: VariableType::Continuous,
80            cost: RB!(1),
81            lower_bound: Some(RB!(0)),
82            upper_bound: None,
83            shift: RB!(0),
84            flipped: false,
85        })
86        .collect::<Vec<_>>();
87    // Note the order of the variables.
88    let variables = [subtraction_amount, column_minimum, column_maximum].concat();
89
90    // Constraint matrix
91    // We collect the constraints row major.
92    let mut row_major_constraints = Vec::new();
93    let mut b = Vec::new();
94    // We'll be choosing a basis column for each row as we construct them.
95    let mut basis_columns = Vec::new();
96
97    // Lowest value in each column
98    let mut nr_upper_bounded_constraints = 0;
99    // Only the constraint containing the variable for the lowest value of a column has an initial
100    // basis column in that same row. For other constraints containing that variable, we use a slack
101    // as a basis column.
102    let mut had_a_min = vec![false; n];
103    // We walk the entire matrix and check whether a `column_maximum_j + x_i >= f(i, j)` constraint
104    // should be added.
105    for j in 0..n {
106        for i in 0..m {
107            if let Some(value) = input_matrix[i][j] {
108                // There is a value there, so we add a constraint.
109                row_major_constraints.push(vec![
110                    (i, RB!(1)),     // The "subtraction amount" variable index
111                    (m + j, RB!(1)), // The "column minimum" variable index (there are `m` of the
112                                     // "subtraction amount" variables)
113                ]);
114                b.push(RB!(value));
115                basis_columns.push(if value == column_min[j] {
116                    // Values that are one of the extreme values in a column "quality" to have their
117                    // constraint contain the "column minimum" variable pivot.
118                    if had_a_min[j] {
119                        // We already have a basis column for the "column minimum" variable, choose
120                        // the constraint slack instead.
121                        m + 2 * n + nr_upper_bounded_constraints
122                    } else {
123                        had_a_min[j] = true;
124                        // We don't yet have a basis column for the "column minimum" variable, so we
125                        // choose its column as the basis column for this constraint.
126                        m + j
127                    }
128                } else {
129                    m + 2 * n + nr_upper_bounded_constraints
130                });
131                nr_upper_bounded_constraints += 1;
132            }
133        }
134    }
135
136    // Largest value in each column
137    let mut nr_lower_bounded_constraints = 0;
138    let mut had_a_max = vec![false; n];
139    for j in 0..n {
140        for i in 0..m {
141            if let Some(value) = input_matrix[i][j] {
142                row_major_constraints.push(vec![(i, RB!(1)), (m + n + j, RB!(1))]);
143                b.push(RB!(value));
144                basis_columns.push(if value == column_max[j] {
145                    if had_a_max[j] {
146                        m + 2 * n + nr_upper_bounded_constraints + nr_lower_bounded_constraints
147                    } else {
148                        had_a_max[j] = true;
149                        m + n + j
150                    }
151                } else {
152                    m + 2 * n + nr_upper_bounded_constraints + nr_lower_bounded_constraints
153                });
154                nr_lower_bounded_constraints += 1;
155            }
156        }
157    }
158
159    // Transposing constraints to column major
160    let mut constraints = vec![vec![]; m + 2 * n];
161    for (row_index, row) in row_major_constraints.into_iter().enumerate() {
162        for (column_index, value) in row {
163            constraints[column_index].push((row_index, value));
164        }
165    }
166    let constraints = ColumnMajor::new(constraints, nr_upper_bounded_constraints + nr_lower_bounded_constraints, m + 2 * n);
167    let b = DenseVector::new(b, nr_upper_bounded_constraints + nr_lower_bounded_constraints);
168
169    // Create the datastructure that will serve as a `MatrixProvider`
170    let matrix = MatrixData::new(
171        &constraints,
172        &b,
173        Vec::with_capacity(0),
174        0, 0, nr_upper_bounded_constraints, nr_lower_bounded_constraints,
175        &variables,
176    );
177
178    // We will maintain a basis inverse explicitly
179    type IM = Carry<RationalBig, BasisInverseRows<RationalBig>>;
180    let inverse_maintainer = IM::from_basis(&basis_columns, &matrix);
181
182    // We create a tableau using the constructed matrix. The basis is initialized using the
183    // specified columns.
184    let mut tableau = Tableau::<_, NonArtificial<_>>::new_with_inverse_maintainer(
185        &matrix, inverse_maintainer, basis_columns.into_iter().collect(),
186    );
187
188    // We apply primal simplex to improve the solution.
189    let result = phase_two::primal::<_, _, FirstProfitable>(&mut tableau);
190    match result {
191        OptimizationResult::FiniteOptimum(vector) => {
192            let solution = matrix.reconstruct_solution(vector);
193            let shifts = (0..m)
194                .map(|i| match solution.get(i) {
195                    None => Zero::zero(),
196                    Some(value) => value.clone(),
197                })
198                .collect::<Vec<_>>();
199            println!("{:?}", shifts);
200        },
201        _ => panic!("We started with a feasible solution, and has at least value 0."),
202    }
203}
column_range/column_range.rs

column_range/
column_range.rs
