Run multiple rounds with decreasing cut violation; remove slack constraints

src/gtsp-subtour.c

@@ -155,7 +155,8 @@ int static add_subtour_cut(
         abort_if(sum <= rhs - LP_EPSILON, "cannot add invalid cut");
     }
 
-    rval = LP_add_rows(lp, 1, newnz, &rhs, &sense, &rmatbeg, rmatind, rmatval);
+    rval = LP_add_rows(lp, 1, newnz, &rhs, &sense, &rmatbeg, rmatind,
+            rmatval);
     abort_if(rval, "LP_add_rows failed");
 
     CLEANUP:
@@ -165,7 +166,7 @@ int static add_subtour_cut(
 }
 
 int find_exact_subtour_cuts(
-        struct LP *lp, struct GTSP *data, int *total_added_cuts)
+        struct LP *lp, struct GTSP *data, int *total_added_cuts, double min_cut_violation)
 {
     int rval = 0;
 
@@ -201,7 +202,7 @@ int find_exact_subtour_cuts(
 
     // Constraints (2.1)
     rval = find_exact_subtour_cuts_cluster_to_cluster(lp, data, &digraph,
-            capacities, &added_cuts_count);
+            capacities, &added_cuts_count, min_cut_violation);
     abort_if(rval, "find_exact_subtour_cuts_cluster_to_cluster failed");
 
     log_debug("Added %d cluster-to-cluster subtour cuts\n", added_cuts_count);
@@ -210,7 +211,7 @@ int find_exact_subtour_cuts(
 
     // Constraints (2.2)
     rval = find_exact_subtour_cuts_node_to_cluster(lp, data, x, &digraph,
-            capacities, &added_cuts_count);
+            capacities, &added_cuts_count, min_cut_violation);
     abort_if(rval, "find_exact_subtour_cuts_node_to_cluster failed");
 
     log_debug("Added %d node-to-cluster subtour cuts\n", added_cuts_count);
@@ -219,7 +220,7 @@ int find_exact_subtour_cuts(
 
     // Constraints (2.3)
     rval = find_exact_subtour_cuts_node_to_node(lp, data, x, &digraph,
-            capacities, &added_cuts_count);
+            capacities, &added_cuts_count, min_cut_violation);
     abort_if(rval, "find_exact_subtour_cuts_node_to_node failed");
 
     log_debug("Added %d node-to-node subtour cuts\n", added_cuts_count);
@@ -238,7 +239,8 @@ int find_exact_subtour_cuts_node_to_node(
         double *x,
         struct Graph *digraph,
         double *capacities,
-        int *added_cuts_count)
+        int *added_cuts_count,
+        double min_cut_violation)
 {
     int rval = 0;
 
@@ -286,7 +288,7 @@ int find_exact_subtour_cuts_node_to_node(
                 &flow_value);
         abort_if(rval, "flow_find_max_flow failed");
 
-        if (flow_value >= 2 * (x[i] + x[j] - 1) - MIN_CUT_VIOLATION)
+        if (flow_value >= 2 * (x[i] + x[j] - 1) - min_cut_violation)
             continue;
 
         log_verbose("Marked nodes:\n");
@@ -323,7 +325,8 @@ int find_exact_subtour_cuts_node_to_cluster(
         double *x,
         struct Graph *digraph,
         double *capacities,
-        int *added_cuts_count)
+        int *added_cuts_count,
+        double min_cut_violation)
 {
     int rval = 0;
 
@@ -367,7 +370,7 @@ int find_exact_subtour_cuts_node_to_cluster(
 
             deactivate_cluster_node(capacities, to);
 
-            if (flow_value + MIN_CUT_VIOLATION >= 2 * x[i])
+            if (flow_value + min_cut_violation >= 2 * x[i])
                 continue;
 
             log_verbose("Marked nodes:\n");
@@ -409,7 +412,8 @@ int find_exact_subtour_cuts_cluster_to_cluster(
         struct GTSP *data,
         struct Graph *digraph,
         double *capacities,
-        int *added_cuts_count)
+        int *added_cuts_count,
+        double min_cut_violation)
 {
     int rval = 0;
 
@@ -455,7 +459,7 @@ int find_exact_subtour_cuts_cluster_to_cluster(
             deactivate_cluster_node(capacities, from);
             deactivate_cluster_node(capacities, to);
 
-            if (flow_value >= 2 - MIN_CUT_VIOLATION) continue;
+            if (flow_value >= 2 - min_cut_violation) continue;
 
             log_verbose("Marked nodes:\n");
             for (int k = 0; k < graph->node_count; k++)

src/gtsp-subtour.h

@@ -7,13 +7,13 @@ void deactivate_cluster_node(double *capacities, struct Node *cluster_node);
 
 void activate_cluster_node(double *capacities, struct Node *cluster_node);
 
-
 int find_exact_subtour_cuts_cluster_to_cluster(
         struct LP *lp,
         struct GTSP *data,
         struct Graph *digraph,
         double *capacities,
-        int *added_cuts_count);
+        int *added_cuts_count,
+        double min_cut_violation);
 
 int find_exact_subtour_cuts_node_to_cluster(
         struct LP *lp,
@@ -21,7 +21,8 @@ int find_exact_subtour_cuts_node_to_cluster(
         double *x,
         struct Graph *digraph,
         double *capacities,
-        int *added_cuts_count);
+        int *added_cuts_count,
+        double min_cut_violation);
 
 int find_exact_subtour_cuts_node_to_node(
         struct LP *lp,
@@ -29,10 +30,13 @@ int find_exact_subtour_cuts_node_to_node(
         double *x,
         struct Graph *digraph,
         double *capacities,
-        int *added_cuts_count);
+        int *added_cuts_count,
+        double min_cut_violation);
 
 int find_exact_subtour_cuts(
-        struct LP *lp, struct GTSP *data, int *total_added_cuts);
-
+        struct LP *lp,
+        struct GTSP *data,
+        int *total_added_cuts,
+        double min_cut_violation);
 
 #endif //_PROJECT_GTSP_SUBTOUR_H_

src/gtsp.c

@@ -190,33 +190,62 @@ int GTSP_add_cutting_planes(struct LP *lp, struct GTSP *data)
     int round = 0;
     int added_cuts_count = 0;
 
+    int violation_total = 3;
+    int violation_current = 0;
+    double violations[] = {1.0, 0.1, LP_EPSILON};
+
     while (1)
     {
         round++;
         int added_cuts_this_round = 0;
 
-        log_debug("Finding subtour cuts, round %d...\n", round);
+        log_debug("Finding subtour cuts, round %d, violation %.4lf...\n", round,
+                violations[violation_current]);
 
-        rval = find_exact_subtour_cuts(lp, data, &added_cuts_this_round);
+        rval = find_exact_subtour_cuts(lp, data, &added_cuts_this_round,
+                violations[violation_current]);
         abort_if(rval, "find_exact_subtour_cuts failed");
 
         if (added_cuts_this_round == 0)
         {
-            log_debug("No more subtour cuts found.\n");
+            ++violation_current;
+            if (violation_current < violation_total)
+            {
+                log_debug("No cuts found. Decreasing minimum cut violation.\n");
+                continue;
+            }
+            else
+            {
+                log_debug("No more 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");
 
         if (is_infeasible) break;
 
@@ -451,7 +480,6 @@ int GTSP_solution_found(struct GTSP *data, double *x)
     return rval;
 }
 
-
 double FLOW_CPU_TIME = 0;
 double LP_CPU_TIME = 0;
 int LP_OPTIMIZE_COUNT = 0;

src/gtsp.h

@@ -19,8 +19,6 @@ struct GTSP
     double *y_coordinates;
 };
 
-static const double MIN_CUT_VIOLATION = 0.5;
-
 int GTSP_create_random_problem(
         int node_count, int cluster_count, int grid_size, struct GTSP *data);
 

src/lp.c

@@ -100,7 +100,6 @@ int LP_add_cols(
 
     rval = CPXaddcols(lp->cplex_env, lp->cplex_lp, newcols, newnz, obj, cmatbeg,
             cmatind, cmatval, lb, ub, (char **) NULL);
-
     abort_if(rval, "CPXaddcols failed");
 
     CLEANUP:
@@ -143,7 +142,7 @@ int LP_optimize(struct LP *lp, int *infeasible)
     }
     else
     {
-        abort_if(solstat != CPX_STAT_OPTIMAL &&
+        abort_if(solstat != CPX_STAT_OPTIMAL && solstat != CPXMIP_OPTIMAL &&
                 solstat != CPX_STAT_OPTIMAL_INFEAS, "Invalid solution status");
     }
 
@@ -151,6 +150,64 @@ int LP_optimize(struct LP *lp, int *infeasible)
     return rval;
 }
 
+int LP_remove_slacks(struct LP *lp, int first_row, double max_slack)
+{
+    int rval = 0;
+
+    double *slacks = 0;
+    int *should_remove = 0;
+
+    int numrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp);
+    if(numrows < 5000) return 0;
+
+    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");
+
+    rval = CPXgetslack(lp->cplex_env, lp->cplex_lp, slacks, 0, numrows - 1);
+    abort_if(rval, "CPXgetslack failed");
+
+    for (int i = 0; i < numrows; i++)
+        should_remove[i] = (slacks[i] < -max_slack);
+    should_remove[numrows] = 0;
+
+    log_debug("Deleting constraints...\n");
+    int start = 0;
+    int end = -1;
+    int count = 0;
+    for (int i = first_row; i < numrows; i++)
+    {
+        if (should_remove[i])
+        {
+            end++;
+        }
+        else
+        {
+            if (end >= start)
+            {
+                rval = CPXdelrows(lp->cplex_env, lp->cplex_lp, start - count,
+                        end - count);
+                abort_if(rval, "CPXdelrows failed");
+                log_verbose("    %d %d (%d)\n", start, end, end-start+1);
+
+                count += end - start + 1;
+            }
+
+            start = i+1;
+            end = i;
+        }
+    }
+
+    log_info("Removed %d of %d constraints\n", count, numrows);
+
+    CLEANUP:
+    if (should_remove) free(should_remove);
+    if (slacks) free(slacks);
+    return rval;
+}
+
 int LP_get_obj_val(struct LP *lp, double *obj)
 {
     int rval = 0;

src/lp.h

@@ -54,4 +54,6 @@ 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);
+
 #endif