pub fn primal<'provider, IM, MP, PR>(
tableau: &mut Tableau<IM, NonArtificial<'provider, MP>>,
) -> OptimizationResult<IM::F>Expand description
Reduces the cost of the basic feasible solution to the minimum.
While calling this method, a number of requirements should be satisfied:
- There should be a valid basis (not necessarily optimal <=> dual feasible <=> c >= 0)
- All constraint values need to be positive (primary feasibility)
TODO(CORRECTNESS): Write debug tests for these requirements
ยงReturn value
An OptimizationResult indicating whether or not the problem has a finite optimum. It cannot be
infeasible, as a feasible solution is needed to start using this method.
Examples found in repository?
examples/column_range.rs (line 189)
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}