diff --git a/src/branch_and_cut.c b/src/branch_and_cut.c index 98cf80b..29465a5 100644 --- a/src/branch_and_cut.c +++ b/src/branch_and_cut.c @@ -6,9 +6,6 @@ #include "util.h" #include "gtsp.h" - -int BNC_NODE_COUNT = 0; - static int BNC_solve_node(struct BNC *bnc, int depth); static int BNC_branch_node(struct BNC *bnc, double *x, int depth); @@ -17,8 +14,6 @@ static int BNC_is_integral(double *x, int num_cols); static int BNC_find_best_branching_var(double *x, int num_cols); -//int optimize_vertex_in_cluster(struct BNC *bnc, double best_val); - int BNC_init(struct BNC *bnc) { int rval = 0; @@ -101,7 +96,7 @@ static int BNC_solve_node(struct BNC *bnc, int depth) if(depth == 1) ROOT_VALUE = objval; - if (ceil(objval) > *best_val + LP_EPSILON) + if (ceil(objval) > *best_val + EPSILON) { log_debug("Branch pruned by bound (%.2lf > %.2lf).\n", objval, *best_val); @@ -127,9 +122,9 @@ static int BNC_solve_node(struct BNC *bnc, int depth) abort_if(rval, "LP_get_obj_val failed"); if(depth == 1) ROOT_VALUE = objval; - abort_if(get_current_time() - INITIAL_TIME >= MAX_TOTAL_TIME, "time limit exceeded"); + abort_if(get_user_time() - INITIAL_TIME >= MAX_TOTAL_TIME, "time limit exceeded"); - if (ceil(objval) > *best_val + LP_EPSILON) + if (ceil(objval) > *best_val + EPSILON) { log_debug("Branch pruned by bound (%.2lf > %.2lf).\n", objval, *best_val); @@ -145,7 +140,7 @@ static int BNC_solve_node(struct BNC *bnc, int depth) { log_debug("Solution is integral\n"); - if (objval + LP_EPSILON < *best_val) + if (objval + EPSILON < *best_val) { if (bnc->best_x) free(bnc->best_x); *best_val = objval; @@ -214,17 +209,11 @@ 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) + if (x[i] > EPSILON && x[i] < 1.0 - EPSILON) return 0; return 1; -#endif } static int BNC_find_best_branching_var(double *x, int num_cols) @@ -244,129 +233,3 @@ static int BNC_find_best_branching_var(double *x, int num_cols) return best_index; } - -/* -int re_optimize_integral(struct BNC *bnc){ - int i = 0 , current_vertex = 0, rval = 0; - struct GTSP* data; - data = bnc->problem_data; - int node_count = data->graph->node_count; - int cluster_count = data->cluster_count; - int edge_count = data->graph->edge_count; - struct Tour * tour = (struct Tour*) NULL; - - //intialize the tour - tour = (struct Tour *) malloc( cluster_count * sizeof(struct Tour)); - for (i = 0; i < edge_count; i++){ - tour[i].vertex = -1; - tour[i].next = -1; - tour[i].prev = -1; - } - - //Constructing the tour with vertices - for (i = 0; i < edge_count; i++){ - if (bnc->best_x[i + node_count] > LP_EPSILON) { - tour[current_vertex].vertex = data->graph->edges[i].from->index; - current_vertex += 1; - printf("From node %d \t", data->graph->edges[i].from->index); - printf("TO node %d \n", data->graph->edges[i].to->index); - } - } - //printf("Edgese in solution %d \n", current_vertex); - - - return rval; - CLEANUP: - if (data) free(data); - -} -*/ -/* -int optimize_vertex_in_cluster(struct BNC *bnc, double best_val) -{ - - int i = 0 , j, current_vertex = 0, rval = 0; - int tour_cost = 0; - struct GTSP* data; - data = bnc->problem_data; - //rval = GTSP_init_data(&data); - //data = bnc->problem_data; - //data = (struct GTSP) malloc(sizeof(struct GTSP)); - - //data = &bnc->problem_data; - int node_count = data->graph->node_count; - int cluster_count = data->cluster_count; - int edge_count = data->graph->edge_count; - int * tour = (int*) NULL; - tour = (int *) malloc( cluster_count * sizeof(int)); - //Constructing the tour with vertices - for (i = 0; i < edge_count; i++) - { //printf(" edge %lf **\n", bnc->best_x[i]); - - if ((bnc->best_x[i] > 1 - LP_EPSILON)){ - //printf(" x[i] = %lf **\n", bnc->best_x[i]); - tour[current_vertex] = (data->graph->edges[i].from)->index; - current_vertex += 1; - //printf(" Edge No = %d **\n", i); - printf(" FROM No = %d **\n", (data->graph->edges[i].from)->index); - printf(" TO No = %d **\n", (data->graph->edges[i].to)->index); - - //printf(" current vertex = %d **\n", current_vertex); - } - } - - //reoptmizing the your with two-opt - //rval = two_opt(cluster_count, tour, data->dist_matrix); - //Optimizing the vertices inside the node_to_cluster - int current_cluster = 0; - int insertion_cost = 0; - - //printf(" o-- val = %.2lf **\n", best_val); - for(i = 1; i < cluster_count - 2; i++){ - //printf(" vertex in tour = %d **\n", tour[current_vertex]); - current_cluster = data->node_to_cluster[tour[i]]; - //printf(" o-- val = %.2lf **\n", best_val); - insertion_cost = data->dist_matrix[tour[i-1]][tour[i]] + - data->dist_matrix[tour[i]][tour[i+1]]; - //printf(" o-- val = %.2lf **\n", best_val); - for(j = 0; j < node_count; j++) - if (current_cluster == data->node_to_cluster[j]) - if (insertion_cost > data->dist_matrix[j][tour[i]] + - data->dist_matrix[j][tour[i+1]]){ - log_info("Optmize vertex in cluster improved the bound\n"); - insertion_cost = data->dist_matrix[j][tour[i]] + - data->dist_matrix[j][tour[i+1]]; - tour[i] = j; - } - } - printf(" o-- val = %.2lf **\n", best_val); - for(i = 0; i< cluster_count ; i++){ - if (i == cluster_count - 1) - tour_cost += data->dist_matrix[tour[i]][tour[0]]; - else - tour_cost += data->dist_matrix[tour[i]][tour[i+1]]; - if(tour_cost < bnc->best_obj_val) - bnc->best_obj_val = tour_cost; - } - return rval; -} -*/ - -/* -static int two_opt(int tour_length, int*tour, int** dist_matrix){ - int rval = 0, i; - for (i = 1; i < tour_length - 2; i++){ - int current_cost = dist_matrix[tour[i-1]][tour[i]] + - dist_matrix[tour[i+1]][tour[i+2]]; - int temp_cost = dist_matrix[tour[i-1]][tour[i+1]] + - dist_matrix[tour[i]][tour[i+2]]; - if(current_cost > temp_cost){ - log_info("Two opt improved the bound\n"); - int temp_vertex = tour[i]; - tour[i] = tour[i+1]; - tour[i+1] = temp_vertex; - } - } - return rval; -} -*/ diff --git a/src/branch_and_cut.h b/src/branch_and_cut.h index 89fbaed..61b5f9e 100644 --- a/src/branch_and_cut.h +++ b/src/branch_and_cut.h @@ -10,8 +10,6 @@ struct BNC double *best_x; double best_obj_val; - double *optimal_x; - int *problem_data; int (*problem_init_lp)(struct LP *, void *); @@ -29,10 +27,6 @@ int BNC_init_lp(struct BNC *bnc); void BNC_free(struct BNC *bnc); -int re_optimize_integral(struct BNC *bnc); - -//int optimize_vertex_in_cluster(struct BNC *bnc, double best_val); - extern int BNC_NODE_COUNT; #endif //_PROJECT_BRANCH_AND_CUT_H_ diff --git a/src/flow.c b/src/flow.c index ac802de..1237b52 100644 --- a/src/flow.c +++ b/src/flow.c @@ -142,7 +142,7 @@ int flow_find_max_flow( struct Edge *e = &digraph->edges[i]; #endif - if (residual_caps[i] < LP_EPSILON) continue; + if (residual_caps[i] < EPSILON) continue; log_verbose("%d %d %.4lf (%d)\n", e->from->index, e->to->index, e->index, residual_caps[e->index]); @@ -206,7 +206,7 @@ int flow_find_augmenting_path( struct Edge *edge = n->adj[i].edge; if (neighbor->mark > 0) continue; - if (residual_caps[edge->index] < LP_EPSILON) continue; + if (residual_caps[edge->index] < EPSILON) continue; parents[neighbor->index] = n; parent_edges[neighbor->index] = edge; @@ -248,112 +248,3 @@ int flow_find_augmenting_path( if (queue) free(queue); return rval; } - -int flow_main(int argc, char **argv) -{ - UNUSED(argc); - UNUSED(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 (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"); - - 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; - double cap; - - rval = fscanf(f, "%d %d %lf ", &from, &to, &cap); - abort_if(rval != 3, "invalid input format (edge specification)"); - - 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[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 (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 f3a9595..2b58e4e 100644 --- a/src/flow.h +++ b/src/flow.h @@ -23,8 +23,4 @@ int flow_find_max_flow( 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/geometry.h b/src/geometry.h index 1d6dfff..68415cd 100644 --- a/src/geometry.h +++ b/src/geometry.h @@ -22,7 +22,9 @@ int get_euclidean_distance( int p2_index); int generate_dist_matrix( - int node_count, + int node_count, double *x_coordinates, - double *y_coordinates, int** dist_matrix); + double *y_coordinates, + int **dist_matrix); + #endif //_PROJECT_GEOMETRY_H_ diff --git a/src/graph.c b/src/graph.c index 1439606..019b091 100644 --- a/src/graph.c +++ b/src/graph.c @@ -1,7 +1,6 @@ #include #include "graph.h" #include "util.h" -#include "lp.h" void graph_init(struct Graph *graph) { @@ -109,50 +108,6 @@ int graph_build( return rval; } -void graph_dfs( - int n, struct Graph *G, double *x, int *island_size, int *island_nodes) -{ - *(island_nodes + (*island_size)) = n; - (*island_size)++; - - struct Node *pn = &G->nodes[n]; - pn->mark = 1; - - for (int i = 0; i < pn->degree; i++) - { - if (x[pn->adj[i].edge_index] > LP_EPSILON) - { - int neighbor = pn->adj[i].neighbor_index; - - if (G->nodes[neighbor].mark == 0) - graph_dfs(neighbor, G, x, island_size, island_nodes); - } - } -} - -void get_delta( - int island_node_count, - int *island_nodes, - int edge_count, - int *edges, - int *delta_count, - int *delta, - int *marks) -{ - for (int i = 0; i < island_node_count; i++) - marks[island_nodes[i]] = 1; - - int k = 0; - for (int i = 0; i < edge_count; i++) - if (marks[edges[2 * i]] + marks[edges[2 * i + 1]] == 1) - delta[k++] = i; - - *delta_count = k; - - for (int i = 0; i < island_node_count; i++) - marks[island_nodes[i]] = 0; -} - int get_cut_edges_from_marks( struct Graph *graph, int *cut_edges_count, struct Edge **cut_edges) { @@ -173,7 +128,7 @@ int get_cut_edges_from_marks( int graph_dump(const struct Graph *graph) { (void) graph; - #if LOG_LEVEL >= LOG_LEVEL_DEBUG +#if LOG_LEVEL >= LOG_LEVEL_DEBUG log_debug("node_count: %d edge_count: %d\n", graph->node_count, graph->edge_count); diff --git a/src/graph.h b/src/graph.h index 9b00587..0b8f8b8 100644 --- a/src/graph.h +++ b/src/graph.h @@ -47,9 +47,6 @@ struct Graph struct Adjacency *adj; }; -void graph_dfs( - int n, struct Graph *G, double *x, int *icount, int *island); - void graph_init(struct Graph *graph); void graph_free(struct Graph *graph); @@ -61,15 +58,6 @@ int graph_build( int is_directed, struct Graph *graph); -void get_delta( - int nsize, - int *nlist, - int ecount, - int *elist, - int *deltacount, - int *delta, - int *marks); - int get_cut_edges_from_marks( struct Graph *graph, int *cut_edges_count, struct Edge **cut_edges); diff --git a/src/gtsp-comb.c b/src/gtsp-comb.c index fb9ec44..de96883 100644 --- a/src/gtsp-comb.c +++ b/src/gtsp-comb.c @@ -2,7 +2,7 @@ #include "util.h" #include -int add_comb_cut( +static int add_comb_cut( struct LP *lp, struct Graph *graph, int current_component, @@ -73,7 +73,7 @@ int add_comb_cut( { for (int i = 0; i < nz; i++) { - if (OPTIMAL_X[rmatind[i]] < LP_EPSILON) continue; + if (OPTIMAL_X[rmatind[i]] < EPSILON) continue; if (rmatind[i] >= node_count) { @@ -98,7 +98,7 @@ int add_comb_cut( for (int i = 0; i < nz; i++) sum += rmatval[i] * OPTIMAL_X[rmatind[i]]; log_verbose("%.2lf >= %.2lf\n", sum, rhs); - abort_if(sum <= rhs - LP_EPSILON, "cannot add invalid cut"); + abort_if(sum <= rhs - EPSILON, "cannot add invalid cut"); } double lhs = 0.0; @@ -107,7 +107,7 @@ int add_comb_cut( log_verbose("Violation: %.4lf >= %.4lf\n", lhs, rhs); - if (lhs + LP_EPSILON > rhs) + if (lhs + EPSILON > rhs) { free(rmatind); free(rmatval); @@ -130,7 +130,7 @@ int add_comb_cut( return rval; } -int find_components( +static int find_components( struct Graph *graph, double *x, int *components, int *component_sizes) { int rval = 0; @@ -168,8 +168,8 @@ int find_components( if (neighbor->mark) continue; double x_e = x[adj->edge->index]; - if (x_e < LP_EPSILON) continue; - if (x_e > 1 - LP_EPSILON) continue; + if (x_e < EPSILON) continue; + if (x_e > 1 - EPSILON) continue; stack[stack_top++] = neighbor; neighbor->mark = 1; @@ -196,7 +196,7 @@ int find_components( return rval; } -int find_teeth( +static int find_teeth( struct Graph *graph, double *x, int current_component, @@ -221,7 +221,7 @@ int find_teeth( struct Node *from = e->from; struct Node *to = e->to; - if (x[e->index] < 1 - LP_EPSILON) continue; + if (x[e->index] < 1 - EPSILON) continue; if (to->mark || from->mark) continue; @@ -242,11 +242,6 @@ int find_teeth( return 0; } -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, @@ -346,10 +341,53 @@ static int shrink_clusters( return rval; } -int find_comb_cuts(struct LP *lp, struct GTSP *data) +static 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] > 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] > 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; +} + +int GTSP_find_comb_cuts(struct LP *lp, struct GTSP *data) { int rval = 0; - double initial_time = get_current_time(); + double initial_time = get_user_time(); double *x = 0; double *shrunken_x = 0; @@ -380,6 +418,8 @@ int find_comb_cuts(struct LP *lp, struct GTSP *data) #if LOG_LEVEL >= LOG_LEVEL_DEBUG rval = write_shrunken_graph(shrunken_x, &shrunken_graph, cluster_count); abort_if(rval, "write_shrunken_graph failed"); +#else + UNUSED(write_shrunken_graph); #endif teeth = (int *) malloc(cluster_count * sizeof(int)); @@ -414,8 +454,8 @@ int find_comb_cuts(struct LP *lp, struct GTSP *data) if (tooth_count % 2 == 0) { - for (int i = 0; i < cluster_count; i++) - if (teeth[i] == tooth_count - 1) teeth[i] = -1; + for (int k = 0; k < cluster_count; k++) + if (teeth[k] == tooth_count - 1) teeth[k] = -1; tooth_count--; } @@ -428,7 +468,7 @@ int find_comb_cuts(struct LP *lp, struct GTSP *data) int added_cuts_count = lp->cut_pool_size - original_cut_pool_size; log_debug(" %d combs\n", added_cuts_count); - COMBS_TIME += get_current_time() - initial_time; + COMBS_TIME += get_user_time() - initial_time; COMBS_COUNT += added_cuts_count; CLEANUP: @@ -440,48 +480,3 @@ int find_comb_cuts(struct LP *lp, struct GTSP *data) if (x) free(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; -} diff --git a/src/gtsp-comb.h b/src/gtsp-comb.h index ee600e8..eb12cce 100644 --- a/src/gtsp-comb.h +++ b/src/gtsp-comb.h @@ -1,6 +1,6 @@ #ifndef PROJECT_GTSP_COMB_H #define PROJECT_GTSP_COMB_H -int find_comb_cuts(struct LP *lp, struct GTSP *data); +int GTSP_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 ffe1adb..1ae2aa6 100644 --- a/src/gtsp-subtour.c +++ b/src/gtsp-subtour.c @@ -1,12 +1,32 @@ #include "gtsp-subtour.h" #include -#include #include "util.h" #include "flow.h" -extern double SUBTOUR_TIME; +static 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; + } +} + +static 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; -int static build_flow_digraph( + capacities[e->index] = 1e10; + capacities[e->reverse->index] = 1e10; + } +} + +static int build_flow_digraph( struct GTSP *data, double *x, struct Graph *digraph, double *capacities) { int rval = 0; @@ -16,9 +36,8 @@ int static build_flow_digraph( 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; + int digraph_node_count = node_count + data->cluster_count; + int digraph_edge_count = 4 * graph->edge_count + 2 * graph->node_count; digraph_edges = (int *) malloc(2 * digraph_edge_count * sizeof(int)); abort_if(!digraph_edges, "could not allocate digraph_edges"); @@ -28,7 +47,7 @@ int static build_flow_digraph( int kc = 0; for (int i = 0; i < graph->edge_count; i++) { - if (x[i] < LP_EPSILON) continue; + if (x[i] < EPSILON) continue; struct Edge *e = &graph->edges[i]; int from = e->from->index; @@ -67,19 +86,6 @@ int static build_flow_digraph( 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); @@ -99,10 +105,8 @@ int static build_flow_digraph( return rval; } -int add_subtour_cut( - struct LP *lp, - struct Edge **cut_edges, - int cut_edges_count) +static int add_subtour_cut( + struct LP *lp, struct Edge **cut_edges, int cut_edges_count) { int rval = 0; @@ -135,7 +139,7 @@ int add_subtour_cut( 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"); + abort_if(sum <= rhs - EPSILON, "cannot add invalid cut"); } struct Row *cut = 0; @@ -155,19 +159,27 @@ int add_subtour_cut( return rval; } -int find_exact_subtour_cuts( - struct LP *lp, struct GTSP *data, double min_cut_violation) +int GTSP_find_exact_subtour_cuts(struct LP *lp, struct GTSP *data) { int rval = 0; double *x = 0; + double *flow = 0; double *capacities = 0; - double initial_time = get_current_time(); + struct Edge **cut_edges = 0; + + double initial_time = get_user_time(); int added_cuts_count = 0; struct Graph *graph = data->graph; + const int edge_count = graph->edge_count; + const int node_count = graph->node_count; + const int cluster_count = data->cluster_count; int num_cols = LP_get_num_cols(lp); + int digraph_edge_count = 4 * edge_count + 2 * node_count; + int original_cut_pool_size = lp->cut_pool_size; + x = (double *) malloc(num_cols * sizeof(double)); abort_if(!x, "could not allocate x"); @@ -181,78 +193,36 @@ 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; - - int original_cut_pool_size = lp->cut_pool_size; + flow = (double *) malloc(digraph_edge_count * sizeof(double)); capacities = (double *) malloc(digraph_edge_count * sizeof(double)); + cut_edges = (struct Edge **) malloc(edge_count * sizeof(struct Edge *)); + + abort_if(!flow, "could not allocate flow"); abort_if(!capacities, "could not allocate capacities"); + abort_if(!cut_edges, "could not allocate cut_edges"); rval = build_flow_digraph(data, x, &digraph, capacities); abort_if(rval, "build_flow_digraph failed"); - rval = find_exact_subtour_cuts_cluster_to_cluster(lp, data, &digraph, - 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; - log_debug(" %d cluster-to-cluster\n", added_cuts_count); - SUBTOUR_CLUSTER_CLUSTER_COUNT += added_cuts_count; - if (added_cuts_count > 0) - goto CLEANUP; - - SUBTOUR_TIME += get_current_time() - initial_time; - - CLEANUP: - graph_free(&digraph); - if (capacities) free(capacities); - if (x) free(x); - return rval; -} - -int find_exact_subtour_cuts_cluster_to_cluster( - struct LP *lp, - struct GTSP *data, - struct Graph *digraph, - double *capacities, - double min_cut_violation) -{ - 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 i = 0; i < cluster_count; i++) { - for (int j = i + 1; j < data->cluster_count; j++) + for (int j = i + 1; j < cluster_count; j++) { - struct Node *from = &digraph->nodes[graph->node_count + i]; - struct Node *to = &digraph->nodes[graph->node_count + j]; + struct Node *from = &digraph.nodes[node_count + i]; + struct Node *to = &digraph.nodes[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, + rval = flow_find_max_flow(&digraph, capacities, from, to, flow, &flow_value); abort_if(rval, "flow_find_max_flow failed"); @@ -262,13 +232,13 @@ int find_exact_subtour_cuts_cluster_to_cluster( deactivate_cluster_node(capacities, from); deactivate_cluster_node(capacities, to); - if (flow_value >= 2 - min_cut_violation) continue; + if (flow_value >= 2 - EPSILON) continue; log_verbose("Marked nodes:\n"); - for (int k = 0; k < graph->node_count; k++) + for (int k = 0; k < node_count; k++) { - graph->nodes[k].mark = digraph->nodes[k].mark; - if (digraph->nodes[k].mark) log_verbose(" %d\n", 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); @@ -286,31 +256,17 @@ int find_exact_subtour_cuts_cluster_to_cluster( } } + added_cuts_count = lp->cut_pool_size - original_cut_pool_size; + log_debug(" %d cluster-to-cluster\n", added_cuts_count); + + SUBTOUR_COUNT += added_cuts_count; + SUBTOUR_TIME += get_user_time() - initial_time; + CLEANUP: + graph_free(&digraph); + if (capacities) free(capacities); if (cut_edges) free(cut_edges); if (flow) free(flow); + if (x) free(x); 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 index 42cfef1..3c61697 100644 --- a/src/gtsp-subtour.h +++ b/src/gtsp-subtour.h @@ -3,36 +3,6 @@ #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, - double min_cut_violation); - -int find_exact_subtour_cuts_node_to_cluster( - struct LP *lp, - struct GTSP *data, - double *x, - struct Graph *digraph, - double *capacities, - double min_cut_violation); - -int find_exact_subtour_cuts_node_to_node( - struct LP *lp, - struct GTSP *data, - double *x, - struct Graph *digraph, - double *capacities, - double min_cut_violation); - -int find_exact_subtour_cuts( - struct LP *lp, - struct GTSP *data, - double min_cut_violation); +int GTSP_find_exact_subtour_cuts(struct LP *lp, struct GTSP *data); #endif //_PROJECT_GTSP_SUBTOUR_H_ diff --git a/src/gtsp.c b/src/gtsp.c index ececd15..ae5b390 100644 --- a/src/gtsp.c +++ b/src/gtsp.c @@ -1,17 +1,16 @@ #include #include -#include -#include #include -#include -#include #include "gtsp.h" #include "geometry.h" #include "util.h" #include "gtsp-subtour.h" #include "gtsp-comb.h" -int large_neighborhood_search(struct Tour *tour, struct GTSP *data, int *tour_cost); +int large_neighborhood_search( + struct Tour *tour, + struct GTSP *data, + int *tour_cost); int build_edge_map(struct GTSP *gtsp, int *edge_map); @@ -231,7 +230,8 @@ int GTSP_add_cutting_planes(struct LP *lp, struct GTSP *data) } current_round++; - abort_if(get_current_time() - INITIAL_TIME >= MAX_TOTAL_TIME, "time limit exceeded"); + abort_if(get_user_time() - INITIAL_TIME >= MAX_TOTAL_TIME, + "time limit exceeded"); int original_cut_pool_size; int added_cuts_count; @@ -239,24 +239,22 @@ int GTSP_add_cutting_planes(struct LP *lp, struct GTSP *data) original_cut_pool_size = lp->cut_pool_size; log_debug("Finding subtour cuts, round %d...\n", current_round); - rval = find_exact_subtour_cuts(lp, data, LP_EPSILON); - abort_if(rval, "find_exact_subtour_cuts failed"); + rval = GTSP_find_exact_subtour_cuts(lp, data); + abort_if(rval, "GTSP_find_exact_subtour_cuts failed"); added_cuts_count = lp->cut_pool_size - original_cut_pool_size; if (added_cuts_count > 0) continue; -#ifdef ENABLE_COMB_INEQUALITIES original_cut_pool_size = lp->cut_pool_size; log_debug("Finding comb cuts, round %d...\n", current_round); - rval = find_comb_cuts(lp, data); - abort_if(rval, "find_comb_cuts failed"); + rval = GTSP_find_comb_cuts(lp, data); + abort_if(rval, "GTSP_find_comb_cuts failed"); added_cuts_count = lp->cut_pool_size - original_cut_pool_size; if (added_cuts_count > 0) continue; -#endif break; } @@ -303,7 +301,7 @@ int GTSP_write_solution(struct GTSP *data, char *filename, double *x) int positive_edge_count = 0; for (int i = 0; i < edge_count; i++) - if (x[i] > LP_EPSILON) + if (x[i] > EPSILON) positive_edge_count++; fprintf(file, "%d %d\n", node_count, edge_count); @@ -311,7 +309,7 @@ int GTSP_write_solution(struct GTSP *data, char *filename, double *x) fprintf(file, "%d\n", positive_edge_count); for (int i = 0; i < edge_count; i++) - if (x[i] > LP_EPSILON) + if (x[i] > EPSILON) fprintf(file, "%d %d %.4lf\n", edges[i].from->index, edges[i].to->index, x[i]); @@ -368,7 +366,7 @@ int GTSP_read_solution(struct GTSP *gtsp, char *filename, double **p_x) for (int i = 0; i < edge_count; i++) { - if (x[i] <= LP_EPSILON) continue; + if (x[i] <= EPSILON) continue; log_debug(" x%-5d = %.6f\n", i, x[i]); } @@ -420,10 +418,10 @@ int GTSP_check_solution(struct GTSP *data, double *x) for (int i = 0; i < edge_count; i++) { - abort_iff(x[i] < 1.0 - LP_EPSILON && x[i] > LP_EPSILON, + abort_iff(x[i] < 1.0 - EPSILON && x[i] > EPSILON, "solution is not integral: x%d = %.4lf", i, x[i]); - abort_iff(x[i] > 1.0 + LP_EPSILON || x[i] < 0.0 - LP_EPSILON, + abort_iff(x[i] > 1.0 + EPSILON || x[i] < 0.0 - EPSILON, "value out of bounds: x%d = %.4lf", i, x[i]); } @@ -435,7 +433,7 @@ int GTSP_check_solution(struct GTSP *data, double *x) int initial; for (initial = 0; initial < edge_count; initial++) - if (x[initial] > 1.0 - LP_EPSILON) break; + if (x[initial] > 1.0 - EPSILON) break; abort_if(initial == edge_count, "no initial node"); @@ -453,7 +451,7 @@ int GTSP_check_solution(struct GTSP *data, double *x) struct Node *neighbor = adj->neighbor; if (neighbor->mark) continue; - if (x[adj->edge->index] < LP_EPSILON) continue; + if (x[adj->edge->index] < EPSILON) continue; stack[stack_top++] = neighbor; neighbor->mark = 1; @@ -479,12 +477,12 @@ int GTSP_solution_found(struct BNC *bnc, struct GTSP *data, double *x) UNUSED(bnc); int rval = 0; - int tour_cost; + int tour_cost; char filename[100]; - double *best_val = &bnc->best_obj_val; - - struct Tour* tour; - tour = (struct Tour*) malloc(data->cluster_count*sizeof(struct Tour)); + double *best_val = &bnc->best_obj_val; + + struct Tour *tour; + tour = (struct Tour *) malloc(data->cluster_count * sizeof(struct Tour)); sprintf(filename, "tmp/gtsp-m%d-n%d-s%d.out", data->cluster_count, data->graph->node_count, SEED); @@ -493,316 +491,25 @@ int GTSP_solution_found(struct BNC *bnc, struct GTSP *data, double *x) rval = GTSP_write_solution(data, filename, x); abort_if(rval, "GTSP_write_solution failed"); - rval = build_tour_from_x(data,tour,x); - abort_if(rval, "build_tour_from_x failed"); - - rval = large_neighborhood_search(tour, data, &tour_cost); - abort_if(rval, "large_neighborhood_search failed"); - - if(tour_cost + LP_EPSILON < *best_val){ - log_info("Local search improve the integral solution\n"); - log_info(" obj val = %f\n",*best_val ); - log_info(" after LNS = %d\n",tour_cost ); - *best_val = tour_cost; - } - - log_info("Checking solution...\n"); - rval = GTSP_check_solution(data, x); - abort_if(rval, "GTSP_check_solution failed"); - - 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 = -1; -static int input_cluster_count = -1; -static int grid_size = 100; - -static char input_x_filename[1000] = {0}; - -static void GTSP_print_usage() -{ - printf("Parameters:\n"); - printf("%4s %-13s %s\n", "-n", "--nodes", "number of nodes"); - 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"); - printf("%4s %-13s %s\n", "-x", "--optimal", - "file containg valid solution (used to assert validity of cuts)"); -} - -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': - strcpy(input_x_filename, optarg); - break; - - case 's': - SEED = (unsigned) atoi(optarg); - break; - - case ':': - fprintf(stderr, "option '-%c' requires an argument\n", optopt); - rval = 1; - goto CLEANUP; - - case '?': - default: - fprintf(stderr, "option '-%c' is invalid\n", optopt); - rval = 1; - goto CLEANUP; - - } - } - - if (input_cluster_count < 0) - { - input_cluster_count = (int) ceil(input_node_count / 5.0); - if (input_cluster_count < 3) input_cluster_count = 3; - } - - if (input_node_count < 0) - { - printf("You must specify the number of nodes.\n"); - rval = 1; - } - - if (input_cluster_count > input_node_count) - { - printf("Number of clusters must be at most number of nodes.\n"); - rval = 1; - } - - if (rval) - { - GTSP_print_usage(); - rval = 1; - } - - CLEANUP: - return rval; -} - -int GTSP_main(int argc, char **argv) -{ - int rval = 0; - - struct BNC bnc; - struct GTSP data; - - double *initial_x = 0; - double initial_time = get_current_time(); - - SEED = (unsigned int) get_real_time() % 1000; - - 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); - - char instance_name[1000]; - sprintf(instance_name, "gtsp-m%d-n%d-s%d", input_cluster_count, - input_node_count, SEED); - - rval = GTSP_create_random_problem(input_node_count, input_cluster_count, - grid_size, &data); - abort_if(rval, "GTSP_create_random_problem failed"); - - char filename[100]; - sprintf(filename, "input/%s.in", instance_name); - log_info("Writing random instance to file %s\n", filename); - rval = GTSP_write_problem(&data, filename); - abort_if(rval, "GTSP_write_problem failed"); - -#if LOG_LEVEL >= LOG_LEVEL_DEBUG - log_info("Writing random instance to file gtsp.in\n"); - rval = GTSP_write_problem(&data, "gtsp.in"); - #endif - - int init_val = 0; -// - initial_x = (double *) malloc( - (data.graph->node_count + data.graph->edge_count) * sizeof(double)); - abort_if(!initial_x, "could not allocate initial_x"); - - rval = inital_tour_value(&data, &init_val, initial_x); - abort_if(rval, "initial_tour_value failed"); - - rval = GTSP_solution_found(&bnc, &data, initial_x); - abort_if(rval, "check_sol failed"); - - bnc.best_x = initial_x; - bnc.best_obj_val = init_val; - 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; - bnc.problem_solution_found = (int (*)( - struct BNC *, void *, double *)) GTSP_solution_found; - - double opt_val = 0.0; - - if (strlen(input_x_filename) == 0) - { - sprintf(input_x_filename, "optimal/%s.out", instance_name); - - FILE *file = fopen(input_x_filename, "r"); - - if (!file) - input_x_filename[0] = 0; - else - fclose(file); - } - - if (strlen(input_x_filename) > 0) - { - rval = GTSP_read_solution(&data, input_x_filename, &OPTIMAL_X); - abort_if(rval, "GTSP_read_solution failed"); - - log_info("Optimal solution is available. Cuts will be checked.\n"); - - for (int i = 0; i < data.graph->edge_count; i++) - { - struct Edge *e = &data.graph->edges[i]; - opt_val += OPTIMAL_X[i] * e->weight; - } - - log_info(" opt = %.2lf\n", opt_val); - } - - 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); - - if (OPTIMAL_X) - { - abort_iff(bnc.best_obj_val - LP_EPSILON > opt_val, - "Solution is not optimal: %.4lf > %.4lf", bnc.best_obj_val, - opt_val); - } - - TOTAL_TIME = get_current_time() - initial_time; - - log_info("Branch-and-bound nodes: %d\n", BNC_NODE_COUNT); - log_info("LP optimize calls: %d\n", LP_SOLVE_COUNT); - log_info("LP solving time: %.2lf\n", LP_SOLVE_TIME); - log_info("LP cut pool management time: %.2lf\n", CUT_POOL_TIME); - - FILE *file = fopen("stats.tab", "a"); - abort_if(!file, "could not open stats.tab"); - + rval = build_tour_from_x(data, tour, x); + abort_if(rval, "build_tour_from_x failed"); - struct stat st; - stat("stats.tab", &st); + rval = large_neighborhood_search(tour, data, &tour_cost); + abort_if(rval, "large_neighborhood_search failed"); - if(st.st_size == 0) + if (tour_cost + EPSILON < *best_val) { - fprintf(file, "%-20s ", "instance"); - fprintf(file, "%-8s ", "time"); - fprintf(file, "%-8s ", "subt-t"); - fprintf(file, "%-8s ", "combs-t"); - fprintf(file, "%-8s ", "pool-t"); - fprintf(file, "%-8s ", "pool-m"); - fprintf(file, "%-8s ", "lp-count"); - fprintf(file, "%-8s ", "lp-time"); - fprintf(file, "%-8s ", "lp-rows"); - fprintf(file, "%-8s ", "lp-cols"); - fprintf(file, "%-8s ", "init-v"); - fprintf(file, "%-8s ", "opt-v"); - fprintf(file, "%-8s ", "root-v"); - fprintf(file, "%-8s ", "nodes"); - fprintf(file, "%-8s ", "subt-cc"); - fprintf(file, "%-8s ", "subt-nc"); - fprintf(file, "%-8s ", "subt-nn"); - fprintf(file, "%-8s ", "combs"); - - fprintf(file, "\n"); + log_info("Local search improve the integral solution\n"); + log_info(" obj val = %f\n", *best_val); + log_info(" after LNS = %d\n", tour_cost); + *best_val = tour_cost; } - fprintf(file, "%-20s ", instance_name); - fprintf(file, "%-8.2lf ", TOTAL_TIME); - fprintf(file, "%-8.2lf ", SUBTOUR_TIME); - fprintf(file, "%-8.2lf ", COMBS_TIME); - fprintf(file, "%-8.2lf ", CUT_POOL_TIME); - fprintf(file, "%-8ld ", CUT_POOL_MAX_MEMORY/1024/1024); - fprintf(file, "%-8d ", LP_SOLVE_COUNT); - fprintf(file, "%-8.2lf ", LP_SOLVE_TIME); - fprintf(file, "%-8d ", LP_MAX_ROWS); - fprintf(file, "%-8d ", LP_MAX_COLS); - fprintf(file, "%-8d ", init_val); - fprintf(file, "%-8.0lf ", bnc.best_obj_val); - fprintf(file, "%-8.0lf ", ROOT_VALUE); - fprintf(file, "%-8d ", BNC_NODE_COUNT); - fprintf(file, "%-8d ", SUBTOUR_CLUSTER_CLUSTER_COUNT); - fprintf(file, "%-8d ", SUBTOUR_NODE_CLUSTER_COUNT); - fprintf(file, "%-8d ", SUBTOUR_NODE_NODE_COUNT); - fprintf(file, "%-8d ", COMBS_COUNT); - fprintf(file, "\n"); - fclose(file); + log_info("Checking solution...\n"); + rval = GTSP_check_solution(data, x); + abort_if(rval, "GTSP_check_solution failed"); CLEANUP: - if (OPTIMAL_X) free(OPTIMAL_X); - GTSP_free(&data); - BNC_free(&bnc); return rval; } @@ -823,11 +530,11 @@ int build_x_from_tour(struct GTSP *data, struct Tour *tour, double *x) for (int i = 0; i < edge_count; i++) x[i] = 0.0; - int next_vertex = tour[0].next; - int current_vertex = tour[0].vertex; + int next_vertex = tour[0].next; + int current_vertex = tour[0].vertex; for (int i = 0; i < data->cluster_count; i++) { - int from = current_vertex; + int from = current_vertex; int to = tour[next_vertex].vertex; current_vertex = tour[next_vertex].vertex; next_vertex = tour[next_vertex].next; @@ -843,7 +550,7 @@ int build_x_from_tour(struct GTSP *data, struct Tour *tour, double *x) int inital_tour_value(struct GTSP *data, int *tour_cost, double *x) { int rval = 0; - + int cluster_count = data->cluster_count; struct Tour *tour = 0; @@ -864,82 +571,94 @@ int inital_tour_value(struct GTSP *data, int *tour_cost, double *x) if (data->node_to_cluster[0] != i) { uncovered_sets[cluster_num] = i; - cluster_num ++; + cluster_num++; } } int new_vertex = 1, cost; - + tour[0].vertex = 0; cluster_in_tour[0] = 1; - + while (new_vertex < data->cluster_count) { int min_dist_vertex = -1; int min_cost = INT_MAX; - + for (int i = 1; i < data->graph->node_count; i++) { if (cluster_in_tour[data->node_to_cluster[i]]) continue; - - for (int k = 0; k < new_vertex; k++) - { - cost = data->dist_matrix[i][tour[k].vertex]; - if (cost < min_cost) - { min_cost = cost; - min_dist_vertex = i; - } - } + + for (int k = 0; k < new_vertex; k++) + { + cost = data->dist_matrix[i][tour[k].vertex]; + if (cost < min_cost) + { + min_cost = cost; + min_dist_vertex = i; + } + } } - + assert(min_dist_vertex >= 0); - int cluster_to_insert = data->node_to_cluster[min_dist_vertex]; - cluster_in_tour[cluster_to_insert] = 1; - - min_cost = INT_MAX; - - int insertion_cost = -1, best_pose = -1, best_vertex = -1; - - for (int k=0 ; k < data->clusters[cluster_to_insert].size; k++){ - for(int j = 0; j < new_vertex; j++){ - int vertex_to_insert = data->clusters[cluster_to_insert].nodes[k]; - if (new_vertex ==1){ - int vertex1 = tour[0].vertex; - insertion_cost = data->dist_matrix[vertex_to_insert][vertex1] - + data->dist_matrix[vertex1][vertex_to_insert]; - }else{ - int vertex1 = tour[j].vertex; - int vertex2 = tour[tour[j].next].vertex; - insertion_cost = data->dist_matrix[vertex1][vertex_to_insert] - + data->dist_matrix[vertex_to_insert][vertex2] - - data->dist_matrix[vertex1][vertex2]; - } - if(insertion_cost < min_cost){ - min_cost = insertion_cost; - best_pose = j; - best_vertex = vertex_to_insert; - } - } - } - - tour[new_vertex].vertex = best_vertex; - tour[new_vertex].prev = best_pose; - - if (new_vertex ==1){ - tour[new_vertex].next = best_pose; - tour[best_pose].prev = new_vertex; - }else{ - int temp_vertex = tour[best_pose].next; - tour[new_vertex].next = temp_vertex; - tour[temp_vertex].prev = new_vertex; - } - tour[best_pose].next = new_vertex; - + int cluster_to_insert = data->node_to_cluster[min_dist_vertex]; + cluster_in_tour[cluster_to_insert] = 1; + + min_cost = INT_MAX; + + int insertion_cost = -1, best_pose = -1, best_vertex = -1; + + for (int k = 0; k < data->clusters[cluster_to_insert].size; k++) + { + for (int j = 0; j < new_vertex; j++) + { + int vertex_to_insert = data->clusters[cluster_to_insert] + .nodes[k]; + if (new_vertex == 1) + { + int vertex1 = tour[0].vertex; + insertion_cost = + data->dist_matrix[vertex_to_insert][vertex1] + + data->dist_matrix[vertex1][vertex_to_insert]; + } + else + { + int vertex1 = tour[j].vertex; + int vertex2 = tour[tour[j].next].vertex; + insertion_cost = + data->dist_matrix[vertex1][vertex_to_insert] + + data->dist_matrix[vertex_to_insert][vertex2] - + data->dist_matrix[vertex1][vertex2]; + } + if (insertion_cost < min_cost) + { + min_cost = insertion_cost; + best_pose = j; + best_vertex = vertex_to_insert; + } + } + } + + tour[new_vertex].vertex = best_vertex; + tour[new_vertex].prev = best_pose; + + if (new_vertex == 1) + { + tour[new_vertex].next = best_pose; + tour[best_pose].prev = new_vertex; + } + else + { + int temp_vertex = tour[best_pose].next; + tour[new_vertex].next = temp_vertex; + tour[temp_vertex].prev = new_vertex; + } + tour[best_pose].next = new_vertex; + new_vertex += 1; } tour[data->cluster_count].vertex = 0; - rval = large_neighborhood_search(tour, data, tour_cost); abort_if(rval, "large_neighborhood_search failed"); @@ -1027,7 +746,10 @@ int two_opt(struct Tour *tour, struct GTSP *data) return 0; } -int large_neighborhood_search(struct Tour *tour, struct GTSP *data, int *tour_cost) +int large_neighborhood_search( + struct Tour *tour, + struct GTSP *data, + int *tour_cost) { int rval = 0; @@ -1035,15 +757,14 @@ int large_neighborhood_search(struct Tour *tour, struct GTSP *data, int *tour_co int *clusters = data->node_to_cluster; int **dist_matrix = data->dist_matrix; struct Cluster *vertex_set = data->clusters; - + //LNS starts for (int iter = 0; iter < 2000; iter++) { //Delete a random vertex - + int delete_vertex = rand() % (cluster_count - 1) + 1; - - + int prev_vertex = tour[delete_vertex].prev; int next_vertex = tour[delete_vertex].next; @@ -1090,19 +811,19 @@ int large_neighborhood_search(struct Tour *tour, struct GTSP *data, int *tour_co tour[delete_vertex].next = next_vertex; tour[best_pose].next = delete_vertex; tour[next_vertex].prev = delete_vertex; - + rval = optimize_vertex_in_cluster(tour, data); abort_if(rval, "optimize_vertex_in_cluster failed"); } - - rval = K_opt(tour, data); + + rval = K_opt(tour, data); abort_if(rval, "two_opt failed"); - + rval = two_opt(tour, data); abort_if(rval, "two_opt failed"); - + *tour_cost = list_length(tour, data); - + CLEANUP: //if (vertex_seq) free(vertex_seq); return rval; @@ -1158,67 +879,69 @@ void print_list(struct Tour *tour, struct GTSP *data) printf("\n"); } -int K_opt(struct Tour* tour, struct GTSP *data) +int K_opt(struct Tour *tour, struct GTSP *data) { - - int cluster_count = data->cluster_count; - int l; - - for(int i = 0 ; i < cluster_count - 3; i++){ - int location_in_path = 0; - - for(int k = 0; k < i; k++) - location_in_path = tour[location_in_path].next; - - int current_vertex = tour[location_in_path].vertex; - - - for(int k = 3; k < cluster_count - i - 2; k++){ - int first_next_location = tour[location_in_path].next; - int first_next_vertex = tour[tour[location_in_path].next].vertex; - int next_location = tour[location_in_path].next; - for(l = 0; l < k; l++){ - next_location = tour[next_location].next; - } - int next_vertex = tour[next_location].vertex; - int next_next_location = tour[next_location].next; - - int next_next_vertex = tour[next_next_location].vertex; - - if(next_next_vertex == current_vertex) break; - - int cost1 = data->dist_matrix[current_vertex][first_next_vertex] + - data->dist_matrix[next_vertex][next_next_vertex]; - int cost2 = data->dist_matrix[current_vertex][next_vertex] + - data->dist_matrix[first_next_vertex][next_next_vertex]; - - if(cost2 < cost1) - { - int tmp_location = next_location; - for(int m =0 ; m < l + 1; m++) - { - int tmp_vertex = tour[tmp_location].next; - tour[tmp_location].next = tour[tmp_location].prev; - tour[tmp_location].prev = tmp_vertex; - tmp_location = tour[tmp_location].next; - } - - tour[location_in_path].next = next_location; - tour[next_location].prev = location_in_path; - - tour[first_next_location].next = next_next_location; - - tour[next_next_location].prev = first_next_location; - - } - } - } - return 0; + + int cluster_count = data->cluster_count; + int l; + + for (int i = 0; i < cluster_count - 3; i++) + { + int location_in_path = 0; + + for (int k = 0; k < i; k++) + location_in_path = tour[location_in_path].next; + + int current_vertex = tour[location_in_path].vertex; + + for (int k = 3; k < cluster_count - i - 2; k++) + { + int first_next_location = tour[location_in_path].next; + int first_next_vertex = tour[tour[location_in_path].next].vertex; + int next_location = tour[location_in_path].next; + for (l = 0; l < k; l++) + { + next_location = tour[next_location].next; + } + int next_vertex = tour[next_location].vertex; + int next_next_location = tour[next_location].next; + + int next_next_vertex = tour[next_next_location].vertex; + + if (next_next_vertex == current_vertex) break; + + int cost1 = data->dist_matrix[current_vertex][first_next_vertex] + + data->dist_matrix[next_vertex][next_next_vertex]; + int cost2 = data->dist_matrix[current_vertex][next_vertex] + + data->dist_matrix[first_next_vertex][next_next_vertex]; + + if (cost2 < cost1) + { + int tmp_location = next_location; + for (int m = 0; m < l + 1; m++) + { + int tmp_vertex = tour[tmp_location].next; + tour[tmp_location].next = tour[tmp_location].prev; + tour[tmp_location].prev = tmp_vertex; + tmp_location = tour[tmp_location].next; + } + + tour[location_in_path].next = next_location; + tour[next_location].prev = location_in_path; + + tour[first_next_location].next = next_next_location; + + tour[next_next_location].prev = first_next_location; + + } + } + } + return 0; } int build_tour_from_x(struct GTSP *data, struct Tour *tour, double *x) { - + int rval = 0; int *cluster_mark = 0; @@ -1228,7 +951,7 @@ int build_tour_from_x(struct GTSP *data, struct Tour *tour, double *x) struct Graph *graph = data->graph; const int node_count = graph->node_count; const int edge_count = graph->edge_count; - + cluster_mark = (int *) malloc(data->cluster_count * sizeof(int)); abort_if(!cluster_mark, "could not allocate cluster_mark"); @@ -1240,10 +963,10 @@ int build_tour_from_x(struct GTSP *data, struct Tour *tour, double *x) for (int i = 0; i < data->cluster_count; i++) cluster_mark[i] = 0; - + int initial; for (initial = 0; initial < edge_count; initial++) - if (x[initial] > 1.0 - LP_EPSILON) break; + if (x[initial] > 1.0 - EPSILON) break; initial = graph->edges[initial].from->index; @@ -1251,10 +974,10 @@ int build_tour_from_x(struct GTSP *data, struct Tour *tour, double *x) stack[stack_top++] = &graph->nodes[initial]; graph->nodes[initial].mark = 1; - - tour[0].vertex = graph->nodes[initial].index; - int next_vertex = 1; - + + tour[0].vertex = graph->nodes[initial].index; + int next_vertex = 1; + while (stack_top > 0) { struct Node *n = stack[--stack_top]; @@ -1266,30 +989,36 @@ int build_tour_from_x(struct GTSP *data, struct Tour *tour, double *x) struct Node *neighbor = adj->neighbor; if (neighbor->mark) continue; - if (x[adj->edge->index] < LP_EPSILON) continue; + if (x[adj->edge->index] < EPSILON) continue; stack[stack_top++] = neighbor; tour[next_vertex].vertex = neighbor->index; - next_vertex++; + next_vertex++; neighbor->mark = 1; - + + } + } + + for (int i = 0; i < data->cluster_count; i++) + { + if (i == 0) + { + tour[i].prev = data->cluster_count - 1; + } + else + { + tour[i].prev = i - 1; + } + if (i == data->cluster_count - 1) + { + tour[i].next = 0; + } + else + { + tour[i].next = i + 1; } } - - - for(int i =0; i < data->cluster_count; i++){ - if(i == 0){ - tour[i].prev = data->cluster_count - 1; - }else{ - tour[i].prev = i-1; - } - if(i == data->cluster_count - 1){ - tour[i].next = 0; - }else{ - tour[i].next = i+1; - } - } - + CLEANUP: if (stack) free(stack); if (cluster_mark) free(cluster_mark); diff --git a/src/gtsp.h b/src/gtsp.h index efa34ba..ec79110 100644 --- a/src/gtsp.h +++ b/src/gtsp.h @@ -46,8 +46,6 @@ int GTSP_write_problem(struct GTSP *data, char *filename); int GTSP_write_solution(struct GTSP *data, char *filename, double *x); -int GTSP_main(int argc, char **argv); - int optimize_vertex_in_cluster(struct Tour * tour, struct GTSP *data); int two_opt(struct Tour * tour, struct GTSP *data); @@ -64,4 +62,10 @@ void print_list(struct Tour * tour, struct GTSP* data); int build_tour_from_x(struct GTSP *data, struct Tour *tour, double *x); +int GTSP_solution_found(struct BNC *bnc, struct GTSP *data, double *x); + +int GTSP_check_solution(struct GTSP *data, double *x); + +int GTSP_read_solution(struct GTSP *gtsp, char *filename, double **p_x); + #endif //_PROJECT_GTSP_H_ diff --git a/src/lp.c b/src/lp.c index df4b65c..906694a 100644 --- a/src/lp.c +++ b/src/lp.c @@ -7,6 +7,170 @@ #include "util.h" #include "main.h" +static int compress_cut_pool(struct LP *lp) +{ + int delete_count = 0; + for (int i = 0; i < lp->cut_pool_size; i++) + { + if (lp->cut_pool[i]->cplex_row_index < 0) + { + free(lp->cut_pool[i]->rmatind); + free(lp->cut_pool[i]->rmatval); + free(lp->cut_pool[i]); + delete_count++; + } + else + { + lp->cut_pool[i - delete_count] = lp->cut_pool[i]; + } + } + + lp->cut_pool_size -= delete_count; + + return 0; +} + +static int remove_old_cuts(struct LP *lp) +{ + int rval = 0; + + int *should_remove = 0; + + int numrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp); + log_verbose(" numrows=%d\n", numrows); + + should_remove = (int *) malloc((numrows + 1) * sizeof(int)); + abort_if(!should_remove, "could not allocate should_remove"); + + for (int i = 0; i < numrows; i++) + should_remove[i] = 0; + should_remove[numrows] = 0; + + log_verbose("Old cplex row index:\n"); + for (int i = 0; i < lp->cut_pool_size; i++) + log_verbose(" %d\n", lp->cut_pool[i]->cplex_row_index); + + log_verbose("Should remove:\n"); + for (int i = 0; i < lp->cut_pool_size; i++) + { + struct Row *cut = lp->cut_pool[i]; + if (cut->age <= MAX_CUT_AGE) continue; + if (cut->cplex_row_index < 0) continue; + + should_remove[cut->cplex_row_index] = 1; + log_verbose(" %d\n", cut->cplex_row_index); + } + + // Update cut->cplex_row_index + int count = 0; + for (int i = 0; i < numrows; i++) + { + if (!should_remove[i]) continue; + + for (int j = 0; j < lp->cut_pool_size; j++) + { + struct Row *cut = lp->cut_pool[j]; + + if (cut->cplex_row_index == i - count) cut->cplex_row_index = -1; + else if (cut->cplex_row_index > i - count) cut->cplex_row_index--; + } + + count++; + } + + log_verbose("New cplex row index:\n"); + for (int i = 0; i < lp->cut_pool_size; i++) + log_verbose(" %d\n", lp->cut_pool[i]->cplex_row_index); + + // Remove from CPLEX, finding the largest intervals of contiguous rows + // that should be removed. + int start = 0; + int end = -1; + count = 0; + for (int i = 0; i < numrows + 1; i++) + { + if (should_remove[i]) + end++; + else + { + if (end >= start) + { + rval = CPXdelrows(lp->cplex_env, lp->cplex_lp, start - count, + end - count); + abort_if(rval, "CPXdelrows failed"); + log_verbose(" %d %d (%d)\n", start, end, end - start + 1); + + count += end - start + 1; + } + + start = i + 1; + end = i; + } + } + + log_debug(" removed %d old cuts\n", count); + + rval = CPXdualopt(lp->cplex_env, lp->cplex_lp); + abort_if(rval, "CPXoptimize failed"); + + log_debug("Compressing cut pool...\n"); + compress_cut_pool(lp); + + + long nz = 0; + long size = 0; + for (int i = 0; i < lp->cut_pool_size; i++) + { + size += sizeof(struct Row); + nz += lp->cut_pool[i]->nz; + size += lp->cut_pool[i]->nz * sizeof(double); + size += lp->cut_pool[i]->nz * sizeof(int); + } + + log_debug(" %ld cuts (%ld nz, %ld MiB)\n", lp->cut_pool_size, nz, size/1024/1024); + + if(size > CUT_POOL_MAX_MEMORY) + CUT_POOL_MAX_MEMORY = size; + + CLEANUP: + if (should_remove) free(should_remove); + return rval; +} + +static int update_cut_ages(struct LP *lp) +{ + int rval = 0; + + double *slacks = 0; + + int numrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp); + + slacks = (double *) malloc(numrows * sizeof(double)); + abort_if(!slacks, "could not allocate slacks"); + + rval = CPXgetslack(lp->cplex_env, lp->cplex_lp, slacks, 0, numrows - 1); + abort_if(rval, "CPXgetslack failed"); + + int count = 0; + for (int i = 0; i < lp->cut_pool_size; i++) + { + struct Row *cut = lp->cut_pool[i]; + if (cut->cplex_row_index < 0) continue; + assert(cut->cplex_row_index < numrows); + + if (slacks[cut->cplex_row_index] < -EPSILON) + cut->age++; + else + cut->age = 0; + + if (cut->age > 5) count++; + } + + CLEANUP: + if (slacks) free(slacks); + return rval; +} + int LP_open(struct LP *lp) { int rval = 0; @@ -50,32 +214,6 @@ void LP_free(struct LP *lp) LP_free_cut_pool(lp); } -int LP_compress_cut_pool(struct LP *lp) -{ - int rval = 0; - - int delete_count = 0; - for (int i = 0; i < lp->cut_pool_size; i++) - { - if (lp->cut_pool[i]->cplex_row_index < 0) - { - free(lp->cut_pool[i]->rmatind); - free(lp->cut_pool[i]->rmatval); - free(lp->cut_pool[i]); - delete_count++; - } - else - { - lp->cut_pool[i - delete_count] = lp->cut_pool[i]; - } - } - - lp->cut_pool_size -= delete_count; - - CLEANUP: - return rval; -} - int LP_create(struct LP *lp, const char *name) { int rval = 0; @@ -98,7 +236,6 @@ int LP_new_row(struct LP *lp, char sense, double rhs) int rval = 0; rval = CPXnewrows(lp->cplex_env, lp->cplex_lp, 1, &rhs, &sense, 0, 0); - abort_if(rval, "CPXnewrows failed"); CLEANUP: @@ -119,7 +256,6 @@ int LP_add_rows( rval = CPXaddrows(lp->cplex_env, lp->cplex_lp, 0, newrows, newnz, rhs, sense, rmatbeg, rmatind, rmatval, 0, 0); - abort_if(rval, "CPXaddrows failed"); CLEANUP: @@ -159,46 +295,6 @@ int LP_change_bound(struct LP *lp, int col, char lower_or_upper, double bnd) return rval; } -extern int LP_SOLVE_COUNT; -extern double LP_SOLVE_TIME; -extern double CUT_POOL_TIME; - -int LP_update_cut_ages(struct LP *lp) -{ - int rval = 0; - - double *slacks = 0; - - int numrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp); - - slacks = (double *) malloc(numrows * sizeof(double)); - abort_if(!slacks, "could not allocate slacks"); - - rval = CPXgetslack(lp->cplex_env, lp->cplex_lp, slacks, 0, numrows - 1); - abort_if(rval, "CPXgetslack failed"); - - int count = 0; - for (int i = 0; i < lp->cut_pool_size; i++) - { - struct Row *cut = lp->cut_pool[i]; - if (cut->cplex_row_index < 0) continue; - assert(cut->cplex_row_index < numrows); - - if (slacks[cut->cplex_row_index] < -LP_EPSILON) - cut->age++; - else - cut->age = 0; - - if (cut->age > 5) count++; - } - - CLEANUP: - if (slacks) free(slacks); - return rval; -} - - - int LP_optimize(struct LP *lp, int *infeasible) { LP_SOLVE_COUNT++; @@ -215,12 +311,12 @@ int LP_optimize(struct LP *lp, int *infeasible) log_debug("Optimizing LP (%d rows %d cols)...\n", numrows, numcols); - double initial_time = get_current_time(); + double initial_time = get_user_time(); rval = CPXdualopt(lp->cplex_env, lp->cplex_lp); abort_if(rval, "CPXdualopt failed"); - LP_SOLVE_TIME += get_current_time() - initial_time; + LP_SOLVE_TIME += get_user_time() - initial_time; solstat = CPXgetstat(lp->cplex_env, lp->cplex_lp); if (solstat == CPX_STAT_INFEASIBLE) @@ -239,131 +335,22 @@ int LP_optimize(struct LP *lp, int *infeasible) abort_if(rval, "LP_get_obj_val failed"); log_debug(" obj val = %.4lf\n", objval); - log_debug(" time = %.4lf\n", get_current_time() - initial_time); + log_debug(" time = %.4lf\n", get_user_time() - initial_time); - initial_time = get_current_time(); - rval = LP_update_cut_ages(lp); - abort_if(rval, "LP_update_cut_ages failed"); + initial_time = get_user_time(); + rval = update_cut_ages(lp); + abort_if(rval, "update_cut_ages failed"); log_debug("Removing old cuts...\n"); - rval = LP_remove_old_cuts(lp); + rval = remove_old_cuts(lp); abort_if(rval, "LP_remove_old_cuts failed"); - CUT_POOL_TIME += get_current_time() - initial_time; + CUT_POOL_TIME += get_user_time() - initial_time; CLEANUP: return rval; } -int LP_remove_old_cuts(struct LP *lp) -{ - int rval = 0; - - int *should_remove = 0; - - int numrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp); - log_verbose(" numrows=%d\n", numrows); - - should_remove = (int *) malloc((numrows + 1) * sizeof(int)); - abort_if(!should_remove, "could not allocate should_remove"); - - for (int i = 0; i < numrows; i++) - should_remove[i] = 0; - should_remove[numrows] = 0; - - log_verbose("Old cplex row index:\n"); - for (int i = 0; i < lp->cut_pool_size; i++) - log_verbose(" %d\n", lp->cut_pool[i]->cplex_row_index); - - log_verbose("Should remove:\n"); - for (int i = 0; i < lp->cut_pool_size; i++) - { - struct Row *cut = lp->cut_pool[i]; - if (cut->age <= MAX_CUT_AGE) continue; - if (cut->cplex_row_index < 0) continue; - - should_remove[cut->cplex_row_index] = 1; - log_verbose(" %d\n", cut->cplex_row_index); - } - - int count = 0; - for (int i = 0; i < numrows; i++) - { - if (!should_remove[i]) continue; - - for (int j = 0; j < lp->cut_pool_size; j++) - { - struct Row *cut = lp->cut_pool[j]; - - if (cut->cplex_row_index == i - count) cut->cplex_row_index = -1; - else if (cut->cplex_row_index > i - count) cut->cplex_row_index--; - } - - count++; - } - - log_verbose("New cplex row index:\n"); - for (int i = 0; i < lp->cut_pool_size; i++) - log_verbose(" %d\n", lp->cut_pool[i]->cplex_row_index); - - int start = 0; - int end = -1; - count = 0; - for (int i = 0; i < numrows + 1; i++) - { - if (should_remove[i]) - { - end++; - } - else - { - if (end >= start) - { - rval = CPXdelrows(lp->cplex_env, lp->cplex_lp, start - count, - end - count); - abort_if(rval, "CPXdelrows failed"); - log_verbose(" %d %d (%d)\n", start, end, end - start + 1); - - count += end - start + 1; - } - - start = i + 1; - end = i; - } - } - - numrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp); - log_verbose(" new numrows=%d\n", numrows); - - log_debug(" removed %d old cuts\n", count); - - rval = CPXdualopt(lp->cplex_env, lp->cplex_lp); - abort_if(rval, "CPXoptimize failed"); - - log_debug("Compressing cut pool...\n"); - LP_compress_cut_pool(lp); - - - long nz = 0; - long size = 0; - for (int i = 0; i < lp->cut_pool_size; i++) - { - size += sizeof(struct Row); - nz += lp->cut_pool[i]->nz; - size += lp->cut_pool[i]->nz * sizeof(double); - size += lp->cut_pool[i]->nz * sizeof(int); - } - - log_debug(" %ld cuts (%ld nz, %ld MiB)\n", lp->cut_pool_size, nz, size/1024/1024); - - if(size > CUT_POOL_MAX_MEMORY) - CUT_POOL_MAX_MEMORY = size; - - CLEANUP: - if (should_remove) free(should_remove); - return rval; -} - int LP_get_obj_val(struct LP *lp, double *obj) { int rval = 0; @@ -414,8 +401,8 @@ int LP_write(struct LP *lp, const char *fname) if((a)>(b)) return 1; #define return_if_neq_epsilon(a, b) \ - if((a+LP_EPSILON)<(b)) return -1; \ - if((a-LP_EPSILON)>(b)) return 1; + if((a+EPSILON)<(b)) return -1; \ + if((a-EPSILON)>(b)) return 1; int compare_cuts(struct Row *cut1, struct Row *cut2) { @@ -430,12 +417,27 @@ int compare_cuts(struct Row *cut1, struct Row *cut2) return 0; } +static int update_hash(struct Row *cut) +{ + unsigned long hash = 0; + + for (int i = 0; i < cut->nz; i++) + { + hash += cut->rmatind[i] * 65521; + hash %= 4294967291; + } + + cut->hash = hash; + + return 0; +} + int LP_add_cut(struct LP *lp, struct Row *cut) { int rval = 0; - double initial_time = get_current_time(); + double initial_time = get_user_time(); - rval = LP_update_hash(cut); + rval = update_hash(cut); abort_if(rval, "LP_update_hash failed"); for (int i = 0; i < lp->cut_pool_size; i++) @@ -465,22 +467,7 @@ int LP_add_cut(struct LP *lp, struct Row *cut) cut->cplex_row_index = CPXgetnumrows(lp->cplex_env, lp->cplex_lp) - 1; cut->age = 0; - CUT_POOL_TIME += get_current_time() - initial_time; + CUT_POOL_TIME += get_user_time() - initial_time; CLEANUP: return rval; } - -int LP_update_hash(struct Row *cut) -{ - unsigned long hash = 0; - - for (int i = 0; i < cut->nz; i++) - { - hash += cut->rmatind[i] * 65521; - hash %= 4294967291; - } - - cut->hash = hash; - - return 0; -} \ No newline at end of file diff --git a/src/lp.h b/src/lp.h index 199f3c4..534f1fd 100644 --- a/src/lp.h +++ b/src/lp.h @@ -69,10 +69,6 @@ int LP_get_x(struct LP *lp, double *x); int LP_get_num_cols(struct LP *lp); -int LP_remove_old_cuts(struct LP *lp); - -int LP_update_hash(struct Row *cut); - int LP_add_cut(struct LP *lp, struct Row *cut); #endif \ No newline at end of file diff --git a/src/main.c b/src/main.c index b5dafbc..6653f85 100644 --- a/src/main.c +++ b/src/main.c @@ -1,14 +1,12 @@ #include +#include +#include +#include #include "main.h" #include "gtsp.h" -#include "tsp.h" -#include "flow.h" +#include "util.h" -char *INPUT_FILENAME = 0; unsigned int SEED = 0; -int GEOMETRIC_DATA = 0; -int NODE_COUNT_RAND = 0; -int GRID_SIZE_RAND = 100; double SUBTOUR_TIME = 0; double COMBS_TIME = 0; @@ -25,53 +23,306 @@ int LP_SOLVE_COUNT = 0; double TOTAL_TIME = 0; double ROOT_VALUE = 0; -int SUBTOUR_CLUSTER_CLUSTER_COUNT = 0; -int SUBTOUR_NODE_CLUSTER_COUNT = 0; -int SUBTOUR_NODE_NODE_COUNT = 0; +int SUBTOUR_COUNT = 0; int COMBS_COUNT = 0; double INITIAL_TIME = 0; -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 BNC_NODE_COUNT = 0; -void GTSP_print_usage() +static int input_node_count = -1; +static int input_cluster_count = -1; +static int grid_size = 100; + +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 char input_x_filename[1000] = {0}; + +static void print_usage(int argc, char **argv) { - printf("wrong usage\n"); + printf("Usage: %s [OPTION]...\n", argv[0]); + printf("Solves the Generalized Traveling Salesman problem for the input graph.\n\n"); + + printf("Parameters:\n"); + printf("%4s %-13s %s\n", "-n", "--nodes", "number of nodes"); + 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"); + printf("%4s %-13s %s\n", "-x", "--optimal", + "file containg valid solution (used to assert validity of cuts)"); +} + +static int 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:h:", 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': + strcpy(input_x_filename, optarg); + break; + + case 's': + SEED = (unsigned) atoi(optarg); + break; + + case 'h': + print_usage(argc, argv); + exit(0); + + case ':': + fprintf(stderr, "%s: option '-%c' requires an argument\n", + argv[0], optopt); + rval = 1; + goto CLEANUP; + + case '?': + default: + fprintf(stderr, "%s: option '-%c' is invalid\n", argv[0], + optopt); + rval = 1; + goto CLEANUP; + + } + } + + if (input_cluster_count < 0) + { + input_cluster_count = (int) ceil(input_node_count / 5.0); + if (input_cluster_count < 3) input_cluster_count = 3; + } + + if (input_node_count < 0) + { + fprintf(stderr, "You must specify the number of nodes.\n"); + rval = 1; + } + else if (input_cluster_count > input_node_count) + { + printf("Number of clusters must be at most number of nodes.\n"); + rval = 1; + } + + CLEANUP: + if (rval) + fprintf(stderr, "Try '%s --help' for more information.\n", argv[0]); + return rval; } int main(int argc, char **argv) { - int c = 0; - int option_index = 0; + int rval = 0; + + struct BNC bnc; + struct GTSP data; + + double *initial_x = 0; + double initial_time = get_user_time(); + + SEED = (unsigned int) get_real_time() % 1000; + + rval = GTSP_init_data(&data); + abort_if(rval, "GTSP_init_data failed"); - c = getopt_long(argc, argv, "htgf", options_tab, &option_index); + rval = BNC_init(&bnc); + abort_if(rval, "BNC_init failed"); - if (c < 0) + rval = 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); + + char instance_name[1000]; + sprintf(instance_name, "gtsp-m%d-n%d-s%d", input_cluster_count, + input_node_count, SEED); + + rval = GTSP_create_random_problem(input_node_count, input_cluster_count, + grid_size, &data); + abort_if(rval, "GTSP_create_random_problem failed"); + + char filename[100]; + sprintf(filename, "input/%s.in", instance_name); + log_info("Writing random instance to file %s\n", filename); + rval = GTSP_write_problem(&data, filename); + abort_if(rval, "GTSP_write_problem failed"); + +#if LOG_LEVEL >= LOG_LEVEL_DEBUG + log_info("Writing random instance to file gtsp.in\n"); + rval = GTSP_write_problem(&data, "gtsp.in"); +#endif + + int init_val = 0; + + initial_x = (double *) malloc( + (data.graph->node_count + data.graph->edge_count) * sizeof(double)); + abort_if(!initial_x, "could not allocate initial_x"); + + rval = inital_tour_value(&data, &init_val, initial_x); + abort_if(rval, "initial_tour_value failed"); + + rval = GTSP_solution_found(&bnc, &data, initial_x); + abort_if(rval, "check_sol failed"); + + bnc.best_x = initial_x; + bnc.best_obj_val = init_val; + 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; + bnc.problem_solution_found = (int (*)( + struct BNC *, void *, double *)) GTSP_solution_found; + + double opt_val = 0.0; + + if (strlen(input_x_filename) == 0) { - GTSP_print_usage(); - return 1; + sprintf(input_x_filename, "optimal/%s.out", instance_name); + + FILE *file = fopen(input_x_filename, "r"); + + if (!file) + input_x_filename[0] = 0; + else + fclose(file); } - switch (c) + if (strlen(input_x_filename) > 0) { - case 'f': - return flow_main(argc, argv); + rval = GTSP_read_solution(&data, input_x_filename, &OPTIMAL_X); + abort_if(rval, "GTSP_read_solution failed"); - case 't': - return TSP_main(argc, argv); + log_info("Optimal solution is available. Cuts will be checked.\n"); - case 'g': - return GTSP_main(argc, argv); + for (int i = 0; i < data.graph->edge_count; i++) + { + struct Edge *e = &data.graph->edges[i]; + opt_val += OPTIMAL_X[i] * e->weight; + } - default: - case 'h': - GTSP_print_usage(); - return 1; + log_info(" opt = %.2lf\n", opt_val); } + 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); + + if (OPTIMAL_X) + { + abort_iff(bnc.best_obj_val - EPSILON > opt_val, + "Solution is not optimal: %.4lf > %.4lf", bnc.best_obj_val, + opt_val); + } + + TOTAL_TIME = get_user_time() - initial_time; + + log_info("Branch-and-bound nodes: %d\n", BNC_NODE_COUNT); + log_info("LP optimize calls: %d\n", LP_SOLVE_COUNT); + log_info("LP solving time: %.2lf\n", LP_SOLVE_TIME); + log_info("LP cut pool management time: %.2lf\n", CUT_POOL_TIME); + + // Write statistics to a file + FILE *file = fopen("stats.tab", "a"); + abort_if(!file, "could not open stats.tab"); + + struct stat st; + stat("stats.tab", &st); + + if (st.st_size == 0) + { + fprintf(file, "%-20s ", "instance"); + fprintf(file, "%-8s ", "time"); + fprintf(file, "%-8s ", "subt-t"); + fprintf(file, "%-8s ", "combs-t"); + fprintf(file, "%-8s ", "pool-t"); + fprintf(file, "%-8s ", "pool-m"); + fprintf(file, "%-8s ", "lp-count"); + fprintf(file, "%-8s ", "lp-time"); + fprintf(file, "%-8s ", "lp-rows"); + fprintf(file, "%-8s ", "lp-cols"); + fprintf(file, "%-8s ", "init-v"); + fprintf(file, "%-8s ", "opt-v"); + fprintf(file, "%-8s ", "root-v"); + fprintf(file, "%-8s ", "nodes"); + fprintf(file, "%-8s ", "subt-cc"); + fprintf(file, "%-8s ", "subt-nc"); + fprintf(file, "%-8s ", "subt-nn"); + fprintf(file, "%-8s ", "combs"); + + fprintf(file, "\n"); + } + + fprintf(file, "%-20s ", instance_name); + fprintf(file, "%-8.2lf ", TOTAL_TIME); + fprintf(file, "%-8.2lf ", SUBTOUR_TIME); + fprintf(file, "%-8.2lf ", COMBS_TIME); + fprintf(file, "%-8.2lf ", CUT_POOL_TIME); + fprintf(file, "%-8ld ", CUT_POOL_MAX_MEMORY / 1024 / 1024); + fprintf(file, "%-8d ", LP_SOLVE_COUNT); + fprintf(file, "%-8.2lf ", LP_SOLVE_TIME); + fprintf(file, "%-8d ", LP_MAX_ROWS); + fprintf(file, "%-8d ", LP_MAX_COLS); + fprintf(file, "%-8d ", init_val); + fprintf(file, "%-8.0lf ", bnc.best_obj_val); + fprintf(file, "%-8.0lf ", ROOT_VALUE); + fprintf(file, "%-8d ", BNC_NODE_COUNT); + fprintf(file, "%-8d ", SUBTOUR_COUNT); + fprintf(file, "%-8d ", COMBS_COUNT); + fprintf(file, "\n"); + fclose(file); + + CLEANUP: + if (OPTIMAL_X) free(OPTIMAL_X); + GTSP_free(&data); + BNC_free(&bnc); + return rval; } + diff --git a/src/main.h b/src/main.h index 7a084b9..091191e 100644 --- a/src/main.h +++ b/src/main.h @@ -1,11 +1,7 @@ #ifndef _PROJECT_MAIN_H_ #define _PROJECT_MAIN_H_ -extern char *INPUT_FILENAME; extern unsigned int SEED; -extern int GEOMETRIC_DATA; -extern int NODE_COUNT_RAND; -extern int GRID_SIZE_RAND; extern double *OPTIMAL_X; extern double SUBTOUR_TIME; @@ -25,10 +21,7 @@ extern double TOTAL_TIME; extern double INITIAL_TIME; extern double ROOT_VALUE; -extern int SUBTOUR_CLUSTER_CLUSTER_COUNT; -extern int SUBTOUR_NODE_CLUSTER_COUNT; -extern int SUBTOUR_NODE_NODE_COUNT; - +extern int SUBTOUR_COUNT; extern int COMBS_COUNT; #endif \ No newline at end of file diff --git a/src/params.h b/src/params.h index 212d42b..052b8da 100644 --- a/src/params.h +++ b/src/params.h @@ -1,17 +1,31 @@ #ifndef PROJECT_PARAMS_H #define PROJECT_PARAMS_H -#define LP_EPSILON 0.000001 +/* + * Error margin for floating point comparisons. + */ +#define EPSILON 0.000001 +/* + * Cut pool settings. + */ #define MAX_CUT_AGE 5 #define MAX_CUT_POOL_SIZE 1000000 +/* + * Available log levels, in decreasing level of verboseness, are: + * + * LOG_LEVEL_VERBOSE + * LOG_LEVEL_DEBUG + * LOG_LEVEL_INFO + * LOG_LEVEL_WARNING + * LOG_LEVEL_ERROR + */ #define LOG_LEVEL LOG_LEVEL_INFO -#define ENABLE_COMB_INEQUALITIES - -#define MAX_TOTAL_TIME 999999 - -//#define ALLOW_FRACTIONAL_SOLUTIONS +/* + * Time limit for the computation (user time, in seconds). + */ +#define MAX_TOTAL_TIME 3600 #endif //PROJECT_PARAMS_H diff --git a/src/tsp.c b/src/tsp.c deleted file mode 100644 index 258d50b..0000000 --- a/src/tsp.c +++ /dev/null @@ -1,577 +0,0 @@ -#include -#include -#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) -{ - data->node_count = 0; - data->edge_count = 0; - data->edge_weights = 0; - data->edge_list = 0; - - return 0; -} - -void TSP_free_data(struct TSPData *data) -{ - if (!data) return; - if (data->edge_list) free(data->edge_list); - if (data->edge_weights) free(data->edge_weights); -} - -int TSP_init_lp(struct LP *lp, struct TSPData *data) -{ - int node_count = data->node_count; - int edge_count = data->edge_count; - int *edge_list = data->edge_list; - int *edge_weights = data->edge_weights; - - int rval = 0; - - /* Build a row for each degree equation */ - for (int i = 0; i < node_count; i++) - { - rval = LP_new_row(lp, 'E', 2.0); - abort_if(rval, "LP_new_row failed"); - } - - /* Build a column for each edge_index of the graph */ - double lb = 0.0; - double ub = 1.0; - int cmatbeg = 0; - double cmatval[] = {1.0, 1.0}; - for (int j = 0; j < edge_count; j++) - { - double obj = (double) edge_weights[j]; - int cmatind[] = {edge_list[2 * j], edge_list[2 * j + 1]}; - - rval = LP_add_cols(lp, 1, 2, &obj, &cmatbeg, cmatind, cmatval, &lb, - &ub); - - abort_if(rval, "LP_add_cols failed"); - } - - CLEANUP: - return rval; -} - -int TSP_find_violated_subtour_elimination_cut( - struct LP *lp, struct TSPData *data) -{ - int ncount = data->node_count; - int edge_count = data->edge_count; - int *edges = data->edge_list; - - int rval = 0; - int is_infeasible = 0; - - double *x = 0; - int *delta = 0; - int *marks = 0; - int *island_nodes = 0; - int *island_start = 0; - int *island_sizes = 0; - - struct Graph G; - graph_init(&G); - - rval = LP_optimize(lp, &is_infeasible); - abort_if(rval, "LP_optimize failed"); - abort_if(is_infeasible, "LP is infeasible"); - - rval = graph_build(ncount, edge_count, edges, 0, &G); - abort_if(rval, "graph_build failed"); - - x = (double *) malloc(edge_count * sizeof(double)); - delta = (int *) malloc(edge_count * sizeof(int)); - marks = (int *) malloc(ncount * sizeof(int)); - abort_if(!x, "Could not allocate memory for x"); - abort_if(!delta, "Could not allocate memory for delta"); - abort_if(!marks, "Could not allocate memory for marks"); - - island_nodes = (int *) malloc(ncount * sizeof(int)); - island_start = (int *) malloc(ncount * sizeof(int)); - island_sizes = (int *) malloc(ncount * sizeof(int)); - abort_if(!island_nodes, "Could not allocate memory for island_nodes"); - abort_if(!island_start, "Could not allocate memory for island_start"); - abort_if(!island_sizes, "Could not allocate memory for island_sizes"); - - for (int i = 0; i < ncount; i++) - marks[i] = 0; - - rval = LP_get_x(lp, x); - abort_if(rval, "LP_get_x failed"); - -#if LOG_LEVEL >= LOG_LEVEL_VERBOSE - int current_round = 0; - #endif - - int delta_count = 0; - int island_count = 0; - - 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++) - { - get_delta(island_sizes[i], island_nodes + island_start[i], - edge_count, edges, &delta_count, delta, marks); - - rval = TSP_add_subtour_elimination_cut(lp, delta_count, delta); - } - - log_verbose("Reoptimizing (round %d)...\n", ++current_round); - abort_if(rval, "TSP_add_subtour_elimination_cut failed"); - - rval = LP_optimize(lp, &is_infeasible); - abort_if(rval, "LP_optimize failed"); - - if (is_infeasible) goto CLEANUP; - - double objval = 0; - 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(" graph is connected\n"); - - CLEANUP: - graph_free(&G); - if (x) free(x); - if (island_nodes) free(island_nodes); - if (delta) free(delta); - if (marks) free(marks); - return rval; -} - -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; - - double *rmatval; - int *rmatind = delta; - - rmatval = (double *) malloc(delta_length * sizeof(double)); - abort_if(!rmatval, "out of memory for rmatval"); - - for (int i = 0; i < delta_length; i++) - rmatval[i] = 1.0; - - rval = LP_add_rows(lp, 1, delta_length, &rhs, &sense, &rmatbeg, rmatind, - rmatval); - - abort_if(rval, "LP_add_rows failed"); - - CLEANUP: - if (rmatval) free(rmatval); - return rval; -} - -int graph_find_connected_components( - struct Graph *G, - double *x, - int *island_count, - int *island_sizes, - int *island_start, - int *island_nodes) -{ - for (int i = 0; i < G->node_count; i++) - { - G->nodes[i].mark = 0; - island_nodes[i] = -1; - } - - int k = 0, current_island = 0; - - for (int i = 0; i < G->node_count; i++) - { - if (G->nodes[i].mark != 0) continue; - - island_sizes[current_island] = 0; - - graph_dfs(i, G, x, island_sizes + current_island, island_nodes + k); - - island_start[current_island] = k; - k += island_sizes[current_island]; - - current_island++; - } - - (*island_count) = current_island; - - return 0; -} - -int TSP_find_closest_neighbor_tour( - int start, - int node_count, - int edge_count, - int *edges, - int *elen, - int *path_length) -{ - int rval; - int current_node = start; - - struct Graph G; - graph_init(&G); - - rval = graph_build(node_count, edge_count, edges, 0, &G); - abort_if(rval, "graph_build failed"); - - for (int j = 0; j < node_count; j++) - G.nodes[j].mark = 0; - - for (int j = 0; j < node_count; j++) - { - if (j == node_count - 1) - G.nodes[start].mark = 0; - - struct Node *pn = &G.nodes[current_node]; - pn->mark = 1; - - int closest_neighbor = -1; - int closest_edge_length = 10000000; - - for (int i = 0; i < pn->degree; i++) - { - int edge = pn->adj[i].edge_index; - int neighbor = pn->adj[i].neighbor_index; - - if (G.nodes[neighbor].mark == 1) continue; - if (elen[edge] > closest_edge_length) continue; - - closest_neighbor = neighbor; - closest_edge_length = elen[edge]; - } - - *path_length += closest_edge_length; - current_node = closest_neighbor; - } - - CLEANUP: - graph_free(&G); - return rval; -} - -int TSP_read_problem(char *filename, struct TSPData *data) -{ - int *p_node_count = &data->node_count; - int *p_edge_count = &data->edge_count; - int **p_edge_list = &data->edge_list; - int **p_edge_weights = &data->edge_weights; - - struct _IO_FILE *f = (struct _IO_FILE *) NULL; - int i, j, end1, end2, w, rval = 0, node_count, edge_count; - int *edge_list = 0, *edge_weights = 0; - double *x = (double *) NULL, *y = (double *) NULL; - - if (filename) - { - if ((f = fopen(filename, "r")) == NULL) - { - fprintf(stderr, "Unable to open %s for input\n", filename); - rval = 1; - goto CLEANUP; - } - } - - if (filename && GEOMETRIC_DATA == 0) - { - if (fscanf(f, "%d %d", &node_count, &edge_count) != 2) - { - fprintf(stderr, "Input file %s has invalid format\n", filename); - rval = 1; - goto CLEANUP; - } - - printf("Nodes: %d Edges: %d\n", node_count, edge_count); - fflush(stdout); - - edge_list = (int *) malloc(2 * edge_count * sizeof(int)); - if (!edge_list) - { - fprintf(stderr, "out of memory for edge_list\n"); - rval = 1; - goto CLEANUP; - } - - edge_weights = (int *) malloc(edge_count * sizeof(int)); - if (!edge_weights) - { - fprintf(stderr, "out of memory for edge_weights\n"); - rval = 1; - goto CLEANUP; - } - - for (i = 0; i < edge_count; i++) - { - if (fscanf(f, "%d %d %d", &end1, &end2, &w) != 3) - { - fprintf(stderr, "%s has invalid input format\n", filename); - rval = 1; - goto CLEANUP; - } - edge_list[2 * i] = end1; - edge_list[2 * i + 1] = end2; - edge_weights[i] = w; - } - } - else - { - if (filename) - { - if (fscanf(f, "%d", &node_count) != 1) - { - fprintf(stderr, "Input file %s has invalid format\n", filename); - rval = 1; - goto CLEANUP; - } - } - else - { - node_count = NODE_COUNT_RAND; - } - - x = (double *) malloc(node_count * sizeof(double)); - y = (double *) malloc(node_count * sizeof(double)); - if (!x || !y) - { - fprintf(stdout, "out of memory for x or y\n"); - rval = 1; - goto CLEANUP; - } - - if (filename) - { - for (i = 0; i < node_count; i++) - { - if (fscanf(f, "%lf %lf", &x[i], &y[i]) != 2) - { - fprintf(stderr, "%s has invalid input format\n", filename); - rval = 1; - goto CLEANUP; - } - } - } - else - { - rval = generate_random_points_2d(node_count, GRID_SIZE_RAND, x, y); - if (rval) - { - fprintf(stderr, "generate_random_points_2d failed\n"); - goto CLEANUP; - } - - printf("%d\n", node_count); - for (i = 0; i < node_count; i++) - { - printf("%.0f %.0f\n", x[i], y[i]); - } - printf("\n"); - } - - edge_count = (node_count * (node_count - 1)) / 2; - log_verbose("Complete graph: %d nodes, %d edges\n", node_count, - edge_count); - - edge_list = (int *) malloc(2 * edge_count * sizeof(int)); - if (!edge_list) - { - fprintf(stderr, "out of memory for edge_list\n"); - rval = 1; - goto CLEANUP; - } - - edge_weights = (int *) malloc(edge_count * sizeof(int)); - if (!edge_weights) - { - fprintf(stderr, "out of memory for edge_weights\n"); - rval = 1; - goto CLEANUP; - } - - edge_count = 0; - for (i = 0; i < node_count; i++) - { - for (j = i + 1; j < node_count; j++) - { - edge_list[2 * edge_count] = i; - edge_list[2 * edge_count + 1] = j; - edge_weights[edge_count] = get_euclidean_distance(x, y, i, j); - edge_count++; - } - } - } - - *p_node_count = node_count; - *p_edge_count = edge_count; - *p_edge_list = edge_list; - *p_edge_weights = edge_weights; - - CLEANUP: - if (f) fclose(f); - if (x) free(x); - if (y) free(y); - return rval; -} - -int TSP_add_cutting_planes(struct LP *lp, struct TSPData *data) -{ - int rval = 0; - - rval = TSP_find_violated_subtour_elimination_cut(lp, data); - abort_if (rval, "TSP_find_violated_subtour_elimination_cut failed\n"); - - CLEANUP: - return rval; -} - -double TSP_find_initial_solution(struct TSPData *data) -{ - double best_val = 1e99; - - log_verbose("Finding closest neighbor_index tour\n"); - for (int i = 0; i < data->node_count; i++) - { - int path_length = 0; - - TSP_find_closest_neighbor_tour(i, data->node_count, data->edge_count, - data->edge_list, data->edge_weights, &path_length); - - if (best_val > path_length) best_val = path_length; - } - - 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 deleted file mode 100644 index 21aaef0..0000000 --- a/src/tsp.h +++ /dev/null @@ -1,50 +0,0 @@ -#ifndef _PROJECT_TSP_H_ -#define _PROJECT_TSP_H_ - -#include "lp.h" -#include "graph.h" - -struct TSPData -{ - int node_count; - int edge_count; - int *edge_list; - int *edge_weights; -}; - -int TSP_init_data(struct TSPData *data); - -void TSP_free_data(struct TSPData *data); - -int TSP_find_violated_subtour_elimination_cut - (struct LP *lp, struct TSPData *data); - -int graph_find_connected_components( - struct Graph *G, - double *x, - int *island_count, - int *island_sizes, - int *island_start, - int *island_nodes); - -int TSP_find_closest_neighbor_tour( - int start, - int node_count, - int edge_count, - int *edges, - int *elen, - int *path_length); - -int TSP_add_subtour_elimination_cut(struct LP *lp, int delta_length, int *delta); - -int TSP_read_problem(char *filename, struct TSPData *data); - -int TSP_add_cutting_planes(struct LP *lp, struct TSPData *data); - -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.c b/src/util.c index 761ebcf..62869c3 100644 --- a/src/util.c +++ b/src/util.c @@ -5,7 +5,7 @@ #include "util.h" #include "main.h" -double get_current_time() +double get_user_time() { struct rusage ru; @@ -23,9 +23,9 @@ double get_real_time() void time_printf(const char *fmt, ...) { if (INITIAL_TIME == 0) - INITIAL_TIME = get_current_time(); + INITIAL_TIME = get_user_time(); - printf("[%10.2lf] ", get_current_time() - INITIAL_TIME); + printf("[%10.2lf] ", get_user_time() - INITIAL_TIME); va_list args; va_start(args, fmt); diff --git a/src/util.h b/src/util.h index b38b173..ba79aa8 100644 --- a/src/util.h +++ b/src/util.h @@ -48,18 +48,11 @@ fprintf(stderr, "%28s:%d " msg "\n", __FILE__, __LINE__, __VA_ARGS__); \ 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) - #define UNUSED(x) (void)(x) void time_printf(const char *fmt, ...); -double get_current_time(void); +double get_user_time(void); double get_real_time();