From ad8e20c62b7fc3945de65aa909c8468a62d0ff4b Mon Sep 17 00:00:00 2001 From: Alinson Xavier Date: Wed, 18 Mar 2015 04:48:43 -0400 Subject: [PATCH] First version of the project --- CMakeLists.txt | 28 +++ input/g10.22.edg | 23 +++ input/pr76.dat | 77 ++++++++ input/r10.dat | 11 ++ input/r15.dat | 16 ++ input/r20.dat | 21 +++ input/r25.dat | 26 +++ input/r30.dat | 31 ++++ input/r40.dat | 41 +++++ input/r50.dat | 51 ++++++ src/branch_and_cut.c | 177 ++++++++++++++++++ src/branch_and_cut.h | 17 ++ src/graph.c | 127 +++++++++++++ src/graph.h | 21 +++ src/lp.c | 175 ++++++++++++++++++ src/lp.h | 40 +++++ src/main.c | 135 ++++++++++++++ src/main.h | 39 ++++ src/tsp.c | 418 +++++++++++++++++++++++++++++++++++++++++++ src/tsp.h | 39 ++++ src/util.c | 115 ++++++++++++ src/util.h | 17 ++ 22 files changed, 1645 insertions(+) create mode 100644 CMakeLists.txt create mode 100644 input/g10.22.edg create mode 100644 input/pr76.dat create mode 100644 input/r10.dat create mode 100644 input/r15.dat create mode 100644 input/r20.dat create mode 100644 input/r25.dat create mode 100644 input/r30.dat create mode 100644 input/r40.dat create mode 100644 input/r50.dat create mode 100644 src/branch_and_cut.c create mode 100644 src/branch_and_cut.h create mode 100644 src/graph.c create mode 100644 src/graph.h create mode 100644 src/lp.c create mode 100644 src/lp.h create mode 100644 src/main.c create mode 100644 src/main.h create mode 100644 src/tsp.c create mode 100644 src/tsp.h create mode 100644 src/util.c create mode 100644 src/util.h diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..0d1f915 --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,28 @@ +cmake_minimum_required(VERSION 3.1) +project(project) + +include_directories(/opt/cplex-12.4/distribution/cplex/include/ilcplex) + +set(CMAKE_BINARY_DIR bin) +set(CMAKE_SOURCE_DIR src) + +set(CMAKE_C_FLAGS "-O3 -g -Wall -pedantic -g --std=c11 -Winline") + + +set(SOURCE_FILES + src/lp.c + src/lp.h + src/main.c + src/main.h + src/util.c + src/util.h + src/graph.c + src/graph.h + src/tsp.h + src/tsp.c + src/branch_and_cut.h + src/branch_and_cut.c) + +add_executable(project ${SOURCE_FILES} ) + +target_link_libraries(project cplex m pthread) \ No newline at end of file diff --git a/input/g10.22.edg b/input/g10.22.edg new file mode 100644 index 0000000..50af5d3 --- /dev/null +++ b/input/g10.22.edg @@ -0,0 +1,23 @@ +10 22 +0 9 1 +0 8 1 +0 1 2 +0 7 4 +0 4 3 +1 8 2 +1 4 3 +1 6 5 +2 5 8 +2 7 4 +2 3 2 +2 6 6 +3 5 6 +3 7 2 +4 7 4 +4 6 2 +5 8 4 +5 9 4 +5 7 5 +6 7 5 +7 9 4 +8 9 1 diff --git a/input/pr76.dat b/input/pr76.dat new file mode 100644 index 0000000..f016742 --- /dev/null +++ b/input/pr76.dat @@ -0,0 +1,77 @@ +76 +3600 2300 +3100 3300 +4700 5750 +5400 5750 +5608 7103 +4493 7102 +3600 6950 +3100 7250 +4700 8450 +5400 8450 +5610 10053 +4492 10052 +3600 10800 +3100 10950 +4700 11650 +5400 11650 +6650 10800 +7300 10950 +7300 7250 +6650 6950 +7300 3300 +6650 2300 +5400 1600 +8350 2300 +7850 3300 +9450 5750 +10150 5750 +10358 7103 +9243 7102 +8350 6950 +7850 7250 +9450 8450 +10150 8450 +10360 10053 +9242 10052 +8350 10800 +7850 10950 +9450 11650 +10150 11650 +11400 10800 +12050 10950 +12050 7250 +11400 6950 +12050 3300 +11400 2300 +10150 1600 +13100 2300 +12600 3300 +14200 5750 +14900 5750 +15108 7103 +13993 7102 +13100 6950 +12600 7250 +14200 8450 +14900 8450 +15110 10053 +13992 10052 +13100 10800 +12600 10950 +14200 11650 +14900 11650 +16150 10800 +16800 10950 +16800 7250 +16150 6950 +16800 3300 +16150 2300 +14900 1600 +19800 800 +19800 10000 +19800 11900 +19800 12200 +200 12200 +200 1100 +200 800 diff --git a/input/r10.dat b/input/r10.dat new file mode 100644 index 0000000..4810717 --- /dev/null +++ b/input/r10.dat @@ -0,0 +1,11 @@ +10 +64 96 +80 39 +69 23 +72 42 +48 67 +58 43 +81 34 +79 17 +30 23 +42 67 diff --git a/input/r15.dat b/input/r15.dat new file mode 100644 index 0000000..b8088d5 --- /dev/null +++ b/input/r15.dat @@ -0,0 +1,16 @@ +15 +95 53 +77 93 +0 28 +40 76 +76 6 +21 70 +26 78 +47 73 +32 5 +77 13 +85 32 +84 32 +43 13 +42 5 +96 29 diff --git a/input/r20.dat b/input/r20.dat new file mode 100644 index 0000000..89e35a2 --- /dev/null +++ b/input/r20.dat @@ -0,0 +1,21 @@ +20 +95 53 +77 93 +0 28 +40 76 +76 6 +21 70 +26 78 +47 73 +32 5 +77 13 +85 32 +84 32 +43 13 +42 5 +96 29 +84 21 +90 18 +98 44 +91 25 +73 70 diff --git a/input/r25.dat b/input/r25.dat new file mode 100644 index 0000000..62f75a0 --- /dev/null +++ b/input/r25.dat @@ -0,0 +1,26 @@ +25 +95 53 +77 93 +0 28 +40 76 +76 6 +21 70 +26 78 +47 73 +32 5 +77 13 +85 32 +84 32 +43 13 +42 5 +96 29 +84 21 +90 18 +98 44 +91 25 +73 70 +38 43 +34 1 +6 17 +75 94 +40 16 diff --git a/input/r30.dat b/input/r30.dat new file mode 100644 index 0000000..6c8c522 --- /dev/null +++ b/input/r30.dat @@ -0,0 +1,31 @@ +30 +95 53 +77 93 +0 28 +40 76 +76 6 +21 70 +26 78 +47 73 +32 5 +77 13 +85 32 +84 32 +43 13 +42 5 +96 29 +84 21 +90 18 +98 44 +91 25 +73 70 +38 43 +34 1 +6 17 +75 94 +40 16 +29 34 +28 68 +5 74 +70 66 +95 56 diff --git a/input/r40.dat b/input/r40.dat new file mode 100644 index 0000000..7df5348 --- /dev/null +++ b/input/r40.dat @@ -0,0 +1,41 @@ +40 +95 53 +77 93 +0 28 +40 76 +76 6 +21 70 +26 78 +47 73 +32 5 +77 13 +85 32 +84 32 +43 13 +42 5 +96 29 +84 21 +90 18 +98 44 +91 25 +73 70 +38 43 +34 1 +6 17 +75 94 +40 16 +29 34 +28 68 +5 74 +70 66 +95 56 +44 22 +10 6 +75 78 +36 32 +79 37 +98 38 +72 61 +91 51 +4 16 +34 76 diff --git a/input/r50.dat b/input/r50.dat new file mode 100644 index 0000000..b12077b --- /dev/null +++ b/input/r50.dat @@ -0,0 +1,51 @@ +50 +95 53 +77 93 +0 28 +40 76 +76 6 +21 70 +26 78 +47 73 +32 5 +77 13 +85 32 +84 32 +43 13 +42 5 +96 29 +84 21 +90 18 +98 44 +91 25 +73 70 +38 43 +34 1 +6 17 +75 94 +40 16 +29 34 +28 68 +5 74 +70 66 +95 56 +44 22 +10 6 +75 78 +36 32 +79 37 +98 38 +72 61 +91 51 +4 16 +34 76 +50 76 +10 40 +85 69 +18 84 +43 20 +15 59 +53 36 +72 0 +71 73 +17 55 diff --git a/src/branch_and_cut.c b/src/branch_and_cut.c new file mode 100644 index 0000000..02d14fb --- /dev/null +++ b/src/branch_and_cut.c @@ -0,0 +1,177 @@ +#include +#include +#include +#include "lp.h" +#include "branch_and_cut.h" +#include "tsp.h" +#include "util.h" + +int bnc_is_integral(double *x, int length); + +int bnc_find_best_branch_var(double *x, int length); + +int bnc_solve_node( + struct LP *lp, double *best_val, int ncount, int ecount, int *elist, + int depth) +{ + int rval = 0; + double *x = (double *) NULL; + + time_printf("Optimizing...\n"); + + int is_infeasible; + rval = lp_optimize(lp, &is_infeasible); + ABORT_IF (rval, "lp_optimize failed\n"); + + if (is_infeasible) + { + time_printf("Branch pruned by infeasibility.\n"); + goto CLEANUP; + } + + double objval; + rval = lp_get_obj_val(lp, &objval); + ABORT_IF (rval, "lp_get_obj_val failed\n"); + + time_printf(" objective value = %.2f\n", objval); + + if (objval > *best_val) + { + time_printf("Branch pruned by bound (%.2lf > %.2lf).\n", objval, + *best_val); + rval = 0; + goto CLEANUP; + } + + x = (double *) malloc(ecount * sizeof(double)); + ABORT_IF(!x, "could not allocate x\n"); + + rval = lp_get_x(lp, x); + ABORT_IF(rval, "lp_get_x failed\n"); + + rval = TSP_add_cutting_planes(ncount, ecount, elist, lp); + ABORT_IF(rval, "TSP_add_cutting_planes failed\n"); + + rval = lp_optimize(lp, &is_infeasible); + ABORT_IF (rval, "lp_optimize failed\n"); + + rval = lp_get_obj_val(lp, &objval); + ABORT_IF(rval, "lp_get_obj_val failed\n"); + + rval = lp_get_x(lp, x); + ABORT_IF(rval, "lp_get_x failed\n"); + + if (bnc_is_integral(x, ecount)) + { + time_printf(" solution is integral\n"); + + if (objval < *best_val) + { + *best_val = objval; + time_printf("Found a better integral solution:\n"); + time_printf(" objval = %.2lf **\n", objval); + } + } else + { + time_printf(" solution is fractional\n"); + rval = bnc_branch_node(lp, x, ncount, ecount, depth, best_val, elist); + ABORT_IF(rval, "bnc_branch_node failed\n"); + } + + CLEANUP: + if (x) free(x); + return rval; +} + +int bnc_branch_node( + struct LP *lp, double *x, int ncount, int ecount, int depth, + double *best_val, int *elist) +{ + int rval = 0; + + int best_index = bnc_find_best_branch_var(x, ecount); + + time_printf("Branching on variable x%d = %.6lf (depth %d)...\n", best_index, + x[best_index], depth); + + time_printf("Fixing variable x%d to one...\n", best_index); + rval = lp_change_bound(lp, best_index, 'L', 1.0); + ABORT_IF(rval, "lp_change_bound failed\n"); + + rval = bnc_solve_node(lp, best_val, ncount, ecount, elist, depth + 1); + ABORT_IF(rval, "bnc_solve_node failed\n"); + + rval = lp_change_bound(lp, best_index, 'L', 0.0); + ABORT_IF(rval, "lp_change_bound failed\n"); + + time_printf("Fixing variable x%d to zero...\n", best_index); + rval = lp_change_bound(lp, best_index, 'U', 0.0); + ABORT_IF(rval, "lp_change_bound failed\n"); + + rval = bnc_solve_node(lp, best_val, ncount, ecount, elist, depth + 1); + ABORT_IF(rval, "bnc_solve_node failed\n"); + + rval = lp_change_bound(lp, best_index, 'U', 1.0); + ABORT_IF(rval, "lp_change_bound failed\n"); + + time_printf("Finished branching on variable %d\n", best_index); + + CLEANUP: + return rval; +} + +int bnc_find_best_branch_var(double *x, int length) +{ + int best_index = 0; + double best_index_frac = 1.0; + + for (int j = 0; j < length; j++) + { + if (fabs(x[j] - 0.5) < best_index_frac) + { + best_index = j; + best_index_frac = fabs(x[j] - 0.5); + } + } + + return best_index; +} + +int bnc_is_integral(double *x, int length) +{ + int all_integral = 1; + + for (int j = 0; j < length; j++) + { + if (x[j] > LP_EPSILON && x[j] < 1.0 - LP_EPSILON) + { + all_integral = 0; + break; + } + } + + return all_integral; +} + +int bnc_init_lp( + struct LP *lp, int node_count, int edge_count, int *edge_list, + int *edge_weights) +{ + int rval = 0; + time_printf("Initializing LP...\n"); + + rval = lp_open(lp); + ABORT_IF(rval, "lp_open failed\n"); + + rval = lp_create(lp, "subtour"); + ABORT_IF(rval, "lp_create failed\n"); + + rval = TSP_init_lp(node_count, lp, edge_count, edge_weights, edge_list); + ABORT_IF(rval, "TSP_init_lp failed\n"); + + rval = lp_write(lp, "subtour.lp"); + ABORT_IF(rval, "lp_write failed\n"); + + CLEANUP: + return rval; +} \ No newline at end of file diff --git a/src/branch_and_cut.h b/src/branch_and_cut.h new file mode 100644 index 0000000..52ff26a --- /dev/null +++ b/src/branch_and_cut.h @@ -0,0 +1,17 @@ +#ifndef _PROJECT_BRANCH_AND_CUT_H_ +#define _PROJECT_BRANCH_AND_CUT_H_ + +#include "lp.h" + +int bnc_solve_node( + struct LP *lp, double *best_val, int ncount, int ecount, int *elist, + int depth); + +int bnc_branch_node( + struct LP *lp, double *x, int ncount, int ecount, int depth, + double *current_val, int *elist); + +int bnc_init_lp( + struct LP *lp, int node_count, int edge_count, int *edge_list, int *edge_weights); + +#endif //_PROJECT_BRANCH_AND_CUT_H_ diff --git a/src/graph.c b/src/graph.c new file mode 100644 index 0000000..62f8b36 --- /dev/null +++ b/src/graph.c @@ -0,0 +1,127 @@ +#include +#include +#include "main.h" +#include "graph.h" +#include "util.h" +#include "lp.h" + +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->node_list[n]; + pn->mark = 1; + + for (int i = 0; i < pn->deg; i++) + { + if (x[pn->adj[i].e] > LP_EPSILON) + { + int neighbor = pn->adj[i].n; + + if (G->node_list[neighbor].mark == 0) + graph_dfs(neighbor, G, x, island_size, island_nodes); + } + } +} + +void graph_init(struct Graph *G) +{ + if (!G) return; + + G->node_list = 0; + G->adj_space = 0; + G->node_count = 0; + G->edge_count = 0; +} + +void graph_free(struct Graph *G) +{ + if (!G) return; + + if (G->node_list) free(G->node_list); + if (G->adj_space) free(G->adj_space); +} + +int graph_build(int node_count, int edge_count, int *edge_list, struct Graph *G) +{ + int rval = 0; + struct Node *n; + struct AdjObj *p; + + G->node_list = (struct Node *) malloc(node_count * sizeof(struct Node)); + G->adj_space = + (struct AdjObj *) malloc(2 * edge_count * sizeof(struct AdjObj)); + + ABORT_IF(!G->node_list, "could not allocate G->node_list\n"); + ABORT_IF(!G->adj_space, "out of memory for node_list or adj_space\n"); + + for (int i = 0; i < node_count; i++) + G->node_list[i].deg = 0; + + for (int i = 0; i < edge_count; i++) + { + int a = edge_list[2 * i]; + int b = edge_list[2 * i + 1]; + G->node_list[a].deg++; + G->node_list[b].deg++; + } + + p = G->adj_space; + for (int i = 0; i < node_count; i++) + { + G->node_list[i].adj = p; + p += G->node_list[i].deg; + G->node_list[i].deg = 0; + } + + for (int i = 0; i < edge_count; i++) + { + int a = edge_list[2 * i]; + int b = edge_list[2 * i + 1]; + n = &G->node_list[a]; + n->adj[n->deg].n = b; + n->adj[n->deg].e = i; + n->deg++; + n = &G->node_list[b]; + n->adj[n->deg].n = a; + n->adj[n->deg].e = i; + n->deg++; + } + + G->node_count = node_count; + G->edge_count = edge_count; + + CLEANUP: + if (rval) + { + if (G->node_list) free(G->node_list); + if (G->adj_space) free(G->adj_space); + } + return rval; +} + +int euclid_edgelen(int i, int j, double *x, double *y) +{ + double t1 = x[i] - x[j], t2 = y[i] - y[j]; + return (int) (sqrt(t1 * t1 + t2 * t2) + 0.5); +} + +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; +} \ No newline at end of file diff --git a/src/graph.h b/src/graph.h new file mode 100644 index 0000000..d5f42ef --- /dev/null +++ b/src/graph.h @@ -0,0 +1,21 @@ +#ifndef __GRAPH_H_ +#define __GRAPH_H_ + +#include "main.h" + +void graph_dfs( + int n, struct Graph *G, double *x, int *icount, int *island); + +void graph_init(struct Graph *G); + +void graph_free(struct Graph *G); + +int graph_build(int node_count, int edge_count, int *edge_list, struct Graph *G); + +int euclid_edgelen(int i, int j, double *x, double *y); + +void get_delta( + int nsize, int *nlist, int ecount, int *elist, int *deltacount, + int *delta, int *marks); + +#endif //___GRAPH_H_ diff --git a/src/lp.c b/src/lp.c new file mode 100644 index 0000000..a5b9bd0 --- /dev/null +++ b/src/lp.c @@ -0,0 +1,175 @@ +#include +#include +#include +#include +#include "lp.h" +#include "util.h" + +#define LP_EPSILON 0.000001 + +int lp_open(struct LP *lp) +{ + int rval = 0; + + lp->cplex_lp = (CPXLPptr)NULL; + lp->cplex_env = CPXopenCPLEX(&rval); + if (rval) + { + fprintf(stderr, "CPXopenCPLEX failed, return code %d\n", rval); + goto CLEANUP; + } + + CLEANUP: + return rval; +} + +void lp_free(struct LP *lp) +{ + if (!lp) return; + if (!lp->cplex_env) return; + + if (lp->cplex_lp) + CPXfreeprob(lp->cplex_env, &(lp->cplex_lp)); + + CPXcloseCPLEX(&lp->cplex_env); + lp->cplex_env = 0; +} + +int lp_create(struct LP *lp, const char *name) +{ + int rval = 0; + char nambuf[MAX_NAME_LENGTH]; + + ABORT_IF(!lp->cplex_env, "cplex_env is null\n"); + + strncpy(nambuf, name, MAX_NAME_LENGTH); + nambuf[MAX_NAME_LENGTH - 1] = '\0'; + + lp->cplex_lp = CPXcreateprob(lp->cplex_env, &rval, nambuf); + ABORT_IF(rval, "CPXcreateprob failed\n"); + + CLEANUP: + return rval; +} + +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, + (double *) NULL, (char **) NULL); + + ABORT_IF(rval, "CPXnewrows failed\n"); + + CLEANUP: + return rval; +} + +int lp_add_rows( + struct LP *lp, int newrows, int newnz, double *rhs, char *sense, + int *rmatbeg, int *rmatind, double *rmatval) +{ + int rval = 0; + + rval = CPXaddrows(lp->cplex_env, lp->cplex_lp, 0, newrows, newnz, rhs, + sense, rmatbeg, rmatind, rmatval, (char **) NULL, (char **) NULL); + + ABORT_IF(rval, "CPXaddrows failed\n"); + + CLEANUP: + return rval; +} + +int lp_add_cols( + struct LP *lp, int newcols, int newnz, double *obj, int *cmatbeg, + int *cmatind, double *cmatval, double *lb, double *ub) +{ + int rval = 0; + + rval = CPXaddcols(lp->cplex_env, lp->cplex_lp, newcols, newnz, obj, cmatbeg, + cmatind, cmatval, lb, ub, (char **) NULL); + + ABORT_IF(rval, "CPXaddcols failed\n"); + + CLEANUP: + return rval; +} + +int lp_change_bound(struct LP *lp, int col, char lower_or_upper, double bnd) +{ + int rval = 0; + + rval = CPXchgbds(lp->cplex_env, lp->cplex_lp, 1, &col, &lower_or_upper, + &bnd); + ABORT_IF(rval, "CPXchgbds failed\n"); + + CLEANUP: + return rval; +} + +int lp_optimize(struct LP *lp, int *infeasible) +{ + int rval = 0, solstat; + + *infeasible = 0; + + rval = CPXdualopt(lp->cplex_env, lp->cplex_lp); + ABORT_IF(rval, "CPXdualopt failed\n"); + + solstat = CPXgetstat(lp->cplex_env, lp->cplex_lp); + if (solstat == CPX_STAT_INFEASIBLE) + { + if (infeasible) + *infeasible = 1; + + } + else + { + ABORT_IF(solstat != CPX_STAT_OPTIMAL + && solstat != CPX_STAT_OPTIMAL_INFEAS, + "Invalid solution status"); + } + + CLEANUP: + return rval; +} + +int lp_get_obj_val(struct LP *lp, double *obj) +{ + int rval = 0; + + rval = CPXgetobjval(lp->cplex_env, lp->cplex_lp, obj); + ABORT_IF(rval, "CPXgetobjval failed\n"); + + CLEANUP: + return rval; +} + +int lp_get_x(struct LP *lp, double *x) +{ + int rval = 0; + + int ncols = CPXgetnumcols(lp->cplex_env, lp->cplex_lp); + ABORT_IF(!ncols, "No columns in LP\n"); + + rval = CPXgetx(lp->cplex_env, lp->cplex_lp, x, 0, ncols - 1); + ABORT_IF(rval, "CPXgetx failed\n"); + + CLEANUP: + return rval; +} + +int lp_write(struct LP *lp, const char *fname) +{ + int rval = 0; + + char nambuf[MAX_NAME_LENGTH]; + strncpy(nambuf, fname, MAX_NAME_LENGTH); + nambuf[MAX_NAME_LENGTH - 1] = '\0'; + + rval = CPXwriteprob(lp->cplex_env, lp->cplex_lp, nambuf, "RLP"); + ABORT_IF(rval, "CPXwriteprob failed\n"); + + CLEANUP: + return rval; +} diff --git a/src/lp.h b/src/lp.h new file mode 100644 index 0000000..d181aef --- /dev/null +++ b/src/lp.h @@ -0,0 +1,40 @@ +#ifndef __MAIN_ +#define __MAIN_ + +#include + +#define LP_EPSILON 0.000001 + +struct LP +{ + CPXENVptr cplex_env; + CPXLPptr cplex_lp; +}; + +static const int MAX_NAME_LENGTH = 100; + +int lp_open(struct LP *lp); + +int lp_create(struct LP *lp, const char *name); + +int lp_write(struct LP *lp, const char *fname); + +void lp_free(struct LP *lp); + +int lp_new_row(struct LP *lp, char sense, double rhs); + +int lp_add_rows(struct LP *lp, int newrows, int newnz, double *rhs, char *sense, + int *rmatbeg, int *rmatind, double *rmatval); + +int lp_add_cols(struct LP *lp, int newcols, int newnz, double *obj, + int *cmatbeg, int *cmatind, double *cmatval, double *lb, double *ub); + +int lp_change_bound(struct LP *lp, int col, char lower_or_upper, double bnd); + +int lp_optimize(struct LP *lp, int *infeasible); + +int lp_get_obj_val(struct LP *lp, double *obj); + +int lp_get_x(struct LP *lp, double *x); + +#endif \ No newline at end of file diff --git a/src/main.c b/src/main.c new file mode 100644 index 0000000..dc436cf --- /dev/null +++ b/src/main.c @@ -0,0 +1,135 @@ +#include +#include +#include +#include "lp.h" +#include "util.h" +#include "main.h" +#include "tsp.h" +#include "branch_and_cut.h" + +char *fname = (char *) NULL; +int seed = 0; +int geometric_data = 0; +int ncount_rand = 0; +int gridsize_rand = 100; +int use_all_subtours = 0; + +int main(int ac, char **av) +{ + int rval = 0; + int *edge_weights = 0; + int *edge_list = 0; + int *tlist = 0; + + seed = (int) util_get_current_time(); + + rval = parseargs(ac, av); + ABORT_IF(rval, "Failed to parse arguments.\n"); + ABORT_IF(!fname && !ncount_rand, "Must specify a problem file or use -k for" + " random prob\n"); + + printf("Seed = %d\n", seed); + srand(seed); + + if (fname) + { + printf("Problem name: %s\n", fname); + if (geometric_data) printf("Geometric data\n"); + } + + int node_count = 0, edge_count = 0; + rval = TSP_read_problem(fname, &node_count, &edge_count, &edge_list, + &edge_weights); + ABORT_IF(rval, "TSP_read_problem failed\n"); + ABORT_IF(use_all_subtours && node_count > 20, "Too many nodes to add all" + " subtours\n"); + + tlist = (int *) malloc((node_count) * sizeof(int)); + + ABORT_IF(!tlist, "out of memory for tlist\n"); + + double best_val = TSP_find_initial_solution(edge_weights, edge_list, + node_count, edge_count); + + initial_time = util_get_current_time(); + + struct LP *lp; + lp = (struct LP *) malloc(sizeof(struct LP *)); + + rval = bnc_init_lp(lp, node_count, edge_count, edge_list, edge_weights); + ABORT_IF(rval, "bnc_init_lp failed\n"); + + rval = bnc_solve_node(lp, &best_val, node_count, edge_count, edge_list, 1); + ABORT_IF(rval, "bnc_solve_node failed\n"); + + time_printf("Optimal integral solution:\n"); + time_printf(" objective value = %.2lf **\n", best_val); + + printf("\nRunning Time: %.2f seconds\n", + util_get_current_time() - initial_time); + fflush(stdout); + + CLEANUP: + if (tlist) free(tlist); + if (edge_list) free(edge_list); + if (edge_weights) free(edge_weights); + return rval; +} + +int parseargs(int ac, char **av) +{ + int c; + + if (ac == 1) + { + usage(av[0]); + return 1; + } + + while ((c = getopt(ac, av, "ab:gk:s:")) != EOF) + { + switch (c) + { + case 'a': + use_all_subtours = 1; + break; + case 'b': + gridsize_rand = atoi(optarg); + break; + case 'g': + geometric_data = 1; + break; + case 'k': + ncount_rand = atoi(optarg); + break; + case 's': + seed = atoi(optarg); + break; + case '?': + default: + usage(av[0]); + return 1; + } + } + + if (optind < ac) fname = av[optind++]; + + if (optind != ac) + { + usage(av[0]); + return 1; + } + + return 0; +} + +void 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); +} + diff --git a/src/main.h b/src/main.h new file mode 100644 index 0000000..9b0eb14 --- /dev/null +++ b/src/main.h @@ -0,0 +1,39 @@ +#ifndef __MAIN_H_ +#define __MAIN_H_ + +struct AdjObj +{ + int n; + /* index of neighbor node */ + int e; /* index of adj joining neighbor */ +}; + +struct Node +{ + int deg; + struct AdjObj *adj; + int mark; +}; + +struct Graph +{ + int node_count; + int edge_count; + struct Node *node_list; + struct AdjObj *adj_space; +}; + +void usage(char *f); + +int parseargs(int ac, char **av); + +extern char *fname; +extern int seed; +extern int geometric_data; +extern int ncount_rand; +extern int gridsize_rand; +extern int use_all_subtours; + +extern double initial_time; + +#endif \ No newline at end of file diff --git a/src/tsp.c b/src/tsp.c new file mode 100644 index 0000000..43d899f --- /dev/null +++ b/src/tsp.c @@ -0,0 +1,418 @@ +#include +#include +#include "main.h" +#include "tsp.h" +#include "util.h" + +int TSP_init_lp( + int node_count, struct LP *lp, int edge_count, int *edge_weights, + int *edge_list) +{ + 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\n"); + } + + /* Build a column for each edge 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\n"); + } + + CLEANUP: + return rval; +} + +int TSP_find_violated_subtour_elimination_cut( + int ncount, int edge_count, int *edges, struct LP *lp) +{ + 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\n"); + ABORT_IF(is_infeasible, "LP is infeasible\n"); + + rval = graph_build(ncount, edge_count, edges, &G); + ABORT_IF(rval, "graph_build failed\n"); + + 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\n"); + + int round = 0; + int delta_count = 0; + int island_count = 0; + + while (!TSP_is_graph_connected(&G, x, &island_count, island_sizes, + island_start, island_nodes)) + { + time_printf("Adding %d bnc_solve_node 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); + } + + time_printf("Reoptimizing (round %d)...\n", ++round); + ABORT_IF(rval, "TSP_add_subtour_elimination_cut failed"); + + rval = lp_optimize(lp, &is_infeasible); + ABORT_IF(rval, "lp_optimize failed\n"); + ABORT_IF(is_infeasible, "LP is infeasible\n"); + + double objval = 0; + rval = lp_get_obj_val(lp, &objval); + ABORT_IF(rval, "lp_get_obj_val failed\n"); + + rval = lp_get_x(lp, x); + ABORT_IF(rval, "lp_get_x failed\n"); + } + + time_printf(" graph is TSP_is_graph_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 deltacount, int *delta) +{ + int rval = 0; + char sense = 'G'; + double rhs = 2.0; + int rmatbeg = 0; + + double *rmatval; + int *rmatind = delta; + + rmatval = (double *) malloc(deltacount * sizeof(double)); + ABORT_IF(!rmatval, "out of memory for rmatval\n"); + + for (int i = 0; i < deltacount; i++) + rmatval[i] = 1.0; + + rval = lp_add_rows(lp, 1, deltacount, &rhs, &sense, &rmatbeg, rmatind, + rmatval); + + ABORT_IF(rval, "lp_add_rows failed"); + + CLEANUP: + if (rmatval) free(rmatval); + return rval; +} + +int TSP_is_graph_connected( + 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->node_list[i].mark = 0; + island_nodes[i] = -1; + } + + int k = 0, current_island = 0; + + for (int i = 0; i < G->node_count; i++) + { + if (G->node_list[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 (*island_count == 1); +} + +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, &G); + ABORT_IF(rval, "graph_build failed\n"); + + for (int j = 0; j < node_count; j++) + G.node_list[j].mark = 0; + + for (int j = 0; j < node_count; j++) + { + if (j == node_count - 1) + G.node_list[start].mark = 0; + + struct Node *pn = &G.node_list[current_node]; + pn->mark = 1; + + int closest_neighbor = -1; + int closest_edge_length = 10000000; + + for (int i = 0; i < pn->deg; i++) + { + int edge = pn->adj[i].e; + int neighbor = pn->adj[i].n; + + if (G.node_list[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, int *p_ncount, int *p_ecount, int **p_elist, + int **p_elen) +{ + struct _IO_FILE *f = (struct _IO_FILE *) NULL; + int i, j, end1, end2, w, rval = 0, ncount, ecount; + int *elist = (int *) NULL, *elen = (int *) NULL; + 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", &ncount, &ecount) != 2) + { + fprintf(stderr, "Input file %s has invalid format\n", filename); + rval = 1; + goto CLEANUP; + } + + printf("Nodes: %d Edges: %d\n", ncount, ecount); + fflush(stdout); + + elist = (int *) malloc(2 * ecount * sizeof(int)); + if (!elist) + { + fprintf(stderr, "out of memory for elist\n"); + rval = 1; + goto CLEANUP; + } + + elen = (int *) malloc(ecount * sizeof(int)); + if (!elen) + { + fprintf(stderr, "out of memory for elen\n"); + rval = 1; + goto CLEANUP; + } + + for (i = 0; i < ecount; i++) + { + if (fscanf(f, "%d %d %d", &end1, &end2, &w) != 3) + { + fprintf(stderr, "%s has invalid input format\n", filename); + rval = 1; + goto CLEANUP; + } + elist[2 * i] = end1; + elist[2 * i + 1] = end2; + elen[i] = w; + } + } + else + { + if (filename) + { + if (fscanf(f, "%d", &ncount) != 1) + { + fprintf(stderr, "Input file %s has invalid format\n", filename); + rval = 1; + goto CLEANUP; + } + } + else + { + ncount = ncount_rand; + } + + x = (double *) malloc(ncount * sizeof(double)); + y = (double *) malloc(ncount * 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 < ncount; 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 = CO759_build_xy(ncount, x, y, gridsize_rand); + if (rval) + { + fprintf(stderr, "CO759_build_xy failed\n"); + goto CLEANUP; + } + + printf("%d\n", ncount); + for (i = 0; i < ncount; i++) + { + printf("%.0f %.0f\n", x[i], y[i]); + } + printf("\n"); + } + + ecount = (ncount * (ncount - 1)) / 2; + time_printf("Complete graph: %d nodes, %d edges\n", ncount, ecount); + + elist = (int *) malloc(2 * ecount * sizeof(int)); + if (!elist) + { + fprintf(stderr, "out of memory for elist\n"); + rval = 1; + goto CLEANUP; + } + + elen = (int *) malloc(ecount * sizeof(int)); + if (!elen) + { + fprintf(stderr, "out of memory for elen\n"); + rval = 1; + goto CLEANUP; + } + + ecount = 0; + for (i = 0; i < ncount; i++) + { + for (j = i + 1; j < ncount; j++) + { + elist[2 * ecount] = i; + elist[2 * ecount + 1] = j; + elen[ecount] = euclid_edgelen(i, j, x, y); + ecount++; + } + } + } + + *p_ncount = ncount; + *p_ecount = ecount; + *p_elist = elist; + *p_elen = elen; + + CLEANUP: + if (f) fclose(f); + if (x) free(x); + if (y) free(y); + return rval; +} + +int TSP_add_cutting_planes(int ncount, int ecount, int *elist, struct LP *lp) +{ + int rval = 0; + + rval = TSP_find_violated_subtour_elimination_cut(ncount, ecount, elist, lp); + ABORT_IF (rval, "TSP_find_violated_subtour_elimination_cut failed\n"); + + CLEANUP: + return rval; +} + +double TSP_find_initial_solution( + int *edge_weights, int *edge_list, int node_count, int edge_count) +{ + double best_val = 1e99; + + time_printf("Finding closest neighbor tour\n"); + for (int i = 0; i < node_count; i++) + { + int path_length = 0; + + TSP_find_closest_neighbor_tour(i, node_count, edge_count, edge_list, + edge_weights, &path_length); + + if (best_val > path_length) best_val = path_length; + } + + time_printf(" length = %lf\n", best_val); + + return best_val; +} \ No newline at end of file diff --git a/src/tsp.h b/src/tsp.h new file mode 100644 index 0000000..da149b2 --- /dev/null +++ b/src/tsp.h @@ -0,0 +1,39 @@ +// +// Created by isoron on 3/17/15. +// + +#ifndef ___TSP_H_ +#define ___TSP_H_ + +#include "lp.h" +#include "graph.h" + +int add_all_subtours(int ncount, int ecount, int *elist, struct LP *lp); + +int TSP_find_violated_subtour_elimination_cut + (int ncount, int ecount, int *elist, struct LP *lp); + +int TSP_is_graph_connected( + 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 deltacount, int *delta); + +int TSP_read_problem( + char *filename, int *p_ncount, int *p_ecount, int **p_elist, + int **p_elen); + +int TSP_add_cutting_planes(int ncount, int ecount, int *elist, struct LP *lp); + +int TSP_init_lp( + int node_count, struct LP *lp, int edge_count, int *edge_weights, + int *edge_list); + +double TSP_find_initial_solution + (int *edge_weights, int *edge_list, int node_count, int edge_count); + +#endif //_PROJECT_TSP_H_ diff --git a/src/util.c b/src/util.c new file mode 100644 index 0000000..7822140 --- /dev/null +++ b/src/util.c @@ -0,0 +1,115 @@ +/****************************************************************************/ +/* */ +/* Utility Functions for CO759 */ +/* */ +/****************************************************************************/ + +#include "main.h" +#include +#include +#include +#include +#include "util.h" + +double util_get_current_time(void) +{ + struct rusage ru; + + getrusage (RUSAGE_SELF, &ru); + + return ((double) ru.ru_utime.tv_sec) + + ((double) ru.ru_utime.tv_usec)/1000000.0; +} + +/* function for creating a random set of points in unit square */ + +int CO759_build_xy (int ncount, double *xlist, double *ylist, int gridsize) +{ + int rval = 0, i, j, winner, x, y; + int **hit = (int **) NULL, *hitcount = (int *) NULL; + + printf ("Random %d point set, gridsize = %d\n", ncount, gridsize); + fflush (stdout); + + hit = (int **) malloc (gridsize * sizeof (int *)); + if (!hit) { + fprintf (stderr, "out of memory for hit\n"); + rval = 1; goto CLEANUP; + } + for (i = 0; i < gridsize; i++) hit[i] = (int *) NULL; + + hitcount = (int *) malloc (gridsize * sizeof (int)); + if (!hitcount) { + fprintf (stderr, "out of memory for hitcount\n"); + rval = 1; goto CLEANUP; + } + for (i = 0; i < gridsize; i++) hitcount[i] = 0; + + for (i = 0; i < ncount; i++) { + winner = 0; + do { + x = (int) (rand () % gridsize); + y = (int) (rand () % gridsize); + + /* check to see if (x,y) is a duplicate point */ + + for (j = 0; j < hitcount[x]; j++) { + if (hit[x][j] == y) break; + } + if (j == hitcount[x]) { + void *tmp_ptr = (void *) hit[x]; + tmp_ptr = realloc (tmp_ptr, (hitcount[x]+1)*sizeof (int)); + if (!tmp_ptr) { + fprintf (stderr, "out of member in realloc of hit\n"); + rval = 1; goto CLEANUP; + } + hit[x] = (int *) tmp_ptr; + hit[x][hitcount[x]] = y; + hitcount[x]++; + winner = 1; + } + if (!winner) { + printf ("X"); fflush (stdout); + } + } while (!winner); + xlist[i] = (double) x; + ylist[i] = (double) y; + } + +CLEANUP: + + printf ("\n"); + + if (hit) { + for (i = 0; i < gridsize; i++) { + if (hit[i]) free (hit[i]); + } + free (hit); + } + if (hitcount) free (hitcount); + return rval; +} + +double initial_time = 0; + +void time_printf(const char *fmt, ...) +{ + if(initial_time == 0) + initial_time = util_get_current_time(); + + printf("[%10.2lf] ", util_get_current_time() - initial_time); + + va_list args; + va_start(args, fmt); + vprintf(fmt, args); + va_end(args); + + fflush(stdout); +} + +void next_set (int sz, int *Set) +{ + int i; + for (i=0; i < sz-1 && Set[i]+1 == Set[i+1]; i++) Set[i] = i; + Set[i] = Set[i]+1; +} \ No newline at end of file diff --git a/src/util.h b/src/util.h new file mode 100644 index 0000000..f871309 --- /dev/null +++ b/src/util.h @@ -0,0 +1,17 @@ +#ifndef __CO759_UTIL_H +#define __CO759_UTIL_H + +#define ABORT_IF(cond, msg) if(cond) { \ + fprintf(stderr, msg); rval = 1; goto CLEANUP; } + +double util_get_current_time(void); + +int CO759_build_xy(int ncount, double *xlist, double *ylist, int gridsize); + +double util_get_current_time(void); + +void time_printf(const char *fmt, ...); + +void next_set(int sz, int *Set); + +#endif /* __CO759_UTIL_H */