Postsolve for infeasibility certificates

.github/workflows/coverage.yml

@@ -38,7 +38,7 @@ jobs:
           --root . \
           --xml-pretty \
           --output coverage.xml \
-          --exclude 'test/*' \
+          --exclude 'tests/*' \
           --exclude 'build/*'
 
     # ---- SEND TO COVERALLS ----

include/PSLP/PSLP_API.h

@@ -133,6 +133,23 @@ extern "C"
     void postsolve(Presolver *presolver, const double *x, const double *y,
                    const double *z);
 
+    /* Postsolve a primal infeasibility ray y of the reduced problem.
+       The function writes the corresponding ray for the original problem
+       to y_orig. It does not check whether y is a valid ray.
+
+       y uses the Farkas sign convention, yi >= 0 when row i is active at its
+       rhs, the opposite of postsolve(). It matches Gurobi's FarkasDual. */
+    void postsolve_primal_infeas_ray(Presolver *presolver, const double *y,
+                                     double *y_orig);
+
+    /* Postsolve a dual infeasibility ray x of the reduced problem.
+       The function writes the corresponding ray for the original problem
+       to x_orig. It does not check whether x is a valid ray.
+
+       x is an unbounded direction with c'x < 0; no sign convention involved. */
+    void postsolve_dual_infeas_ray(Presolver *presolver, const double *x,
+                                   double *x_orig);
+
 #ifdef __cplusplus
 }
 #endif

include/core/Postsolver.h

@@ -24,6 +24,11 @@
 
 #include "PSLP_sol.h"
 #include "Tags.h"
+#include "glbopts.h"
+
+#define COL_NOT_RETRIEVED INF
+#define ROW_NOT_RETRIEVED INF
+#define DUMMY_VALUE -382749
 
 struct u16Vec;
 struct dVec;
@@ -82,6 +87,15 @@ void postsolver_update(PostsolveInfo *info, size_t n_cols_reduced,
                        const int *row_map);
 void postsolver_run(const PostsolveInfo *info, Solution *sol, const double *x,
                     const double *y, const double *z);
+void postsolver_run_primal_infeas_ray(const PostsolveInfo *info, Solution *sol,
+                                      const double *y, const double *z);
+void postsolver_run_dual_infeas_ray(const PostsolveInfo *info, Solution *sol,
+                                    const double *x);
+
+void retrieve_deleted_row(Solution *sol, int row, double val);
+void retrieve_added_row(Solution *sol, const int *rows, const double *vals);
+void retrieve_added_rows(Solution *sol, int i, const int *rows, const double *vals,
+                         int len, double aik);
 
 /* Saves the information required to retrieve variable xk that was fixed
    to val. To recover the dual variable we need zk = ck - ak^T y

src/core/Postsolver.c

@@ -31,9 +31,6 @@
 #include <assert.h>
 
 #define INIT_FRAC_POSTSOLVE 0.3
-#define COL_NOT_RETRIEVED INF
-#define ROW_NOT_RETRIEVED INF
-#define DUMMY_VALUE -382749
 
 PostsolveInfo *postsolve_info_new(size_t n_rows, size_t n_cols)
 {

src/core/Presolver.c

@@ -768,6 +768,86 @@ void postsolve(Presolver *presolver, const double *x, const double *y,
     }
 }
 
+static void compute_primal_infeas_ray_z(double *z, const double *Ax, const int *Ai,
+                                        const int *Ap, const double *y,
+                                        size_t n_rows, size_t n_cols)
+{
+    for (size_t j = 0; j < n_cols; ++j)
+    {
+        z[j] = 0.0;
+    }
+
+    for (size_t i = 0; i < n_rows; ++i)
+    {
+        for (int p = Ap[i]; p < Ap[i + 1]; ++p)
+        {
+            z[Ai[p]] -= Ax[p] * y[i];
+        }
+    }
+}
+
+void postsolve_primal_infeas_ray(Presolver *presolver, const double *y,
+                                 double *y_orig)
+{
+    Timer timer;
+    Solution *sol = presolver->sol;
+    PresolveStats *stats = presolver->stats;
+    State *data = presolver->prob->constraints->state;
+    PostsolveInfo *postsolve_info = data->postsolve_info;
+    PresolvedProblem *reduced_prob = presolver->reduced_prob;
+    size_t len_z = MAX((size_t) 1, stats->n_cols_reduced);
+    double *z;
+    assert(reduced_prob != NULL);
+    z = (double *) ps_malloc(len_z, sizeof(double));
+    assert(z != NULL);
+
+    clock_gettime(CLOCK_MONOTONIC, &timer.start);
+    compute_primal_infeas_ray_z(z, reduced_prob->Ax, reduced_prob->Ai,
+                                reduced_prob->Ap, y, stats->n_rows_reduced,
+                                stats->n_cols_reduced);
+    postsolver_update(postsolve_info, stats->n_cols_reduced, stats->n_rows_reduced,
+                      data->work->mappings->cols, data->work->mappings->rows);
+    postsolver_run_primal_infeas_ray(postsolve_info, sol, y, z);
+    for (int i = 0; i < sol->dim_y; ++i)
+    {
+        y_orig[i] = sol->y[i];
+    }
+    PS_FREE(z);
+    clock_gettime(CLOCK_MONOTONIC, &timer.end);
+    stats->time_postsolve = GET_ELAPSED_SECONDS(timer);
+
+    if (presolver->stgs->verbose)
+    {
+        printf("PSLP postsolve time: %.4f seconds\n", stats->time_postsolve);
+    }
+}
+
+void postsolve_dual_infeas_ray(Presolver *presolver, const double *x, double *x_orig)
+{
+    Timer timer;
+    Solution *sol = presolver->sol;
+    PresolveStats *stats = presolver->stats;
+    State *data = presolver->prob->constraints->state;
+    PostsolveInfo *postsolve_info = data->postsolve_info;
+    assert(presolver->reduced_prob != NULL);
+
+    clock_gettime(CLOCK_MONOTONIC, &timer.start);
+    postsolver_update(postsolve_info, stats->n_cols_reduced, stats->n_rows_reduced,
+                      data->work->mappings->cols, data->work->mappings->rows);
+    postsolver_run_dual_infeas_ray(postsolve_info, sol, x);
+    for (int i = 0; i < sol->dim_x; ++i)
+    {
+        x_orig[i] = sol->x[i];
+    }
+    clock_gettime(CLOCK_MONOTONIC, &timer.end);
+    stats->time_postsolve = GET_ELAPSED_SECONDS(timer);
+
+    if (presolver->stgs->verbose)
+    {
+        printf("PSLP postsolve time: %.4f seconds\n", stats->time_postsolve);
+    }
+}
+
 void free_presolver(Presolver *presolver)
 {
     if (presolver == NULL)