better numerics stoncols + parallel cols postsolve for non-vertex sol… (#49)

* better numerics stoncols + parallel cols postsolve for non-vertex solution

* test

* formatter

Author: dance858

.github/workflows/formatting.yml

@@ -28,11 +28,9 @@ jobs:
       # Check formatting
       - name: Check C/C++ formatting
         run: |
-          misformatted=$(find . -name '*.c' -o -name '*.h' -print0 | xargs -0 clang-format -style=file -output-replacements-xml | grep "<replacement " || true)
-          if [ -n "$misformatted" ]; then
-            echo "ERROR: Some files are not properly formatted. Run clang-format -i."
+          if ! git ls-files '*.c' '*.h' | xargs clang-format --dry-run -Werror; then
+            echo "ERROR: Some files are not properly formatted."
+            echo "Fix with: git ls-files '*.c' '*.h' | xargs clang-format -i"
             exit 1
-          else
-            echo "All files are properly formatted."
-            exit 0
-          fi
+          fi
+          echo "All files are properly formatted."

include/core/Numerics.h

@@ -24,6 +24,7 @@
 #define FEAS_TOL 1e-6
 #define BOUND_MARGINAL 0.5 * FEAS_TOL
 #define ZERO_TOL 1e-10
+#define CANCEL_TOL_REL 1e-12
 #define ZERO_TOL_DUAL_POSTSOLVE 1e-6
 #define HUGE_VAL_PS 1e7
 

src/core/Postsolver.c

@@ -270,46 +270,7 @@ static void retrieve_parallel_col(Solution *sol, const int *indices,
     sol->x[k] = xk_val;
 
     // dual postsolve
-    bool is_xj_at_bound =
-        (!HAS_TAG(cTag_j, C_TAG_LB_INF) && IS_EQUAL_FEAS_TOL(xj_val, lb_j)) ||
-        (!HAS_TAG(cTag_j, C_TAG_UB_INF) && IS_EQUAL_FEAS_TOL(xj_val, ub_j));
-    bool is_xk_at_bound =
-        (!HAS_TAG(cTag_k, C_TAG_LB_INF) && IS_EQUAL_FEAS_TOL(xk_val, lb_k)) ||
-        (!HAS_TAG(cTag_k, C_TAG_UB_INF) && IS_EQUAL_FEAS_TOL(xk_val, ub_k));
-
-    if (is_xj_at_bound && is_xk_at_bound)
-    {
-        sol->z[k] = ratio * sol->z[j];
-    }
-    else
-    {
-        sol->z[k] = 0.0;
-    }
-
-#ifndef NDEBUG
-    SET_ZERO_IF_SMALL_DUAL_POSTSOLVE(sol->z[j]);
-
-    // if the multiplier is positive the variable should be at its lower bound
-    if (sol->z[j] > 0)
-    {
-        assert(!HAS_TAG(cTag_j, C_TAG_LB_INF) && IS_EQUAL_FEAS_TOL(sol->x[j], lb_j));
-    }
-    // if the multiplier is negative the variable should be at its upper bound
-    else if (sol->z[j] < 0)
-    {
-        assert(!HAS_TAG(cTag_j, C_TAG_UB_INF) && IS_EQUAL_FEAS_TOL(sol->x[j], ub_j));
-    }
-
-    // similar checks for variable k
-    if (sol->z[k] > 0)
-    {
-        assert(!HAS_TAG(cTag_k, C_TAG_LB_INF) && IS_EQUAL_FEAS_TOL(sol->x[k], lb_k));
-    }
-    else if (sol->z[k] < 0)
-    {
-        assert(!HAS_TAG(cTag_k, C_TAG_UB_INF) && IS_EQUAL_FEAS_TOL(sol->x[k], ub_k));
-    }
-#endif
+    sol->z[k] = ratio * sol->z[j];
 }
 
 void retrieve_deleted_row(Solution *sol, int row, double val)

src/core/Problem.c

@@ -20,6 +20,7 @@
 #include "Constraints.h"
 #include "Matrix.h"
 #include "Memory_wrapper.h"
+#include "Numerics.h"
 #include "State.h"
 #include "Workspace.h"
 #include "utils.h"
@@ -103,7 +104,17 @@ void sub_var_in_obj(Objective *obj, const double *vals, const int *cols, int len
     double ratio = obj->c[k] / aik;
     for (int i = 0; i < len; ++i)
     {
-        obj->c[cols[i]] -= ratio * vals[i];
+        double upd = ratio * vals[i];
+        double new_c = obj->c[cols[i]] - upd;
+
+        /* Exact cancellation leaves a residue whose sign is rounding noise;
+           downstream sign tests (e.g. the unboundedness check in
+           process_colston_ineq) must never see it. */
+        if (ABS(new_c) <= CANCEL_TOL_REL * MAX(ABS(obj->c[cols[i]]), ABS(upd)))
+        {
+            new_c = 0.0;
+        }
+        obj->c[cols[i]] = new_c;
     }
 
     obj->offset += rhs * ratio;

src/core/radix_sort.c

@@ -140,8 +140,7 @@ void radix_sort_rows(int *rows, size_t n, const int *sparsity_IDs,
 
 // --- Single-key radix sort ---
 
-static void insertion_sort_by_key(int *indices, size_t n,
-                                  const int *keys)
+static void insertion_sort_by_key(int *indices, size_t n, const int *keys)
 {
     for (size_t i = 1; i < n; i++)
     {
@@ -157,8 +156,7 @@ static void insertion_sort_by_key(int *indices, size_t n,
     }
 }
 
-void radix_sort_by_key(int *indices, size_t n,
-                       const int *keys, int *aux)
+void radix_sort_by_key(int *indices, size_t n, const int *keys, int *aux)
 {
     if (n < 256)
     {
@@ -177,8 +175,7 @@ void radix_sort_by_key(int *indices, size_t n,
         memset(counts, 0, 256 * sizeof(size_t));
         for (size_t i = 0; i < n; i++)
         {
-            unsigned byte =
-                ((uint32_t) keys[src[i]] >> shift) & 0xFF;
+            unsigned byte = ((uint32_t) keys[src[i]] >> shift) & 0xFF;
             counts[byte]++;
         }
 
@@ -212,8 +209,7 @@ void radix_sort_by_key(int *indices, size_t n,
         // scatter (forward — stable)
         for (size_t i = 0; i < n; i++)
         {
-            unsigned byte =
-                ((uint32_t) keys[src[i]] >> shift) & 0xFF;
+            unsigned byte = ((uint32_t) keys[src[i]] >> shift) & 0xFF;
             dst[counts[byte]++] = src[i];
         }
 

tests/test_postsolve.h

@@ -993,6 +993,54 @@ static char *test_fix_col_inf()
     return 0;
 }
 
+/*  Parallel-column dual postsolve: z_k = ratio * z_j is an exact identity.
+
+    min. [1  2  0.1] x
+    s.t. -20 <= [1  2  1] x <= 20
+         -40 <= [3  6  1] x <= 40
+    x1 in [0, 3], x2 in [0, 5], x3 in [-10, 10] */
+static char *test_parallel_col_dual_identity()
+{
+    double Ax[] = {1, 2, 1, 3, 6, 1};
+    int Ai[] = {0, 1, 2, 0, 1, 2};
+    int Ap[] = {0, 3, 6};
+    int nnz = 6;
+    int n_rows = 2;
+    int n_cols = 3;
+
+    double lhs[] = {-20, -40};
+    double rhs[] = {20, 40};
+    double lbs[] = {0, 0, -10};
+    double ubs[] = {3, 5, 10};
+    double c[] = {1, 2, 0.1};
+
+    Settings *stgs = default_settings();
+    set_settings_false(stgs);
+    stgs->parallel_cols = true;
+    Presolver *presolver =
+        new_presolver(Ax, Ai, Ap, n_rows, n_cols, nnz, lhs, rhs, lbs, ubs, c, stgs);
+
+    run_presolver(presolver);
+
+    // reduced problem: the merged column (box [0, 13]) and x3
+    double x[] = {7, 10};
+    double y[] = {0.3, 0.1};
+    double z[] = {0.4, -0.3};
+    postsolve(presolver, x, y, z);
+
+    double correct_x[] = {0, 3.5, 10};
+    double correct_y[] = {0.3, 0.1};
+    double correct_z[] = {0.4, 0.8, -0.3};
+
+    mu_assert("parallel col dual identity error",
+              is_solution_correct(presolver->sol->x, correct_x, presolver->sol->y,
+                                  correct_y, presolver->sol->z, correct_z, n_rows,
+                                  n_cols, POSTSOLVE_TOL_FEAS));
+    PS_FREE(stgs);
+    free_presolver(presolver);
+    return 0;
+}
+
 static const char *all_tests_postsolve()
 {
     mu_run_test(test_0_postsolve, counter_postsolve);
@@ -1017,6 +1065,7 @@ static const char *all_tests_postsolve()
     mu_run_test(test_pathological_ston_one, counter_postsolve);
     mu_run_test(test_pathological_ston_two, counter_postsolve);
     mu_run_test(test_fix_col_inf, counter_postsolve);
+    mu_run_test(test_parallel_col_dual_identity, counter_postsolve);
     //        all tests above pass
 
     return 0;

tests/test_ston.h

@@ -1200,6 +1200,45 @@ static char *test_16_ston()
     return 0;
 }
 
+/* Test 17: exact cost cancellation must snap to zero, not flip the
+   unboundedness check.
+
+    min. -0.1125*x1 + 0.0375*x2 + x3
+    s.t. [-0.3   0.1   0] x  = 0.6
+         [   0   0.1   1] x >= 0
+    x1 free, x2 >= 0, x3 >= 0 */
+static char *test_17_ston()
+{
+    double Ax[] = {-0.3, 0.1, 0.1, 1};
+    int Ai[] = {0, 1, 1, 2};
+    int Ap[] = {0, 2, 4};
+    int nnz = 4;
+    int n_rows = 2;
+    int n_cols = 3;
+
+    double lhs[] = {0.6, 0};
+    double rhs[] = {0.6, INF};
+    double lbs[] = {-INF, 0, 0};
+    double ubs[] = {INF, INF, INF};
+    double c[] = {-0.1125, 0.0375, 1};
+
+    Settings *stgs = default_settings();
+    Presolver *presolver =
+        new_presolver(Ax, Ai, Ap, n_rows, n_cols, nnz, lhs, rhs, lbs, ubs, c, stgs);
+
+    Problem *prob = presolver->prob;
+
+    PresolveStatus status = remove_ston_cols(prob);
+
+    mu_assert("cancellation residue misread as unbounded", status != UNBNDORINFEAS);
+    mu_assert("partner cost must snap to exact zero", prob->obj->c[1] == 0.0);
+
+    PS_FREE(stgs);
+    free_presolver(presolver);
+
+    return 0;
+}
+
 static const char *all_tests_ston()
 {
     mu_run_test(test_01_ston, counter_ston); // (✓)
@@ -1217,6 +1256,7 @@ static const char *all_tests_ston()
     mu_run_test(test_14_ston, counter_ston); // (✓)
     mu_run_test(test_15_ston, counter_ston); // (✓)
     mu_run_test(test_16_ston, counter_ston); // (✓)
+    mu_run_test(test_17_ston, counter_ston); // (✓)
     return 0;
 }