1 /*++
2 Copyright (c) 2017 Microsoft Corporation
3 
4 Module Name:
5 
6     <name>
7 
8 Abstract:
9 
10     <abstract>
11 
12 Author:
13 
14     Lev Nachmanson (levnach)
15 
16 Revision History:
17 
18 
19 --*/
20 #include "math/lp/lp_dual_simplex.h"
21 namespace lp{
22 
decide_on_status_after_stage1()23 template <typename T, typename X> void lp_dual_simplex<T, X>::decide_on_status_after_stage1() {
24     switch (m_core_solver->get_status()) {
25     case lp_status::OPTIMAL:
26         if (this->m_settings.abs_val_is_smaller_than_artificial_tolerance(m_core_solver->get_cost())) {
27             this->m_status = lp_status::FEASIBLE;
28         } else {
29             this->m_status = lp_status::UNBOUNDED;
30         }
31         break;
32     case lp_status::DUAL_UNBOUNDED:
33         lp_unreachable();
34     case lp_status::ITERATIONS_EXHAUSTED:
35         this->m_status = lp_status::ITERATIONS_EXHAUSTED;
36         break;
37     case lp_status::TIME_EXHAUSTED:
38         this->m_status = lp_status::TIME_EXHAUSTED;
39         break;
40     case lp_status::FLOATING_POINT_ERROR:
41         this->m_status = lp_status::FLOATING_POINT_ERROR;
42         break;
43     default:
44         lp_unreachable();
45     }
46 }
47 
fix_logical_for_stage2(unsigned j)48 template <typename T, typename X> void lp_dual_simplex<T, X>::fix_logical_for_stage2(unsigned j) {
49     lp_assert(j >= this->number_of_core_structurals());
50     switch (m_column_types_of_logicals[j - this->number_of_core_structurals()]) {
51     case column_type::lower_bound:
52         m_lower_bounds[j] = numeric_traits<T>::zero();
53         m_column_types_of_core_solver[j] = column_type::lower_bound;
54         m_can_enter_basis[j] = true;
55         break;
56     case column_type::fixed:
57         this->m_upper_bounds[j] = m_lower_bounds[j] = numeric_traits<T>::zero();
58         m_column_types_of_core_solver[j] = column_type::fixed;
59         m_can_enter_basis[j] = false;
60         break;
61     default:
62         lp_unreachable();
63     }
64 }
65 
fix_structural_for_stage2(unsigned j)66 template <typename T, typename X> void lp_dual_simplex<T, X>::fix_structural_for_stage2(unsigned j) {
67     column_info<T> * ci = this->m_map_from_var_index_to_column_info[this->m_core_solver_columns_to_external_columns[j]];
68     switch (ci->get_column_type()) {
69     case column_type::lower_bound:
70         m_lower_bounds[j] = numeric_traits<T>::zero();
71         m_column_types_of_core_solver[j] = column_type::lower_bound;
72         m_can_enter_basis[j] = true;
73         break;
74     case column_type::fixed:
75     case column_type::upper_bound:
76         lp_unreachable();
77     case column_type::boxed:
78         this->m_upper_bounds[j] = ci->get_adjusted_upper_bound() / this->m_column_scale[j];
79         m_lower_bounds[j] = numeric_traits<T>::zero();
80         m_column_types_of_core_solver[j] = column_type::boxed;
81         m_can_enter_basis[j] = true;
82         break;
83     case column_type::free_column:
84         m_can_enter_basis[j] = true;
85         m_column_types_of_core_solver[j] = column_type::free_column;
86         break;
87     default:
88         lp_unreachable();
89     }
90     //    T cost_was = this->m_costs[j];
91     this->set_scaled_cost(j);
92 }
93 
unmark_boxed_and_fixed_columns_and_fix_structural_costs()94 template <typename T, typename X> void lp_dual_simplex<T, X>::unmark_boxed_and_fixed_columns_and_fix_structural_costs() {
95     unsigned j = this->m_A->column_count();
96     while (j-- > this->number_of_core_structurals()) {
97         fix_logical_for_stage2(j);
98     }
99     j = this->number_of_core_structurals();
100     while (j--) {
101         fix_structural_for_stage2(j);
102     }
103 }
104 
restore_right_sides()105 template <typename T, typename X> void lp_dual_simplex<T, X>::restore_right_sides() {
106     unsigned i = this->m_A->row_count();
107     while (i--) {
108         this->m_b[i] = m_b_copy[i];
109     }
110 }
111 
solve_for_stage2()112 template <typename T, typename X> void lp_dual_simplex<T, X>::solve_for_stage2() {
113     m_core_solver->restore_non_basis();
114     m_core_solver->solve_yB(m_core_solver->m_y);
115     m_core_solver->fill_reduced_costs_from_m_y_by_rows();
116     m_core_solver->start_with_initial_basis_and_make_it_dual_feasible();
117     m_core_solver->set_status(lp_status::FEASIBLE);
118     m_core_solver->solve();
119     switch (m_core_solver->get_status()) {
120     case lp_status::OPTIMAL:
121         this->m_status = lp_status::OPTIMAL;
122         break;
123     case lp_status::DUAL_UNBOUNDED:
124         this->m_status = lp_status::INFEASIBLE;
125         break;
126     case lp_status::TIME_EXHAUSTED:
127         this->m_status = lp_status::TIME_EXHAUSTED;
128         break;
129     case lp_status::FLOATING_POINT_ERROR:
130         this->m_status = lp_status::FLOATING_POINT_ERROR;
131         break;
132     default:
133         lp_unreachable();
134     }
135     this->m_second_stage_iterations = m_core_solver->total_iterations();
136     this->m_total_iterations = (this->m_first_stage_iterations + this->m_second_stage_iterations);
137 }
138 
fill_x_with_zeros()139 template <typename T, typename X> void lp_dual_simplex<T, X>::fill_x_with_zeros() {
140     unsigned j = this->m_A->column_count();
141     while (j--) {
142         this->m_x[j] = numeric_traits<T>::zero();
143     }
144 }
145 
stage1()146 template <typename T, typename X> void lp_dual_simplex<T, X>::stage1() {
147     lp_assert(m_core_solver == nullptr);
148     this->m_x.resize(this->m_A->column_count(), numeric_traits<T>::zero());
149     if (this->m_settings.get_message_ostream() != nullptr)
150         this->print_statistics_on_A(*this->m_settings.get_message_ostream());
151     m_core_solver = new lp_dual_core_solver<T, X>(
152                                                   *this->m_A,
153                                                   m_can_enter_basis,
154                                                   this->m_b, // the right side vector
155                                                   this->m_x,
156                                                   this->m_basis,
157                                                   this->m_nbasis,
158                                                   this->m_heading,
159                                                   this->m_costs,
160                                                   this->m_column_types_of_core_solver,
161                                                   this->m_lower_bounds,
162                                                   this->m_upper_bounds,
163                                                   this->m_settings,
164                                                   *this);
165     m_core_solver->fill_reduced_costs_from_m_y_by_rows();
166     m_core_solver->start_with_initial_basis_and_make_it_dual_feasible();
167     if (this->m_settings.abs_val_is_smaller_than_artificial_tolerance(m_core_solver->get_cost())) {
168         // skipping stage 1
169         m_core_solver->set_status(lp_status::OPTIMAL);
170         m_core_solver->set_total_iterations(0);
171     } else {
172         m_core_solver->solve();
173     }
174     decide_on_status_after_stage1();
175     this->m_first_stage_iterations = m_core_solver->total_iterations();
176 }
177 
stage2()178 template <typename T, typename X> void lp_dual_simplex<T, X>::stage2() {
179     unmark_boxed_and_fixed_columns_and_fix_structural_costs();
180     restore_right_sides();
181     solve_for_stage2();
182 }
183 
fill_first_stage_solver_fields()184 template <typename T, typename X> void lp_dual_simplex<T, X>::fill_first_stage_solver_fields() {
185     unsigned slack_var = this->number_of_core_structurals();
186     unsigned artificial = this->number_of_core_structurals() + this->m_slacks;
187 
188     for (unsigned row = 0; row < this->row_count(); row++) {
189         fill_first_stage_solver_fields_for_row_slack_and_artificial(row, slack_var, artificial);
190     }
191     fill_costs_and_bounds_and_column_types_for_the_first_stage_solver();
192 }
193 
get_column_type(unsigned j)194 template <typename T, typename X> column_type lp_dual_simplex<T, X>::get_column_type(unsigned j) {
195     lp_assert(j < this->m_A->column_count());
196     if (j >= this->number_of_core_structurals()) {
197         return m_column_types_of_logicals[j - this->number_of_core_structurals()];
198     }
199     return this->m_map_from_var_index_to_column_info[this->m_core_solver_columns_to_external_columns[j]]->get_column_type();
200 }
201 
fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_structural_column(unsigned j)202 template <typename T, typename X> void lp_dual_simplex<T, X>::fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_structural_column(unsigned j) {
203     // see 4.7 in the dissertation of Achim Koberstein
204     lp_assert(this->m_core_solver_columns_to_external_columns.find(j) !=
205                 this->m_core_solver_columns_to_external_columns.end());
206 
207     T free_bound = T(1e4); // see 4.8
208     unsigned jj = this->m_core_solver_columns_to_external_columns[j];
209     lp_assert(this->m_map_from_var_index_to_column_info.find(jj) != this->m_map_from_var_index_to_column_info.end());
210     column_info<T> * ci = this->m_map_from_var_index_to_column_info[jj];
211     switch (ci->get_column_type()) {
212     case column_type::upper_bound: {
213         std::stringstream s;
214         s << "unexpected bound type " << j << " "
215           << column_type_to_string(get_column_type(j));
216         throw_exception(s.str());
217         break;
218     }
219     case column_type::lower_bound: {
220         m_can_enter_basis[j] = true;
221         this->set_scaled_cost(j);
222         this->m_lower_bounds[j] = numeric_traits<T>::zero();
223         this->m_upper_bounds[j] =numeric_traits<T>::one();
224         break;
225     }
226     case column_type::free_column: {
227         m_can_enter_basis[j] = true;
228         this->set_scaled_cost(j);
229         this->m_upper_bounds[j] = free_bound;
230         this->m_lower_bounds[j] =  -free_bound;
231         break;
232     }
233     case column_type::boxed:
234         m_can_enter_basis[j] = false;
235         this->m_costs[j] = numeric_traits<T>::zero();
236         this->m_upper_bounds[j] = this->m_lower_bounds[j] =  numeric_traits<T>::zero(); // is it needed?
237         break;
238     default:
239         lp_unreachable();
240     }
241     m_column_types_of_core_solver[j] = column_type::boxed;
242 }
243 
fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_logical_column(unsigned j)244 template <typename T, typename X> void lp_dual_simplex<T, X>::fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_logical_column(unsigned j) {
245     this->m_costs[j] = 0;
246     lp_assert(get_column_type(j) != column_type::upper_bound);
247     if ((m_can_enter_basis[j] = (get_column_type(j) == column_type::lower_bound))) {
248         m_column_types_of_core_solver[j] = column_type::boxed;
249         this->m_lower_bounds[j] = numeric_traits<T>::zero();
250         this->m_upper_bounds[j] = numeric_traits<T>::one();
251     } else {
252         m_column_types_of_core_solver[j] = column_type::fixed;
253         this->m_lower_bounds[j] = numeric_traits<T>::zero();
254         this->m_upper_bounds[j] = numeric_traits<T>::zero();
255     }
256 }
257 
fill_costs_and_bounds_and_column_types_for_the_first_stage_solver()258 template <typename T, typename X> void lp_dual_simplex<T, X>::fill_costs_and_bounds_and_column_types_for_the_first_stage_solver() {
259     unsigned j = this->m_A->column_count();
260     while (j-- > this->number_of_core_structurals()) { // go over logicals here
261         fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_logical_column(j);
262     }
263     j = this->number_of_core_structurals();
264     while (j--) {
265         fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_structural_column(j);
266     }
267 }
268 
fill_first_stage_solver_fields_for_row_slack_and_artificial(unsigned row,unsigned & slack_var,unsigned & artificial)269 template <typename T, typename X> void lp_dual_simplex<T, X>::fill_first_stage_solver_fields_for_row_slack_and_artificial(unsigned row,
270                                                                                                                           unsigned & slack_var,
271                                                                                                                           unsigned & artificial) {
272     lp_assert(row < this->row_count());
273     auto & constraint = this->m_constraints[this->m_core_solver_rows_to_external_rows[row]];
274     // we need to bring the program to the form Ax = b
275     T rs = this->m_b[row];
276     switch (constraint.m_relation) {
277     case Equal: // no slack variable here
278         set_type_for_logical(artificial, column_type::fixed);
279         this->m_basis[row] = artificial;
280         this->m_costs[artificial] = numeric_traits<T>::zero();
281         (*this->m_A)(row, artificial) = numeric_traits<T>::one();
282         artificial++;
283         break;
284 
285     case Greater_or_equal:
286         set_type_for_logical(slack_var, column_type::lower_bound);
287         (*this->m_A)(row, slack_var) = - numeric_traits<T>::one();
288         if (rs > 0) {
289             // adding one artificial
290             set_type_for_logical(artificial, column_type::fixed);
291             (*this->m_A)(row, artificial) = numeric_traits<T>::one();
292             this->m_basis[row] = artificial;
293             this->m_costs[artificial] = numeric_traits<T>::zero();
294             artificial++;
295         } else {
296             // we can put a slack_var into the basis, and avoid adding an artificial variable
297             this->m_basis[row] = slack_var;
298             this->m_costs[slack_var] = numeric_traits<T>::zero();
299         }
300         slack_var++;
301         break;
302     case Less_or_equal:
303         // introduce a non-negative slack variable
304         set_type_for_logical(slack_var, column_type::lower_bound);
305         (*this->m_A)(row, slack_var) = numeric_traits<T>::one();
306         if (rs < 0) {
307             // adding one artificial
308             set_type_for_logical(artificial, column_type::fixed);
309             (*this->m_A)(row, artificial) = - numeric_traits<T>::one();
310             this->m_basis[row] = artificial;
311             this->m_costs[artificial] = numeric_traits<T>::zero();
312             artificial++;
313         } else {
314             // we can put slack_var into the basis, and avoid adding an artificial variable
315             this->m_basis[row] = slack_var;
316             this->m_costs[slack_var] = numeric_traits<T>::zero();
317         }
318         slack_var++;
319         break;
320     }
321 }
322 
augment_matrix_A_and_fill_x_and_allocate_some_fields()323 template <typename T, typename X> void lp_dual_simplex<T, X>::augment_matrix_A_and_fill_x_and_allocate_some_fields() {
324     this->count_slacks_and_artificials();
325     this->m_A->add_columns_at_the_end(this->m_slacks + this->m_artificials);
326     unsigned n = this->m_A->column_count();
327     this->m_column_types_of_core_solver.resize(n);
328     m_column_types_of_logicals.resize(this->m_slacks + this->m_artificials);
329     this->m_costs.resize(n);
330     this->m_upper_bounds.resize(n);
331     this->m_lower_bounds.resize(n);
332     m_can_enter_basis.resize(n);
333     this->m_basis.resize(this->m_A->row_count());
334 }
335 
336 
337 
copy_m_b_aside_and_set_it_to_zeros()338 template <typename T, typename X> void lp_dual_simplex<T, X>::copy_m_b_aside_and_set_it_to_zeros() {
339     for (unsigned i = 0; i < this->m_b.size(); i++) {
340         m_b_copy.push_back(this->m_b[i]);
341         this->m_b[i] = numeric_traits<T>::zero(); // preparing for the first stage
342     }
343 }
344 
find_maximal_solution()345 template <typename T, typename X> void lp_dual_simplex<T, X>::find_maximal_solution(){
346     if (this->problem_is_empty()) {
347         this->m_status = lp_status::EMPTY;
348         return;
349     }
350 
351     this->flip_costs(); // do it for now, todo ( remove the flipping)
352 
353     this->cleanup();
354     if (this->m_status == lp_status::INFEASIBLE) {
355         return;
356     }
357     this->fill_matrix_A_and_init_right_side();
358     this->fill_m_b();
359     this->scale();
360     augment_matrix_A_and_fill_x_and_allocate_some_fields();
361     fill_first_stage_solver_fields();
362     copy_m_b_aside_and_set_it_to_zeros();
363     stage1();
364     if (this->m_status == lp_status::FEASIBLE) {
365         stage2();
366     }
367 }
368 
369 
get_current_cost()370 template <typename T, typename X> T lp_dual_simplex<T, X>::get_current_cost() const {
371     T ret = numeric_traits<T>::zero();
372     for (auto it : this->m_map_from_var_index_to_column_info) {
373         ret += this->get_column_cost_value(it.first, it.second);
374     }
375     return -ret; // we flip costs for now
376 }
377 }
378