From f82002aa11568d2a7f1cd0210c4383329a54d50d Mon Sep 17 00:00:00 2001 From: Alinson Xavier Date: Mon, 30 Mar 2015 11:43:40 -0400 Subject: [PATCH] Implement comb cuts separation --- src/branch_and_cut.c | 12 +- src/graph.c | 6 +- src/graph.h | 3 + src/gtsp-comb.c | 476 +++++++++++++++++++++++++++++++++++++++++++ src/gtsp-comb.h | 6 + src/gtsp-subtour.c | 42 ++-- src/gtsp.c | 123 +++++------ src/gtsp.h | 7 - src/lp.c | 6 +- src/params.h | 10 + src/util.h | 3 +- 11 files changed, 600 insertions(+), 94 deletions(-) create mode 100644 src/gtsp-comb.c create mode 100644 src/gtsp-comb.h create mode 100644 src/params.h diff --git a/src/branch_and_cut.c b/src/branch_and_cut.c index ef4baad..82d7520 100644 --- a/src/branch_and_cut.c +++ b/src/branch_and_cut.c @@ -100,7 +100,7 @@ static int BNC_solve_node(struct BNC *bnc, int depth) if (ceil(objval) > *best_val + LP_EPSILON) { log_debug("Branch pruned by bound (%.2lf > %.2lf).\n", objval, - *best_val); + *best_val); rval = 0; goto CLEANUP; } @@ -125,7 +125,7 @@ static int BNC_solve_node(struct BNC *bnc, int depth) if (ceil(objval) > *best_val + LP_EPSILON) { log_debug("Branch pruned by bound (%.2lf > %.2lf).\n", objval, - *best_val); + *best_val); rval = 0; goto CLEANUP; } @@ -174,7 +174,7 @@ static int BNC_branch_node(struct BNC *bnc, double *x, int depth) int best_branch_var = BNC_find_best_branching_var(x, num_cols); log_debug("Branching on variable x%d = %.6lf (depth %d)...\n", - best_branch_var, x[best_branch_var], depth); + best_branch_var, x[best_branch_var], depth); log_debug("Fixing variable x%d to one...\n", best_branch_var); rval = LP_change_bound(lp, best_branch_var, 'L', 1.0); @@ -204,11 +204,17 @@ static int BNC_branch_node(struct BNC *bnc, double *x, int depth) static int BNC_is_integral(double *x, int num_cols) { +#ifdef ALLOW_FRACTIONAL_SOLUTIONS + UNUSED(num_cols); + UNUSED(x); + return 1; +#else for (int i = 0; i < num_cols; i++) if (x[i] > LP_EPSILON && x[i] < 1.0 - LP_EPSILON) return 0; return 1; +#endif } static int BNC_find_best_branching_var(double *x, int num_cols) diff --git a/src/graph.c b/src/graph.c index 2ae1f9f..1439606 100644 --- a/src/graph.c +++ b/src/graph.c @@ -11,6 +11,8 @@ void graph_init(struct Graph *graph) graph->adj = 0; graph->node_count = 0; graph->edge_count = 0; + graph->x_coordinates = 0; + graph->y_coordinates = 0; } void graph_free(struct Graph *graph) @@ -20,6 +22,8 @@ void graph_free(struct Graph *graph) if (graph->edges) free(graph->edges); if (graph->nodes) free(graph->nodes); if (graph->adj) free(graph->adj); + if (graph->x_coordinates) free(graph->x_coordinates); + if (graph->y_coordinates) free(graph->y_coordinates); } int graph_build( @@ -86,7 +90,7 @@ int graph_build( n = &graph->nodes[b]; n->adj[n->degree].neighbor_index = a; n->adj[n->degree].edge_index = i; - n->adj[n->degree].neighbor = &graph->nodes[b]; + n->adj[n->degree].neighbor = &graph->nodes[a]; n->adj[n->degree].edge = &graph->edges[i]; n->degree++; } diff --git a/src/graph.h b/src/graph.h index 67eed3d..9b00587 100644 --- a/src/graph.h +++ b/src/graph.h @@ -41,6 +41,9 @@ struct Graph struct Edge *edges; struct Node *nodes; + double *x_coordinates; + double *y_coordinates; + struct Adjacency *adj; }; diff --git a/src/gtsp-comb.c b/src/gtsp-comb.c new file mode 100644 index 0000000..31e4dcf --- /dev/null +++ b/src/gtsp-comb.c @@ -0,0 +1,476 @@ +#include +#include "gtsp.h" +#include "util.h" + +int add_comb_cut( + struct LP *lp, + struct Graph *graph, + int current_component, + int *clusters, + int *components, + int *component_sizes, + int *teeth, + int tooth_count, + double *x) +{ + int rval = 0; + + char sense = 'G'; + + const int node_count = graph->node_count; + const int edge_count = graph->edge_count; + + struct Row *cut = 0; + int *rmatind = 0; + double *rmatval = 0; + + rmatind = (int *) malloc((node_count + edge_count) * sizeof(int)); + rmatval = (double *) malloc((node_count + edge_count) * sizeof(double)); + abort_if(!rmatind, "could not allocate rmatind"); + abort_if(!rmatval, "could not allocate rmatval"); + + double rhs = -component_sizes[current_component] - tooth_count + + (tooth_count + 1) / 2; + + int nz = 0; + + // Edges inside handle + for (int i = 0; i < edge_count; i++) + { + struct Edge *e = &graph->edges[i]; + if (components[clusters[e->from->index]] != current_component) continue; + if (components[clusters[e->to->index]] != current_component) continue; + + rmatind[nz] = node_count + e->index; + rmatval[nz] = -1.0; + nz++; + + log_verbose(" handle (%d %d)\n", e->from->index, e->to->index); + } + + // Edges inside each tooth + for (int i = 0; i < edge_count; i++) + { + struct Edge *e = &graph->edges[i]; + struct Node *from = e->from; + struct Node *to = e->to; + + if (teeth[clusters[from->index]] < 0) continue; + if (teeth[clusters[to->index]] < 0) continue; + if (teeth[clusters[from->index]] != teeth[clusters[to->index]]) + continue; + + log_verbose(" tooth (%d %d)\n", e->from->index, e->to->index); + + rmatind[nz] = node_count + e->index; + rmatval[nz] = -1.0; + nz++; + } + + // Lifting of the nodes + for (int i = 0; i < node_count; i++) + { + double val; + struct Node *n = &graph->nodes[i]; + int c = clusters[n->index]; + + if (components[c] == current_component) + val = (teeth[c] < 0 ? 1.0 : 2.0); + else + val = (teeth[c] < 0 ? 0.0 : 1.0); + + if (val == 0.0) continue; + + rmatind[nz] = n->index; + rmatval[nz] = val; + nz++; + + rhs += val; + } + + log_debug("Generated cut:\n"); + for (int i = 0; i < nz; i++) + log_debug(" %.2lf x%d (%.4lf)\n", rmatval[i], rmatind[i], + x[rmatind[i]]); + log_debug(" %c %.2lf\n", sense, rhs); + + if (OPTIMAL_X) + { + double sum = 0; + for (int i = 0; i < nz; i++) + sum += rmatval[i] * OPTIMAL_X[rmatind[i]]; + abort_if(sum <= rhs - LP_EPSILON, "cannot add invalid cut"); + } + + double lhs = 0.0; + for (int i = 0; i < nz; i++) + lhs += rmatval[i] * x[rmatind[i]]; + + log_verbose("Violation: %.4lf >= %.4lf\n", lhs, rhs); + + if (lhs + LP_EPSILON > rhs) goto CLEANUP; + + cut = (struct Row *) malloc(sizeof(struct Row)); + abort_if(!cut, "could not allocate cut"); + + cut->nz = nz; + cut->sense = sense; + cut->rhs = rhs; + cut->rmatval = rmatval; + cut->rmatind = rmatind; + + rval = LP_add_cut(lp, cut); + abort_if(rval, "LP_add_cut failed"); + + CLEANUP: + if (rmatind) free(rmatind); + if (rmatval) free(rmatval); + return rval; +} + +int find_components( + struct Graph *graph, double *x, int *components, int *component_sizes) +{ + int rval = 0; + struct Node **stack = 0; + + const int node_count = graph->node_count; + + for (int i = 0; i < node_count; i++) + { + components[i] = -1; + graph->nodes[i].mark = 0; + } + + int stack_top = 0; + stack = (struct Node **) malloc(node_count * sizeof(struct Node *)); + abort_if(!stack, "could not allocate stack"); + + for (int i = 0; i < node_count; i++) + { + struct Node *root = &graph->nodes[i]; + if (root->mark) continue; + + stack[stack_top++] = root; + + while (stack_top > 0) + { + struct Node *n = stack[--stack_top]; + components[n->index] = i; + + for (int j = 0; j < n->degree; j++) + { + struct Adjacency *adj = &n->adj[j]; + struct Node *neighbor = adj->neighbor; + + if (neighbor->mark) continue; + + double x_e = x[adj->edge->index]; + if (x_e < LP_EPSILON) continue; + if (x_e > 1 - LP_EPSILON) continue; + + stack[stack_top++] = neighbor; + neighbor->mark = 1; + } + } + } + + for (int i = 0; i < node_count; i++) + component_sizes[i] = 0; + + for (int i = 0; i < node_count; i++) + component_sizes[components[i]]++; + + log_verbose("Components:\n"); + for (int i = 0; i < graph->node_count; i++) + log_verbose(" %d %d\n", i, components[i]); + + log_verbose("Component sizes:\n"); + for (int i = 0; i < graph->node_count; i++) + log_verbose(" %d %d\n", i, component_sizes[i]); + + CLEANUP: + if (stack) free(stack); + return rval; +} + +int find_teeth( + struct Graph *graph, + double *x, + int current_component, + int *components, + int *teeth, + int *tooth_count) +{ + int rval = 0; + + const int node_count = graph->node_count; + const int edge_count = graph->edge_count; + + for (int i = 0; i < node_count; i++) + { + graph->nodes[i].mark = 0; + teeth[i] = -1; + } + + *tooth_count = 0; + + for (int i = 0; i < edge_count; i++) + { + struct Edge *e = &graph->edges[i]; + struct Node *from = e->from; + struct Node *to = e->to; + + if (x[e->index] < 1 - LP_EPSILON) continue; + + if (to->mark || from->mark) continue; + + int z = 0; + if (components[from->index] == current_component) z++; + if (components[to->index] == current_component) z++; + if (z != 1) continue; + + to->mark = 1; + from->mark = 1; + + teeth[to->index] = *tooth_count; + teeth[from->index] = *tooth_count; + + (*tooth_count)++; + } + + CLEANUP: + return rval; +} + +int write_shrunken_graph( + double *shrunken_x, + struct Graph *shrunken_graph, + int const cluster_count); + +static int shrink_clusters( + const struct GTSP *data, + double *x, + struct Graph *shrunken_graph, + double *shrunken_x) +{ + int rval = 0; + + double *x_coords = 0; + double *y_coords = 0; + int *cluster_sizes = 0; + + const int *clusters = data->clusters; + const int cluster_count = data->cluster_count; + const struct Graph *graph = data->graph; + + int *edges = 0; + int *edge_map = 0; + int edge_count = (cluster_count * (cluster_count - 1)) / 2; + + edge_map = (int *) malloc(cluster_count * cluster_count * sizeof(int)); + abort_if(!edge_map, "could not allocate edge_map"); + + edges = (int *) malloc(2 * edge_count * sizeof(int)); + abort_if(!edges, "could not allocate edges"); + + cluster_sizes = (int *) malloc(cluster_count * sizeof(int)); + x_coords = (double *) malloc(cluster_count * sizeof(double)); + y_coords = (double *) malloc(cluster_count * sizeof(double)); + + abort_if(!cluster_sizes, "could not allocate cluster_sizes"); + abort_if(!x_coords, "could not allocate x_coords"); + abort_if(!y_coords, "could not allocate y_coords"); + + for (int i = 0; i < cluster_count; i++) + { + x_coords[i] = 0.0; + y_coords[i] = 0.0; + cluster_sizes[i] = 0; + } + + for (int i = 0; i < graph->node_count; i++) + { + struct Node *n = &graph->nodes[i]; + int c = clusters[n->index]; + + cluster_sizes[c]++; + x_coords[c] += graph->x_coordinates[n->index]; + y_coords[c] += graph->y_coordinates[n->index]; + } + + for (int i = 0; i < cluster_count; i++) + { + x_coords[i] = x_coords[i] / cluster_sizes[i]; + y_coords[i] = y_coords[i] / cluster_sizes[i]; + } + + shrunken_graph->x_coordinates = x_coords; + shrunken_graph->y_coordinates = y_coords; + + int curr_edge = 0; + for (int i = 0; i < cluster_count; i++) + { + for (int j = i + 1; j < cluster_count; j++) + { + edges[curr_edge * 2] = i; + edges[curr_edge * 2 + 1] = j; + edge_map[i * cluster_count + j] = curr_edge; + edge_map[j * cluster_count + i] = curr_edge; + curr_edge++; + } + } + + assert(curr_edge == edge_count); + + rval = graph_build(cluster_count, edge_count, edges, 0, shrunken_graph); + abort_if(rval, "graph_build failed"); + + for (int i = 0; i < edge_count; i++) + shrunken_x[i] = 0.0; + + for (int i = 0; i < graph->edge_count; i++) + { + struct Edge *e = &graph->edges[i]; + + int from = clusters[e->from->index]; + int to = clusters[e->to->index]; + int shunk_e_index = edge_map[from * cluster_count + to]; + + shrunken_x[shunk_e_index] += x[graph->node_count + e->index]; + } + + CLEANUP: + if (edges) free(edges); + if (edge_map) free(edge_map); + if (cluster_sizes) free(cluster_sizes); + return rval; +} + +int find_comb_cuts(struct LP *lp, struct GTSP *data) +{ + int rval = 0; + + double *x = 0; + double *shrunken_x = 0; + + int *teeth = 0; + int *components = 0; + int *component_sizes = 0; + + int num_cols = LP_get_num_cols(lp); + x = (double *) malloc(num_cols * sizeof(double)); + abort_if(!x, "could not allocate x"); + + rval = LP_get_x(lp, x); + abort_if(rval, "LP_get_x failed"); + + struct Graph shrunken_graph; + graph_init(&shrunken_graph); + + const int cluster_count = data->cluster_count; + const int shrunken_edge_count = (cluster_count * (cluster_count - 1)) / 2; + + shrunken_x = (double *) malloc(shrunken_edge_count * sizeof(double)); + abort_if(!shrunken_x, "could not allocate shrunken_x"); + + rval = shrink_clusters(data, x, &shrunken_graph, shrunken_x); + abort_if(rval, "shrink_clusters failed"); + + rval = write_shrunken_graph(shrunken_x, &shrunken_graph, cluster_count); + abort_if(rval, "write_shrunken_graph failed"); + + teeth = (int *) malloc(cluster_count * sizeof(int)); + components = (int *) malloc(cluster_count * sizeof(int)); + component_sizes = (int *) malloc(cluster_count * sizeof(int)); + + abort_if(!teeth, "could not allocate teeth"); + abort_if(!components, "could not allocate components"); + abort_if(!component_sizes, "could not allocate component_sizes"); + + rval = find_components(&shrunken_graph, shrunken_x, components, + component_sizes); + abort_if(rval, "find_components failed"); + + int original_cut_pool_size = lp->cut_pool_size; + for (int i = 0; i < cluster_count; i++) + { + if (component_sizes[i] < 3) continue; + + int tooth_count; + rval = find_teeth(&shrunken_graph, shrunken_x, i, components, teeth, + &tooth_count); + abort_if(rval, "find_teeth failed"); + + log_verbose("Component %d has %d teeth:\n", i, tooth_count); + for (int j = 0; j < cluster_count; j++) + { + if (teeth[j] < 0) continue; + log_verbose(" %d %d\n", j, teeth[j]); + } + + if (tooth_count % 2 == 0) continue; + + rval = add_comb_cut(lp, data->graph, i, data->clusters, components, + component_sizes, teeth, tooth_count, x); + abort_if(rval, "add_comb_cut failed"); + + } + int added_cuts_count = lp->cut_pool_size - original_cut_pool_size; + + if(added_cuts_count > 0) + log_debug(" %d combs\n", added_cuts_count); + + CLEANUP: + graph_free(&shrunken_graph); + if (teeth) free(teeth); + if (components) free(components); + if (component_sizes) free(component_sizes); + if (shrunken_x) free(shrunken_x); + return rval; +} + +int write_shrunken_graph( + double *shrunken_x, + struct Graph *shrunken_graph, + int const cluster_count) +{ + int rval = 0; + + FILE *file = 0; + + file = fopen("gtsp-shrunken.in", "w"); + abort_if(!file, "could not open file"); + + fprintf(file, "%d %d\n", (*shrunken_graph).node_count, cluster_count); + for (int i = 0; i < (*shrunken_graph).node_count; i++) + { + fprintf(file, "%.2lf %.2lf %d\n", (*shrunken_graph).x_coordinates[i], + (*shrunken_graph).y_coordinates[i], i); + } + + fclose(file); + + file = fopen("gtsp-shrunken.out", "w"); + abort_if(!file, "could not open file"); + + int positive_edge_count = 0; + for (int i = 0; i < (*shrunken_graph).edge_count; i++) + if (shrunken_x[i] > LP_EPSILON) + positive_edge_count++; + + fprintf(file, "%d %d\n", (*shrunken_graph).node_count, + (*shrunken_graph).edge_count); + + fprintf(file, "%d\n", positive_edge_count); + + for (int i = 0; i < (*shrunken_graph).edge_count; i++) + if (shrunken_x[i] > LP_EPSILON) + fprintf(file, "%d %d %.4lf\n", + (*shrunken_graph).edges[i].from->index, + (*shrunken_graph).edges[i].to->index, shrunken_x[i]); + fclose(file); + + CLEANUP: + return rval; +} \ No newline at end of file diff --git a/src/gtsp-comb.h b/src/gtsp-comb.h new file mode 100644 index 0000000..ee600e8 --- /dev/null +++ b/src/gtsp-comb.h @@ -0,0 +1,6 @@ +#ifndef PROJECT_GTSP_COMB_H +#define PROJECT_GTSP_COMB_H + +int find_comb_cuts(struct LP *lp, struct GTSP *data); + +#endif //PROJECT_GTSP_COMB_H diff --git a/src/gtsp-subtour.c b/src/gtsp-subtour.c index fa9d1ea..e60b2d3 100644 --- a/src/gtsp-subtour.c +++ b/src/gtsp-subtour.c @@ -18,7 +18,7 @@ int static build_flow_digraph( int digraph_node_count = node_count + data->cluster_count + 1; int digraph_edge_count = 4 * graph->edge_count + 2 * graph->node_count + - 2 * data->cluster_count; + 2 * data->cluster_count; digraph_edges = (int *) malloc(2 * digraph_edge_count * sizeof(int)); abort_if(!digraph_edges, "could not allocate digraph_edges"); @@ -143,7 +143,7 @@ int static add_subtour_cut( log_verbose("Generated cut:\n"); for (int i = 0; i < newnz; i++) - log_verbose("%8.2f x%d\n", rmatval[i], rmatind[i]); + log_verbose(" %8.2f x%d\n", rmatval[i], rmatind[i]); log_verbose(" %c %.2lf\n", sense, rhs); if (OPTIMAL_X) @@ -197,7 +197,7 @@ int find_exact_subtour_cuts( struct Graph digraph; graph_init(&digraph); int digraph_edge_count = 4 * graph->edge_count + 2 * graph->node_count + - 2 * data->cluster_count; + 2 * data->cluster_count; int original_cut_pool_size = lp->cut_pool_size; @@ -209,40 +209,42 @@ int find_exact_subtour_cuts( // Constraints (2.1) rval = find_exact_subtour_cuts_cluster_to_cluster(lp, data, &digraph, - capacities, min_cut_violation); + capacities, + min_cut_violation); abort_if(rval, "find_exact_subtour_cuts_cluster_to_cluster failed"); added_cuts_count = lp->cut_pool_size - original_cut_pool_size; if (added_cuts_count > 0) { - log_debug("Added %d cluster-to-cluster subtour cuts\n", - added_cuts_count); + log_debug(" %d cluster-to-cluster\n", + added_cuts_count); goto CLEANUP; } // Constraints (2.2) original_cut_pool_size = lp->cut_pool_size; rval = find_exact_subtour_cuts_node_to_cluster(lp, data, x, &digraph, - capacities, min_cut_violation); + capacities, + min_cut_violation); abort_if(rval, "find_exact_subtour_cuts_node_to_cluster failed"); added_cuts_count = lp->cut_pool_size - original_cut_pool_size; if (added_cuts_count > 0) { - log_debug("Added %d node-to-cluster subtour cuts\n", added_cuts_count); + log_debug(" %d node-to-cluster\n", added_cuts_count); goto CLEANUP; } // Constraints (2.3) original_cut_pool_size = lp->cut_pool_size; rval = find_exact_subtour_cuts_node_to_node(lp, data, x, &digraph, - capacities, min_cut_violation); + capacities, min_cut_violation); abort_if(rval, "find_exact_subtour_cuts_node_to_node failed"); added_cuts_count = lp->cut_pool_size - original_cut_pool_size; if (added_cuts_count > 0) { - log_debug("Added %d node-to-node subtour cuts\n", added_cuts_count); + log_debug(" %d node-to-node\n", added_cuts_count); goto CLEANUP; } @@ -304,7 +306,7 @@ int find_exact_subtour_cuts_node_to_node( double flow_value; rval = flow_find_max_flow(digraph, capacities, from, to, flow, - &flow_value); + &flow_value); abort_if(rval, "flow_find_max_flow failed"); if (flow_value >= 2 * (x[i] + x[j] - 1) - min_cut_violation) @@ -323,10 +325,10 @@ int find_exact_subtour_cuts_node_to_node( log_verbose("Cut edges:\n"); for (int k = 0; k < cut_edges_count; k++) log_verbose(" %d %d (%d)\n", cut_edges[k]->from->index, - cut_edges[k]->to->index, cut_edges[k]->index); + cut_edges[k]->to->index, cut_edges[k]->index); rval = add_subtour_cut(lp, graph, from, to, cut_edges, cut_edges_count, - 2); + 2); abort_if(rval, "add_subtour_cut failed"); } @@ -379,7 +381,7 @@ int find_exact_subtour_cuts_node_to_cluster( int cut_edges_count; rval = flow_find_max_flow(digraph, capacities, from, to, flow, - &flow_value); + &flow_value); abort_if(rval, "flow_find_max_flow failed"); log_verbose(" flow value = %.4lf\n", flow_value); @@ -405,11 +407,11 @@ int find_exact_subtour_cuts_node_to_cluster( struct Edge *e = cut_edges[k]; assert(e->from->mark != e->to->mark); log_verbose(" %d (%d) %d (%d) [%d]\n", e->from->index, - e->from->mark, e->to->index, e->to->mark, e->index); + e->from->mark, e->to->index, e->to->mark, e->index); } rval = add_subtour_cut(lp, graph, from, 0, cut_edges, - cut_edges_count, 1); + cut_edges_count, 1); abort_if(rval, "add_subtour_cut failed"); cuts_count++; @@ -445,7 +447,7 @@ int find_exact_subtour_cuts_cluster_to_cluster( abort_if(!flow, "could not allocate flow"); struct Node *root_node = &digraph->nodes[graph->node_count + - data->cluster_count]; + data->cluster_count]; for (int i = 0; i < data->cluster_count; i++) { @@ -464,7 +466,7 @@ int find_exact_subtour_cuts_cluster_to_cluster( log_verbose("Sending flow from cluster %d to cluster %d\n", i, j); rval = flow_find_max_flow(digraph, capacities, from, to, flow, - &flow_value); + &flow_value); abort_if(rval, "flow_find_max_flow failed"); @@ -488,10 +490,10 @@ int find_exact_subtour_cuts_cluster_to_cluster( log_verbose("Cut edges:\n"); for (int k = 0; k < cut_edges_count; k++) log_verbose(" %d %d (%d)\n", cut_edges[k]->from->index, - cut_edges[k]->to->index, cut_edges[k]->index); + cut_edges[k]->to->index, cut_edges[k]->index); rval = add_subtour_cut(lp, graph, 0, 0, cut_edges, cut_edges_count, - 0); + 0); abort_if(rval, "add_subtour_cut failed"); cuts_count++; diff --git a/src/gtsp.c b/src/gtsp.c index c98ba56..1966b31 100644 --- a/src/gtsp.c +++ b/src/gtsp.c @@ -2,13 +2,14 @@ #include #include #include +#include #include "gtsp.h" #include "geometry.h" #include "util.h" #include "flow.h" #include "branch_and_cut.h" -#include "math.h" #include "gtsp-subtour.h" +#include "gtsp-comb.h" double *OPTIMAL_X = 0; @@ -30,8 +31,6 @@ int GTSP_init_data(struct GTSP *data) data->clusters = 0; data->cluster_count = 0; - data->x_coordinates = 0; - data->y_coordinates = 0; data->graph = (struct Graph *) malloc(sizeof(struct Graph)); abort_if(!data->graph, "could not allocate data->graph"); @@ -50,8 +49,6 @@ void GTSP_free(struct GTSP *data) free(data->graph); if (data->clusters) free(data->clusters); - if (data->x_coordinates) free(data->x_coordinates); - if (data->y_coordinates) free(data->y_coordinates); } int GTSP_create_random_problem( @@ -90,7 +87,7 @@ int GTSP_create_random_problem( abort_if(!y_coords, "could not allocate y_coords\n"); rval = generate_random_clusters_2d(node_count, cluster_count, grid_size, - x_coords, y_coords, clusters); + x_coords, y_coords, clusters); abort_if(rval, "generate_random_clusters_2d failed"); int curr_edge = 0; @@ -102,7 +99,7 @@ int GTSP_create_random_problem( edges[curr_edge * 2] = i; edges[curr_edge * 2 + 1] = j; weights[curr_edge] = get_euclidean_distance(x_coords, y_coords, i, - j); + j); curr_edge++; } @@ -118,8 +115,8 @@ int GTSP_create_random_problem( data->graph = graph; data->clusters = clusters; data->cluster_count = cluster_count; - data->x_coordinates = x_coords; - data->y_coordinates = y_coords; + graph->x_coordinates = x_coords; + graph->y_coordinates = y_coords; CLEANUP: if (weights) free(weights); @@ -164,7 +161,7 @@ int GTSP_init_lp(struct LP *lp, struct GTSP *data) int cmatind[] = {i, node_count + clusters[i]}; rval = LP_add_cols(lp, 1, 2, &obj, &cmatbeg, cmatind, cmatval, &lb, - &ub); + &ub); abort_if(rval, "LP_add_cols failed"); } @@ -175,7 +172,7 @@ int GTSP_init_lp(struct LP *lp, struct GTSP *data) int cmatind[] = {edges[i].from->index, edges[i].to->index}; rval = LP_add_cols(lp, 1, 2, &obj, &cmatbeg, cmatind, cmatval, &lb, - &ub); + &ub); abort_if(rval, "LP_add_cols failed"); } @@ -186,43 +183,48 @@ int GTSP_init_lp(struct LP *lp, struct GTSP *data) int GTSP_add_cutting_planes(struct LP *lp, struct GTSP *data) { int rval = 0; - - int round = 0; - - int violation_total = 3; - int violation_current = 0; - double violations[] = {1.0, 0.1, LP_EPSILON}; + int current_round = 0; while (1) { - round++; + if (current_round > 0) + { + int is_infeasible; + rval = LP_optimize(lp, &is_infeasible); + abort_if(rval, "LP_optimize failed"); + + if (is_infeasible) break; + } + + current_round++; + + int original_cut_pool_size; + int added_cuts_count; - log_debug("Finding subtour cuts, round %d, violation %.4lf...\n", round, - violations[violation_current]); + original_cut_pool_size = lp->cut_pool_size; + log_debug("Finding subtour cuts, round %d...\n", current_round); - int original_cut_pool_size = lp->cut_pool_size; - rval = find_exact_subtour_cuts(lp, data, violations[violation_current]); + rval = find_exact_subtour_cuts(lp, data, LP_EPSILON); abort_if(rval, "find_exact_subtour_cuts failed"); - if (lp->cut_pool_size - original_cut_pool_size == 0) - { - if (++violation_current < violation_total) - { - log_debug("No cuts found. Decreasing minimum cut violation.\n"); - continue; - } - else - { - log_debug("No additional cuts found.\n"); - break; - } - } + added_cuts_count = lp->cut_pool_size - original_cut_pool_size; + if (added_cuts_count > 0) + continue; - int is_infeasible; - rval = LP_optimize(lp, &is_infeasible); - abort_if(rval, "LP_optimize failed"); +#ifdef ENABLE_COMB_INEQUALITIES + original_cut_pool_size = lp->cut_pool_size; + log_debug("Finding comb cuts, round %d...\n", current_round); - if (is_infeasible) break; + rval = find_comb_cuts(lp, data); + abort_if(rval, "find_comb_cuts failed"); + + added_cuts_count = lp->cut_pool_size - original_cut_pool_size; + if (added_cuts_count > 0) + continue; +#endif + + log_debug("No additional cuts found.\n"); + break; } CLEANUP: @@ -238,12 +240,14 @@ int GTSP_write_problem(struct GTSP *data, char *filename) file = fopen(filename, "w"); abort_if(!file, "could not open file"); - fprintf(file, "%d %d\n", data->graph->node_count, data->cluster_count); + const struct Graph *graph = data->graph; + + fprintf(file, "%d %d\n", graph->node_count, data->cluster_count); - for (int i = 0; i < data->graph->node_count; i++) + for (int i = 0; i < graph->node_count; i++) { - fprintf(file, "%.2lf %.2lf %d\n", data->x_coordinates[i], - data->y_coordinates[i], data->clusters[i]); + fprintf(file, "%.2lf %.2lf %d\n", graph->x_coordinates[i], + graph->y_coordinates[i], data->clusters[i]); } CLEANUP: @@ -314,23 +318,24 @@ int GTSP_read_solution(char *filename, double **p_x) for (int i = 0; i < edge_count; i++) { int from, to, edge; - rval = fscanf(file, "%d %d", &from, &to); - abort_if(rval != 2, "invalid input format (edge endpoints)"); + double val; + rval = fscanf(file, "%d %d %lf", &from, &to, &val); + abort_if(rval != 3, "invalid input format (edge endpoints)"); if (from > to) swap(from, to); edge = get_edge_num(node_count, from, to); abort_if(edge > num_cols, "invalid edge"); - x[from] += 0.5; - x[to] += 0.5; - x[edge] = 1; + x[from] += val/2; + x[to] += val/2; + x[edge] = val; } for (int i = 0; i < num_cols; i++) { if (x[i] <= LP_EPSILON) continue; - log_debug(" x%-3d = %.2f\n", i, x[i]); + log_debug(" x%-5d = %.2f\n", i, x[i]); } *p_x = x; @@ -340,13 +345,13 @@ int GTSP_read_solution(char *filename, double **p_x) return rval; } -static const struct option options_tab[] = {{"help", no_argument, 0, 'h'}, - {"nodes", required_argument, 0, 'n'}, - {"clusters", required_argument, 0, 'm'}, - {"grid-size", required_argument, 0, 'g'}, - {"optimal", required_argument, 0, 'x'}, - {"seed", required_argument, 0, 's'}, - {(char *) 0, (int) 0, (int *) 0, (int) 0}}; +static const struct option options_tab[] = {{"help", no_argument, 0, 'h'}, + {"nodes", required_argument, 0, 'n'}, + {"clusters", required_argument, 0, 'm'}, + {"grid-size", required_argument, 0, 'g'}, + {"optimal", required_argument, 0, 'x'}, + {"seed", required_argument, 0, 's'}, + {(char *) 0, (int) 0, (int *) 0, (int) 0}}; static int input_node_count = -1; static int input_cluster_count = -1; @@ -359,9 +364,9 @@ static void GTSP_print_usage() printf("%4s %-13s %s\n", "-m", "--clusters", "number of clusters"); printf("%4s %-13s %s\n", "-s", "--seed", "random seed"); printf("%4s %-13s %s\n", "-g", "--grid-size", - "size of the box used for generating random points"); + "size of the box used for generating random points"); printf("%4s %-13s %s\n", "-x", "--optimal", - "file containg valid solution (used to assert validity of cuts)"); + "file containg valid solution (used to assert validity of cuts)"); } static int GTSP_parse_args(int argc, char **argv) @@ -485,7 +490,7 @@ int GTSP_main(int argc, char **argv) log_info(" grid_size = %d\n", grid_size); rval = GTSP_create_random_problem(input_node_count, input_cluster_count, - grid_size, &data); + grid_size, &data); abort_if(rval, "GTSP_create_random_problem failed"); log_info("Writing random instance to file gtsp.in\n"); diff --git a/src/gtsp.h b/src/gtsp.h index 96397ab..2d3d049 100644 --- a/src/gtsp.h +++ b/src/gtsp.h @@ -1,7 +1,3 @@ -// -// Created by isoron on 18/03/15. -// - #ifndef _PROJECT_GTSP_H_ #define _PROJECT_GTSP_H_ @@ -14,9 +10,6 @@ struct GTSP int *clusters; int cluster_count; - - double *x_coordinates; - double *y_coordinates; }; int GTSP_create_random_problem( diff --git a/src/lp.c b/src/lp.c index f2cd0dd..2d05148 100644 --- a/src/lp.c +++ b/src/lp.c @@ -165,10 +165,11 @@ int LP_optimize(struct LP *lp, int *infeasible) *infeasible = 0; + #if LOG_LEVEL >= LOG_LEVEL_DEBUG 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); + #endif double time_before = get_current_time(); rval = CPXdualopt(lp->cplex_env, lp->cplex_lp); @@ -298,7 +299,7 @@ int LP_remove_old_cuts(struct LP *lp) if (count > 0) { - log_info("Found and removed %d old cuts\n", count); + log_debug("Found and removed %d old cuts\n", count); rval = CPXdualopt(lp->cplex_env, lp->cplex_lp); abort_if(rval, "CPXoptimize failed"); } @@ -364,6 +365,7 @@ int LP_write(struct LP *lp, const char *fname) int compare_cuts(struct Row *cut1, struct Row *cut2) { return_if_neq(cut1->nz, cut2->nz); + assert(cut1->nz == cut2->nz); for (int i = 0; i < cut1->nz; i++) { diff --git a/src/params.h b/src/params.h new file mode 100644 index 0000000..31fb664 --- /dev/null +++ b/src/params.h @@ -0,0 +1,10 @@ +#ifndef PROJECT_PARAMS_H +#define PROJECT_PARAMS_H + +#define LOG_LEVEL LOG_LEVEL_INFO + +#define ENABLE_COMB_INEQUALITIES + +//#define ALLOW_FRACTIONAL_SOLUTIONS + +#endif //PROJECT_PARAMS_H diff --git a/src/util.h b/src/util.h index 61552cb..ff34696 100644 --- a/src/util.h +++ b/src/util.h @@ -2,6 +2,7 @@ #define _PROJECT_UTIL_H_ #include +#include "params.h" #define LOG_LEVEL_ERROR 10 #define LOG_LEVEL_WARNING 20 @@ -9,8 +10,6 @@ #define LOG_LEVEL_DEBUG 40 #define LOG_LEVEL_VERBOSE 50 -#define LOG_LEVEL LOG_LEVEL_INFO - #if LOG_LEVEL < LOG_LEVEL_VERBOSE #define log_verbose(...) #else