diff --git a/scripts/draw-solution.sage b/scripts/draw-solution.sage index aa2d2eb..78edabf 100755 --- a/scripts/draw-solution.sage +++ b/scripts/draw-solution.sage @@ -100,9 +100,9 @@ else: c = white.blend(red, 0.1 + 0.9 * edges[k][2]) plot = plot + line([all_points[edges[k][0]], all_points[edges[k][1]]], color=c) - print ('Drawing labels...') - for i in range(node_count): - plot = plot + text(str(i), all_points[i] + text_offset, color='gray') + #print ('Drawing labels...') + #for i in range(node_count): + # plot = plot + text(str(i), all_points[i] + text_offset, color='gray') print ('Drawing clusters...') for i in range(cluster_count): diff --git a/src/branch-and-cut.c b/src/branch-and-cut.c index cc8d91f..aabc912 100644 --- a/src/branch-and-cut.c +++ b/src/branch-and-cut.c @@ -22,13 +22,13 @@ #include "util.h" #include "gtsp.h" -static int BNC_solve_node(struct BNC *bnc, int depth); +static int solve_node(struct BNC *bnc, int depth); -static int BNC_branch_node(struct BNC *bnc, double *x, int depth); +static int branch_node(struct BNC *bnc, double *x, int depth); -static int BNC_is_integral(double *x, int num_cols); +static int is_integral(double *x, int num_cols); -static int BNC_find_best_branching_var(double *x, int num_cols); +static int find_best_branching_var(double *x, int num_cols); int BNC_init(struct BNC *bnc) { @@ -86,10 +86,10 @@ int BNC_init_lp(struct BNC *bnc) int BNC_solve(struct BNC *bnc) { - return BNC_solve_node(bnc, 1); + return solve_node(bnc, 1); } -static int BNC_solve_node(struct BNC *bnc, int depth) +static int solve_node(struct BNC *bnc, int depth) { struct LP *lp = bnc->lp; double *best_val = &bnc->best_obj_val; @@ -151,11 +151,17 @@ static int BNC_solve_node(struct BNC *bnc, int depth) goto CLEANUP; } + free(x); + num_cols = LP_get_num_cols(lp); + + x = (double *) malloc(num_cols * sizeof(double)); + abort_if(!x, "could not allocate x"); + rval = LP_get_x(lp, x); abort_if(rval, "LP_get_x failed"); } - if (BNC_is_integral(x, num_cols)) + if (is_integral(x, num_cols)) { log_debug("Solution is integral\n"); @@ -179,8 +185,8 @@ static int BNC_solve_node(struct BNC *bnc, int depth) else { log_debug("Solution is fractional\n"); - rval = BNC_branch_node(bnc, x, depth); - abort_if(rval, "BNC_branch_node failed"); + rval = branch_node(bnc, x, depth); + abort_if(rval, "branch_node failed"); } CLEANUP: @@ -188,14 +194,14 @@ static int BNC_solve_node(struct BNC *bnc, int depth) return rval; } -static int BNC_branch_node(struct BNC *bnc, double *x, int depth) +static int branch_node(struct BNC *bnc, double *x, int depth) { int rval = 0; struct LP *lp = bnc->lp; int num_cols = LP_get_num_cols(lp); - int best_branch_var = BNC_find_best_branching_var(x, num_cols); + int best_branch_var = find_best_branching_var(x, num_cols); log_debug("Branching on variable x%d = %.6lf (depth %d)...\n", best_branch_var, x[best_branch_var], depth); @@ -204,8 +210,8 @@ static int BNC_branch_node(struct BNC *bnc, double *x, int depth) rval = LP_change_bound(lp, best_branch_var, 'L', 1.0); abort_if(rval, "LP_change_bound failed"); - rval = BNC_solve_node(bnc, depth + 1); - abort_if(rval, "BNC_solve_node failed"); + rval = solve_node(bnc, depth + 1); + abort_if(rval, "solve_node failed"); rval = LP_change_bound(lp, best_branch_var, 'L', 0.0); abort_if(rval, "LP_change_bound failed"); @@ -214,8 +220,8 @@ static int BNC_branch_node(struct BNC *bnc, double *x, int depth) rval = LP_change_bound(lp, best_branch_var, 'U', 0.0); abort_if(rval, "LP_change_bound failed"); - rval = BNC_solve_node(bnc, depth + 1); - abort_if(rval, "BNC_solve_node failed"); + rval = solve_node(bnc, depth + 1); + abort_if(rval, "solve_node failed"); rval = LP_change_bound(lp, best_branch_var, 'U', 1.0); abort_if(rval, "LP_change_bound failed"); @@ -226,7 +232,7 @@ static int BNC_branch_node(struct BNC *bnc, double *x, int depth) return rval; } -static int BNC_is_integral(double *x, int num_cols) +static int is_integral(double *x, int num_cols) { for (int i = 0; i < num_cols; i++) if (x[i] > EPSILON && x[i] < 1.0 - EPSILON) @@ -235,7 +241,7 @@ static int BNC_is_integral(double *x, int num_cols) return 1; } -static int BNC_find_best_branching_var(double *x, int num_cols) +static int find_best_branching_var(double *x, int num_cols) { int best_index = 0; double best_index_frac = 1.0; diff --git a/src/branch-and-cut.h b/src/branch-and-cut.h index 3a995c8..9c590d7 100644 --- a/src/branch-and-cut.h +++ b/src/branch-and-cut.h @@ -21,18 +21,44 @@ struct BNC { + struct LP *lp; + /* + * Best integral solution found so far. + */ double *best_x; + + /* + * Value of the best integral solution found so far. + */ double best_obj_val; - int *problem_data; + /* + * Pointer to problem-specific data. This pointer is passed to other + * problem-specific functions, such as problem_init_lp and + * problem_add_cutting_planes. + */ + void *problem_data; - int (*problem_init_lp)(struct LP *, void *); + /* + * This callback is called at the beginning of the branch-and-cut search. + * It should initialize the LP, create the initial rows and columns. + */ + int (*problem_init_lp)(struct LP *lp, void *data); - int (*problem_add_cutting_planes)(struct LP *, void *); + /* + * This callback is called after an LP relaxation is solved and before + * checking whether the solution is integral. It should find and add + * any problem-specific cutting planes to the LP. + */ + int (*problem_add_cutting_planes)(struct LP *lp, void *data); - int (*problem_solution_found)(struct BNC*, void *data, double *x); + /* + * This callback is called after a better integral solution has been found. + * This solution satisfies all the problem-specific cuts. + */ + int (*problem_solution_found)(struct BNC *, void *data, double *x); }; int BNC_init(struct BNC *bnc); @@ -43,6 +69,4 @@ int BNC_init_lp(struct BNC *bnc); void BNC_free(struct BNC *bnc); -extern int BNC_NODE_COUNT; - #endif //_PROJECT_BRANCH_AND_CUT_H_ diff --git a/src/flow.c b/src/flow.c index 5f393dd..8e01a48 100644 --- a/src/flow.c +++ b/src/flow.c @@ -22,53 +22,17 @@ int FLOW_MAX_FLOW_COUNT = 0; -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 = 0; - - 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]); +static int mark_reachable_nodes( + const struct Graph *graph, double *residual_caps, struct Node *from); - CLEANUP: - if(parents) free(parents); - if (stack) free(stack); - return rval; -} +static int find_augmenting_path( + const struct Graph *graph, + const double *residual_caps, + struct Node *from, + struct Node *to, + int *path_length, + struct Edge **path_edges, + double *path_capacity); int flow_find_max_flow( const struct Graph *digraph, @@ -85,27 +49,18 @@ int flow_find_max_flow( for (int i = 0; i < digraph->node_count; i++) digraph->nodes[i].mark = 0; - log_verbose("Input graph:\n"); - - #if LOG_LEVEL >= LOG_LEVEL_VERBOSE - graph_dump(digraph); - #endif - 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; + log_verbose("%d %d %.4lf\n", digraph->edges[i].from->index, + digraph->edges[i].to->index, capacities[i]); + int path_length = 0; + double path_capacity = 0; double *residual_caps = 0; + struct Edge **path_edges = 0; residual_caps = (double *) malloc(digraph->edge_count * sizeof(double)); abort_if(!residual_caps, "could not allocate residual_caps"); @@ -128,8 +83,8 @@ int flow_find_max_flow( while (1) { - flow_find_augmenting_path(digraph, residual_caps, from, to, - &path_length, path_edges, &path_capacity); + find_augmenting_path(digraph, residual_caps, from, to, &path_length, + path_edges, &path_capacity); if (path_length == 0) break; @@ -142,7 +97,8 @@ int flow_find_max_flow( { struct Edge *e = &digraph->edges[path_edges[i]->index]; - log_verbose(" %d %d (%d)\n", e->from->index, e->to->index, e->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; @@ -151,24 +107,24 @@ int flow_find_max_flow( flow[e->reverse->index] -= path_capacity; } +#if LOG_LEVEL >= LOG_LEVEL_VERBOSE log_verbose("New residual capacities:\n"); for (int i = 0; i < digraph->edge_count; i++) { - #if LOG_LEVEL >= LOG_LEVEL_VERBOSE struct Edge *e = &digraph->edges[i]; - #endif 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]); + log_verbose("%d %d %.4lf (%d)\n", e->from->index, e->to->index, + e->index, residual_caps[e->index]); } +#endif } log_verbose("No more paths found.\n"); - rval = flow_mark_reachable_nodes(digraph, residual_caps, from); - abort_if(rval, "flow_mark_reachable_nodes failed"); + rval = mark_reachable_nodes(digraph, residual_caps, from); + abort_if(rval, "mark_reachable_nodes failed"); CLEANUP: if (path_edges) free(path_edges); @@ -176,7 +132,56 @@ int flow_find_max_flow( return rval; } -int flow_find_augmenting_path( +int static 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 = 0; + + 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 static find_augmenting_path( const struct Graph *graph, const double *residual_caps, struct Node *from, diff --git a/src/flow.h b/src/flow.h index 99c7d26..427de39 100644 --- a/src/flow.h +++ b/src/flow.h @@ -19,15 +19,6 @@ #include "graph.h" -int flow_find_augmenting_path( - const struct Graph *graph, - const double *residual_caps, - struct Node *from, - struct Node *to, - int *path_length, - struct Edge **path_edges, - double *path_capacity); - int flow_find_max_flow( const struct Graph *digraph, const double *capacities, @@ -36,7 +27,4 @@ int flow_find_max_flow( double *flow, double *value); -int flow_mark_reachable_nodes( - const struct Graph *graph, double *residual_caps, struct Node *from); - #endif //_PROJECT_FLOW_H_ diff --git a/src/geometry.c b/src/geometry.c index cc58205..f252dee 100644 --- a/src/geometry.c +++ b/src/geometry.c @@ -21,8 +21,6 @@ #include "geometry.h" #include "util.h" -/* function for creating a random set of points in unit square */ - int generate_random_points_2d( int node_count, int grid_size, @@ -110,8 +108,8 @@ int generate_random_clusters_2d( for (int j = 0; j < cluster_count; j++) { - int distance = - get_euclidean_distance(x_coordinates, y_coordinates, i, j); + int distance = get_euclidean_distance(x_coordinates, y_coordinates, + i, j); if (distance < closest_distance) { @@ -128,18 +126,21 @@ int generate_random_clusters_2d( } 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) { - int i,j; - for (i = 0; i < node_count; i++){ - for (j = 0; j < node_count; j++){ - dist_matrix[i][j] = - get_euclidean_distance(x_coordinates, y_coordinates, i, j); - } - } - return 0; + int i, j; + for (i = 0; i < node_count; i++) + { + for (j = 0; j < node_count; j++) + { + dist_matrix[i][j] = get_euclidean_distance(x_coordinates, + y_coordinates, i, j); + } + } + return 0; } int get_euclidean_distance( diff --git a/src/graph.c b/src/graph.c index 8261f0c..50fe985 100644 --- a/src/graph.c +++ b/src/graph.c @@ -78,6 +78,7 @@ int graph_build( graph->edges[i].index = i; graph->edges[i].from = &graph->nodes[a]; graph->edges[i].to = &graph->nodes[b]; + graph->edges[i].column = -1; } p = graph->adj; @@ -136,29 +137,3 @@ int get_cut_edges_from_marks( return 0; } - -int graph_dump(const struct Graph *graph) -{ - (void) graph; -#if LOG_LEVEL >= LOG_LEVEL_DEBUG - 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 (e->reverse) printf("reverse: %d ", e->reverse->index); - printf("\n"); - - } - #endif - return 0; -} diff --git a/src/graph.h b/src/graph.h index 407c514..bfe5280 100644 --- a/src/graph.h +++ b/src/graph.h @@ -37,12 +37,30 @@ struct Node struct Edge { + /* + * Index of the edge. Each edge is numbered from 0 to graph->edge_count. + */ int index; + + /* + * Weight of the edge. + */ int weight; + /* + * If this edge corresponds to a column of the LP, this field contains + * the index of that column. Otherwise, this field contains a negative + * value. + */ + int column; + struct Node *from; struct Node *to; + /* + * Pointer to an edge that points in the opposite direction of this one. + * Used by flow algorithms. + */ struct Edge *reverse; }; @@ -71,9 +89,11 @@ int graph_build( int is_directed, struct Graph *graph); +/* + * Returns the list of edges e=uv such that either u or v (but not both) are + * marked nodes. + */ int get_cut_edges_from_marks( struct Graph *graph, int *cut_edges_count, struct Edge **cut_edges); -int graph_dump(const struct Graph *graph); - #endif diff --git a/src/gtsp-cols.c b/src/gtsp-cols.c new file mode 100644 index 0000000..ddeaf0a --- /dev/null +++ b/src/gtsp-cols.c @@ -0,0 +1,111 @@ +/* Copyright (c) 2015 Alinson Xavier + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include "gtsp-cols.h" +#include "util.h" + +int GTSP_find_columns(struct LP *lp, struct GTSP *data) +{ + int rval = 0; + + double *y = 0; + double initial_time = get_user_time(); + int num_rows = LP_get_num_rows(lp); + + y = (double *) malloc(num_rows * sizeof(double)); + abort_if(!y, "could not allocate y"); + + rval = LP_get_y(lp, y); + abort_if(rval, "LP_get_y failed"); + + int violated_count = 0; + int edge_count = data->graph->edge_count; + + log_debug("Finding new columns...\n"); + + for (int i = 0; i < edge_count; i++) + { + struct Edge *e = &data->graph->edges[i]; + if (e->column >= 0) continue; + + double sum = 0; + + for (int j = 0; j < lp->cut_pool_size; j++) + { + struct Row *cut = lp->cut_pool[j]; + if (!cut->edges[e->index]) continue; + + sum += y[cut->cplex_row_index]; + } + + if (sum - EPSILON < e->weight) continue; + + rval = GTSP_add_column(lp, e); + abort_if(rval, "GTSP_add_column failed"); + + violated_count++; + } + + log_debug(" %d of %d edges are violated\n", violated_count, edge_count); + + COLUMNS_TIME += get_user_time() - initial_time; + + CLEANUP: + if (y) free(y); + return rval; +} + +int GTSP_add_column(struct LP *lp, struct Edge *e) +{ + int rval = 0; + + int *cmatind = 0; + double *cmatval = 0; + + int num_rows = LP_get_num_rows(lp); + int num_cols = LP_get_num_cols(lp); + + cmatind = (int *) malloc(num_rows * sizeof(int)); + cmatval = (double *) malloc(num_rows * sizeof(double)); + + abort_if(!cmatind, "could not allocate cmatbeg"); + abort_if(!cmatval, "could not allocate cmatval"); + + int nz = 0; + int cmatbeg = 0; + double lb = 0.0; + double ub = 1.0; + double obj = e->weight; + + for (int j = 0; j < lp->cut_pool_size; j++) + { + struct Row *cut = lp->cut_pool[j]; + if (!cut->edges[e->index]) continue; + + cmatind[nz] = cut->cplex_row_index; + cmatval[nz] = cut->edge_val; + nz++; + } + + e->column = num_cols; + LP_add_cols(lp, 1, nz, &obj, &cmatbeg, cmatind, cmatval, &lb, &ub); + + CLEANUP: + if (cmatind) free(cmatind); + if (cmatval) free(cmatval); + return rval; +} + diff --git a/src/gtsp-cols.h b/src/gtsp-cols.h new file mode 100644 index 0000000..6ba2c8c --- /dev/null +++ b/src/gtsp-cols.h @@ -0,0 +1,27 @@ +/* Copyright (c) 2015 Alinson Xavier + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef PROJECT_GTSP_COLS_H +#define PROJECT_GTSP_COLS_H + +#include "lp.h" +#include "gtsp.h" + +int GTSP_find_columns(struct LP *lp, struct GTSP *data); + +int GTSP_add_column(struct LP *lp, struct Edge *e); + +#endif //PROJECT_GTSP_COLS_H diff --git a/src/gtsp-comb.c b/src/gtsp-comb.c index 80a6455..21a5071 100644 --- a/src/gtsp-comb.c +++ b/src/gtsp-comb.c @@ -45,6 +45,15 @@ static int add_comb_cut( abort_if(!rmatind, "could not allocate rmatind"); abort_if(!rmatval, "could not allocate rmatval"); + cut = (struct Row *) malloc(sizeof(struct Row)); + cut->edges = (char *) malloc(graph->edge_count * sizeof(char)); + + abort_if(!cut, "could not allocate cut"); + abort_if(!cut->edges, "could not allocate cut->edges"); + + for (int i = 0; i < graph->edge_count; i++) + cut->edges[i] = 0; + double rhs = -component_sizes[current_component] - tooth_count + (tooth_count + 1) / 2; @@ -57,7 +66,10 @@ static int add_comb_cut( if (components[clusters[e->from->index]] != current_component) continue; if (components[clusters[e->to->index]] != current_component) continue; - rmatind[nz] = e->index; + cut->edges[e->index] = 1; + if (e->column < 0) continue; + + rmatind[nz] = e->column; rmatval[nz] = -1.0; nz++; @@ -76,11 +88,16 @@ static int add_comb_cut( if (teeth[clusters[from->index]] != teeth[clusters[to->index]]) continue; + cut->edges[e->index] = 1; + if (e->column < 0) continue; + log_verbose(" tooth (%d %d)\n", e->from->index, e->to->index); - rmatind[nz] = e->index; + rmatind[nz] = e->column; rmatval[nz] = -1.0; nz++; + + } #if LOG_LEVEL >= LOG_LEVEL_VERBOSE @@ -123,21 +140,20 @@ static int add_comb_cut( log_verbose("Violation: %.4lf >= %.4lf\n", lhs, rhs); - if (lhs + EPSILON > rhs) - { - free(rmatind); - free(rmatval); - goto CLEANUP; - } - - cut = (struct Row *) malloc(sizeof(struct Row)); - abort_if(!cut, "could not allocate cut"); +// if (lhs + EPSILON > rhs) +// { +// free(rmatind); +// free(rmatval); +// goto CLEANUP; +// } cut->nz = nz; cut->sense = sense; cut->rhs = rhs; cut->rmatval = rmatval; cut->rmatind = rmatind; + cut->edge_val = -1.0; + cut->edge_count = graph->edge_count; rval = LP_add_cut(lp, cut); abort_if(rval, "LP_add_cut failed"); @@ -184,6 +200,7 @@ static int find_components( if (neighbor->mark) continue; double x_e = x[adj->edge->index]; + if (x_e < EPSILON) continue; if (x_e > 1 - EPSILON) continue; @@ -238,7 +255,6 @@ static int find_teeth( struct Node *to = e->to; if (x[e->index] < 1 - EPSILON) continue; - if (to->mark || from->mark) continue; int z = 0; @@ -342,12 +358,13 @@ static int shrink_clusters( for (int i = 0; i < graph->edge_count; i++) { struct Edge *e = &graph->edges[i]; + if(e->column < 0) continue; int from = clusters[e->from->index]; int to = clusters[e->to->index]; int shunk_e_index = edge_map[from * cluster_count + to]; - shrunken_x[shunk_e_index] += x[e->index]; + shrunken_x[shunk_e_index] += x[e->column]; } CLEANUP: diff --git a/src/gtsp-heur.c b/src/gtsp-heur.c new file mode 100644 index 0000000..4b091ac --- /dev/null +++ b/src/gtsp-heur.c @@ -0,0 +1,565 @@ +/* Copyright (c) 2015 Armin Sadeghi + * Copyright (c) 2015 Alinson Xavier + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#include +#include "util.h" +#include "gtsp.h" +#include "gtsp-heur.h" + +static int large_neighborhood_search( + struct Tour *tour, struct GTSP *data, int *tour_cost); + +static int build_tour_from_x(struct GTSP *data, struct Tour *tour, double *x); + +static int build_x_from_tour(struct GTSP *data, struct Tour *tour, double *x); + +static int optimize_vertex_in_cluster(struct Tour *tour, struct GTSP *data); + +static int two_opt(struct Tour *tour, struct GTSP *data); + +//static int tour_length(int *tour, struct GTSP *data); + +static int get_tour_length(struct Tour *tour, struct GTSP *data); + +static int K_opt(struct Tour *tour, struct GTSP *data); + +int GTSP_find_initial_tour(struct GTSP *data, int *tour_cost, double *x) +{ + int rval = 0; + + int cluster_count = data->cluster_count; + + struct Tour *tour = 0; + int *uncovered_sets = 0; + int *cluster_in_tour = 0; + + tour = (struct Tour *) malloc((cluster_count + 1) * sizeof(struct Tour)); + uncovered_sets = (int *) malloc((cluster_count - 1) * sizeof(int)); + cluster_in_tour = (int *) malloc(cluster_count * sizeof(int)); + abort_if(!tour, "could not allocate tour"); + abort_if(!uncovered_sets, "could not allocate uncovered_sets"); + abort_if(!cluster_in_tour, "could not allocate cluster_in_tour"); + + int cluster_num = 0; + for (int i = 0; i < cluster_count; i++) + { + cluster_in_tour[i] = 0; + if (data->node_to_cluster[0] != i) + { + uncovered_sets[cluster_num] = i; + 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; + } + } + } + + assert(min_dist_vertex >= 0); + assert(min_dist_vertex < data->graph->node_count); + + 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"); + + log_info("Initial upper-bound: %d \n", *tour_cost); + + rval = build_x_from_tour(data, tour, x); + abort_if(rval, "build_x_from_tour failed"); + + CLEANUP: + if (tour) free(tour); + if (cluster_in_tour) free(cluster_in_tour); + if (uncovered_sets) free(uncovered_sets); + return rval; +} + +int GTSP_improve_solution(struct BNC *bnc, struct GTSP *data, double *x) +{ + int rval = 0; + + int tour_cost; + double *best_val = &bnc->best_obj_val; + double **best_x = &bnc->best_x; + + struct Tour *tour; + tour = (struct Tour *) malloc( + (data->cluster_count + 1) * sizeof(struct Tour)); + + rval = build_tour_from_x(data, tour, x); + abort_if(rval, "build_tour_from_x failed"); + + for (int i = 0; i < data->cluster_count; i++) + { + log_debug(" %d\n", tour[i]); + } + + rval = large_neighborhood_search(tour, data, &tour_cost); + abort_if(rval, "large_neighborhood_search failed"); + + if (tour_cost + EPSILON < *best_val) + { + log_info("Local search improved the integral solution\n"); + log_info(" before = %f\n", *best_val); + log_info(" after = %f\n", tour_cost); + + build_x_from_tour(data, tour, x); + + *best_val = tour_cost; + *best_x = x; + } + + CLEANUP: + return rval; +} + +int static optimize_vertex_in_cluster(struct Tour *tour, struct GTSP *data) +{ + int current_cluster; + int insertion_cost; + + int **dist_matrix = data->dist_matrix; + int cluster_count = data->cluster_count; + struct Cluster *vertex_set = data->clusters; + + for (int i = 0; i < cluster_count; i++) + { + int vertex = tour[i].vertex; + int prev_vertex = tour[tour[i].prev].vertex; + int next_vertex = tour[tour[i].next].vertex; + + current_cluster = data->node_to_cluster[vertex]; + + insertion_cost = dist_matrix[prev_vertex][vertex] + + dist_matrix[vertex][next_vertex]; + + for (int j = 0; j < vertex_set[current_cluster].size; j++) + { + int vertex_in_cluster = vertex_set[current_cluster].nodes[j]; + int cost = dist_matrix[vertex_in_cluster][prev_vertex] + + dist_matrix[vertex_in_cluster][next_vertex]; + if (insertion_cost > cost) + { + insertion_cost = cost; + tour[i].vertex = vertex_in_cluster; + } + } + } + + return 0; +} + +int static two_opt(struct Tour *tour, struct GTSP *data) +{ + int **dist_matrix = data->dist_matrix; + + for (int i = 0; i < data->cluster_count; i++) + { + int v1 = tour[i].vertex; + int v2 = tour[tour[i].prev].vertex; + int v3 = tour[tour[i].next].vertex; + int v4 = tour[tour[tour[i].next].next].vertex; + + int current_cost = dist_matrix[v2][v1] + dist_matrix[v3][v4]; + int temp_cost = dist_matrix[v2][v3] + dist_matrix[v1][v4]; + + if (current_cost > temp_cost) + { + int temp_next = tour[i].next; + int temp_prev = tour[i].prev; + + tour[i].next = tour[temp_next].next; + tour[i].prev = temp_next; + + tour[tour[temp_next].next].prev = i; + + tour[temp_next].next = i; + tour[temp_next].prev = temp_prev; + + tour[temp_prev].next = temp_next; + + } + } + + return 0; +} + +static int large_neighborhood_search( + struct Tour *tour, struct GTSP *data, int *tour_cost) +{ + int rval = 0; + + int cluster_count = data->cluster_count; + int *node_to_cluster = 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; + + tour[prev_vertex].next = next_vertex; + tour[next_vertex].prev = prev_vertex; + + int cluster_to_insert = node_to_cluster[tour[delete_vertex].vertex]; + + int best_pose = -1; + int best_vertex = -1; + int min_cost = INT_MAX; + + for (int i = 0; i < vertex_set[cluster_to_insert].size; i++) + { + int vertex_to_insert = vertex_set[cluster_to_insert].nodes[i]; + + int next_edge = tour[0].next; + for (int j = 1; j < cluster_count; j++) + { + int vertex1 = tour[next_edge].vertex; + int vertex2 = tour[tour[next_edge].next].vertex; + + int insert_cost = dist_matrix[vertex1][vertex_to_insert] + + dist_matrix[vertex_to_insert][vertex2] - + dist_matrix[vertex1][vertex2]; + + if (insert_cost < min_cost) + { + min_cost = insert_cost; + best_pose = next_edge; + best_vertex = vertex_to_insert; + } + + next_edge = tour[next_edge].next; + } + } + + assert(best_pose >= 0); + assert(best_vertex >= 0); + + next_vertex = tour[best_pose].next; + tour[delete_vertex].prev = best_pose; + tour[delete_vertex].vertex = best_vertex; + 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); + abort_if(rval, "two_opt failed"); + + rval = two_opt(tour, data); + abort_if(rval, "two_opt failed"); + + *tour_cost = get_tour_length(tour, data); + + CLEANUP: + //if (vertex_seq) free(vertex_seq); + return rval; +} + +int static get_tour_length(struct Tour *tour, struct GTSP *data) +{ + int tour_cost = 0; + for (int i = 0; i < data->cluster_count; i++) + { + int vertex1 = tour[i].vertex; + int vertex2 = tour[tour[i].next].vertex; + tour_cost += data->dist_matrix[vertex1][vertex2]; + } + return tour_cost; +} + +int static 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 static build_tour_from_x(struct GTSP *data, struct Tour *tour, double *x) +{ + int rval = 0; + int *cluster_mark = 0; + + struct Node **stack = 0; + int stack_top = 0; + + 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"); + + stack = (struct Node **) malloc(graph->node_count * sizeof(struct Node *)); + abort_if(!stack, "could not allocate stack"); + + for (int i = 0; i < node_count; i++) + graph->nodes[i].mark = 0; + + for (int i = 0; i < data->cluster_count; i++) + cluster_mark[i] = 0; + + int initial = -1; + + for (int i = 0; i < edge_count; i++) + { + int col = graph->edges[i].column; + if (col < 0) continue; + + if (x[col] > 1.0 - EPSILON) + { + initial = i; + break; + } + } + + initial = graph->edges[initial].from->index; + + abort_if(initial < 0, "no initial node"); + + stack[stack_top++] = &graph->nodes[initial]; + graph->nodes[initial].mark = 1; + + tour[0].vertex = graph->nodes[initial].index; + int next_vertex = 1; + + while (stack_top > 0) + { + struct Node *n = stack[--stack_top]; + cluster_mark[data->node_to_cluster[n->index]]++; + + for (int i = 0; i < n->degree; i++) + { + struct Adjacency *adj = &n->adj[i]; + struct Node *neighbor = adj->neighbor; + + if (neighbor->mark) continue; + + if (adj->edge->column < 0) continue; + if (x[adj->edge->column] < EPSILON) continue; + + stack[stack_top++] = neighbor; + tour[next_vertex].vertex = neighbor->index; + 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; + } + } + + CLEANUP: + if (stack) free(stack); + if (cluster_mark) free(cluster_mark); + return rval; +} + +static int build_x_from_tour(struct GTSP *data, struct Tour *tour, double *x) +{ + int rval = 0; + int *edge_map = 0; + + int node_count = data->graph->node_count; + int edge_count = data->graph->edge_count; + + edge_map = (int *) malloc(node_count * node_count * sizeof(int)); + abort_if(!edge_map, "could not allocate edge_map"); + + rval = GTSP_build_edge_map(data, edge_map); + abort_if(rval, "GTSP_build_edge_map failed"); + + for (int i = 0; i < edge_count; i++) + x[i] = 0.0; + + int k = 0; + 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 to = tour[next_vertex].vertex; + current_vertex = tour[next_vertex].vertex; + next_vertex = tour[next_vertex].next; + struct Edge *e = &data->graph->edges[edge_map[from * node_count + to]]; + + if (e->column < 0) e->column = k++; + x[e->column] = 1.0; + } + + CLEANUP: + if (edge_map) free(edge_map); + return rval; +} + diff --git a/src/gtsp-heur.h b/src/gtsp-heur.h new file mode 100644 index 0000000..c75744e --- /dev/null +++ b/src/gtsp-heur.h @@ -0,0 +1,27 @@ +/* Copyright (c) 2015 Armin Sadeghi + * Copyright (c) 2015 Alinson Xavier + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ + +#ifndef PROJECT_GTSP_HEUR_H +#define PROJECT_GTSP_HEUR_H + +#include "gtsp.h" + +int GTSP_find_initial_tour(struct GTSP *data, int *value, double *x); + +int GTSP_improve_solution(struct BNC *bnc, struct GTSP *data, double *x); + +#endif //PROJECT_GTSP_HEUR_H diff --git a/src/gtsp-subtour.c b/src/gtsp-subtour.c index b632b9a..48b3fa1 100644 --- a/src/gtsp-subtour.c +++ b/src/gtsp-subtour.c @@ -19,7 +19,8 @@ #include "util.h" #include "flow.h" -static void deactivate_cluster_node(double *capacities, struct Node *cluster_node) +static void deactivate_cluster_node( + double *capacities, struct Node *cluster_node) { for (int i = 0; i < cluster_node->degree; i++) { @@ -63,7 +64,10 @@ static int build_flow_digraph( int kc = 0; for (int i = 0; i < graph->edge_count; i++) { - if (x[i] < EPSILON) continue; + int col = graph->edges[i].column; + + if (col < 0) continue; + if (x[col] < EPSILON) continue; struct Edge *e = &graph->edges[i]; int from = e->from->index; @@ -71,7 +75,7 @@ static int build_flow_digraph( digraph_edges[ke++] = from; digraph_edges[ke++] = to; - capacities[kc++] = x[i]; + capacities[kc++] = x[col]; digraph_edges[ke++] = to; digraph_edges[ke++] = from; @@ -79,7 +83,7 @@ static int build_flow_digraph( digraph_edges[ke++] = to; digraph_edges[ke++] = from; - capacities[kc++] = x[i]; + capacities[kc++] = x[col]; digraph_edges[ke++] = from; digraph_edges[ke++] = to; @@ -122,51 +126,68 @@ static int build_flow_digraph( } static int add_subtour_cut( - struct LP *lp, struct Edge **cut_edges, int cut_edges_count) + struct LP *lp, + struct Edge **cut_edges, + int cut_edges_count, + struct Graph *graph) { int rval = 0; char sense = 'G'; double rhs = 2.0; - int newnz = cut_edges_count; + int nz = cut_edges_count; int *rmatind = 0; double *rmatval = 0; + struct Row *cut = 0; - rmatind = (int *) malloc(newnz * sizeof(int)); + rmatind = (int *) malloc(nz * sizeof(int)); abort_if(!rmatind, "could not allocate rmatind"); - rmatval = (double *) malloc(newnz * sizeof(double)); + rmatval = (double *) malloc(nz * sizeof(double)); abort_if(!rmatval, "could not allocate rmatval"); + nz = 0; for (int i = 0; i < cut_edges_count; i++) { - rmatind[i] = cut_edges[i]->index; - rmatval[i] = 1.0; + if (cut_edges[i]->column < 0) continue; + rmatind[nz] = cut_edges[i]->column; + rmatval[nz] = 1.0; + nz++; } log_verbose("Generated cut:\n"); - for (int i = 0; i < newnz; i++) + for (int i = 0; i < nz; i++) log_verbose(" %8.2f x%d\n", rmatval[i], rmatind[i]); log_verbose(" %c %.2lf\n", sense, rhs); if (OPTIMAL_X) { double sum = 0; - for (int i = 0; i < newnz; i++) + for (int i = 0; i < nz; i++) sum += rmatval[i] * OPTIMAL_X[rmatind[i]]; abort_if(sum <= rhs - EPSILON, "cannot add invalid cut"); } - struct Row *cut = 0; cut = (struct Row *) malloc(sizeof(struct Row)); abort_if(!cut, "could not allocate cut"); - cut->nz = newnz; + cut->edges = (char *) malloc(graph->edge_count * sizeof(char)); + abort_if(!cut->edges, "could not allocate cut->edges"); + + for (int i = 0; i < graph->edge_count; i++) + cut->edges[i] = 0; + + for (int i = 0; i < cut_edges_count; i++) + cut->edges[cut_edges[i]->index] = 1; + + cut->nz = nz; cut->sense = sense; cut->rhs = rhs; cut->rmatval = rmatval; cut->rmatind = rmatind; + cut->edge_count = graph->edge_count; + cut->edge_val = 1.0; rval = LP_add_cut(lp, cut); abort_if(rval, "LP_add_cut failed"); @@ -203,7 +224,7 @@ int GTSP_find_exact_subtour_cuts(struct LP *lp, struct GTSP *data) abort_if(rval, "LP_get_x failed"); #if LOG_LEVEL >= LOG_LEVEL_DEBUG - if(strlen(FRAC_SOLUTION_FILENAME) > 0) + if (strlen(FRAC_SOLUTION_FILENAME) > 0) { rval = GTSP_write_solution(data, FRAC_SOLUTION_FILENAME, x); abort_if(rval, "GTSP_write_solution failed"); @@ -268,7 +289,7 @@ int GTSP_find_exact_subtour_cuts(struct LP *lp, struct GTSP *data) log_verbose(" %d %d (%d)\n", cut_edges[k]->from->index, cut_edges[k]->to->index, cut_edges[k]->index); - rval = add_subtour_cut(lp, cut_edges, cut_edges_count); + rval = add_subtour_cut(lp, cut_edges, cut_edges_count, graph); abort_if(rval, "add_subtour_cut failed"); cuts_count++; diff --git a/src/gtsp.c b/src/gtsp.c index 1ac0865..bc57d10 100644 --- a/src/gtsp.c +++ b/src/gtsp.c @@ -17,19 +17,12 @@ #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 build_edge_map(struct GTSP *gtsp, int *edge_map); - -double *OPTIMAL_X = 0; +#include "gtsp-cols.h" int GTSP_init_data(struct GTSP *data) { @@ -76,6 +69,7 @@ int GTSP_create_random_problem( int *weights = 0; int *node_to_cluster = 0; struct Cluster *clusters = 0; + int *edge_map = 0; int **dist_matrix = 0; @@ -180,6 +174,7 @@ int GTSP_create_random_problem( data->clusters = clusters; CLEANUP: + if (edge_map) free(edge_map); if (weights) free(weights); if (edges) free(edges); if (rval) @@ -196,31 +191,22 @@ int GTSP_init_lp(struct LP *lp, struct GTSP *data) int rval = 0; int edge_count = data->graph->edge_count; - int cluster_count = data->cluster_count; - int *clusters = data->node_to_cluster; struct Edge *edges = data->graph->edges; - for (int i = 0; i < cluster_count; i++) - { - rval = LP_new_row(lp, 'E', 2.0); - abort_if(rval, "LP_new_row failed"); - } - double lb = 0.0; double ub = 1.0; - int cmatbeg = 0; + int empty = 0; + double emptyd = 0.0; + int k = 0; for (int i = 0; i < edge_count; i++) { - struct Node *from = edges[i].from; - struct Node *to = edges[i].to; + if (edges[i].column < 0) continue; + edges[i].column = k++; double obj = (double) edges[i].weight; - double cmatval[] = {1.0, 1.0}; - int cmatind[] = {clusters[from->index], clusters[to->index]}; - rval = LP_add_cols(lp, 1, 2, &obj, &cmatbeg, cmatind, cmatval, &lb, - &ub); + rval = LP_add_cols(lp, 1, 0, &obj, &empty, &empty, &emptyd, &lb, &ub); abort_if(rval, "LP_add_cols failed"); } @@ -251,6 +237,7 @@ int GTSP_add_cutting_planes(struct LP *lp, struct GTSP *data) int original_cut_pool_size; int added_cuts_count; + // Generalized subtour cuts original_cut_pool_size = lp->cut_pool_size; log_debug("Finding subtour cuts, round %d...\n", current_round); @@ -261,6 +248,7 @@ int GTSP_add_cutting_planes(struct LP *lp, struct GTSP *data) if (added_cuts_count > 0) continue; + // Generalized comb cuts original_cut_pool_size = lp->cut_pool_size; log_debug("Finding comb cuts, round %d...\n", current_round); @@ -271,6 +259,17 @@ int GTSP_add_cutting_planes(struct LP *lp, struct GTSP *data) if (added_cuts_count > 0) continue; + // Column generation + int original_cols_count = LP_get_num_cols(lp); + + rval = GTSP_find_columns(lp, data); + abort_if(rval, "GTSP_find_columns failed"); + + int added_cols_count = LP_get_num_cols(lp) - original_cols_count; + + if (added_cols_count > 0) + continue; + break; } @@ -309,9 +308,13 @@ int GTSP_print_solution(struct GTSP *data, double *x) int edge_count = data->graph->edge_count; for (int i = 0; i < edge_count; i++) - if (x[i] > EPSILON) + { + int col = edges[i].column; + if (col < 0) continue; + if (x[col] > EPSILON) log_info(" %-3d %-3d %8.4lf\n", edges[i].from->index, - edges[i].to->index, x[i]); + edges[i].to->index, x[col]); + } return 0; } @@ -330,15 +333,27 @@ 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] > EPSILON) - positive_edge_count++; + { + struct Edge *e = &edges[i]; + + if (e->column < 0) continue; + if (x[e->column] < EPSILON) continue; + + positive_edge_count++; + } fprintf(file, "%d\n", positive_edge_count); for (int i = 0; i < edge_count; i++) - if (x[i] > EPSILON) - fprintf(file, "%d %d %.4lf\n", edges[i].from->index, - edges[i].to->index, x[i]); + { + struct Edge *e = &edges[i]; + + if (e->column < 0) continue; + if (x[e->column] < EPSILON) continue; + + fprintf(file, "%d %d %.4lf\n", edges[i].from->index, edges[i].to->index, + x[e->column]); + } CLEANUP: if (file) fclose(file); @@ -376,8 +391,8 @@ int GTSP_read_solution(struct GTSP *gtsp, char *filename, double **p_x) edge_map = (int *) malloc(node_count * node_count * sizeof(int)); abort_if(!edge_map, "could not allocate edge_map"); - rval = build_edge_map(gtsp, edge_map); - abort_if(rval, "build_edge_map failed"); + rval = GTSP_build_edge_map(gtsp, edge_map); + abort_if(rval, "GTSP_build_edge_map failed"); for (int i = 0; i < edge_count; i++) { @@ -406,25 +421,6 @@ int GTSP_read_solution(struct GTSP *gtsp, char *filename, double **p_x) return rval; } -int build_edge_map(struct GTSP *gtsp, int *edge_map) -{ - int node_count = gtsp->graph->node_count; - - int k = 0; - for (int i = 0; i < node_count; i++) - { - for (int j = i + 1; j < node_count; j++) - { - if (gtsp->node_to_cluster[i] == gtsp->node_to_cluster[j]) continue; - edge_map[i * node_count + j] = k; - edge_map[j * node_count + i] = k; - k++; - } - } - - return 0; -} - int GTSP_check_solution(struct GTSP *data, double *x) { int rval = 0; @@ -445,11 +441,14 @@ int GTSP_check_solution(struct GTSP *data, double *x) for (int i = 0; i < edge_count; i++) { - abort_iff(x[i] < 1.0 - EPSILON && x[i] > EPSILON, - "solution is not integral: x%d = %.4lf", i, x[i]); + int col = graph->edges[i].column; + if (col < 0) continue; - abort_iff(x[i] > 1.0 + EPSILON || x[i] < 0.0 - EPSILON, - "value out of bounds: x%d = %.4lf", i, x[i]); + abort_iff(x[col] < 1.0 - EPSILON && x[col] > EPSILON, + "solution is not integral: x%d = %.4lf", i, x[col]); + + abort_iff(x[col] > 1.0 + EPSILON || x[col] < 0.0 - EPSILON, + "value out of bounds: x%d = %.4lf", i, x[col]); } for (int i = 0; i < node_count; i++) @@ -464,8 +463,15 @@ int GTSP_check_solution(struct GTSP *data, double *x) abort_if(initial == edge_count, "no initial node"); - stack[stack_top++] = graph->edges[initial].from; - graph->edges[initial].from->mark = 1; + for (int i = 0; i < edge_count; i++) + { + if (graph->edges[i].column == initial) + { + stack[stack_top++] = graph->edges[i].from; + graph->edges[i].from->mark = 1; + break; + } + } while (stack_top > 0) { @@ -478,7 +484,9 @@ int GTSP_check_solution(struct GTSP *data, double *x) struct Node *neighbor = adj->neighbor; if (neighbor->mark) continue; - if (x[adj->edge->index] < EPSILON) continue; + + if (adj->edge->column < 0) continue; + if (x[adj->edge->column] < EPSILON) continue; stack[stack_top++] = neighbor; neighbor->mark = 1; @@ -486,10 +494,7 @@ int GTSP_check_solution(struct GTSP *data, double *x) } for (int i = 0; i < data->cluster_count; i++) - { - abort_iff(cluster_mark[i] > 1, "cluster %d visited multiple times", i); abort_iff(cluster_mark[i] < 1, "cluster %d not visited", i); - } log_info(" solution is valid\n"); @@ -501,14 +506,7 @@ int GTSP_check_solution(struct GTSP *data, double *x) int GTSP_solution_found(struct BNC *bnc, struct GTSP *data, double *x) { - UNUSED(bnc); - int rval = 0; - int tour_cost; - double *best_val = &bnc->best_obj_val; - - struct Tour *tour; - tour = (struct Tour *) malloc(data->cluster_count * sizeof(struct Tour)); if (strlen(SOLUTION_FILENAME) > 0) { @@ -517,20 +515,6 @@ int GTSP_solution_found(struct BNC *bnc, struct GTSP *data, double *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 + 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"); @@ -539,513 +523,21 @@ int GTSP_solution_found(struct BNC *bnc, struct GTSP *data, double *x) return rval; } -int build_x_from_tour(struct GTSP *data, struct Tour *tour, double *x) +int GTSP_build_edge_map(struct GTSP *gtsp, int *edge_map) { - int rval = 0; - int *edge_map = 0; - - int node_count = data->graph->node_count; - int edge_count = data->graph->edge_count; - - edge_map = (int *) malloc(node_count * node_count * sizeof(int)); - abort_if(!edge_map, "could not allocate edge_map"); - - rval = build_edge_map(data, edge_map); - abort_if(rval, "build_edge_map failed"); - - for (int i = 0; i < edge_count; i++) - x[i] = 0.0; - - 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 to = tour[next_vertex].vertex; - current_vertex = tour[next_vertex].vertex; - next_vertex = tour[next_vertex].next; - - x[edge_map[from * node_count + to]] = 1.0; - } - - CLEANUP: - if (edge_map) free(edge_map); - return rval; -} - -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; - int *uncovered_sets = 0; - int *cluster_in_tour = 0; - - tour = (struct Tour *) malloc((cluster_count + 1) * sizeof(struct Tour)); - uncovered_sets = (int *) malloc((cluster_count - 1) * sizeof(int)); - cluster_in_tour = (int *) malloc(cluster_count * sizeof(int)); - abort_if(!tour, "could not allocate tour"); - abort_if(!uncovered_sets, "could not allocate uncovered_sets"); - abort_if(!cluster_in_tour, "could not allocate cluster_in_tour"); - - int cluster_num = 0; - for (int i = 0; i < cluster_count; i++) - { - cluster_in_tour[i] = 0; - if (data->node_to_cluster[0] != i) - { - uncovered_sets[cluster_num] = i; - 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; - } - } - } - - 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; - - new_vertex += 1; - } - - tour[data->cluster_count].vertex = 0; - - rval = large_neighborhood_search(tour, data, tour_cost); - abort_if(rval, "large_neighborhood_search failed"); - - log_info("Initial upper-bound: %d \n", *tour_cost); - - rval = build_x_from_tour(data, tour, x); - abort_if(rval, "build_x_from_tour failed"); - - CLEANUP: - if (tour) free(tour); - if (cluster_in_tour) free(cluster_in_tour); - if (uncovered_sets) free(uncovered_sets); - return rval; -} - -int optimize_vertex_in_cluster(struct Tour *tour, struct GTSP *data) -{ - int current_cluster; - int insertion_cost; - - int **dist_matrix = data->dist_matrix; - int cluster_count = data->cluster_count; - struct Cluster *vertex_set = data->clusters; - - for (int i = 0; i < cluster_count; i++) - { - int vertex = tour[i].vertex; - int prev_vertex = tour[tour[i].prev].vertex; - int next_vertex = tour[tour[i].next].vertex; - - current_cluster = data->node_to_cluster[vertex]; - - insertion_cost = dist_matrix[prev_vertex][vertex] + - dist_matrix[vertex][next_vertex]; - - for (int j = 0; j < vertex_set[current_cluster].size; j++) - { - int vertex_in_cluster = vertex_set[current_cluster].nodes[j]; - int cost = dist_matrix[vertex_in_cluster][prev_vertex] + - dist_matrix[vertex_in_cluster][next_vertex]; - if (insertion_cost > cost) - { - insertion_cost = cost; - tour[i].vertex = vertex_in_cluster; - } - } - } - - return 0; -} - -int two_opt(struct Tour *tour, struct GTSP *data) -{ - int **dist_matrix = data->dist_matrix; - - for (int i = 0; i < data->cluster_count; i++) - { - int v1 = tour[i].vertex; - int v2 = tour[tour[i].prev].vertex; - int v3 = tour[tour[i].next].vertex; - int v4 = tour[tour[tour[i].next].next].vertex; - - int current_cost = dist_matrix[v2][v1] + dist_matrix[v3][v4]; - int temp_cost = dist_matrix[v2][v3] + dist_matrix[v1][v4]; - - if (current_cost > temp_cost) - { - int temp_next = tour[i].next; - int temp_prev = tour[i].prev; - - tour[i].next = tour[temp_next].next; - tour[i].prev = temp_next; - - tour[tour[temp_next].next].prev = i; - - tour[temp_next].next = i; - tour[temp_next].prev = temp_prev; - - tour[temp_prev].next = temp_next; - - } - } - - return 0; -} - -int large_neighborhood_search( - struct Tour *tour, struct GTSP *data, int *tour_cost) -{ - int rval = 0; - - int cluster_count = data->cluster_count; - 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; - - tour[prev_vertex].next = next_vertex; - tour[next_vertex].prev = prev_vertex; - - int cluster_to_insert = clusters[tour[delete_vertex].vertex]; - - int best_pose = -1; - int best_vertex = -1; - int min_cost = INT_MAX; - - for (int i = 0; i < vertex_set[cluster_to_insert].size; i++) - { - int vertex_to_insert = vertex_set[cluster_to_insert].nodes[i]; - - int next_edge = tour[0].next; - for (int j = 1; j < cluster_count; j++) - { - int vertex1 = tour[next_edge].vertex; - int vertex2 = tour[tour[next_edge].next].vertex; - - int insert_cost = dist_matrix[vertex1][vertex_to_insert] + - dist_matrix[vertex_to_insert][vertex2] - - dist_matrix[vertex1][vertex2]; - - if (insert_cost < min_cost) - { - min_cost = insert_cost; - best_pose = next_edge; - best_vertex = vertex_to_insert; - } - - next_edge = tour[next_edge].next; - } - } - - assert(best_pose >= 0); - assert(best_vertex >= 0); - - next_vertex = tour[best_pose].next; - tour[delete_vertex].prev = best_pose; - tour[delete_vertex].vertex = best_vertex; - 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); - 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; -} - -int tour_length(int *tour, struct GTSP *data) -{ - int tour_cost = 0; - for (int i = 0; i < data->cluster_count; i++) - { - if (i == data->cluster_count - 1) - tour_cost += data->dist_matrix[tour[i]][tour[0]]; - else - tour_cost += data->dist_matrix[tour[i]][tour[i + 1]]; - - } - return tour_cost; -} - -int list_length(struct Tour *tour, struct GTSP *data) -{ - int tour_cost = 0; - for (int i = 0; i < data->cluster_count; i++) - { - int vertex1 = tour[i].vertex; - int vertex2 = tour[tour[i].next].vertex; - tour_cost += data->dist_matrix[vertex1][vertex2]; - } - return tour_cost; -} - -void print_tour(int *tour, struct GTSP *data) -{ - for (int i = 0; i < data->cluster_count; i++) - { - printf("%d\t", tour[i]); - } - - printf("\n"); -} - -void print_list(struct Tour *tour, struct GTSP *data) -{ - printf("%d\t", tour[0].vertex); - int vertex_next = tour[0].next; - - for (int i = 1; i < data->cluster_count; i++) - { - printf("%d\t", tour[vertex_next].vertex); - vertex_next = tour[vertex_next].next; - } - - printf("\n"); -} - -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 build_tour_from_x(struct GTSP *data, struct Tour *tour, double *x) -{ - - int rval = 0; - int *cluster_mark = 0; - - struct Node **stack = 0; - int stack_top = 0; - - 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"); - - stack = (struct Node **) malloc(graph->node_count * sizeof(struct Node *)); - abort_if(!stack, "could not allocate stack"); + int node_count = gtsp->graph->node_count; + int k = 0; for (int i = 0; i < node_count; i++) - graph->nodes[i].mark = 0; - - 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 - EPSILON) break; - - initial = graph->edges[initial].from->index; - - abort_if(initial == edge_count, "no initial node"); - - stack[stack_top++] = &graph->nodes[initial]; - graph->nodes[initial].mark = 1; - - tour[0].vertex = graph->nodes[initial].index; - int next_vertex = 1; - - while (stack_top > 0) { - struct Node *n = stack[--stack_top]; - cluster_mark[data->node_to_cluster[n->index]]++; - - for (int i = 0; i < n->degree; i++) - { - struct Adjacency *adj = &n->adj[i]; - struct Node *neighbor = adj->neighbor; - - if (neighbor->mark) continue; - if (x[adj->edge->index] < EPSILON) continue; - - stack[stack_top++] = neighbor; - tour[next_vertex].vertex = neighbor->index; - 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 + for (int j = i + 1; j < node_count; j++) { - tour[i].next = i + 1; + if (gtsp->node_to_cluster[i] == gtsp->node_to_cluster[j]) continue; + edge_map[i * node_count + j] = k; + edge_map[j * node_count + i] = k; + k++; } } - CLEANUP: - if (stack) free(stack); - if (cluster_mark) free(cluster_mark); - return rval; + return 0; } - diff --git a/src/gtsp.h b/src/gtsp.h index f9110e8..1119f61 100644 --- a/src/gtsp.h +++ b/src/gtsp.h @@ -31,26 +31,36 @@ struct Tour struct Cluster { - int size; - int* nodes; + /* + * Number of nodes inside the cluster. + */ + int size; + + /* + * List of nodes inside the cluster. + */ + int *nodes; }; struct GTSP { struct Graph *graph; - int** dist_matrix; + int **dist_matrix; - int cluster_count; + /* + * Mapping between a node and the cluster which contains it. If a node + * has index i and is contained in cluster k, then node_to_cluster[i] = k. + */ int *node_to_cluster; + + int cluster_count; struct Cluster *clusters; }; int GTSP_create_random_problem( int node_count, int cluster_count, int grid_size, struct GTSP *data); -int inital_tour_value(struct GTSP *data, int *value, double *x); - void GTSP_free(struct GTSP *data); int GTSP_init_data(struct GTSP *data); @@ -63,30 +73,15 @@ int GTSP_write_problem(struct GTSP *data, char *filename); int GTSP_write_solution(struct GTSP *data, char *filename, double *x); -int optimize_vertex_in_cluster(struct Tour * tour, struct GTSP *data); - -int two_opt(struct Tour * tour, struct GTSP *data); - -int K_opt(struct Tour* tour, struct GTSP *data); - -int tour_length(int* tour, struct GTSP* data); - -void print_tour(int* tour, struct GTSP* data); - -int list_length(struct Tour * tour, struct GTSP* data); - -void print_list(struct Tour * tour, struct GTSP* data); - int GTSP_print_solution(struct GTSP *data, double *x); - - -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); +int GTSP_build_edge_map(struct GTSP *gtsp, int *edge_map); + + #endif //_PROJECT_GTSP_H_ diff --git a/src/lp.c b/src/lp.c index 1ee162d..aa2d1fd 100644 --- a/src/lp.c +++ b/src/lp.c @@ -21,18 +21,19 @@ #include #include "lp.h" #include "util.h" -#include "main.h" +#include "gtsp-subtour.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) + struct Row *cut = lp->cut_pool[i]; + + if (cut->cplex_row_index < 0) { - free(lp->cut_pool[i]->rmatind); - free(lp->cut_pool[i]->rmatval); - free(lp->cut_pool[i]); + free(cut->edges); + free(cut); delete_count++; } else @@ -87,8 +88,10 @@ static int remove_old_cuts(struct LP *lp) { 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--; + 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++; @@ -132,20 +135,23 @@ static int remove_old_cuts(struct LP *lp) 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); + struct Row *cut = lp->cut_pool[i]; + + nz += cut->nz; + size += cut->nz * sizeof(double); + size += cut->nz * sizeof(int); + size += cut->edge_count * sizeof(char); } - log_debug(" %ld cuts (%ld nz, %ld MiB)\n", lp->cut_pool_size, nz, size/1024/1024); + log_debug(" %ld cuts (%ld nz, %ld MiB)\n", lp->cut_pool_size, nz, + size / 1024 / 1024); - if(size > CUT_POOL_MAX_MEMORY) + if (size > CUT_POOL_MAX_MEMORY) CUT_POOL_MAX_MEMORY = size; CLEANUP: @@ -167,7 +173,6 @@ static int update_cut_ages(struct LP *lp) 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]; @@ -178,8 +183,6 @@ static int update_cut_ages(struct LP *lp) cut->age++; else cut->age = 0; - - if (cut->age > 5) count++; } CLEANUP: @@ -208,9 +211,9 @@ void LP_free_cut_pool(struct LP *lp) { for (int i = 0; i < lp->cut_pool_size; i++) { - free(lp->cut_pool[i]->rmatind); - free(lp->cut_pool[i]->rmatval); - free(lp->cut_pool[i]); + struct Row *cut = lp->cut_pool[i]; + free(cut->edges); + free(cut); } if (lp->cut_pool) free(lp->cut_pool); @@ -322,8 +325,8 @@ int LP_optimize(struct LP *lp, int *infeasible) int numrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp); int numcols = CPXgetnumcols(lp->cplex_env, lp->cplex_lp); - if(numrows > LP_MAX_ROWS) LP_MAX_ROWS = numrows; - if(numrows > LP_MAX_COLS) LP_MAX_COLS = numcols; + if (numrows > LP_MAX_ROWS) LP_MAX_ROWS = numrows; + if (numrows > LP_MAX_COLS) LP_MAX_COLS = numcols; log_debug("Optimizing LP (%d rows %d cols)...\n", numrows, numcols); @@ -337,6 +340,8 @@ int LP_optimize(struct LP *lp, int *infeasible) solstat = CPXgetstat(lp->cplex_env, lp->cplex_lp); if (solstat == CPX_STAT_INFEASIBLE) { + log_debug(" infeasible\n"); + *infeasible = 1; goto CLEANUP; } @@ -357,9 +362,12 @@ int LP_optimize(struct LP *lp, int *infeasible) rval = update_cut_ages(lp); abort_if(rval, "update_cut_ages failed"); - log_debug("Removing old cuts...\n"); - rval = remove_old_cuts(lp); - abort_if(rval, "LP_remove_old_cuts failed"); + if (LP_SOLVE_COUNT % MAX_CUT_AGE == 0) + { + log_debug("Removing old cuts...\n"); + rval = remove_old_cuts(lp); + abort_if(rval, "LP_remove_old_cuts failed"); + } CUT_POOL_TIME += get_user_time() - initial_time; @@ -392,11 +400,30 @@ int LP_get_x(struct LP *lp, double *x) return rval; } +int LP_get_y(struct LP *lp, double *y) +{ + int rval = 0; + + int nrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp); + abort_if(!nrows, "No rows in LP"); + + rval = CPXgetpi(lp->cplex_env, lp->cplex_lp, y, 0, nrows - 1); + abort_iff(rval, "CPXgetpi failed (errno = %d)", rval); + + CLEANUP: + return rval; +} + int LP_get_num_cols(struct LP *lp) { return CPXgetnumcols(lp->cplex_env, lp->cplex_lp); } +int LP_get_num_rows(struct LP *lp) +{ + return CPXgetnumrows(lp->cplex_env, lp->cplex_lp); +} + int LP_write(struct LP *lp, const char *fname) { int rval = 0; @@ -427,12 +454,10 @@ int LP_write(struct LP *lp, const char *fname) int compare_cuts(struct Row *cut1, struct Row *cut2) { - return_if_neq(cut1->nz, cut2->nz); - for (int i = 0; i < cut1->nz; i++) + for (int i = 0; i < cut1->edge_count; i++) { - return_if_neq(cut1->rmatind[i], cut2->rmatind[i]); - return_if_neq_epsilon(cut1->rmatval[i], cut2->rmatval[i]); + return_if_neq(cut1->edges[i], cut2->edges[i]); } return 0; @@ -442,9 +467,9 @@ static int update_hash(struct Row *cut) { unsigned long hash = 0; - for (int i = 0; i < cut->nz; i++) + for (int i = 0; i < cut->edge_count; i++) { - hash += cut->rmatind[i] * 65521; + hash += cut->edges[i] * 65521; hash %= 4294967291; } @@ -472,6 +497,7 @@ int LP_add_cut(struct LP *lp, struct Row *cut) free(cut->rmatval); free(cut->rmatind); + free(cut->edges); free(cut); return 0; } @@ -488,6 +514,11 @@ 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; + free(cut->rmatval); + free(cut->rmatind); + cut->rmatind = 0; + cut->rmatval = 0; + CUT_POOL_TIME += get_user_time() - initial_time; CLEANUP: return rval; diff --git a/src/lp.h b/src/lp.h index b5d522d..e8e2c5d 100644 --- a/src/lp.h +++ b/src/lp.h @@ -40,6 +40,10 @@ struct Row double rhs; double *rmatval; int *rmatind; + + int edge_count; + double edge_val; + char *edges; }; static const int MAX_NAME_LENGTH = 100; @@ -83,8 +87,12 @@ int LP_get_obj_val(struct LP *lp, double *obj); int LP_get_x(struct LP *lp, double *x); +int LP_get_y(struct LP *lp, double *y); + int LP_get_num_cols(struct LP *lp); +int LP_get_num_rows(struct LP *lp); + 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 434a640..55e1aaf 100644 --- a/src/main.c +++ b/src/main.c @@ -21,11 +21,11 @@ #include "main.h" #include "gtsp.h" #include "util.h" - -unsigned int SEED = 0; +#include "gtsp-heur.h" double SUBTOUR_TIME = 0; double COMBS_TIME = 0; +double COLUMNS_TIME = 0; double CUT_POOL_TIME = 0; long CUT_POOL_MAX_MEMORY = 0; @@ -52,19 +52,24 @@ char FRAC_SOLUTION_FILENAME[100] = {0}; char WRITE_GRAPH_FILENAME[100] = {0}; char STATS_FILENAME[100] = {0}; +double *OPTIMAL_X = 0; +unsigned int SEED = 0; + 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'}, {"out", required_argument, 0, 'o'}, - {"stats", required_argument, 0, 't'}, {"lp", required_argument, 0, 'l'}, - {"write-graph", required_argument, 0, 'w'}, - {(char *) 0, (int) 0, (int *) 0, (int) 0}}; +static const struct option options_tab[] = {{"help", no_argument, 0, 'h'}, + {"nodes", required_argument, 0, 'n'}, + {"clusters", required_argument, 0, 'm'}, + {"grid-size", required_argument, 0, 'g'}, + {"optimal", required_argument, 0, 'x'}, + {"seed", required_argument, 0, 's'}, + {"out", required_argument, 0, 'o'}, + {"stats", required_argument, 0, 't'}, + {"lp", required_argument, 0, 'l'}, + {"write-graph", required_argument, 0, 'w'}, + {(char *) 0, (int) 0, (int *) 0, (int) 0}}; static char input_x_filename[1000] = {0}; @@ -103,7 +108,7 @@ static int parse_args(int argc, char **argv) { int c = 0; int option_index = 0; - c = getopt_long(argc, argv, "n:m:g:x:s:h:o:t:l:w:", options_tab, + c = getopt_long(argc, argv, "n:m:g:x:s:h:o:t:l:w:f:", options_tab, &option_index); if (c < 0) break; @@ -146,6 +151,10 @@ static int parse_args(int argc, char **argv) strcpy(WRITE_GRAPH_FILENAME, optarg); break; + case 'f': + strcpy(FRAC_SOLUTION_FILENAME, optarg); + break; + case 'h': print_usage(argv); exit(0); @@ -235,11 +244,10 @@ int main(int argc, char **argv) int init_val = 0; - initial_x = (double *) malloc( - (data.graph->node_count + data.graph->edge_count) * sizeof(double)); + initial_x = (double *) malloc((data.graph->edge_count) * sizeof(double)); abort_if(!initial_x, "could not allocate initial_x"); - rval = inital_tour_value(&data, &init_val, initial_x); + rval = GTSP_find_initial_tour(&data, &init_val, initial_x); abort_if(rval, "initial_tour_value failed"); rval = GTSP_solution_found(&bnc, &data, initial_x); @@ -289,7 +297,7 @@ int main(int argc, char **argv) log_info("Starting branch-and-cut solver...\n"); rval = BNC_solve(&bnc); - abort_if(rval, "BNC_solve_node failed"); + abort_if(rval, "solve_node failed"); abort_if(!bnc.best_x, "problem has no feasible solution"); @@ -317,6 +325,7 @@ int main(int argc, char **argv) fprintf(file, "%-8s ", "time"); fprintf(file, "%-8s ", "subt-t"); fprintf(file, "%-8s ", "combs-t"); + fprintf(file, "%-8s ", "cols-t"); fprintf(file, "%-8s ", "pool-t"); fprintf(file, "%-8s ", "pool-m"); fprintf(file, "%-8s ", "lp-count"); @@ -337,6 +346,7 @@ int main(int argc, char **argv) fprintf(file, "%-8.2lf ", TOTAL_TIME); fprintf(file, "%-8.2lf ", SUBTOUR_TIME); fprintf(file, "%-8.2lf ", COMBS_TIME); + fprintf(file, "%-8.2lf ", COLUMNS_TIME); fprintf(file, "%-8.2lf ", CUT_POOL_TIME); fprintf(file, "%-8ld ", CUT_POOL_MAX_MEMORY / 1024 / 1024); fprintf(file, "%-8d ", LP_SOLVE_COUNT); @@ -353,7 +363,7 @@ int main(int argc, char **argv) fclose(file); } - if(strlen(SOLUTION_FILENAME) == 0) + if (strlen(SOLUTION_FILENAME) == 0) { log_info("Optimal solution:\n"); rval = GTSP_print_solution(&data, bnc.best_x); @@ -363,7 +373,6 @@ int main(int argc, char **argv) log_info("Optimal solution value:\n"); log_info(" %.4lf\n", bnc.best_obj_val); - CLEANUP: if (OPTIMAL_X) free(OPTIMAL_X); GTSP_free(&data); diff --git a/src/main.h b/src/main.h index c556369..ae443b3 100644 --- a/src/main.h +++ b/src/main.h @@ -22,6 +22,7 @@ extern unsigned int SEED; extern double *OPTIMAL_X; extern double SUBTOUR_TIME; extern double COMBS_TIME; +extern double COLUMNS_TIME; extern double LP_SOLVE_TIME; extern int LP_SOLVE_COUNT; @@ -44,4 +45,6 @@ extern char LP_FILENAME[100]; extern char SOLUTION_FILENAME[100]; extern char FRAC_SOLUTION_FILENAME[100]; +extern int BNC_NODE_COUNT; + #endif \ No newline at end of file diff --git a/src/params.h b/src/params.h index e96fd1c..cad93bf 100644 --- a/src/params.h +++ b/src/params.h @@ -42,6 +42,6 @@ /* * Time limit for the computation (user time, in seconds). */ -#define MAX_TOTAL_TIME 3600 +#define MAX_TOTAL_TIME (24*3600) #endif //PROJECT_PARAMS_H