Remove cuts according to age

src/branch_and_cut.c

@@ -97,9 +97,7 @@ static int BNC_solve_node(struct BNC *bnc, int depth)
     rval = LP_get_obj_val(lp, &objval);
     abort_if(rval, "LP_get_obj_val failed\n");
 
-    log_debug("    obj value = %.2f\n", objval);
-
-    if (objval > *best_val + LP_EPSILON)
+    if (ceil(objval) > *best_val + LP_EPSILON)
     {
         log_debug("Branch pruned by bound (%.2lf > %.2lf).\n", objval,
                 *best_val);
@@ -120,20 +118,25 @@ static int BNC_solve_node(struct BNC *bnc, int depth)
         log_debug("Adding problem cutting planes...\n");
         rval = bnc->problem_add_cutting_planes(lp, bnc->problem_data);
         abort_if(rval, "problem_add_cutting_planes failed");
-    }
 
-    rval = LP_optimize(lp, &is_infeasible);
-    abort_if(rval, "LP_optimize failed\n");
+        rval = LP_get_obj_val(lp, &objval);
+        abort_if(rval, "LP_get_obj_val failed");
 
-    rval = LP_get_obj_val(lp, &objval);
-    abort_if(rval, "LP_get_obj_val failed");
+        if (ceil(objval) > *best_val + LP_EPSILON)
+        {
+            log_debug("Branch pruned by bound (%.2lf > %.2lf).\n", objval,
+                    *best_val);
+            rval = 0;
+            goto CLEANUP;
+        }
 
-    rval = LP_get_x(lp, x);
-    abort_if(rval, "LP_get_x failed");
+        rval = LP_get_x(lp, x);
+        abort_if(rval, "LP_get_x failed");
+    }
 
     if (BNC_is_integral(x, num_cols))
     {
-        log_debug("    solution is integral\n");
+        log_debug("Solution is integral\n");
 
         if (objval + LP_EPSILON < *best_val)
         {
@@ -148,9 +151,10 @@ static int BNC_solve_node(struct BNC *bnc, int depth)
             if (bnc->problem_solution_found)
                 bnc->problem_solution_found(bnc->problem_data, bnc->best_x);
         }
-    } else
+    }
+    else
     {
-        log_debug("    solution is fractional\n");
+        log_debug("Solution is fractional\n");
         rval = BNC_branch_node(bnc, x, depth);
         abort_if(rval, "BNC_branch_node failed");
     }

src/gtsp-subtour.c

@@ -172,9 +172,7 @@ int static add_subtour_cut(
 }
 
 int find_exact_subtour_cuts(
-        struct LP *lp,
-        struct GTSP *data,
-        double min_cut_violation)
+        struct LP *lp, struct GTSP *data, double min_cut_violation)
 {
     int rval = 0;
 
@@ -215,8 +213,12 @@ int find_exact_subtour_cuts(
     abort_if(rval, "find_exact_subtour_cuts_cluster_to_cluster failed");
 
     added_cuts_count = lp->cut_pool_size - original_cut_pool_size;
-    log_debug("Added %d cluster-to-cluster subtour cuts\n", added_cuts_count);
-    if (added_cuts_count > 0) goto CLEANUP;
+    if (added_cuts_count > 0)
+    {
+        log_debug("Added %d cluster-to-cluster subtour cuts\n",
+                added_cuts_count);
+        goto CLEANUP;
+    }
 
     // Constraints (2.2)
     original_cut_pool_size = lp->cut_pool_size;
@@ -225,8 +227,11 @@ int find_exact_subtour_cuts(
     abort_if(rval, "find_exact_subtour_cuts_node_to_cluster failed");
 
     added_cuts_count = lp->cut_pool_size - original_cut_pool_size;
-    log_debug("Added %d node-to-cluster subtour cuts\n", added_cuts_count);
-    if (added_cuts_count > 0) goto CLEANUP;
+    if (added_cuts_count > 0)
+    {
+        log_debug("Added %d node-to-cluster subtour cuts\n", added_cuts_count);
+        goto CLEANUP;
+    }
 
     // Constraints (2.3)
     original_cut_pool_size = lp->cut_pool_size;
@@ -235,7 +240,11 @@ int find_exact_subtour_cuts(
     abort_if(rval, "find_exact_subtour_cuts_node_to_node failed");
 
     added_cuts_count = lp->cut_pool_size - original_cut_pool_size;
-    log_debug("Added %d node-to-node subtour cuts\n", added_cuts_count);
+    if (added_cuts_count > 0)
+    {
+        log_debug("Added %d node-to-node subtour cuts\n", added_cuts_count);
+        goto CLEANUP;
+    }
 
     CLEANUP:
     graph_free(&digraph);

src/gtsp.c

@@ -213,34 +213,12 @@ int GTSP_add_cutting_planes(struct LP *lp, struct GTSP *data)
             }
             else
             {
-                log_debug("No more cuts found.\n");
+                log_debug("No additional cuts found.\n");
                 break;
             }
         }
 
-        log_debug("Reoptimizing...\n");
-
-        double time_before_lp = get_current_time();
         int is_infeasible;
-        rval = LP_optimize(lp, &is_infeasible);
-        abort_if(rval, "LP_optimize failed");
-        double time_after_lp = get_current_time();
-
-        double obj_val;
-        rval = LP_get_obj_val(lp, &obj_val);
-        abort_if(rval, "LP_get_obj_val failed");
-
-        log_debug("    obj val = %.4lf\n", obj_val);
-        log_debug("    time = %.2lf\n", time_after_lp - time_before_lp);
-
-        if (time_after_lp - time_before_lp > 10.0)
-        {
-            log_debug("LP is too slow. Removing slack constraints...\n");
-            int start = data->graph->node_count + data->cluster_count;
-            rval = LP_remove_slacks(lp, start, LP_EPSILON);
-            abort_if(rval, "LP_remove_slacks failed");
-        }
-
         rval = LP_optimize(lp, &is_infeasible);
         abort_if(rval, "LP_optimize failed");
 
@@ -476,7 +454,8 @@ int GTSP_solution_found(struct GTSP *data, double *x)
 }
 
 double FLOW_CPU_TIME = 0;
-double LP_CPU_TIME = 0;
+double LP_SOLVE_TIME = 0;
+double LP_CUT_POOL_TIME = 0;
 int LP_OPTIMIZE_COUNT = 0;
 
 int GTSP_main(int argc, char **argv)
@@ -556,7 +535,8 @@ int GTSP_main(int argc, char **argv)
     log_info("Max-flow calls: %d\n", FLOW_MAX_FLOW_COUNT);
     log_info("Max-flow computation time: %.2lf\n", FLOW_CPU_TIME);
     log_info("LP optimize calls: %d\n", LP_OPTIMIZE_COUNT);
-    log_info("LP solving time: %.2lf\n", LP_CPU_TIME);
+    log_info("LP solving time: %.2lf\n", LP_SOLVE_TIME);
+    log_info("LP cut pool management time: %.2lf\n", LP_CUT_POOL_TIME);
 
     CLEANUP:
     GTSP_free(&data);

src/lp.c

@@ -2,12 +2,10 @@
 #include <stdlib.h>
 #include <string.h>
 #include <cplex.h>
+#include <assert.h>
 #include "lp.h"
 #include "util.h"
 
-#define LP_EPSILON 0.000001
-#define MAX_CUT_POOL_SIZE 100000
-
 int LP_open(struct LP *lp)
 {
     int rval = 0;
@@ -16,6 +14,7 @@ int LP_open(struct LP *lp)
     lp->cplex_env = CPXopenCPLEX(&rval);
     abort_if(rval, "CPXopenCPLEX failed");
 
+    lp->cut_pool_size = 0;
     lp->cut_pool = (struct Row **) malloc(
             MAX_CUT_POOL_SIZE * sizeof(struct Row *));
     abort_if(!lp->cut_pool, "could not allocate cut_pool");
@@ -121,7 +120,42 @@ int LP_change_bound(struct LP *lp, int col, char lower_or_upper, double bnd)
 }
 
 extern int LP_OPTIMIZE_COUNT;
-extern double LP_CPU_TIME;
+extern double LP_SOLVE_TIME;
+extern double LP_CUT_POOL_TIME;
+
+int LP_update_cut_ages(struct LP *lp)
+{
+    int rval = 0;
+
+    double *slacks = 0;
+
+    int numrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp);
+
+    slacks = (double *) malloc(numrows * sizeof(double));
+    abort_if(!slacks, "could not allocate slacks");
+
+    rval = CPXgetslack(lp->cplex_env, lp->cplex_lp, slacks, 0, numrows - 1);
+    abort_if(rval, "CPXgetslack failed");
+
+    int count = 0;
+    for (int i = 0; i < lp->cut_pool_size; i++)
+    {
+        struct Row *cut = lp->cut_pool[i];
+        if (cut->cplex_row_index < 0) continue;
+        assert(cut->cplex_row_index < numrows);
+
+        if (slacks[cut->cplex_row_index] < -LP_EPSILON)
+            cut->age++;
+        else
+            cut->age = 0;
+
+        if (cut->age > 5) count++;
+    }
+
+    CLEANUP:
+    if (slacks) free(slacks);
+    return rval;
+}
 
 int LP_optimize(struct LP *lp, int *infeasible)
 {
@@ -131,10 +165,17 @@ int LP_optimize(struct LP *lp, int *infeasible)
 
     *infeasible = 0;
 
-    double current = get_current_time();
+    int numrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp);
+    int numcols = CPXgetnumcols(lp->cplex_env, lp->cplex_lp);
+
+    log_debug("Optimizing LP (%d rows %d cols)...\n", numrows, numcols);
+
+    double time_before = get_current_time();
     rval = CPXdualopt(lp->cplex_env, lp->cplex_lp);
     abort_if(rval, "CPXdualopt failed");
-    LP_CPU_TIME += get_current_time() - current;
+
+    double time_after = get_current_time();
+    LP_SOLVE_TIME += time_after - time_before;
 
     solstat = CPXgetstat(lp->cplex_env, lp->cplex_lp);
     if (solstat == CPX_STAT_INFEASIBLE)
@@ -148,38 +189,82 @@ int LP_optimize(struct LP *lp, int *infeasible)
                 solstat != CPX_STAT_OPTIMAL_INFEAS, "Invalid solution status");
     }
 
+    double objval;
+    rval = LP_get_obj_val(lp, &objval);
+    abort_if(rval, "LP_get_obj_val failed");
+
+    log_debug("    obj val = %.4lf\n", objval);
+    log_debug("    time = %.4lf\n", time_after - time_before);
+
+    time_before = get_current_time();
+    rval = LP_update_cut_ages(lp);
+    abort_if(rval, "LP_update_cut_ages failed");
+
+    rval = LP_remove_old_cuts(lp);
+    abort_if(rval, "LP_remove_old_cuts failed");
+
+    time_after = get_current_time();
+    LP_CUT_POOL_TIME += time_after - time_before;
+
     CLEANUP:
     return rval;
 }
 
-int LP_remove_slacks(struct LP *lp, int first_row, double max_slack)
+int LP_remove_old_cuts(struct LP *lp)
 {
     int rval = 0;
 
-    double *slacks = 0;
     int *should_remove = 0;
 
     int numrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp);
-    if (numrows < 5000) return 0;
+    log_verbose("    numrows=%d\n", numrows);
 
     should_remove = (int *) malloc((numrows + 1) * sizeof(int));
     abort_if(!should_remove, "could not allocate should_remove");
 
-    slacks = (double *) malloc(numrows * sizeof(double));
-    abort_if(!slacks, "could not allocate slacks");
+    for (int i = 0; i < numrows; i++)
+        should_remove[i] = 0;
+    should_remove[numrows] = 0;
 
-    rval = CPXgetslack(lp->cplex_env, lp->cplex_lp, slacks, 0, numrows - 1);
-    abort_if(rval, "CPXgetslack failed");
+    log_verbose("Old cplex row index:\n");
+    for (int i = 0; i < lp->cut_pool_size; i++)
+        log_verbose("    %d\n", lp->cut_pool[i]->cplex_row_index);
+
+    log_verbose("Should remove:\n");
+    for (int i = 0; i < lp->cut_pool_size; i++)
+    {
+        struct Row *cut = lp->cut_pool[i];
+        if (cut->age <= MAX_CUT_AGE) continue;
+        if (cut->cplex_row_index < 0) continue;
 
+        should_remove[cut->cplex_row_index] = 1;
+        log_verbose("    %d\n", cut->cplex_row_index);
+    }
+
+    int count = 0;
     for (int i = 0; i < numrows; i++)
-        should_remove[i] = (slacks[i] < -max_slack);
-    should_remove[numrows] = 0;
+    {
+        if (!should_remove[i]) continue;
+
+        for (int j = 0; j < lp->cut_pool_size; j++)
+        {
+            struct Row *cut = lp->cut_pool[j];
+
+            if (cut->cplex_row_index == i - count) cut->cplex_row_index = 0;
+            else if (cut->cplex_row_index > i - count) cut->cplex_row_index--;
+        }
+
+        count++;
+    }
+
+    log_verbose("New cplex row index:\n");
+    for (int i = 0; i < lp->cut_pool_size; i++)
+        log_verbose("    %d\n", lp->cut_pool[i]->cplex_row_index);
 
-    log_debug("Deleting constraints...\n");
     int start = 0;
     int end = -1;
-    int count = 0;
-    for (int i = first_row; i < numrows; i++)
+    count = 0;
+    for (int i = 0; i < numrows + 1; i++)
     {
         if (should_remove[i])
         {
@@ -202,11 +287,24 @@ int LP_remove_slacks(struct LP *lp, int first_row, double max_slack)
         }
     }
 
-    log_info("Removed %d of %d constraints\n", count, numrows);
+    numrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp);
+    log_verbose("    new numrows=%d\n", numrows);
+
+    for (int i = 0; i < lp->cut_pool_size; i++)
+    {
+        struct Row *cut = lp->cut_pool[i];
+        assert(cut->cplex_row_index < numrows);
+    }
+
+    if (count > 0)
+    {
+        log_info("Found and removed %d old cuts\n", count);
+        rval = CPXdualopt(lp->cplex_env, lp->cplex_lp);
+        abort_if(rval, "CPXoptimize failed");
+    }
 
     CLEANUP:
     if (should_remove) free(should_remove);
-    if (slacks) free(slacks);
     return rval;
 }
 
@@ -280,11 +378,13 @@ int LP_add_cut(struct LP *lp, struct Row *cut)
 {
     int rval = 0;
 
+    double time_before = get_current_time();
     rval = LP_update_hash(cut);
     abort_if(rval, "LP_update_hash failed");
 
     for (int i = 0; i < lp->cut_pool_size; i++)
     {
+        if (lp->cut_pool[i]->cplex_row_index < 0) continue;
         if (lp->cut_pool[i]->hash != cut->hash) continue;
         if (!compare_cuts(lp->cut_pool[i], cut))
         {
@@ -295,6 +395,7 @@ int LP_add_cut(struct LP *lp, struct Row *cut)
         }
     }
 
+    abort_if(lp->cut_pool_size > MAX_CUT_POOL_SIZE, "Cut pool is too large");
     lp->cut_pool[lp->cut_pool_size++] = cut;
 
     int rmatbeg = 0;
@@ -302,6 +403,12 @@ int LP_add_cut(struct LP *lp, struct Row *cut)
             cut->rmatind, cut->rmatval);
     abort_if(rval, "LP_add_rows failed");
 
+    cut->cplex_row_index = CPXgetnumrows(lp->cplex_env, lp->cplex_lp) - 1;
+    cut->age = 0;
+
+    double time_after = get_current_time();
+    LP_CUT_POOL_TIME += time_after - time_before;
+
     CLEANUP:
     return rval;
 }

src/lp.h

@@ -4,6 +4,8 @@
 #include <cplex.h>
 
 #define LP_EPSILON 0.000001
+#define MAX_CUT_POOL_SIZE 1000000
+#define MAX_CUT_AGE 10
 
 struct LP
 {
@@ -17,7 +19,8 @@ struct LP
 struct Row
 {
     unsigned long hash;
-//    int cplex_row_index;
+    int cplex_row_index;
+    int age;
 
     int nz;
     char sense;
@@ -69,7 +72,7 @@ int LP_get_x(struct LP *lp, double *x);
 
 int LP_get_num_cols(struct LP *lp);
 
-int LP_remove_slacks(struct LP *lp, int start, double max_slack);
+int LP_remove_old_cuts(struct LP *lp);
 
 int LP_update_hash(struct Row *cut);