First version of the project

This commit is contained in:
2015-03-18 04:48:43 -04:00
commit ad8e20c62b
22 changed files with 1645 additions and 0 deletions

177
src/branch_and_cut.c Normal file
View File

@@ -0,0 +1,177 @@
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#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;
}

17
src/branch_and_cut.h Normal file
View File

@@ -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_

127
src/graph.c Normal file
View File

@@ -0,0 +1,127 @@
#include <malloc.h>
#include <math.h>
#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;
}

21
src/graph.h Normal file
View File

@@ -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_

175
src/lp.c Normal file
View File

@@ -0,0 +1,175 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cplex.h>
#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;
}

40
src/lp.h Normal file
View File

@@ -0,0 +1,40 @@
#ifndef __MAIN_
#define __MAIN_
#include <cplex.h>
#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

135
src/main.c Normal file
View File

@@ -0,0 +1,135 @@
#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
#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);
}

39
src/main.h Normal file
View File

@@ -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

418
src/tsp.c Normal file
View File

@@ -0,0 +1,418 @@
#include <stdio.h>
#include <stdlib.h>
#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;
}

39
src/tsp.h Normal file
View File

@@ -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_

115
src/util.c Normal file
View File

@@ -0,0 +1,115 @@
/****************************************************************************/
/* */
/* Utility Functions for CO759 */
/* */
/****************************************************************************/
#include "main.h"
#include <stdio.h>
#include <stdlib.h>
#include <sys/resource.h>
#include <stdarg.h>
#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;
}

17
src/util.h Normal file
View File

@@ -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 */