From 5fe42908b40954d80df9338625cb75e5c99892c5 Mon Sep 17 00:00:00 2001 From: Alinson Xavier Date: Tue, 24 Mar 2015 10:00:45 -0400 Subject: [PATCH] Fix subtour cuts, create separate main functions --- scripts/draw_solution.sage | 8 +- src/branch_and_cut.c | 21 +- src/branch_and_cut.h | 2 + src/flow.c | 194 ++++++++++++-- src/flow.h | 11 +- src/graph.c | 26 ++ src/graph.h | 2 + src/gtsp.c | 505 ++++++++++++++++++++++++++++++++++++- src/gtsp.h | 6 +- src/main.c | 287 ++------------------- src/tsp.c | 117 ++++++++- src/tsp.h | 2 + src/util.h | 11 +- 13 files changed, 884 insertions(+), 308 deletions(-) diff --git a/scripts/draw_solution.sage b/scripts/draw_solution.sage index b6147aa..b2ee462 100644 --- a/scripts/draw_solution.sage +++ b/scripts/draw_solution.sage @@ -56,9 +56,10 @@ for i in range(cluster_count): for i in range(node_count): (x,y,cluster) = [ int(float(x)) for x in raw_input().split(' ') ] points[cluster].append((x,y)) - all_points.append((x,y)) + all_points.append(vector([x,y])) # solutions file +(node_count, edge_count) = [ int(x) for x in raw_input().split(' ') ] edges_count = int(raw_input()) edges = [] for i in range(edges_count): @@ -67,6 +68,11 @@ for i in range(edges_count): plot = list_plot([], xmax=100, xmin=0, ymax=100, ymin=0) +#plot = plot + sum([text(str(i), points[i]) for i in range(len(points))]) + +for i in range(node_count): + plot = plot + text(str(i), all_points[i] + vector([0,-2]), color='gray') + for i in range(cluster_count): plot = plot + list_plot(points[i], color='gray', figsize=FIGURE_SIZE, pointsize=POINT_SIZE) diff --git a/src/branch_and_cut.c b/src/branch_and_cut.c index 8b43c4a..751f276 100644 --- a/src/branch_and_cut.c +++ b/src/branch_and_cut.c @@ -45,7 +45,6 @@ void BNC_free(struct BNC *bnc) int BNC_init_lp(struct BNC *bnc) { int rval = 0; - log_verbose("Initializing LP...\n"); rval = LP_open(bnc->lp); abort_if(rval, "LP_open failed"); @@ -76,7 +75,7 @@ static int BNC_solve_node(struct BNC *bnc, int depth) int rval = 0; double *x = (double *) NULL; - log_verbose("Optimizing...\n"); + log_debug("Optimizing...\n"); int is_infeasible; rval = LP_optimize(lp, &is_infeasible); @@ -84,7 +83,7 @@ static int BNC_solve_node(struct BNC *bnc, int depth) if (is_infeasible) { - log_verbose("Branch pruned by infeasibility.\n"); + log_debug("Branch pruned by infeasibility.\n"); goto CLEANUP; } @@ -92,11 +91,11 @@ 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_verbose(" objective value = %.2f\n", objval); + log_debug(" objective value = %.2f\n", objval); if (objval > *best_val) { - log_verbose("Branch pruned by bound (%.2lf > %.2lf).\n", objval, + log_debug("Branch pruned by bound (%.2lf > %.2lf).\n", objval, *best_val); rval = 0; goto CLEANUP; @@ -127,7 +126,7 @@ static int BNC_solve_node(struct BNC *bnc, int depth) if (BNC_is_integral(x, num_cols)) { - log_verbose(" solution is integral\n"); + log_debug(" solution is integral\n"); if (objval < *best_val) { @@ -141,7 +140,7 @@ static int BNC_solve_node(struct BNC *bnc, int depth) } else { - log_verbose(" solution is fractional\n"); + log_debug(" solution is fractional\n"); rval = BNC_branch_node(bnc, x, depth); abort_if(rval, "BNC_branch_node failed"); } @@ -160,10 +159,10 @@ static int BNC_branch_node(struct BNC *bnc, double *x, int depth) int num_cols = LP_get_num_cols(lp); int best_branch_var = BNC_find_best_branching_var(x, num_cols); - log_verbose("Branching on variable x%d = %.6lf (depth %d)...\n", + log_debug("Branching on variable x%d = %.6lf (depth %d)...\n", best_branch_var, x[best_branch_var], depth); - log_verbose("Fixing variable x%d to one...\n", best_branch_var); + log_debug("Fixing variable x%d to one...\n", best_branch_var); rval = LP_change_bound(lp, best_branch_var, 'L', 1.0); abort_if(rval, "LP_change_bound failed"); @@ -173,7 +172,7 @@ static int BNC_branch_node(struct BNC *bnc, double *x, int depth) rval = LP_change_bound(lp, best_branch_var, 'L', 0.0); abort_if(rval, "LP_change_bound failed"); - log_verbose("Fixing variable x%d to zero...\n", best_branch_var); + log_debug("Fixing variable x%d to zero...\n", best_branch_var); rval = LP_change_bound(lp, best_branch_var, 'U', 0.0); abort_if(rval, "LP_change_bound failed"); @@ -183,7 +182,7 @@ static int BNC_branch_node(struct BNC *bnc, double *x, int depth) rval = LP_change_bound(lp, best_branch_var, 'U', 1.0); abort_if(rval, "LP_change_bound failed"); - log_verbose("Finished branching on variable %d\n", best_branch_var); + log_debug("Finished branching on variable %d\n", best_branch_var); CLEANUP: return rval; diff --git a/src/branch_and_cut.h b/src/branch_and_cut.h index f6e811e..2be9b00 100644 --- a/src/branch_and_cut.h +++ b/src/branch_and_cut.h @@ -10,6 +10,8 @@ struct BNC double *best_x; double best_obj_val; + double *optimal_x; + int *problem_data; int (*problem_init_lp)(struct LP *, void *); diff --git a/src/flow.c b/src/flow.c index 6082dfa..57f2905 100644 --- a/src/flow.c +++ b/src/flow.c @@ -1,8 +1,57 @@ #include #include +#include "flow.h" #include "gtsp.h" #include "util.h" +int flow_mark_reachable_nodes( + const struct Graph *graph, double *residual_caps, struct Node *from) +{ + int rval = 0; + + struct Node **stack; + int stack_top = 0; + int *parents; + + stack = (struct Node **) malloc(graph->node_count * sizeof(struct Node *)); + abort_if(!stack, "could not allocate stack"); + + parents = (int *) malloc(graph->node_count * sizeof(int )); + abort_if(!parents, "could not allocate parents"); + + stack[stack_top++] = from; + from->mark = 1; + + while (stack_top > 0) + { + struct Node *n = stack[--stack_top]; + + for (int j = 0; j < n->degree; j++) + { + struct Edge *e = n->adj[j].edge; + struct Node *neighbor = n->adj[j].neighbor; + + if (neighbor->mark) continue; + if (residual_caps[e->index] <= 0) continue; + + stack[stack_top++] = neighbor; + neighbor->mark = 1; + parents[neighbor->index] = n->index; + } + + } + + log_verbose("Reachable nodes:\n"); + for (int i = 0; i < graph->node_count; i++) + if (graph->nodes[i].mark) + log_verbose(" %d from %d\n", graph->nodes[i].index, parents[i]); + + CLEANUP: + if(parents) free(parents); + if (stack) free(stack); + return rval; +} + int flow_find_max_flow( const struct Graph *digraph, const double *capacities, @@ -13,6 +62,22 @@ int flow_find_max_flow( { int rval = 0; + for (int i = 0; i < digraph->node_count; i++) + digraph->nodes[i].mark = 0; + + log_verbose("Input graph:\n"); + graph_dump(digraph); + + log_verbose("Solving flow problem:\n"); + + log_verbose("%d %d\n", digraph->node_count, digraph->edge_count); + log_verbose("%d %d\n", from->index, to->index); + for (int i = 0; i < digraph->edge_count; i++) + { + log_verbose("%d %d %.4lf\n", digraph->edges[i].from->index, + digraph->edges[i].to->index, capacities[i]); + } + int path_length; struct Edge **path_edges = 0; double path_capacity; @@ -32,6 +97,8 @@ int flow_find_max_flow( residual_caps[i] = capacities[i]; abort_if(!digraph->edges[i].reverse, "digraph must have reverse edge information"); + abort_if(digraph->edges[i].reverse->reverse != &digraph->edges[i], + "invalid reverse edge"); } *value = 0; @@ -43,20 +110,37 @@ int flow_find_max_flow( if (path_length == 0) break; + log_verbose("Found augmenting path of capacity %.4lf:\n", + path_capacity); + (*value) += path_capacity; for (int i = 0; i < path_length; i++) { struct Edge *e = &digraph->edges[path_edges[i]->index]; + log_verbose(" %d %d (%d)\n", e->from->index, e->to->index, e->index); + residual_caps[e->index] -= path_capacity; residual_caps[e->reverse->index] += path_capacity; flow[e->index] += path_capacity; flow[e->reverse->index] -= path_capacity; } + + log_verbose("New residual capacities:\n"); + for (int i = 0; i < digraph->edge_count; i++) + { + struct Edge *e = &digraph->edges[i]; + if (residual_caps[i] < LP_EPSILON) continue; + + log_verbose("%d %d %.4lf (%d)\n", e->from->index, e->to->index, e->index, + residual_caps[e->index]); + } } + log_verbose("No more paths found.\n"); + rval = flow_mark_reachable_nodes(digraph, residual_caps, from); abort_if(rval, "flow_mark_reachable_nodes failed"); @@ -155,42 +239,108 @@ int flow_find_augmenting_path( return rval; } -int flow_mark_reachable_nodes( - struct Graph *graph, double *residual_caps, struct Node *from) +int flow_main(int argc, char **argv) { int rval = 0; - struct Node **stack; - int stack_top = 0; + int *edges = 0; + double *capacities = 0; + double *flow = 0; + double flow_value; - stack = (struct Node**) malloc(graph->node_count * sizeof(struct Node*)); - abort_if(!stack, "could not allocate stack"); + struct Edge **cut_edges = 0; - for (int i = 0; i < graph->node_count; i++) - graph->nodes[i].mark = 0; + FILE *f = fopen("tmp/flow.in", "r"); + abort_if(!f, "could not open input file"); - from->mark = 1; - stack[stack_top++] = from; + struct Graph graph; + graph_init(&graph); + + int node_count, edge_count; + + rval = fscanf(f, "%d %d ", &node_count, &edge_count); + abort_if(rval != 2, "invalid input format (node count, edge count)"); + + int sink, source; + rval = fscanf(f, "%d %d", &source, &sink); + abort_if(rval != 2, "invalid input format (source, sink)"); + + edges = (int *) malloc(4 * edge_count * sizeof(int)); + abort_if(!edges, "could not allocate edges\n"); - while(stack_top > 0) + capacities = (double *) malloc(2 * edge_count * sizeof(double)); + abort_if(!capacities, "could not allocate capacities"); + + for (int i = 0; i < edge_count; i++) { - struct Node *n = stack[--stack_top]; + int from, to; + double cap; - for (int j = 0; j < n->degree; j++) - { - struct Edge *e = n->adj[j].edge; - struct Node *neighbor = n->adj[j].neighbor; + rval = fscanf(f, "%d %d %lf ", &from, &to, &cap); + abort_if(rval != 3, "invalid input format (edge specification)"); - if(neighbor->mark) continue; - if(residual_caps[e->index] <= 0) continue; + edges[i * 4] = edges[i * 4 + 3] = from; + edges[i * 4 + 1] = edges[i * 4 + 2] = to; + capacities[2 * i] = cap; + capacities[2 * i + 1] = 0; + } - stack[stack_top++] = neighbor; - neighbor->mark = 1; - } + rval = graph_build(node_count, 2 * edge_count, edges, 1, &graph); + abort_if(rval, "graph_build failed"); + + for (int i = 0; i < edge_count; i++) + { + graph.edges[2 * i].reverse = &graph.edges[2 * i + 1]; + graph.edges[2 * i + 1].reverse = &graph.edges[2 * i]; + } + + flow = (double *) malloc(graph.edge_count * sizeof(double)); + abort_if(!flow, "could not allocate flow"); + + struct Node *from = &graph.nodes[source]; + struct Node *to = &graph.nodes[sink]; + + rval = flow_find_max_flow(&graph, capacities, from, to, flow, &flow_value); + abort_if(rval, "flow_find_max_flow failed"); + log_info("Optimal flow has value %f\n", flow_value); + for (int i = 0; i < graph.edge_count; i++) + { + struct Edge *e = &graph.edges[i]; + if (flow[e->index] <= 0) continue; + + log_info(" %d %d %6.2f / %6.2f\n", e->from->index, e->to->index, + flow[e->index], capacities[e->index]); + } + + log_info("Nodes reachable from origin on residual graph:\n"); + for (int i = 0; i < graph.node_count; i++) + { + struct Node *n = &graph.nodes[i]; + if (n->mark) + log_info(" %d\n", n->index); + } + + int cut_edges_count = 0; + cut_edges = + (struct Edge **) malloc(graph.edge_count * sizeof(struct Edge *)); + abort_if(!cut_edges, "could not allocate cut_edges"); + + rval = get_cut_edges_from_marks(&graph, &cut_edges_count, cut_edges); + abort_if(rval, "get_cut_edges_from_marks failed"); + + log_info("Min cut edges:\n"); + for (int i = 0; i < cut_edges_count; i++) + { + struct Edge *e = cut_edges[i]; + if (capacities[e->index] <= 0) continue; + log_info(" %d %d\n", e->from->index, e->to->index); } CLEANUP: - if(stack) free(stack); + if (cut_edges) free(cut_edges); + if (capacities) free(capacities); + if (edges) free(edges); + if (flow) free(flow); return rval; } \ No newline at end of file diff --git a/src/flow.h b/src/flow.h index 80d829c..f3a9595 100644 --- a/src/flow.h +++ b/src/flow.h @@ -1,10 +1,8 @@ -// -// Created by isoron on 19/03/15. -// - #ifndef _PROJECT_FLOW_H_ #define _PROJECT_FLOW_H_ +#include "graph.h" + int flow_find_augmenting_path( const struct Graph *graph, const double *residual_caps, @@ -22,6 +20,11 @@ int flow_find_max_flow( double *flow, double *value); +int flow_mark_reachable_nodes( + const struct Graph *graph, double *residual_caps, struct Node *from); + +int flow_main(int argc, char **argv); + #include "graph.h" #endif //_PROJECT_FLOW_H_ diff --git a/src/graph.c b/src/graph.c index 6862fac..a7656e8 100644 --- a/src/graph.c +++ b/src/graph.c @@ -196,3 +196,29 @@ int get_cut_edges_from_marks( return 0; } + +int graph_dump(struct Graph *graph) +{ + int rval = 0; + + log_debug("node_count: %d edge_count: %d\n", graph->node_count, graph->edge_count); + + for (int i = 0; i < graph->node_count; i++) + { + struct Node *n = &graph->nodes[i]; + log_debug("%3d degree: %d mark: %d\n", n->index, n->degree, n->mark); + } + + for (int i = 0; i < graph->edge_count; i++) + { + struct Edge *e = &graph->edges[i]; + log_debug("%3d (%d, %d) weight: %d ", e->index, e->from->index, e->to->index, e->weight); + #if LOG_LEVEL >= LOG_LEVEL_VERBOSE + if (e->reverse) printf("reverse: %d ", e->reverse->index); + printf("\n"); + #endif + } + + CLEANUP: + return rval; +} diff --git a/src/graph.h b/src/graph.h index b5849d7..b5a2c3e 100644 --- a/src/graph.h +++ b/src/graph.h @@ -73,4 +73,6 @@ int graph_build_directed_from_undirected int get_cut_edges_from_marks( struct Graph *graph, int *cut_edges_count, struct Edge **cut_edges); +int graph_dump(struct Graph *graph); + #endif diff --git a/src/gtsp.c b/src/gtsp.c index 164d0f8..53c75cf 100644 --- a/src/gtsp.c +++ b/src/gtsp.c @@ -1,24 +1,42 @@ #include #include +#include +#include +#include #include "gtsp.h" #include "geometry.h" #include "util.h" +#include "flow.h" +#include "branch_and_cut.h" + +static double *OPTIMAL_X = 0; int GTSP_init_data(struct GTSP *data) { + int rval = 0; + 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"); + graph_init(data->graph); - return 0; + CLEANUP: + return rval; } void GTSP_free(struct GTSP *data) { if (!data) return; - if (data->graph) graph_free(data->graph); + if (data->graph) + { + graph_free(data->graph); + 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); @@ -40,7 +58,7 @@ int GTSP_create_random_problem( int edge_count = (node_count * (node_count - 1)) / 2; - graph = (struct Graph*) malloc(sizeof(struct Graph)); + graph = (struct Graph *) malloc(sizeof(struct Graph)); abort_if(!graph, "could not allocate graph\n"); graph_init(graph); @@ -78,6 +96,9 @@ int GTSP_create_random_problem( rval = graph_build(node_count, edge_count, edges, 0, graph); abort_if(rval, "graph_build failed"); + for (int i = 0; i < edge_count; i++) + graph->edges[i].weight = weights[i]; + data->graph = graph; data->clusters = clusters; data->cluster_count = cluster_count; @@ -145,7 +166,271 @@ int GTSP_init_lp(struct LP *lp, struct GTSP *data) return rval; } -int GTSP_write_data(struct GTSP *data, char *filename) +int GTSP_add_subtour_elimination_cut( + struct LP *lp, + struct Graph *graph, + struct Node *from, + struct Node *to, + struct Edge **cut_edges, + int cut_edges_count) +{ + int rval = 0; + + char sense = 'G'; + double rhs = -2.0; + int newnz = cut_edges_count + 2; + + 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; + } + + rmatind[cut_edges_count] = from->index; + rmatval[cut_edges_count] = -2.0; + + rmatind[cut_edges_count + 1] = to->index; + rmatval[cut_edges_count + 1] = -2.0; + + log_debug("Generated cut:\n"); + for (int i = 0; i < newnz; i++) + log_debug("%8.2f x%d\n", rmatval[i], rmatind[i]); + log_debug(" %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 GTSP_find_exact_subtour_elimination_cuts( + struct LP *lp, struct GTSP *data, int *added_cuts_count) +{ + int rval = 0; + + int *clusters = data->clusters; + + double *x = 0; + double *capacities = 0; + + double *flow = 0; + struct Edge **cut_edges = 0; + + int *digraph_edges = 0; + + struct Graph *graph = data->graph; + int node_count = graph->node_count; + + 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 digraph; + graph_init(&digraph); + + digraph_edges = (int *) malloc(8 * graph->edge_count * sizeof(int)); + flow = (double *) malloc(4 * graph->edge_count * sizeof(double)); + capacities = (double *) malloc(4 * graph->edge_count * sizeof(double)); + cut_edges = + (struct Edge **) malloc( + 4 * graph->edge_count * sizeof(struct Edge *)); + + abort_if(!digraph_edges, "could not allocate digraph_edges"); + abort_if(!flow, "could not allocate flow"); + abort_if(!capacities, "could not allocate capacities"); + abort_if(!cut_edges, "could not allocate cut_edges"); + + // Create four directed edges for each edge of the original graph. + for (int i = 0; i < graph->edge_count; i++) + { + assert(node_count + i < num_cols); + + struct Edge *e = &graph->edges[i]; + int from = e->from->index; + int to = e->to->index; + + digraph_edges[8 * i] = from; + digraph_edges[8 * i + 1] = to; + capacities[4 * i] = x[node_count + i]; + + digraph_edges[8 * i + 2] = to; + digraph_edges[8 * i + 3] = from; + capacities[4 * i + 1] = 0; + + digraph_edges[8 * i + 4] = to; + digraph_edges[8 * i + 5] = from; + capacities[4 * i + 2] = x[node_count + i]; + + digraph_edges[8 * i + 6] = from; + digraph_edges[8 * i + 7] = to; + capacities[4 * i + 3] = 0; + } + + rval = graph_build(node_count, 4 * graph->edge_count, digraph_edges, 1, &digraph); + abort_if(rval, "graph_build failed"); + + for (int i = 0; i < graph->edge_count; i++) + { + digraph.edges[4 * i].reverse = &digraph.edges[4 * i + 1]; + digraph.edges[4 * i + 1].reverse = &digraph.edges[4 * i]; + + digraph.edges[4 * i + 2].reverse = &digraph.edges[4 * i + 3]; + digraph.edges[4 * i + 3].reverse = &digraph.edges[4 * i + 2]; + } + + 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 < digraph.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]; + + log_verbose("Calculating max flow from %d to %to\n", from->index, + to->index); + double flow_value; + rval = flow_find_max_flow(&digraph, capacities, from, to, flow, + &flow_value); + abort_if(rval, "flow_find_max_flow failed"); + + log_verbose(" %.2lf\n", flow_value); + + if (flow_value >= 2 * (x[i] + x[j] - 1) - LP_EPSILON) continue; + + + log_verbose("violation: %.2lf >= %.2lf\n", flow_value, + 2 * (x[i] + x[j] - 1)); + + int cut_edges_count; + rval = get_cut_edges_from_marks(&digraph, &cut_edges_count, cut_edges); + abort_if(rval, "get_cut_edges_from_marks failed"); + + log_verbose("Adding cut for i=%d j=%d, cut edges:\n", i, j); + for (int k = 0; k < cut_edges_count/2; k++) + { + cut_edges[k] = &graph->edges[cut_edges[k*2]->index / 4]; + log_verbose(" %d %d\n", cut_edges[k*2]->from->index, + cut_edges[k*2]->to->index); + } + + rval = GTSP_add_subtour_elimination_cut(lp, graph, from, to, cut_edges, + cut_edges_count/2); + abort_if(rval, "GTSP_add_subtour_elimination_cut failed"); + + (*added_cuts_count)++; + } + + CLEANUP: + if (digraph_edges) free(digraph_edges); + if (flow) free(flow); + if (cut_edges) free(cut_edges); + if (capacities) free(capacities); + if (x) free(x); + return rval; +} + +int GTSP_add_cutting_planes(struct LP *lp, struct GTSP *data) +{ + int rval = 0; + double *x = 0; + + int num_cols = LP_get_num_cols(lp); + x = (double *) malloc(num_cols * sizeof(double)); + abort_if(!x, "could not allocate x"); + + while (1) + { + int added_cuts_count = 0; + + rval = GTSP_find_exact_subtour_elimination_cuts(lp, data, + &added_cuts_count); + abort_if(rval, "GTSP_find_exact_subtour_elimination_cuts failed"); + + log_verbose("Reoptimizing...\n"); + int is_infeasible; + rval = LP_optimize(lp, &is_infeasible); + abort_if(rval, "LP_optimize failed"); + + if (is_infeasible) break; + + double objval; + rval = LP_get_obj_val(lp, &objval); + abort_if(rval, "LP_get_obj_val failed"); + + rval = LP_get_x(lp, x); + abort_if(rval, "LP_get_x failed"); + + log_verbose("Current solution:\n"); + for (int i = 0; i < data->graph->node_count; i++) + if (x[i] > LP_EPSILON) log_verbose(" node %d = %.2f\n", i, x[i]); + + for (int i = 0; i < data->graph->edge_count; i++) + { + struct Edge *e = &data->graph->edges[i]; + int idx = e->index + data->graph->node_count; + if (x[idx] > LP_EPSILON) + { + log_verbose(" edge (%d, %d) = %.2f\n", e->from->index, + e->to->index, x[idx]); + } + } + + log_debug(" obj val = %f\n", objval); + + if (added_cuts_count > 0) + { + log_debug("Found %d subtour elimination cuts using exact " + "separation\n", added_cuts_count); + } + else break; + } + + CLEANUP: + if (x) free(x); + return rval; +} + +int GTSP_write_input_data(struct GTSP *data, char *filename) { int rval = 0; @@ -184,6 +469,8 @@ int GTSP_write_solution(struct GTSP *data, char *filename, double *x) if (x[i + node_count] > 0.5) positive_edge_count++; + fprintf(file, "%d %d\n", node_count, edge_count); + fprintf(file, "%d\n", positive_edge_count); for (int i = 0; i < edge_count; i++) @@ -195,3 +482,213 @@ int GTSP_write_solution(struct GTSP *data, char *filename, double *x) return rval; } +int get_edge_num(int node_count, int from, int to) +{ + int idx = node_count; + + for (int k = 0; k < from; k++) + idx += node_count - k - 1; + + idx += to - from - 1; + + return idx; +} + +int GTSP_read_x(char *filename, double **p_x) +{ + int rval = 0; + + int node_count; + int edge_count; + + double *x; + + FILE *file; + + log_info("Reading optimal solution from file %s\n", filename); + + file = fopen(filename, "r"); + abort_if(!file, "could not open file"); + + rval = fscanf(file, "%d %d", &node_count, &edge_count); + abort_if(rval != 2, "invalid input format (node and edge count)"); + + int num_cols = node_count + edge_count; + + x = (double *) malloc(num_cols * sizeof(double)); + abort_if(!x, "could not allocate x"); + + for (int i = 0; i < node_count + edge_count; i++) x[i] = 0.0; + + rval = fscanf(file, "%d", &edge_count); + abort_if(rval != 1, "invalid input format (positive edge count)"); + + 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)"); + + if (from > to) swap(from, to); + + edge = get_edge_num(node_count, from, to); + abort_if(edge > num_cols, "invalid edge"); + + x[from] = x[ + to] = 1.0; + x[edge] = 1; + } + + for (int i = 0; i < num_cols; i++) + { + if (x[i] <= LP_EPSILON) continue; + log_verbose(" x%-3d = %.2f\n", i, x[i]); + } + + *p_x = x; + rval = 0; + + CLEANUP: + 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 int input_node_count = 20; +static int input_cluster_count = 5; +static int grid_size = 100; + +static void GTSP_print_usage(char **argv) +{ + printf("wrong usage\n"); +} + +static int GTSP_parse_args(int argc, char **argv) +{ + int rval = 0; + + opterr = 0; + + while (1) + { + int c = 0; + int option_index = 0; + c = getopt_long(argc, argv, "n:m:g:x:s:", options_tab, &option_index); + + if (c < 0) break; + + switch (c) + { + case 'n': + input_node_count = atoi(optarg); + break; + + case 'm': + input_cluster_count = atoi(optarg); + break; + + case 'g': + grid_size = atoi(optarg); + break; + + case 'x': + rval = GTSP_read_x(optarg, &OPTIMAL_X); + abort_if(rval, "GTSP_read_x failed"); + break; + + case 's': + SEED = atoi(optarg); + break; + + case ':': + fprintf(stderr, "option '-%c' requires an argument\n", optopt); + return 1; + + case '?': + default: + fprintf(stderr, "option '-%c' is invalid\n", optopt); + return 1; + + } + } + + CLEANUP: + return rval; +} + +int GTSP_main(int argc, char **argv) +{ + int rval = 0; + + struct BNC bnc; + struct GTSP data; + + SEED = (unsigned int) get_real_time() % 10000; + + rval = GTSP_init_data(&data); + abort_if(rval, "GTSP_init_data failed"); + + rval = BNC_init(&bnc); + abort_if(rval, "BNC_init failed"); + + rval = GTSP_parse_args(argc, argv); + if (rval) return 1; + + srand(SEED); + + log_info("Generating random GTSP instance...\n"); + log_info(" seed = %d\n", SEED); + log_info(" input_node_count = %d\n", input_node_count); + log_info(" input_cluster_count = %d\n", input_cluster_count); + log_info(" grid_size = %d\n", grid_size); + + rval = GTSP_create_random_problem(input_node_count, input_cluster_count, + grid_size, + &data); + abort_if(rval, "GTSP_create_random_problem failed"); + + log_info("Writing random instance to file gtsp.in\n"); + rval = GTSP_write_input_data(&data, "gtsp.in"); + abort_if(rval, "GTSP_write_problem failed"); + + bnc.best_obj_val = DBL_MAX; + bnc.problem_data = (void *) &data; + bnc.problem_init_lp = (int (*)(struct LP *, void *)) GTSP_init_lp; + bnc.problem_add_cutting_planes = + (int (*)(struct LP *, void *)) GTSP_add_cutting_planes; + + if (OPTIMAL_X) + log_info("Optimal solution is available. Cuts will be checked.\n"); + + log_info("Initializing LP...\n"); + rval = BNC_init_lp(&bnc); + abort_if(rval, "BNC_init_lp failed"); + + log_info("Writing LP to file gtsp.lp...\n"); + rval = LP_write(bnc.lp, "gtsp.lp"); + abort_if(rval, "LP_write failed"); + + log_info("Starting branch-and-cut solver...\n"); + rval = BNC_solve(&bnc); + abort_if(rval, "BNC_solve_node failed"); + + abort_if(!bnc.best_x, "problem has no feasible solution"); + + log_info("Optimal integral solution:\n"); + log_info(" obj value = %.2lf **\n", bnc.best_obj_val); + + rval = GTSP_write_solution(&data, "gtsp.out", bnc.best_x); + abort_if(rval, "GTSP_write_solution failed"); + + CLEANUP: + GTSP_free(&data); + BNC_free(&bnc); + return rval; +} \ No newline at end of file diff --git a/src/gtsp.h b/src/gtsp.h index 0198cce..b69e1ce 100644 --- a/src/gtsp.h +++ b/src/gtsp.h @@ -28,8 +28,12 @@ int GTSP_init_data(struct GTSP *data); int GTSP_init_lp(struct LP *lp, struct GTSP *data); -int GTSP_write_data(struct GTSP *data, char *filename); +int GTSP_write_input_data(struct GTSP *data, char *filename); int GTSP_write_solution(struct GTSP *data, char *filename, double *x); +int GTSP_add_cutting_planes(struct LP *lp, struct GTSP *data); + +int GTSP_main(int argc, char **argv); + #endif //_PROJECT_GTSP_H_ diff --git a/src/main.c b/src/main.c index dee42ed..77451c9 100644 --- a/src/main.c +++ b/src/main.c @@ -1,13 +1,7 @@ -#include -#include #include -#include -#include "lp.h" -#include "util.h" #include "main.h" -#include "tsp.h" -#include "branch_and_cut.h" #include "gtsp.h" +#include "tsp.h" #include "flow.h" char *INPUT_FILENAME = 0; @@ -16,276 +10,45 @@ int GEOMETRIC_DATA = 0; int NODE_COUNT_RAND = 0; int GRID_SIZE_RAND = 100; -static int parse_arguments_tsp(int ac, char **av); - -static void print_usage_tsp(char *f); +static const struct option options_tab[] = { + {"help", no_argument, 0, 'h'}, {"tsp", no_argument, 0, 't'}, + {"gtsp", no_argument, 0, 'g'}, {"flow", no_argument, 0, 'f'}, + {(char *) 0, (int) 0, (int *) 0, (int) 0} +}; -int test_max_flow() +void GTSP_print_usage(char **argv) { - int rval = 0; - - int *edges = 0; - double *capacities = 0; - double *flow = 0; - double flow_value; - - struct Edge **cut_edges = 0; - - FILE *f = fopen("tmp/flow.in", "r"); - abort_if(!f, "could not open input file"); - - struct Graph graph; - graph_init(&graph); - - int node_count, edge_count; - - rval = fscanf(f, "%d %d ", &node_count, &edge_count); - abort_if(rval != 2, "invalid input format"); - - edges = (int *) malloc(4 * edge_count * sizeof(int)); - abort_if(!edges, "could not allocate edges\n"); - - capacities = (double *) malloc(2 * edge_count * sizeof(double)); - abort_if(!capacities, "could not allocate capacities"); - - for (int i = 0; i < edge_count; i++) - { - int from, to, cap; - - rval = fscanf(f, "%d %d %d ", &from, &to, &cap); - abort_if(rval != 3, "invalid input format"); - - edges[i * 4] = edges[i * 4 + 3] = from; - edges[i * 4 + 1] = edges[i * 4 + 2] = to; - capacities[2 * i] = cap; - capacities[2 * i + 1] = 0; - } - - rval = graph_build(node_count, 2 * edge_count, edges, 1, &graph); - abort_if(rval, "graph_build failed"); - - for (int i = 0; i < edge_count; i++) - { - graph.edges[2 * i].reverse = &graph.edges[2 * i + 1]; - graph.edges[2 * i + 1].reverse = &graph.edges[2 * i]; - } - - flow = (double *) malloc(graph.edge_count * sizeof(double)); - abort_if(!flow, "could not allocate flow"); - - struct Node *from = &graph.nodes[0]; - struct Node *to = &graph.nodes[graph.node_count - 1]; - - rval = flow_find_max_flow(&graph, capacities, from, to, flow, &flow_value); - abort_if(rval, "flow_find_max_flow failed"); - - log_info("Optimal flow has value %f\n", flow_value); - for (int i = 0; i < graph.edge_count; i++) - { - struct Edge *e = &graph.edges[i]; - if (flow[e->index] <= 0) continue; - - log_info(" %d %d %6.2f / %6.2f\n", e->from->index, e->to->index, - flow[e->index], capacities[e->index]); - } - - log_info("Nodes reachable from origin on residual graph:\n"); - for (int i = 0; i < graph.node_count; i++) - { - struct Node *n = &graph.nodes[i]; - if(n->mark) - log_info(" %d\n", n->index); - } - - int cut_edges_count = 0; - cut_edges = (struct Edge**) malloc(graph.edge_count * sizeof(struct Edge*)); - abort_if(!cut_edges, "could not allocate cut_edges"); - - rval = get_cut_edges_from_marks(&graph, &cut_edges_count, cut_edges); - abort_if(rval, "get_cut_edges_from_marks failed"); - - log_info("Min cut edges:\n"); - for (int i = 0; i < cut_edges_count; i++) - { - struct Edge *e = cut_edges[i]; - if(capacities[e->index] <= 0) continue; - log_info(" %d %d\n", e->from->index, e->to->index); - } - - - CLEANUP: - if(cut_edges) free(cut_edges); - if (capacities) free(capacities); - if (edges) free(edges); - if (flow) free(flow); - return rval; + printf("wrong usage\n"); } -int main_tsp(int ac, char **av) +int main(int argc, char **argv) { - int rval = 0; - - SEED = (unsigned int) get_real_time(); - - struct BNC bnc; - struct TSPData data; - - rval = TSP_init_data(&data); - abort_if(rval, "TSP_init_data failed"); - - rval = BNC_init(&bnc); - abort_if(rval, "BNC_init failed"); + int c = 0; + int option_index = 0; - rval = parse_arguments_tsp(ac, av); - abort_if(rval, "Failed to parse arguments."); + c = getopt_long(argc, argv, "htgf", options_tab, &option_index); - printf("Seed = %d\n", SEED); - srand(SEED); - - rval = TSP_read_problem(INPUT_FILENAME, &data); - abort_if(rval, "TSP_read_problem failed"); - - 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; - - rval = BNC_init_lp(&bnc); - abort_if(rval, "BNC_init_lp failed"); - - rval = BNC_solve(&bnc); - abort_if(rval, "BNC_solve_node failed"); - - log_info("Optimal integral solution:\n"); - log_info(" obj value = %.2lf **\n", bnc.best_obj_val); - - CLEANUP: - BNC_free(&bnc); - TSP_free_data(&data); - return rval; -} - -int main_gtsp(int ac, char **av) -{ - int rval = 0; - - SEED = (unsigned int) get_real_time(); - srand(SEED); - - int node_count = 50; - int cluster_count = node_count / 5; - int grid_size = 100; - - struct BNC bnc; - struct GTSP data; - - rval = GTSP_init_data(&data); - abort_if(rval, "GTSP_init_data failed"); - - rval = GTSP_create_random_problem(node_count, cluster_count, grid_size, - &data); - abort_if(rval, "GTSP_create_random_problem failed"); - - rval = GTSP_write_data(&data, "gtsp.in"); - abort_if(rval, "GTSP_write_problem failed"); - - rval = BNC_init(&bnc); - abort_if(rval, "BNC_init failed"); - - log_info("Setting seed = %d\n", SEED); - srand(SEED); - - bnc.best_obj_val = DBL_MAX; - bnc.problem_data = (void *) &data; - bnc.problem_init_lp = (int (*)(struct LP *, void *)) GTSP_init_lp; - - rval = BNC_init_lp(&bnc); - abort_if(rval, "BNC_init_lp failed"); - - log_info("Starting branch-and-cut solver...\n"); - rval = BNC_solve(&bnc); - abort_if(rval, "BNC_solve_node failed"); - - log_info("Optimal integral solution:\n"); - log_info(" obj value = %.2lf **\n", bnc.best_obj_val); - - rval = GTSP_write_solution(&data, "gtsp.out", bnc.best_x); - abort_if(rval, "GTSP_write_solution failed"); - - CLEANUP: - GTSP_free(&data); - BNC_free(&bnc); - return rval; - -} - -static int parse_arguments_tsp(int ac, char **av) -{ - int rval = 0; - - int c; - - if (ac == 1) + if (c < 0) { - print_usage_tsp(av[0]); + GTSP_print_usage(argv); return 1; } - while ((c = getopt(ac, av, "ab:gk:s:")) != EOF) - { - switch (c) - { - case 'a':; - break; - case 'b': - GRID_SIZE_RAND = atoi(optarg); - break; - case 'g': - GEOMETRIC_DATA = 1; - break; - case 'k': - NODE_COUNT_RAND = atoi(optarg); - break; - case 's': - SEED = (unsigned) atoi(optarg); - break; - case '?': - default: - print_usage_tsp(av[0]); - return 1; - } - } - - if (optind < ac) INPUT_FILENAME = av[optind++]; - - if (optind != ac) + switch (c) { - print_usage_tsp(av[0]); - return 1; - } + case 'h': + GTSP_print_usage(argv); + return 1; - abort_if(!INPUT_FILENAME && !NODE_COUNT_RAND, - "Must specify an input file or use -k for random problem\n"); + case 'f': + return flow_main(argc, argv); - CLEANUP: - return rval; -} + case 't': + return TSP_main(argc, argv); -int main(int ac, char **av) -{ - return test_max_flow(); -// return main_gtsp(ac, av); -// return main_tsp(ac, av); -} + case 'g': + return GTSP_main(argc, argv); + } -static void print_usage_tsp(char *f) -{ - fprintf(stderr, "Usage: %s [-see below-] [prob_file]\n" - " -a add all subtours cuts at once\n" - " -b d gridsize d for random problems\n" - " -g prob_file has x-y coordinates\n" - " -k d generate problem with d cities\n" - " -s d random SEED\n", f); } diff --git a/src/tsp.c b/src/tsp.c index 33052e0..b668ebf 100644 --- a/src/tsp.c +++ b/src/tsp.c @@ -1,9 +1,14 @@ #include #include -#include "main.h" +#include #include "tsp.h" #include "util.h" #include "geometry.h" +#include "branch_and_cut.h" + +static int TSP_parse_arguments(int argc, char **argv); + +static void TSP_print_usage(char *f); int TSP_init_data(struct TSPData *data) { @@ -113,7 +118,7 @@ int TSP_find_violated_subtour_elimination_cut( while (!TSP_is_graph_connected(&G, x, &island_count, island_sizes, island_start, island_nodes)) { - log_verbose("Adding %d BNC_solve_node inequalities...\n", island_count); + log_verbose("Adding %d subtour inequalities...\n", island_count); for (int i = 0; i < island_count; i++) { get_delta(island_sizes[i], island_nodes + island_start[i], @@ -152,6 +157,7 @@ int TSP_find_violated_subtour_elimination_cut( int TSP_add_subtour_elimination_cut(struct LP *lp, int delta_length, int *delta) { int rval = 0; + char sense = 'G'; double rhs = 2.0; int rmatbeg = 0; @@ -453,4 +459,111 @@ double TSP_find_initial_solution(struct TSPData *data) log_verbose(" length = %lf\n", best_val); return best_val; +} + +int TSP_main(int argc, char **argv) +{ + int rval = 0; + + SEED = (unsigned int) get_real_time(); + + struct BNC bnc; + struct TSPData data; + + rval = TSP_init_data(&data); + abort_if(rval, "TSP_init_data failed"); + + rval = BNC_init(&bnc); + abort_if(rval, "BNC_init failed"); + + rval = TSP_parse_arguments(argc, argv); + abort_if(rval, "Failed to parse arguments."); + + printf("Seed = %d\n", SEED); + srand(SEED); + + rval = TSP_read_problem(INPUT_FILENAME, &data); + abort_if(rval, "TSP_read_problem failed"); + + 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; + + rval = BNC_init_lp(&bnc); + abort_if(rval, "BNC_init_lp failed"); + + rval = BNC_solve(&bnc); + abort_if(rval, "BNC_solve_node failed"); + + log_info("Optimal integral solution:\n"); + log_info(" obj value = %.2lf **\n", bnc.best_obj_val); + + CLEANUP: + BNC_free(&bnc); + TSP_free_data(&data); + return rval; +} + +static int TSP_parse_arguments(int argc, char **argv) +{ + int rval = 0; + + int c; + + if (argc == 1) + { + TSP_print_usage(argv[0]); + return 1; + } + + while ((c = getopt(argc, argv, "ab:gk:s:")) != EOF) + { + switch (c) + { + case 'a':; + break; + case 'b': + GRID_SIZE_RAND = atoi(optarg); + break; + case 'g': + GEOMETRIC_DATA = 1; + break; + case 'k': + NODE_COUNT_RAND = atoi(optarg); + break; + case 's': + SEED = (unsigned) atoi(optarg); + break; + case '?': + default: + TSP_print_usage(argv[0]); + return 1; + } + } + + if (optind < argc) INPUT_FILENAME = argv[optind++]; + + if (optind != argc) + { + TSP_print_usage(argv[0]); + return 1; + } + + abort_if(!INPUT_FILENAME && !NODE_COUNT_RAND, + "Must specify an input file or use -k for random problem\n"); + + CLEANUP: + return rval; +} + +static void TSP_print_usage(char *f) +{ + fprintf(stderr, "Usage: %s [-see below-] [prob_file]\n" + " -a add all subtours cuts at once\n" + " -b d gridsize d for random problems\n" + " -g prob_file has x-y coordinates\n" + " -k d generate problem with d cities\n" + " -s d random SEED\n", f); } \ No newline at end of file diff --git a/src/tsp.h b/src/tsp.h index 2c4160e..3361669 100644 --- a/src/tsp.h +++ b/src/tsp.h @@ -45,4 +45,6 @@ int TSP_init_lp(struct LP *lp, struct TSPData *data); double TSP_find_initial_solution(struct TSPData *data); +int TSP_main(int argc, char **argv); + #endif diff --git a/src/util.h b/src/util.h index c99bf8e..43f3bb3 100644 --- a/src/util.h +++ b/src/util.h @@ -1,13 +1,15 @@ #ifndef _PROJECT_UTIL_H_ #define _PROJECT_UTIL_H_ +#include + #define LOG_LEVEL_ERROR 10 #define LOG_LEVEL_WARNING 20 #define LOG_LEVEL_INFO 30 #define LOG_LEVEL_DEBUG 40 #define LOG_LEVEL_VERBOSE 50 -#define LOG_LEVEL LOG_LEVEL_DEBUG +#define LOG_LEVEL LOG_LEVEL_INFO #if LOG_LEVEL < LOG_LEVEL_DEBUG #define log_debug(...) @@ -43,6 +45,13 @@ fprintf(stderr, "%20s:%d " msg "\n", __FILE__, __LINE__); \ rval = 1; goto CLEANUP; } +#define swap(x, y) do \ + { unsigned char swap_temp[sizeof(x)]; \ + memcpy(swap_temp,&y,sizeof(x)); \ + memcpy(&y,&x, sizeof(x)); \ + memcpy(&x,swap_temp,sizeof(x)); \ + } while(0) + void time_printf(const char *fmt, ...); double get_current_time(void);