Fix bug in scaling of dual slacks and sign of dual variables for >= constraints (#191)

Fixes a bug in how dual slacks (reduced costs) were unscaled. 

Fixes #183 : dual simplex produces the wrong sign for dual variables with >= constraints

## Issue
closes #183

Authors:
  - Chris Maes (https://github.com/chris-maes)

Approvers:
  - Rajesh Gandham (https://github.com/rg20)

URL: #191

Author: chris-maes

cpp/src/dual_simplex/presolve.cpp

@@ -842,6 +842,30 @@ void uncrush_primal_solution(const user_problem_t<i_t, f_t>& user_problem,
   std::copy(solution.begin(), solution.begin() + user_problem.num_cols, user_solution.data());
 }
 
+template <typename i_t, typename f_t>
+void uncrush_dual_solution(const user_problem_t<i_t, f_t>& user_problem,
+                           const lp_problem_t<i_t, f_t>& problem,
+                           const std::vector<f_t>& y,
+                           const std::vector<f_t>& z,
+                           std::vector<f_t>& user_y,
+                           std::vector<f_t>& user_z)
+{
+  // Reduced costs are uncrushed just like the primal solution
+  uncrush_primal_solution(user_problem, problem, z, user_z);
+
+  // Adjust the sign of the dual variables y
+  // We should have A^T y + z = c
+  // In convert_user_problem, we converted >= to <=, so we need to adjust the sign of the dual
+  // variables
+  for (i_t i = 0; i < user_problem.num_rows; i++) {
+    if (user_problem.row_sense[i] == 'G') {
+      user_y[i] = -y[i];
+    } else {
+      user_y[i] = y[i];
+    }
+  }
+}
+
 template <typename i_t, typename f_t>
 void uncrush_solution(const presolve_info_t<i_t, f_t>& presolve_info,
                       const std::vector<f_t>& crushed_x,
@@ -903,6 +927,13 @@ template void uncrush_primal_solution<int, double>(const user_problem_t<int, dou
                                                    const std::vector<double>& solution,
                                                    std::vector<double>& user_solution);
 
+template void uncrush_dual_solution<int, double>(const user_problem_t<int, double>& user_problem,
+                                                 const lp_problem_t<int, double>& problem,
+                                                 const std::vector<double>& y,
+                                                 const std::vector<double>& z,
+                                                 std::vector<double>& user_y,
+                                                 std::vector<double>& user_z);
+
 template void uncrush_solution<int, double>(const presolve_info_t<int, double>& presolve_info,
                                             const std::vector<double>& crushed_x,
                                             const std::vector<double>& crushed_z,

cpp/src/dual_simplex/presolve.hpp

@@ -115,6 +115,14 @@ void uncrush_primal_solution(const user_problem_t<i_t, f_t>& user_problem,
                              const std::vector<f_t>& solution,
                              std::vector<f_t>& user_solution);
 
+template <typename i_t, typename f_t>
+void uncrush_dual_solution(const user_problem_t<i_t, f_t>& user_problem,
+                           const lp_problem_t<i_t, f_t>& problem,
+                           const std::vector<f_t>& y,
+                           const std::vector<f_t>& z,
+                           std::vector<f_t>& user_y,
+                           std::vector<f_t>& user_z);
+
 template <typename i_t, typename f_t>
 void uncrush_solution(const presolve_info_t<i_t, f_t>& presolve_info,
                       const std::vector<f_t>& crushed_x,

cpp/src/dual_simplex/scaling.cpp

@@ -88,7 +88,7 @@ void unscale_solution(const std::vector<f_t>& column_scaling,
   unscaled_z.resize(n);
   for (i_t j = 0; j < n; ++j) {
     unscaled_x[j] = scaled_x[j] / column_scaling[j];
-    unscaled_z[j] = scaled_z[j] / column_scaling[j];
+    unscaled_z[j] = scaled_z[j] * column_scaling[j];
   }
 }
 

cpp/src/dual_simplex/solve.cpp

@@ -252,8 +252,8 @@ lp_status_t solve_linear_program(const user_problem_t<i_t, f_t>& user_problem,
   lp_status_t status = solve_linear_program_advanced(
     original_lp, start_time, settings, lp_solution, vstatus, edge_norms);
   uncrush_primal_solution(user_problem, original_lp, lp_solution.x, solution.x);
-  uncrush_primal_solution(user_problem, original_lp, lp_solution.z, solution.z);
-  solution.y                  = lp_solution.y;
+  uncrush_dual_solution(
+    user_problem, original_lp, lp_solution.y, lp_solution.z, solution.y, solution.z);
   solution.objective          = lp_solution.objective;
   solution.user_objective     = lp_solution.user_objective;
   solution.iterations         = lp_solution.iterations;

cpp/tests/dual_simplex/unit_tests/solve.cpp

@@ -262,4 +262,74 @@ TEST(dual_simplex, empty_columns)
   EXPECT_NEAR(solution.x[8], 0, 1e-6);
 }
 
+TEST(dual_simplex, dual_variable_greater_than)
+{
+  // minimize   3*x0 + 2 * x1
+  // subject to  x0 + x1  >= 1
+  //             x0 + 2x1 >= 3
+  //             x0, x1 >= 0
+
+  cuopt::linear_programming::dual_simplex::user_problem_t<int, double> user_problem;
+  constexpr int m  = 2;
+  constexpr int n  = 2;
+  constexpr int nz = 4;
+
+  user_problem.num_rows = m;
+  user_problem.num_cols = n;
+  user_problem.objective.resize(n);
+  user_problem.objective[0] = 3.0;
+  user_problem.objective[1] = 2.0;
+  user_problem.A.m          = m;
+  user_problem.A.n          = n;
+  user_problem.A.nz_max     = nz;
+  user_problem.A.reallocate(nz);
+  user_problem.A.col_start.resize(n + 1);
+  user_problem.A.col_start[0] = 0;  // x0 start
+  user_problem.A.col_start[1] = 2;
+  user_problem.A.col_start[2] = 4;
+
+  int nnz                 = 0;
+  user_problem.A.i[nnz]   = 0;
+  user_problem.A.x[nnz++] = 1.0;
+  user_problem.A.i[nnz]   = 1;
+  user_problem.A.x[nnz++] = 1.0;
+  user_problem.A.i[nnz]   = 0;
+  user_problem.A.x[nnz++] = 1.0;
+  user_problem.A.i[nnz]   = 1;
+  user_problem.A.x[nnz++] = 2.0;
+  user_problem.A.print_matrix();
+  EXPECT_EQ(nnz, nz);
+
+  user_problem.rhs.resize(m);
+  user_problem.rhs[0] = 1.0;
+  user_problem.rhs[1] = 3.0;
+
+  user_problem.row_sense.resize(m);
+  user_problem.row_sense[0] = 'G';
+  user_problem.row_sense[1] = 'G';
+
+  user_problem.lower.resize(n);
+  user_problem.lower[0] = 0.0;
+  user_problem.lower[1] = 0.0;
+
+  user_problem.upper.resize(n);
+  user_problem.upper[0] = dual_simplex::inf;
+  user_problem.upper[1] = dual_simplex::inf;
+
+  user_problem.num_range_rows = 0;
+  user_problem.problem_name   = "dual_variable_greater_than";
+
+  dual_simplex::simplex_solver_settings_t<int, double> settings;
+  dual_simplex::lp_solution_t<int, double> solution(user_problem.num_rows, user_problem.num_cols);
+  EXPECT_EQ((dual_simplex::solve_linear_program(user_problem, settings, solution)),
+            dual_simplex::lp_status_t::OPTIMAL);
+  EXPECT_NEAR(solution.objective, 3.0, 1e-6);
+  EXPECT_NEAR(solution.x[0], 0.0, 1e-6);
+  EXPECT_NEAR(solution.x[1], 1.5, 1e-6);
+  EXPECT_NEAR(solution.y[0], 0.0, 1e-6);
+  EXPECT_NEAR(solution.y[1], 1.0, 1e-6);
+  EXPECT_NEAR(solution.z[0], 2.0, 1e-6);
+  EXPECT_NEAR(solution.z[1], 0.0, 1e-6);
+}
+
 }  // namespace cuopt::linear_programming::dual_simplex::test