From 6e4d69e289704c17f3ba98234f4d58b3015e620b Mon Sep 17 00:00:00 2001 From: Alinson Xavier Date: Thu, 26 Mar 2015 18:21:28 -0400 Subject: [PATCH] Split gtsp.c --- src/flow.c | 5 + src/gtsp-subtour.c | 513 ++++++++++++++++++++++++++++++++++++++++++++ src/gtsp-subtour.h | 38 ++++ src/gtsp.c | 515 +-------------------------------------------- src/gtsp.h | 5 + src/lp.c | 13 +- src/tsp.c | 27 ++- src/tsp.h | 2 +- 8 files changed, 599 insertions(+), 519 deletions(-) create mode 100644 src/gtsp-subtour.c create mode 100644 src/gtsp-subtour.h diff --git a/src/flow.c b/src/flow.c index ac802de..0faf6e4 100644 --- a/src/flow.c +++ b/src/flow.c @@ -54,6 +54,8 @@ int flow_mark_reachable_nodes( return rval; } +extern double FLOW_CPU_TIME; + int flow_find_max_flow( const struct Graph *digraph, const double *capacities, @@ -65,6 +67,7 @@ int flow_find_max_flow( int rval = 0; FLOW_MAX_FLOW_COUNT++; + double initial_time = get_current_time(); for (int i = 0; i < digraph->node_count; i++) digraph->nodes[i].mark = 0; @@ -154,6 +157,8 @@ int flow_find_max_flow( rval = flow_mark_reachable_nodes(digraph, residual_caps, from); abort_if(rval, "flow_mark_reachable_nodes failed"); + FLOW_CPU_TIME += get_current_time() - initial_time; + CLEANUP: if (path_edges) free(path_edges); if (residual_caps) free(residual_caps); diff --git a/src/gtsp-subtour.c b/src/gtsp-subtour.c new file mode 100644 index 0000000..0fa4c6f --- /dev/null +++ b/src/gtsp-subtour.c @@ -0,0 +1,513 @@ +#include "gtsp-subtour.h" +#include +#include +#include "util.h" +#include "flow.h" + +extern double FLOW_CPU_TIME; + +int static build_flow_digraph( + struct GTSP *data, double *x, struct Graph *digraph, double *capacities) +{ + int rval = 0; + + int *digraph_edges = 0; + + int node_count = data->graph->node_count; + struct Graph *graph = data->graph; + + 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; + + digraph_edges = (int *) malloc(2 * digraph_edge_count * sizeof(int)); + abort_if(!digraph_edges, "could not allocate digraph_edges"); + + // Create four directed edges for each edge of the original graph + int ke = 0; + int kc = 0; + for (int i = 0; i < graph->edge_count; i++) + { + if (x[node_count + i] < LP_EPSILON) continue; + + struct Edge *e = &graph->edges[i]; + int from = e->from->index; + int to = e->to->index; + + digraph_edges[ke++] = from; + digraph_edges[ke++] = to; + capacities[kc++] = x[node_count + i]; + + digraph_edges[ke++] = to; + digraph_edges[ke++] = from; + capacities[kc++] = 0; + + digraph_edges[ke++] = to; + digraph_edges[ke++] = from; + capacities[kc++] = x[node_count + i]; + + digraph_edges[ke++] = from; + digraph_edges[ke++] = to; + capacities[kc++] = 0; + } + + // Create an extra node for each cluster and connect it to the vertices + // of the cluster through some edge with very high capacity + for (int i = 0; i < node_count; i++) + { + struct Node *n = &graph->nodes[i]; + int cl = data->clusters[n->index]; + + digraph_edges[ke++] = n->index; + digraph_edges[ke++] = node_count + cl; + capacities[kc++] = 0; + + digraph_edges[ke++] = node_count + cl; + digraph_edges[ke++] = n->index; + capacities[kc++] = 0; + } + + // Create an extra root node and connect it to each cluster node through + // some edge with zero capacity + for (int i = 0; i < data->cluster_count; i++) + { + digraph_edges[ke++] = node_count + i; + digraph_edges[ke++] = node_count + data->cluster_count; + capacities[kc++] = 0; + + digraph_edges[ke++] = node_count + data->cluster_count; + digraph_edges[ke++] = node_count + i; + capacities[kc++] = 0; + } + + assert(ke <= 2 * digraph_edge_count); + assert(kc <= digraph_edge_count); + + digraph_edge_count = kc; + + rval = graph_build(digraph_node_count, kc, digraph_edges, 1, digraph); + abort_if(rval, "graph_build failed"); + + for (int i = 0; i < digraph_edge_count; i += 2) + { + digraph->edges[i].reverse = &digraph->edges[i + 1]; + digraph->edges[i + 1].reverse = &digraph->edges[i]; + } + + CLEANUP: + if (digraph_edges) free(digraph_edges); + return rval; +} + +int static add_subtour_cut( + struct LP *lp, + struct Graph *graph, + struct Node *from, + struct Node *to, + struct Edge **cut_edges, + int cut_edges_count, + int type) +{ + int rval = 0; + + char sense = 'G'; + double rhs = 2.0 - 2.0 * type; + int newnz = cut_edges_count + type; + + int rmatbeg = 0; + int *rmatind = 0; + double *rmatval = 0; + + rmatind = (int *) malloc(newnz * sizeof(int)); + abort_if(!rmatind, "could not allocate rmatind"); + + rmatval = (double *) malloc(newnz * sizeof(double)); + abort_if(!rmatval, "could not allocate rmatval"); + + for (int i = 0; i < cut_edges_count; i++) + { + rmatind[i] = cut_edges[i]->index + graph->node_count; + rmatval[i] = 1.0; + } + + if (type >= 1) + { + rmatind[cut_edges_count] = from->index; + rmatval[cut_edges_count] = -2.0; + } + + if (type >= 2) + { + rmatind[cut_edges_count + 1] = to->index; + rmatval[cut_edges_count + 1] = -2.0; + } + + 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(" %c %.2lf\n", sense, rhs); + + if (OPTIMAL_X) + { + double sum = 0; + for (int i = 0; i < newnz; i++) + sum += rmatval[i] * OPTIMAL_X[rmatind[i]]; + abort_if(sum <= rhs - LP_EPSILON, "cannot add invalid cut"); + } + + rval = LP_add_rows(lp, 1, newnz, &rhs, &sense, &rmatbeg, rmatind, rmatval); + abort_if(rval, "LP_add_rows failed"); + + CLEANUP: + if (rmatval) free(rmatval); + if (rmatind) free(rmatind); + return rval; +} + +int find_exact_subtour_cuts( + struct LP *lp, struct GTSP *data, int *total_added_cuts) +{ + int rval = 0; + + double *x = 0; + double *capacities = 0; + + int added_cuts_count = 0; + struct Graph *graph = data->graph; + + 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"); + +#if LOG_LEVEL >= LOG_LEVEL_DEBUG + rval = GTSP_write_solution(data, "gtsp-frac.out", x); + abort_if(rval, "GTSP_write_solution failed"); +#endif + + struct Graph digraph; + graph_init(&digraph); + + int digraph_edge_count = 4 * graph->edge_count + 2 * graph->node_count + + 2 * data->cluster_count; + + capacities = (double *) malloc(digraph_edge_count * sizeof(double)); + abort_if(!capacities, "could not allocate capacities"); + + rval = build_flow_digraph(data, x, &digraph, capacities); + abort_if(rval, "build_flow_digraph failed"); + + // Constraints (2.1) + rval = find_exact_subtour_cuts_cluster_to_cluster(lp, data, &digraph, + capacities, &added_cuts_count); + abort_if(rval, "find_exact_subtour_cuts_cluster_to_cluster failed"); + + log_debug("Added %d cluster-to-cluster subtour cuts\n", added_cuts_count); + (*total_added_cuts) += added_cuts_count; + if (added_cuts_count > 0) goto CLEANUP; + + // Constraints (2.2) + rval = find_exact_subtour_cuts_node_to_cluster(lp, data, x, &digraph, + capacities, &added_cuts_count); + abort_if(rval, "find_exact_subtour_cuts_node_to_cluster failed"); + + log_debug("Added %d node-to-cluster subtour cuts\n", added_cuts_count); + (*total_added_cuts) += added_cuts_count; + if (added_cuts_count > 0) goto CLEANUP; + + // Constraints (2.3) + rval = find_exact_subtour_cuts_node_to_node(lp, data, x, &digraph, + capacities, &added_cuts_count); + abort_if(rval, "find_exact_subtour_cuts_node_to_node failed"); + + log_debug("Added %d node-to-node subtour cuts\n", added_cuts_count); + (*total_added_cuts) += added_cuts_count; + + CLEANUP: + graph_free(&digraph); + if (capacities) free(capacities); + if (x) free(x); + return rval; +} + +int find_exact_subtour_cuts_node_to_node( + struct LP *lp, + struct GTSP *data, + double *x, + struct Graph *digraph, + double *capacities, + int *added_cuts_count) +{ + int rval = 0; + + struct Edge **cut_edges = 0; + double *flow = 0; + + struct Graph *graph = data->graph; + int *clusters = data->clusters; + + cut_edges = (struct Edge **) malloc( + graph->edge_count * sizeof(struct Edge *)); + flow = (double *) malloc(digraph->edge_count * sizeof(double)); + + abort_if(!cut_edges, "could not allocate cut_edges"); + abort_if(!flow, "could not allocate flow"); + + int max_x_index = 0; + double max_x = DBL_MIN; + + for (int i = 0; i < graph->node_count; i++) + { + struct Node *n = &graph->nodes[i]; + if (x[n->index] > max_x) + { + max_x = x[n->index]; + max_x_index = i; + } + } + + int i = max_x_index; + + for (int j = 0; j < graph->node_count; j++) + { + if (i == j) continue; + if (clusters[i] == clusters[j]) continue; + if (x[i] + x[j] - 1 <= LP_EPSILON) continue; + + struct Node *from = &digraph->nodes[i]; + struct Node *to = &digraph->nodes[j]; + + int cut_edges_count; + double flow_value; + + rval = flow_find_max_flow(digraph, capacities, from, to, flow, + &flow_value); + abort_if(rval, "flow_find_max_flow failed"); + + if (flow_value >= 2 * (x[i] + x[j] - 1) - MIN_CUT_VIOLATION) + continue; + + log_verbose("Marked nodes:\n"); + for (int k = 0; k < graph->node_count; k++) + { + graph->nodes[k].mark = digraph->nodes[k].mark; + if (digraph->nodes[k].mark) log_verbose(" %d\n", k); + } + + rval = get_cut_edges_from_marks(graph, &cut_edges_count, cut_edges); + abort_if(rval, "get_cut_edges_from_marks failed"); + + 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); + + rval = add_subtour_cut(lp, graph, from, to, cut_edges, cut_edges_count, + 2); + abort_if(rval, "add_subtour_cut failed"); + + (*added_cuts_count)++; + } + + CLEANUP: + if (flow) free(flow); + if (cut_edges) free(cut_edges); + return rval; +} + +int find_exact_subtour_cuts_node_to_cluster( + struct LP *lp, + struct GTSP *data, + double *x, + struct Graph *digraph, + double *capacities, + int *added_cuts_count) +{ + int rval = 0; + + int cuts_count = 0; + struct Edge **cut_edges = 0; + double *flow = 0; + + struct Graph *graph = data->graph; + int *clusters = data->clusters; + + cut_edges = (struct Edge **) malloc( + graph->edge_count * sizeof(struct Edge *)); + flow = (double *) malloc(digraph->edge_count * sizeof(double)); + abort_if(!cut_edges, "could not allocate cut_edges"); + abort_if(!flow, "could not allocate flow"); + + for (int i = 0; i < graph->node_count; i++) + { + for (int j = 0; j < data->cluster_count; j++) + { + if (clusters[i] == j) continue; + if (x[i] < LP_EPSILON) continue; + + struct Node *from = &digraph->nodes[i]; + struct Node *to = &digraph->nodes[graph->node_count + j]; + + log_verbose( + "Sending flow from node %d to cluster %d (must be >= %.4lf)\n", + i, j, 2 * x[i]); + + activate_cluster_node(capacities, to); + + double flow_value; + int cut_edges_count; + + rval = flow_find_max_flow(digraph, capacities, from, to, flow, + &flow_value); + abort_if(rval, "flow_find_max_flow failed"); + + log_verbose(" flow value = %.4lf\n", flow_value); + + deactivate_cluster_node(capacities, to); + + if (flow_value + MIN_CUT_VIOLATION >= 2 * x[i]) + continue; + + log_verbose("Marked nodes:\n"); + for (int k = 0; k < graph->node_count; k++) + { + graph->nodes[k].mark = digraph->nodes[k].mark; + if (graph->nodes[k].mark) log_verbose(" %d\n", k); + } + + rval = get_cut_edges_from_marks(graph, &cut_edges_count, cut_edges); + abort_if(rval, "get_cut_edges_from_marks failed"); + + log_verbose("Cut edges:\n"); + for (int k = 0; k < cut_edges_count; k++) + { + 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); + } + + rval = add_subtour_cut(lp, graph, from, 0, cut_edges, + cut_edges_count, 1); + abort_if(rval, "add_subtour_cut failed"); + + (*added_cuts_count)++; + cuts_count++; + } + } + + CLEANUP: + if (cut_edges) free(cut_edges); + if (flow) free(flow); + return rval; +} + +int find_exact_subtour_cuts_cluster_to_cluster( + struct LP *lp, + struct GTSP *data, + struct Graph *digraph, + double *capacities, + int *added_cuts_count) +{ + int rval = 0; + + double *flow = 0; + struct Edge **cut_edges = 0; + + int cuts_count = 0; + + struct Graph *graph = data->graph; + + cut_edges = (struct Edge **) malloc( + graph->edge_count * sizeof(struct Edge *)); + flow = (double *) malloc(digraph->edge_count * sizeof(double)); + abort_if(!cut_edges, "could not allocate cut_edges"); + abort_if(!flow, "could not allocate flow"); + + struct Node *root_node = &digraph->nodes[graph->node_count + + data->cluster_count]; + + for (int i = 0; i < data->cluster_count; i++) + { + for (int j = i + 1; j < data->cluster_count; j++) + { + struct Node *from = &digraph->nodes[graph->node_count + i]; + struct Node *to = &digraph->nodes[graph->node_count + j]; + + double flow_value; + int cut_edges_count; + + activate_cluster_node(capacities, from); + activate_cluster_node(capacities, to); + deactivate_cluster_node(capacities, root_node); + + 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); + + abort_if(rval, "flow_find_max_flow failed"); + + log_verbose(" flow value = %.4lf\n", flow_value); + + deactivate_cluster_node(capacities, from); + deactivate_cluster_node(capacities, to); + + if (flow_value >= 2 - MIN_CUT_VIOLATION) continue; + + log_verbose("Marked nodes:\n"); + for (int k = 0; k < graph->node_count; k++) + { + graph->nodes[k].mark = digraph->nodes[k].mark; + if (digraph->nodes[k].mark) log_verbose(" %d\n", k); + } + + rval = get_cut_edges_from_marks(graph, &cut_edges_count, cut_edges); + abort_if(rval, "get_cut_edges_from_marks failed"); + + 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); + + rval = add_subtour_cut(lp, graph, 0, 0, cut_edges, cut_edges_count, + 0); + abort_if(rval, "add_subtour_cut failed"); + + (*added_cuts_count)++; + cuts_count++; + } + +// if(cuts_count > 0) break; + } + + CLEANUP: + if (cut_edges) free(cut_edges); + if (flow) free(flow); + return rval; +} + +void deactivate_cluster_node(double *capacities, struct Node *cluster_node) +{ + for (int i = 0; i < cluster_node->degree; i++) + { + struct Adjacency *adj = &cluster_node->adj[i]; + struct Edge *e = adj->edge; + + capacities[e->index] = 0; + } +} + +void activate_cluster_node(double *capacities, struct Node *cluster_node) +{ + for (int i = 0; i < cluster_node->degree; i++) + { + struct Adjacency *adj = &cluster_node->adj[i]; + struct Edge *e = adj->edge; + + capacities[e->index] = 1e10; + capacities[e->reverse->index] = 1e10; + } +} \ No newline at end of file diff --git a/src/gtsp-subtour.h b/src/gtsp-subtour.h new file mode 100644 index 0000000..9f991a8 --- /dev/null +++ b/src/gtsp-subtour.h @@ -0,0 +1,38 @@ +#ifndef _PROJECT_GTSP_SUBTOUR_H_ +#define _PROJECT_GTSP_SUBTOUR_H_ + +#include "gtsp.h" + +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 find_exact_subtour_cuts_node_to_cluster( + struct LP *lp, + struct GTSP *data, + double *x, + struct Graph *digraph, + double *capacities, + int *added_cuts_count); + +int find_exact_subtour_cuts_node_to_node( + struct LP *lp, + struct GTSP *data, + double *x, + struct Graph *digraph, + double *capacities, + int *added_cuts_count); + +int find_exact_subtour_cuts( + struct LP *lp, struct GTSP *data, int *total_added_cuts); + + +#endif //_PROJECT_GTSP_SUBTOUR_H_ diff --git a/src/gtsp.c b/src/gtsp.c index 30c415d..65c9f21 100644 --- a/src/gtsp.c +++ b/src/gtsp.c @@ -1,7 +1,6 @@ #include #include #include -#include #include #include "gtsp.h" #include "geometry.h" @@ -9,8 +8,9 @@ #include "flow.h" #include "branch_and_cut.h" #include "math.h" +#include "gtsp-subtour.h" -static double *OPTIMAL_X = 0; +double *OPTIMAL_X = 0; static int get_edge_num(int node_count, int from, int to) { @@ -183,507 +183,6 @@ int GTSP_init_lp(struct LP *lp, struct GTSP *data) return rval; } -static int add_subtour_cut( - struct LP *lp, - struct Graph *graph, - struct Node *from, - struct Node *to, - struct Edge **cut_edges, - int cut_edges_count, - int type) -{ - int rval = 0; - - char sense = 'G'; - double rhs = 2.0 - 2.0 * type; - int newnz = cut_edges_count + type; - - int rmatbeg = 0; - int *rmatind = 0; - double *rmatval = 0; - - rmatind = (int *) malloc(newnz * sizeof(int)); - abort_if(!rmatind, "could not allocate rmatind"); - - rmatval = (double *) malloc(newnz * sizeof(double)); - abort_if(!rmatval, "could not allocate rmatval"); - - for (int i = 0; i < cut_edges_count; i++) - { - rmatind[i] = cut_edges[i]->index + graph->node_count; - rmatval[i] = 1.0; - } - - if (type >= 1) - { - rmatind[cut_edges_count] = from->index; - rmatval[cut_edges_count] = -2.0; - } - - if (type >= 2) - { - rmatind[cut_edges_count + 1] = to->index; - rmatval[cut_edges_count + 1] = -2.0; - } - - 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(" %c %.2lf\n", sense, rhs); - - if (OPTIMAL_X) - { - double sum = 0; - for (int i = 0; i < newnz; i++) - sum += rmatval[i] * OPTIMAL_X[rmatind[i]]; - abort_if(sum <= rhs - LP_EPSILON, "cannot add invalid cut"); - } - - rval = LP_add_rows(lp, 1, newnz, &rhs, &sense, &rmatbeg, rmatind, rmatval); - abort_if(rval, "LP_add_rows failed"); - - CLEANUP: - if (rmatval) free(rmatval); - if (rmatind) free(rmatind); - return rval; -} - -static int build_flow_digraph( - struct GTSP *data, double *x, struct Graph *digraph, double *capacities) -{ - int rval = 0; - - int *digraph_edges = 0; - - int node_count = data->graph->node_count; - struct Graph *graph = data->graph; - - 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; - - digraph_edges = (int *) malloc(2 * digraph_edge_count * sizeof(int)); - abort_if(!digraph_edges, "could not allocate digraph_edges"); - - // Create four directed edges for each edge of the original graph - int ke = 0; - int kc = 0; - for (int i = 0; i < graph->edge_count; i++) - { -// if (x[node_count + i] < LP_EPSILON) continue; - - struct Edge *e = &graph->edges[i]; - int from = e->from->index; - int to = e->to->index; - - digraph_edges[ke++] = from; - digraph_edges[ke++] = to; - capacities[kc++] = x[node_count + i]; - - digraph_edges[ke++] = to; - digraph_edges[ke++] = from; - capacities[kc++] = 0; - - digraph_edges[ke++] = to; - digraph_edges[ke++] = from; - capacities[kc++] = x[node_count + i]; - - digraph_edges[ke++] = from; - digraph_edges[ke++] = to; - capacities[kc++] = 0; - } - - // Create an extra node for each cluster and connect it to the vertices - // of the cluster through some edge with very high capacity - for (int i = 0; i < node_count; i++) - { - struct Node *n = &graph->nodes[i]; - int cl = data->clusters[n->index]; - - digraph_edges[ke++] = n->index; - digraph_edges[ke++] = node_count + cl; - capacities[kc++] = 0; - - digraph_edges[ke++] = node_count + cl; - digraph_edges[ke++] = n->index; - capacities[kc++] = 0; - } - - // Create an extra root node and connect it to each cluster node through - // some edge with zero capacity - for (int i = 0; i < data->cluster_count; i++) - { - digraph_edges[ke++] = node_count + i; - digraph_edges[ke++] = node_count + data->cluster_count; - capacities[kc++] = 0; - - digraph_edges[ke++] = node_count + data->cluster_count; - digraph_edges[ke++] = node_count + i; - capacities[kc++] = 0; - } - - assert(ke <= 2 * digraph_edge_count); - assert(kc <= digraph_edge_count); - -// digraph_edge_count = kc; - - rval = graph_build(digraph_node_count, kc, digraph_edges, 1, digraph); - abort_if(rval, "graph_build failed"); - - for (int i = 0; i < digraph_edge_count; i += 2) - { - digraph->edges[i].reverse = &digraph->edges[i + 1]; - digraph->edges[i + 1].reverse = &digraph->edges[i]; - } - - CLEANUP: - if (digraph_edges) free(digraph_edges); - return rval; -} - -void deactivate_cluster_node(double *capacities, struct Node *cluster_node) -{ - for (int i = 0; i < cluster_node->degree; i++) - { - struct Adjacency *adj = &cluster_node->adj[i]; - struct Edge *e = adj->edge; - - capacities[e->index] = 0; - } -} - -void activate_cluster_node(double *capacities, struct Node *cluster_node) -{ - for (int i = 0; i < cluster_node->degree; i++) - { - struct Adjacency *adj = &cluster_node->adj[i]; - struct Edge *e = adj->edge; - - capacities[e->index] = 1e10; - capacities[e->reverse->index] = 1e10; - } -} - -int find_exact_subtour_cuts_cluster_to_cluster( - struct LP *lp, - struct GTSP *data, - struct Graph *digraph, - double *capacities, - int *added_cuts_count) -{ - int rval = 0; - - double *flow = 0; - struct Edge **cut_edges = 0; - - int cuts_count = 0; - - struct Graph *graph = data->graph; - - cut_edges = (struct Edge **) malloc( - graph->edge_count * sizeof(struct Edge *)); - flow = (double *) malloc(digraph->edge_count * sizeof(double)); - abort_if(!cut_edges, "could not allocate cut_edges"); - abort_if(!flow, "could not allocate flow"); - - struct Node *root_node = &digraph->nodes[graph->node_count + - data->cluster_count]; - - for (int i = 0; i < data->cluster_count; i++) - { - for (int j = i + 1; j < data->cluster_count; j++) - { - struct Node *from = &digraph->nodes[graph->node_count + i]; - struct Node *to = &digraph->nodes[graph->node_count + j]; - - double flow_value; - int cut_edges_count; - - activate_cluster_node(capacities, from); - activate_cluster_node(capacities, to); - deactivate_cluster_node(capacities, root_node); - - 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); - abort_if(rval, "flow_find_max_flow failed"); - - log_verbose(" flow value = %.4lf\n", flow_value); - - deactivate_cluster_node(capacities, from); - deactivate_cluster_node(capacities, to); - - if (flow_value >= 2 - LP_EPSILON) continue; - - log_verbose("Marked nodes:\n"); - for (int k = 0; k < graph->node_count; k++) - { - graph->nodes[k].mark = digraph->nodes[k].mark; - if (digraph->nodes[k].mark) log_verbose(" %d\n", k); - } - - rval = get_cut_edges_from_marks(graph, &cut_edges_count, cut_edges); - abort_if(rval, "get_cut_edges_from_marks failed"); - - 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); - - rval = add_subtour_cut(lp, graph, 0, 0, cut_edges, cut_edges_count, - 0); - abort_if(rval, "add_subtour_cut failed"); - - (*added_cuts_count)++; - cuts_count++; - } - } - - CLEANUP: - if (cut_edges) free(cut_edges); - if (flow) free(flow); - return rval; -} - -int find_exact_subtour_cuts_node_to_cluster( - struct LP *lp, - struct GTSP *data, - double *x, - struct Graph *digraph, - double *capacities, - int *added_cuts_count) -{ - int rval = 0; - - int cuts_count = 0; - struct Edge **cut_edges = 0; - double *flow = 0; - - struct Graph *graph = data->graph; - int *clusters = data->clusters; - - cut_edges = (struct Edge **) malloc( - digraph->edge_count * sizeof(struct Edge *)); - flow = (double *) malloc(digraph->edge_count * sizeof(double)); - abort_if(!cut_edges, "could not allocate cut_edges"); - abort_if(!flow, "could not allocate flow"); - - for (int i = 0; i < graph->node_count; i++) - { - for (int j = 0; j < data->cluster_count; j++) - { - if (clusters[i] == j) continue; - if (x[i] < LP_EPSILON) continue; - - struct Node *from = &digraph->nodes[i]; - struct Node *to = &digraph->nodes[graph->node_count + j]; - - log_verbose("Sending flow from node %d to cluster %d (must be >= %.4lf)\n", i, j, 2*x[i]); - - activate_cluster_node(capacities, to); - - double flow_value; - int cut_edges_count; - - rval = flow_find_max_flow(digraph, capacities, from, to, flow, - &flow_value); - abort_if(rval, "flow_find_max_flow failed"); - - log_verbose(" flow value = %.4lf\n", flow_value); - - deactivate_cluster_node(capacities, to); - - if (flow_value + LP_EPSILON >= 2 * x[i]) - continue; - - log_verbose("Marked nodes:\n"); - for (int k = 0; k < graph->node_count; k++) - { - graph->nodes[k].mark = digraph->nodes[k].mark; - if (graph->nodes[k].mark) log_verbose(" %d\n", k); - } - - rval = get_cut_edges_from_marks(graph, &cut_edges_count, cut_edges); - abort_if(rval, "get_cut_edges_from_marks failed"); - - log_verbose("Cut edges:\n"); - for (int k = 0; k < cut_edges_count; k++) - { - 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); - } - - rval = add_subtour_cut(lp, graph, from, 0, cut_edges, - cut_edges_count, 1); - abort_if(rval, "add_subtour_cut failed"); - - (*added_cuts_count)++; - cuts_count++; - } - } - - CLEANUP: - if (cut_edges) free(cut_edges); - if (flow) free(flow); - return rval; -} - -int find_exact_subtour_cuts_node_to_node( - struct LP *lp, - struct GTSP *data, - double *x, - struct Graph *digraph, - double *capacities, - int *added_cuts_count) -{ - int rval = 0; - - struct Edge **cut_edges = 0; - double *flow = 0; - - struct Graph *graph = data->graph; - int *clusters = data->clusters; - - cut_edges = (struct Edge **) malloc( - digraph->edge_count * sizeof(struct Edge *)); - flow = (double *) malloc(digraph->edge_count * sizeof(double)); - - abort_if(!cut_edges, "could not allocate cut_edges"); - abort_if(!flow, "could not allocate flow"); - - int max_x_index = 0; - double max_x = DBL_MIN; - - for (int i = 0; i < graph->node_count; i++) - { - struct Node *n = &graph->nodes[i]; - if (x[n->index] > max_x) - { - max_x = x[n->index]; - max_x_index = i; - } - } - - int i = max_x_index; - - for (int j = 0; j < graph->node_count; j++) - { - if (i == j) continue; - if (clusters[i] == clusters[j]) continue; - if (x[i] + x[j] - 1 <= LP_EPSILON) continue; - - struct Node *from = &digraph->nodes[i]; - struct Node *to = &digraph->nodes[j]; - - int cut_edges_count; - double flow_value; - - rval = flow_find_max_flow(digraph, capacities, from, to, flow, - &flow_value); - abort_if(rval, "flow_find_max_flow failed"); - - if (flow_value >= 2 * (x[i] + x[j] - 1) - LP_EPSILON) - continue; - - log_verbose("Marked nodes:\n"); - for (int k = 0; k < graph->node_count; k++) - { - graph->nodes[k].mark = digraph->nodes[k].mark; - if (digraph->nodes[k].mark) log_verbose(" %d\n", k); - } - - rval = get_cut_edges_from_marks(graph, &cut_edges_count, cut_edges); - abort_if(rval, "get_cut_edges_from_marks failed"); - - 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); - - rval = add_subtour_cut(lp, graph, from, to, cut_edges, cut_edges_count, - 2); - abort_if(rval, "add_subtour_cut failed"); - - (*added_cuts_count)++; - } - - CLEANUP: - if (flow) free(flow); - if (cut_edges) free(cut_edges); - return rval; -} - -int find_exact_subtour_cuts( - struct LP *lp, struct GTSP *data, int *total_added_cuts) -{ - int rval = 0; - - double *x = 0; - double *capacities = 0; - - int added_cuts_count = 0; - struct Graph *graph = data->graph; - - 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"); - -#if LOG_LEVEL >= LOG_LEVEL_DEBUG - rval = GTSP_write_solution(data, "gtsp-frac.out", x); - abort_if(rval, "GTSP_write_solution failed"); -#endif - - struct Graph digraph; - graph_init(&digraph); - - int digraph_edge_count = 4 * graph->edge_count + 2 * graph->node_count + - 2 * data->cluster_count; - - capacities = (double *) malloc(digraph_edge_count * sizeof(double)); - abort_if(!capacities, "could not allocate capacities"); - - rval = build_flow_digraph(data, x, &digraph, capacities); - abort_if(rval, "build_flow_digraph failed"); - - // Constraints (2.1) - rval = find_exact_subtour_cuts_cluster_to_cluster(lp, data, &digraph, - capacities, &added_cuts_count); - abort_if(rval, "find_exact_subtour_cuts_cluster_to_cluster failed"); - - log_debug("Added %d cluster-to-cluster subtour cuts\n", added_cuts_count); - (*total_added_cuts) += added_cuts_count; - - if (added_cuts_count > 0) goto CLEANUP; - - // Constraints (2.2) - rval = find_exact_subtour_cuts_node_to_cluster(lp, data, x, &digraph, - capacities, &added_cuts_count); - abort_if(rval, "find_exact_subtour_cuts_node_to_cluster failed"); - - log_debug("Added %d node-to-cluster subtour cuts\n", added_cuts_count); - (*total_added_cuts) += added_cuts_count; - - // Constraints (2.3) - rval = find_exact_subtour_cuts_node_to_node(lp, data, x, &digraph, - capacities, &added_cuts_count); - abort_if(rval, "find_exact_subtour_cuts_node_to_node failed"); - - log_debug("Added %d node-to-node subtour cuts\n", added_cuts_count); - (*total_added_cuts) += added_cuts_count; - - CLEANUP: - graph_free(&digraph); - if (capacities) free(capacities); - if (x) free(x); - return rval; -} - int GTSP_add_cutting_planes(struct LP *lp, struct GTSP *data) { int rval = 0; @@ -952,6 +451,11 @@ 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; + int GTSP_main(int argc, char **argv) { int rval = 0; @@ -1026,7 +530,10 @@ int GTSP_main(int argc, char **argv) log_info(" obj value = %.2lf **\n", bnc.best_obj_val); log_info("Branch-and-bound nodes: %d\n", BNC_NODE_COUNT); - log_info("Max-flow computations: %d\n", FLOW_MAX_FLOW_COUNT); + 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); CLEANUP: GTSP_free(&data); diff --git a/src/gtsp.h b/src/gtsp.h index e73f898..b76f653 100644 --- a/src/gtsp.h +++ b/src/gtsp.h @@ -19,6 +19,8 @@ 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); @@ -36,4 +38,7 @@ int GTSP_write_solution(struct GTSP *data, char *filename, double *x); int GTSP_main(int argc, char **argv); +extern double *OPTIMAL_X; +extern double FLOW_CPU_TIME; + #endif //_PROJECT_GTSP_H_ diff --git a/src/lp.c b/src/lp.c index ea9dfca..50162db 100644 --- a/src/lp.c +++ b/src/lp.c @@ -119,25 +119,32 @@ int LP_change_bound(struct LP *lp, int col, char lower_or_upper, double bnd) return rval; } +extern int LP_OPTIMIZE_COUNT; +extern double LP_CPU_TIME; + int LP_optimize(struct LP *lp, int *infeasible) { + LP_OPTIMIZE_COUNT++; + int rval = 0, solstat; *infeasible = 0; + double current = get_current_time(); rval = CPXdualopt(lp->cplex_env, lp->cplex_lp); abort_if(rval, "CPXdualopt failed"); + LP_CPU_TIME += get_current_time() - current; solstat = CPXgetstat(lp->cplex_env, lp->cplex_lp); if (solstat == CPX_STAT_INFEASIBLE) { *infeasible = 1; + goto CLEANUP; } else { - abort_if(solstat != CPX_STAT_OPTIMAL - && solstat != CPX_STAT_OPTIMAL_INFEAS, - "Invalid solution status"); + abort_if(solstat != CPX_STAT_OPTIMAL && + solstat != CPX_STAT_OPTIMAL_INFEAS, "Invalid solution status"); } CLEANUP: diff --git a/src/tsp.c b/src/tsp.c index 0dddb44..258d50b 100644 --- a/src/tsp.c +++ b/src/tsp.c @@ -22,7 +22,7 @@ int TSP_init_data(struct TSPData *data) void TSP_free_data(struct TSPData *data) { - if(!data) return; + if (!data) return; if (data->edge_list) free(data->edge_list); if (data->edge_weights) free(data->edge_weights); } @@ -64,8 +64,7 @@ int TSP_init_lp(struct LP *lp, struct TSPData *data) } int TSP_find_violated_subtour_elimination_cut( - struct LP *lp, - struct TSPData *data) + struct LP *lp, struct TSPData *data) { int ncount = data->node_count; int edge_count = data->edge_count; @@ -111,16 +110,21 @@ int TSP_find_violated_subtour_elimination_cut( rval = LP_get_x(lp, x); abort_if(rval, "LP_get_x failed"); - #if LOG_LEVEL >= LOG_LEVEL_VERBOSE +#if LOG_LEVEL >= LOG_LEVEL_VERBOSE int current_round = 0; #endif int delta_count = 0; int island_count = 0; - while (!TSP_is_graph_connected(&G, x, &island_count, island_sizes, - island_start, island_nodes)) + while (1) { + rval = graph_find_connected_components(&G, x, &island_count, + island_sizes, island_start, island_nodes); + abort_if(rval, "graph_find_connected_components failed"); + + if(island_count == 1) break; + log_verbose("Adding %d subtour inequalities...\n", island_count); for (int i = 0; i < island_count; i++) { @@ -136,7 +140,7 @@ int TSP_find_violated_subtour_elimination_cut( rval = LP_optimize(lp, &is_infeasible); abort_if(rval, "LP_optimize failed"); - if(is_infeasible) goto CLEANUP; + if (is_infeasible) goto CLEANUP; double objval = 0; rval = LP_get_obj_val(lp, &objval); @@ -184,7 +188,7 @@ int TSP_add_subtour_elimination_cut(struct LP *lp, int delta_length, int *delta) return rval; } -int TSP_is_graph_connected( +int graph_find_connected_components( struct Graph *G, double *x, int *island_count, @@ -216,7 +220,7 @@ int TSP_is_graph_connected( (*island_count) = current_island; - return (*island_count == 1); + return 0; } int TSP_find_closest_neighbor_tour( @@ -491,8 +495,9 @@ int TSP_main(int argc, char **argv) bnc.best_obj_val = TSP_find_initial_solution(&data); bnc.problem_data = (void *) &data; bnc.problem_init_lp = (int (*)(struct LP *, void *)) TSP_init_lp; - bnc.problem_add_cutting_planes = - (int (*)(struct LP *, void *)) TSP_add_cutting_planes; + bnc.problem_add_cutting_planes = (int (*)( + struct LP *, + void *)) TSP_add_cutting_planes; rval = BNC_init_lp(&bnc); abort_if(rval, "BNC_init_lp failed"); diff --git a/src/tsp.h b/src/tsp.h index 3361669..21aaef0 100644 --- a/src/tsp.h +++ b/src/tsp.h @@ -19,7 +19,7 @@ void TSP_free_data(struct TSPData *data); int TSP_find_violated_subtour_elimination_cut (struct LP *lp, struct TSPData *data); -int TSP_is_graph_connected( +int graph_find_connected_components( struct Graph *G, double *x, int *island_count,