/* main.c */
/* Written by Andrew Makhorin <mao@gnu.org>, October 2015. */
/***********************************************************************
* This program is a stand-alone solver intended for solving Symmetric
* Traveling Salesman Problem (TSP) with the branch-and-bound method.
*
* Please note that this program is only an illustrative example. It is
* *not* a state-of-the-art code, so only small TSP instances (perhaps,
* having up to 150-200 nodes) can be solved with this code.
*
* To run this program use the following command:
*
* tspsol tsp-file
*
* where tsp-file specifies an input text file containing TSP data in
* TSPLIB 95 format.
*
* Detailed description of the input format recognized by this program
* is given in the report: Gerhard Reinelt, "TSPLIB 95". This report as
* well as TSPLIB, a library of sample TSP instances (and other related
* problems), are freely available for research purposes at the webpage
* <http://www.iwr.uni-heidelberg.de/groups/comopt/software/TSPLIB95/>.
*
* Symmetric Traveling Salesman Problem
* ------------------------------------
* Let a complete undirected graph be given:
*
* K = (V, E), (1)
*
* where V = {1, ..., n} is a set of nodes, E = V cross V is a set of
* edges. Let also each edge e = (i,j) be assigned a positive number
* c[i,j], which is the length of e. The Symmetric Traveling Salesman
* Problem (TSP) is to find a tour in K of minimal length.
*
* Integer programming formulation of TSP
* --------------------------------------
* For a set of nodes W within V introduce the following notation:
*
* d(W) = {(i,j):i in W and j not in W or i not in W and j in W}, (2)
*
* i.e. d(W) is the set of edges which have exactly one endnode in W.
* If W = {v}, i.e. W consists of the only node, we write simply d(v).
*
* The integer programming formulation of TSP is the following:
*
* minimize sum c[i,j] * x[i,j] (3)
* i,j
*
* subject to sum x[i,j] = 2 for all v in V (4)
* (i,j) in d(v)
*
* sum x[i,j] >= 2 for all W within V, (5)
* (i,j) in d(W) W != empty, W != V
*
* x[i,j] in {0, 1} for all i, j (6)
*
* The binary variables x[i,j] have conventional meaning: if x[i,j] = 1,
* the edge (i,j) is included in the tour, otherwise, if x[i,j] = 0, the
* edge is not included in the tour.
*
* The constraints (4) are called degree constraints. They require that
* for each node v in V there must be exactly two edges included in the
* tour which are incident to v.
*
* The constraints (5) are called subtour elimination constraints. They
* are used to forbid subtours. Note that the number of the subtour
* elimination constraints grows exponentially on the size of the TSP
* instance, so these constraints are not included explicitly in the
* IP, but generated dynamically during the B&B search.
***********************************************************************/
#include <errno.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <glpk.h>
#include "maxflow.h"
#include "mincut.h"
#include "misc.h"
#include "tsplib.h"
int n;
/* number of nodes in the problem, n >= 2 */
int *c; /* int c[1+n*(n-1)/2]; */
/* upper triangle (without diagonal entries) of the (symmetric) matrix
* C = (c[i,j]) in row-wise format, where c[i,j] specifies a length of
* edge e = (i,j), 1 <= i < j <= n */
int *tour; /* int tour[1+n]; */
/* solution to TSP, which is a tour specified by the list of node
* numbers tour[1] -> ... -> tour[nn] -> tour[1] in the order the nodes
* are visited; note that any tour is a permutation of node numbers */
glp_prob *P;
/* integer programming problem object */
/***********************************************************************
* loc - determine reduced index of element of symmetric matrix
*
* Given indices i and j of an element of a symmetric nxn-matrix,
* 1 <= i, j <= n, i != j, this routine returns the index of that
* element in an array, which is the upper triangle (without diagonal
* entries) of the matrix in row-wise format. */
int loc(int i, int j)
{ xassert(1 <= i && i <= n);
xassert(1 <= j && j <= n);
xassert(i != j);
if (i < j)
return ((n - 1) + (n - i + 1)) * (i - 1) / 2 + (j - i);
else
return loc(j, i);
}
/***********************************************************************
* read_data - read TSP data
*
* This routine reads TSP data from a specified text file in TSPLIB 95
* format. */
void read_data(const char *fname)
{ TSP *tsp;
int i, j;
tsp = tsp_read_data(fname);
if (tsp == NULL)
{ xprintf("TSP data file processing error\n");
exit(EXIT_FAILURE);
}
if (tsp->type != TSP_TSP)
{ xprintf("Invalid TSP data type\n");
exit(EXIT_FAILURE);
}
n = tsp->dimension;
xassert(n >= 2);
if (n > 32768)
{ xprintf("TSP instance too large\n");
exit(EXIT_FAILURE);
}
c = xalloc(1+loc(n-1, n), sizeof(int));
for (i = 1; i <= n; i++)
{ for (j = i+1; j <= n; j++)
c[loc(i, j)] = tsp_distance(tsp, i, j);
}
tsp_free_data(tsp);
return;
}
/***********************************************************************
* build_prob - build initial integer programming problem
*
* This routine builds the initial integer programming problem, which
* includes all variables (6), objective (3) and all degree constraints
* (4). Subtour elimination constraints (5) are considered "lazy" and
* not included in the initial problem. */
void build_prob(void)
{ int i, j, k, *ind;
double *val;
char name[50];
/* create problem object */
P = glp_create_prob();
/* add all binary variables (6) */
for (i = 1; i <= n; i++)
{ for (j = i+1; j <= n; j++)
{ k = glp_add_cols(P, 1);
xassert(k == loc(i,j));
sprintf(name, "x[%d,%d]", i, j);
glp_set_col_name(P, k, name);
glp_set_col_kind(P, k, GLP_BV);
/* set objective coefficient (3) */
glp_set_obj_coef(P, k, c[k]);
}
}
/* add all degree constraints (4) */
ind = xalloc(1+n, sizeof(int));
val = xalloc(1+n, sizeof(double));
for (i = 1; i <= n; i++)
{ k = glp_add_rows(P, 1);
xassert(k == i);
sprintf(name, "v[%d]", i);
glp_set_row_name(P, i, name);
glp_set_row_bnds(P, i, GLP_FX, 2, 2);
k = 0;
for (j = 1; j <= n; j++)
{ if (i != j)
k++, ind[k] = loc(i,j), val[k] = 1;
}
xassert(k == n-1);
glp_set_mat_row(P, i, n-1, ind, val);
}
xfree(ind);
xfree(val);
return;
}
/***********************************************************************
* build_tour - build tour for corresponding solution to IP
*
* Given a feasible solution to IP (3)-(6) this routine builds the
* corresponding solution to TSP, which is a tour specified by the list
* of node numbers tour[1] -> ... -> tour[nn] -> tour[1] in the order
* the nodes are to be visited */
void build_tour(void)
{ int i, j, k, kk, *beg, *end;
/* solution to MIP should be feasible */
switch (glp_mip_status(P))
{ case GLP_FEAS:
case GLP_OPT:
break;
default:
xassert(P != P);
}
/* build the list of edges included in the tour */
beg = xalloc(1+n, sizeof(int));
end = xalloc(1+n, sizeof(int));
k = 0;
for (i = 1; i <= n; i++)
{ for (j = i+1; j <= n; j++)
{ double x;
x = glp_mip_col_val(P, loc(i,j));
xassert(x == 0 || x == 1);
if (x)
{ k++;
xassert(k <= n);
beg[k] = i, end[k] = j;
}
}
}
xassert(k == n);
/* reorder edges in the list as they follow in the tour */
for (k = 1; k <= n; k++)
{ /* find k-th edge of the tour */
j = (k == 1 ? 1 : end[k-1]);
for (kk = k; kk <= n; kk++)
{ if (beg[kk] == j)
break;
if (end[kk] == j)
{ end[kk] = beg[kk], beg[kk] = j;
break;
}
}
xassert(kk <= n);
/* put the edge to k-th position in the list */
i = beg[k], beg[k] = beg[kk], beg[kk] = i;
j = end[k], end[k] = end[kk], end[kk] = j;
}
/* build the tour starting from node 1 */
xassert(beg[1] == 1);
for (k = 1; k <= n; k++)
{ if (k > 1)
xassert(end[k-1] == beg[k]);
tour[k] = beg[k];
}
xassert(end[n] == 1);
xfree(beg);
xfree(end);
return;
}
/***********************************************************************
* tour_length - calculate tour length
*
* This routine calculates the length of the specified tour, which is
* the sum of corresponding edge lengths. */
int tour_length(const int tour[/*1+n*/])
{ int i, j, sum;
sum = 0;
for (i = 1; i <= n; i++)
{ j = (i < n ? i+1 : 1);
sum += c[loc(tour[i], tour[j])];
}
return sum;
}
/***********************************************************************
* write_tour - write tour to text file in TSPLIB format
*
* This routine writes the specified tour to a text file in TSPLIB
* format. */
void write_tour(const char *fname, const int tour[/*1+n*/])
{ FILE *fp;
int i;
xprintf("Writing TSP solution to '%s'...\n", fname);
fp = fopen(fname, "w");
if (fp == NULL)
{ xprintf("Unable to create '%s' - %s\n", fname,
strerror(errno));
return;
}
fprintf(fp, "NAME : %s\n", fname);
fprintf(fp, "COMMENT : Tour length is %d\n", tour_length(tour));
fprintf(fp, "TYPE : TOUR\n");
fprintf(fp, "DIMENSION : %d\n", n);
fprintf(fp, "TOUR_SECTION\n");
for (i = 1; i <= n; i++)
fprintf(fp, "%d\n", tour[i]);
fprintf(fp, "-1\n");
fprintf(fp, "EOF\n");
fclose(fp);
return;
}
/***********************************************************************
* gen_subt_row - generate violated subtour elimination constraint
*
* This routine is called from the MIP solver to generate a violated
* subtour elimination constraint (5).
*
* Constraints of this class has the form:
*
* sum x[i,j] >= 2, i in W, j in V \ W,
*
* for all W, where W is a proper nonempty subset of V, V is the set of
* nodes of the given graph.
*
* In order to find a violated constraint of this class this routine
* finds a min cut in a capacitated network, which has the same sets of
* nodes and edges as the original graph, and where capacities of edges
* are values of variables x[i,j] in a basic solution to LP relaxation
* of the current subproblem. */
void gen_subt(glp_tree *T)
{ int i, j, ne, nz, *beg, *end, *cap, *cut, *ind;
double sum, *val;
/* MIP preprocessor should not be used */
xassert(glp_ios_get_prob(T) == P);
/* if some variable x[i,j] is zero in basic solution, then the
* capacity of corresponding edge in the associated network is
* zero, so we may not include such edge in the network */
/* count number of edges having non-zero capacity */
ne = 0;
for (i = 1; i <= n; i++)
{ for (j = i+1; j <= n; j++)
{ if (glp_get_col_prim(P, loc(i,j)) >= .001)
ne++;
}
}
/* build the capacitated network */
beg = xalloc(1+ne, sizeof(int));
end = xalloc(1+ne, sizeof(int));
cap = xalloc(1+ne, sizeof(int));
nz = 0;
for (i = 1; i <= n; i++)
{ for (j = i+1; j <= n; j++)
{ if (glp_get_col_prim(P, loc(i,j)) >= .001)
{ nz++;
xassert(nz <= ne);
beg[nz] = i, end[nz] = j;
/* scale all edge capacities to make them integral */
cap[nz] = ceil(1000 * glp_get_col_prim(P, loc(i,j)));
}
}
}
xassert(nz == ne);
/* find minimal cut in the capacitated network */
cut = xalloc(1+n, sizeof(int));
min_cut(n, ne, beg, end, cap, cut);
/* determine the number of non-zero coefficients in the subtour
* elimination constraint and calculate its left-hand side which
* is the (unscaled) capacity of corresponding min cut */
ne = 0, sum = 0;
for (i = 1; i <= n; i++)
{ for (j = i+1; j <= n; j++)
{ if (cut[i] && !cut[j] || !cut[i] && cut[j])
{ ne++;
sum += glp_get_col_prim(P, loc(i,j));
}
}
}
/* if the (unscaled) capacity of min cut is less than 2, the
* corresponding subtour elimination constraint is violated */
if (sum <= 1.999)
{ /* build the list of non-zero coefficients */
ind = xalloc(1+ne, sizeof(int));
val = xalloc(1+ne, sizeof(double));
nz = 0;
for (i = 1; i <= n; i++)
{ for (j = i+1; j <= n; j++)
{ if (cut[i] && !cut[j] || !cut[i] && cut[j])
{ nz++;
xassert(nz <= ne);
ind[nz] = loc(i,j);
val[nz] = 1;
}
}
}
xassert(nz == ne);
/* add violated tour elimination constraint to the current
* subproblem */
i = glp_add_rows(P, 1);
glp_set_row_bnds(P, i, GLP_LO, 2, 0);
glp_set_mat_row(P, i, nz, ind, val);
xfree(ind);
xfree(val);
}
/* free working arrays */
xfree(beg);
xfree(end);
xfree(cap);
xfree(cut);
return;
}
/***********************************************************************
* cb_func - application callback routine
*
* This routine is called from the MIP solver at various points of
* the branch-and-cut algorithm. */
void cb_func(glp_tree *T, void *info)
{ xassert(info == info);
switch (glp_ios_reason(T))
{ case GLP_IROWGEN:
/* generate one violated subtour elimination constraint */
gen_subt(T);
break;
}
return;
}
/***********************************************************************
* main - TSP solver main program
*
* This main program parses command-line arguments, reads specified TSP
* instance from a text file, and calls the MIP solver to solve it. */
int main(int argc, char *argv[])
{ int j;
char *in_file = NULL, *out_file = NULL;
time_t start;
glp_iocp iocp;
/* parse command-line arguments */
# define p(str) (strcmp(argv[j], str) == 0)
for (j = 1; j < argc; j++)
{ if (p("--output") || p("-o"))
{ j++;
if (j == argc || argv[j][0] == '\0' || argv[j][0] == '-')
{ xprintf("No solution output file specified\n");
exit(EXIT_FAILURE);
}
if (out_file != NULL)
{ xprintf("Only one solution output file allowed\n");
exit(EXIT_FAILURE);
}
out_file = argv[j];
}
else if (p("--help") || p("-h"))
{ xprintf("Usage: %s [options...] tsp-file\n", argv[0]);
xprintf("\n");
xprintf("Options:\n");
xprintf(" -o filename, --output filename\n");
xprintf(" write solution to filename\n")
;
xprintf(" -h, --help display this help information"
" and exit\n");
exit(EXIT_SUCCESS);
}
else if (argv[j][0] == '-' ||
(argv[j][0] == '-' && argv[j][1] == '-'))
{ xprintf("Invalid option '%s'; try %s --help\n", argv[j],
argv[0]);
exit(EXIT_FAILURE);
}
else
{ if (in_file != NULL)
{ xprintf("Only one input file allowed\n");
exit(EXIT_FAILURE);
}
in_file = argv[j];
}
}
if (in_file == NULL)
{ xprintf("No input file specified; try %s --help\n", argv[0]);
exit(EXIT_FAILURE);
}
# undef p
/* display program banner */
xprintf("TSP Solver for GLPK %s\n", glp_version());
/* remove output solution file specified in command-line */
if (out_file != NULL)
remove(out_file);
/* read TSP instance from input data file */
read_data(in_file);
/* build initial IP problem */
start = time(NULL);
build_prob();
tour = xalloc(1+n, sizeof(int));
/* solve LP relaxation of initial IP problem */
xprintf("Solving initial LP relaxation...\n");
xassert(glp_simplex(P, NULL) == 0);
xassert(glp_get_status(P) == GLP_OPT);
/* solve IP problem with "lazy" constraints */
glp_init_iocp(&iocp);
iocp.br_tech = GLP_BR_MFV; /* most fractional variable */
iocp.bt_tech = GLP_BT_BLB; /* best local bound */
iocp.sr_heur = GLP_OFF; /* disable simple rounding heuristic */
iocp.gmi_cuts = GLP_ON; /* enable Gomory cuts */
iocp.cb_func = cb_func;
glp_intopt(P, &iocp);
build_tour();
/* display some statistics */
xprintf("Time used: %.1f secs\n", difftime(time(NULL), start));
{ size_t tpeak;
glp_mem_usage(NULL, NULL, NULL, &tpeak);
xprintf("Memory used: %.1f Mb (%.0f bytes)\n",
(double)tpeak / 1048576.0, (double)tpeak);
}
/* write solution to output file, if required */
if (out_file != NULL)
write_tour(out_file, tour);
/* deallocate working objects */
xfree(c);
xfree(tour);
glp_delete_prob(P);
/* check that no memory blocks are still allocated */
{ int count;
size_t total;
glp_mem_usage(&count, NULL, &total, NULL);
if (count != 0)
xerror("Error: %d memory block(s) were lost\n", count);
xassert(total == 0);
}
return 0;
}
/* eof */