isoron/gtsp
Archived
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Refactoring; remove unused code

Browse Source
This commit is contained in:
Alinson S. Xavier
2015-04-04 06:57:34 -04:00
parent b802256426
commit 2bc1b897bb
22 changed files with 850 additions and 1900 deletions
Download Patch File Download Diff File Expand all files Collapse all files

147
src/branch_and_cut.c
View File

@@ -6,9 +6,6 @@
#include "util.h" #include "util.h"
#include "gtsp.h" #include "gtsp.h"
int BNC_NODE_COUNT = 0;
static int BNC_solve_node(struct BNC *bnc, int depth); static int BNC_solve_node(struct BNC *bnc, int depth);
static int BNC_branch_node(struct BNC *bnc, double *x, int depth); static int BNC_branch_node(struct BNC *bnc, double *x, int depth);
@@ -17,8 +14,6 @@ static int BNC_is_integral(double *x, int num_cols);
static int BNC_find_best_branching_var(double *x, int num_cols); static int BNC_find_best_branching_var(double *x, int num_cols);
//int optimize_vertex_in_cluster(struct BNC *bnc, double best_val);
int BNC_init(struct BNC *bnc) int BNC_init(struct BNC *bnc)
{ {
int rval = 0; int rval = 0;
@@ -101,7 +96,7 @@ static int BNC_solve_node(struct BNC *bnc, int depth)
if(depth == 1) ROOT_VALUE = objval; if(depth == 1) ROOT_VALUE = objval;
if (ceil(objval) > *best_val + LP_EPSILON) if (ceil(objval) > *best_val + EPSILON)
{ {
log_debug("Branch pruned by bound (%.2lf > %.2lf).\n", objval, log_debug("Branch pruned by bound (%.2lf > %.2lf).\n", objval,
*best_val); *best_val);
@@ -127,9 +122,9 @@ static int BNC_solve_node(struct BNC *bnc, int depth)
abort_if(rval, "LP_get_obj_val failed"); abort_if(rval, "LP_get_obj_val failed");
if(depth == 1) ROOT_VALUE = objval; if(depth == 1) ROOT_VALUE = objval;
abort_if(get_current_time() - INITIAL_TIME >= MAX_TOTAL_TIME, "time limit exceeded"); abort_if(get_user_time() - INITIAL_TIME >= MAX_TOTAL_TIME, "time limit exceeded");
if (ceil(objval) > *best_val + LP_EPSILON) if (ceil(objval) > *best_val + EPSILON)
{ {
log_debug("Branch pruned by bound (%.2lf > %.2lf).\n", objval, log_debug("Branch pruned by bound (%.2lf > %.2lf).\n", objval,
*best_val); *best_val);
@@ -145,7 +140,7 @@ static int BNC_solve_node(struct BNC *bnc, int depth)
{ {
log_debug("Solution is integral\n"); log_debug("Solution is integral\n");
if (objval + LP_EPSILON < *best_val) if (objval + EPSILON < *best_val)
{ {
if (bnc->best_x) free(bnc->best_x); if (bnc->best_x) free(bnc->best_x);
*best_val = objval; *best_val = objval;
@@ -214,17 +209,11 @@ static int BNC_branch_node(struct BNC *bnc, double *x, int depth)
static int BNC_is_integral(double *x, int num_cols) static int BNC_is_integral(double *x, int num_cols)
{ {
#ifdef ALLOW_FRACTIONAL_SOLUTIONS
UNUSED(num_cols);
UNUSED(x);
return 1;
#else
for (int i = 0; i < num_cols; i++) for (int i = 0; i < num_cols; i++)
if (x[i] > LP_EPSILON && x[i] < 1.0 - LP_EPSILON) if (x[i] > EPSILON && x[i] < 1.0 - EPSILON)
return 0; return 0;
return 1; return 1;
#endif
} }
static int BNC_find_best_branching_var(double *x, int num_cols) static int BNC_find_best_branching_var(double *x, int num_cols)
@@ -244,129 +233,3 @@ static int BNC_find_best_branching_var(double *x, int num_cols)
return best_index; return best_index;
} }
/*
int re_optimize_integral(struct BNC *bnc){
int i = 0 , current_vertex = 0, rval = 0;
struct GTSP* data;
data = bnc->problem_data;
int node_count = data->graph->node_count;
int cluster_count = data->cluster_count;
int edge_count = data->graph->edge_count;
struct Tour * tour = (struct Tour*) NULL;
//intialize the tour
tour = (struct Tour *) malloc( cluster_count * sizeof(struct Tour));
for (i = 0; i < edge_count; i++){
tour[i].vertex = -1;
tour[i].next = -1;
tour[i].prev = -1;
}
//Constructing the tour with vertices
for (i = 0; i < edge_count; i++){
if (bnc->best_x[i + node_count] > LP_EPSILON) {
tour[current_vertex].vertex = data->graph->edges[i].from->index;
current_vertex += 1;
printf("From node %d \t", data->graph->edges[i].from->index);
printf("TO node %d \n", data->graph->edges[i].to->index);
}
}
//printf("Edgese in solution %d \n", current_vertex);
return rval;
CLEANUP:
if (data) free(data);
}
*/
/*
int optimize_vertex_in_cluster(struct BNC *bnc, double best_val)
{
int i = 0 , j, current_vertex = 0, rval = 0;
int tour_cost = 0;
struct GTSP* data;
data = bnc->problem_data;
//rval = GTSP_init_data(&data);
//data = bnc->problem_data;
//data = (struct GTSP) malloc(sizeof(struct GTSP));
//data = &bnc->problem_data;
int node_count = data->graph->node_count;
int cluster_count = data->cluster_count;
int edge_count = data->graph->edge_count;
int * tour = (int*) NULL;
tour = (int *) malloc( cluster_count * sizeof(int));
//Constructing the tour with vertices
for (i = 0; i < edge_count; i++)
{ //printf(" edge %lf **\n", bnc->best_x[i]);
if ((bnc->best_x[i] > 1 - LP_EPSILON)){
//printf(" x[i] = %lf **\n", bnc->best_x[i]);
tour[current_vertex] = (data->graph->edges[i].from)->index;
current_vertex += 1;
//printf(" Edge No = %d **\n", i);
printf(" FROM No = %d **\n", (data->graph->edges[i].from)->index);
printf(" TO No = %d **\n", (data->graph->edges[i].to)->index);
//printf(" current vertex = %d **\n", current_vertex);
}
}
//reoptmizing the your with two-opt
//rval = two_opt(cluster_count, tour, data->dist_matrix);
//Optimizing the vertices inside the node_to_cluster
int current_cluster = 0;
int insertion_cost = 0;
//printf(" o-- val = %.2lf **\n", best_val);
for(i = 1; i < cluster_count - 2; i++){
//printf(" vertex in tour = %d **\n", tour[current_vertex]);
current_cluster = data->node_to_cluster[tour[i]];
//printf(" o-- val = %.2lf **\n", best_val);
insertion_cost = data->dist_matrix[tour[i-1]][tour[i]] +
data->dist_matrix[tour[i]][tour[i+1]];
//printf(" o-- val = %.2lf **\n", best_val);
for(j = 0; j < node_count; j++)
if (current_cluster == data->node_to_cluster[j])
if (insertion_cost > data->dist_matrix[j][tour[i]] +
data->dist_matrix[j][tour[i+1]]){
log_info("Optmize vertex in cluster improved the bound\n");
insertion_cost = data->dist_matrix[j][tour[i]] +
data->dist_matrix[j][tour[i+1]];
tour[i] = j;
}
}
printf(" o-- val = %.2lf **\n", best_val);
for(i = 0; i< cluster_count ; i++){
if (i == cluster_count - 1)
tour_cost += data->dist_matrix[tour[i]][tour[0]];
else
tour_cost += data->dist_matrix[tour[i]][tour[i+1]];
if(tour_cost < bnc->best_obj_val)
bnc->best_obj_val = tour_cost;
}
return rval;
}
*/
/*
static int two_opt(int tour_length, int*tour, int** dist_matrix){
int rval = 0, i;
for (i = 1; i < tour_length - 2; i++){
int current_cost = dist_matrix[tour[i-1]][tour[i]] +
dist_matrix[tour[i+1]][tour[i+2]];
int temp_cost = dist_matrix[tour[i-1]][tour[i+1]] +
dist_matrix[tour[i]][tour[i+2]];
if(current_cost > temp_cost){
log_info("Two opt improved the bound\n");
int temp_vertex = tour[i];
tour[i] = tour[i+1];
tour[i+1] = temp_vertex;
}
}
return rval;
}
*/

6
src/branch_and_cut.h
View File

@@ -10,8 +10,6 @@ struct BNC
double *best_x; double *best_x;
double best_obj_val; double best_obj_val;
double *optimal_x;
int *problem_data; int *problem_data;
int (*problem_init_lp)(struct LP *, void *); int (*problem_init_lp)(struct LP *, void *);
@@ -29,10 +27,6 @@ int BNC_init_lp(struct BNC *bnc);
void BNC_free(struct BNC *bnc); void BNC_free(struct BNC *bnc);
int re_optimize_integral(struct BNC *bnc);
//int optimize_vertex_in_cluster(struct BNC *bnc, double best_val);
extern int BNC_NODE_COUNT; extern int BNC_NODE_COUNT;
#endif //_PROJECT_BRANCH_AND_CUT_H_ #endif //_PROJECT_BRANCH_AND_CUT_H_

113
src/flow.c
View File

@@ -142,7 +142,7 @@ int flow_find_max_flow(
struct Edge *e = &digraph->edges[i]; struct Edge *e = &digraph->edges[i];
#endif #endif
if (residual_caps[i] < LP_EPSILON) continue; if (residual_caps[i] < EPSILON) continue;
log_verbose("%d %d %.4lf (%d)\n", e->from->index, e->to->index, e->index, log_verbose("%d %d %.4lf (%d)\n", e->from->index, e->to->index, e->index,
residual_caps[e->index]); residual_caps[e->index]);
@@ -206,7 +206,7 @@ int flow_find_augmenting_path(
struct Edge *edge = n->adj[i].edge; struct Edge *edge = n->adj[i].edge;
if (neighbor->mark > 0) continue; if (neighbor->mark > 0) continue;
if (residual_caps[edge->index] < LP_EPSILON) continue; if (residual_caps[edge->index] < EPSILON) continue;
parents[neighbor->index] = n; parents[neighbor->index] = n;
parent_edges[neighbor->index] = edge; parent_edges[neighbor->index] = edge;
@@ -248,112 +248,3 @@ int flow_find_augmenting_path(
if (queue) free(queue); if (queue) free(queue);
return rval; return rval;
} }
int flow_main(int argc, char **argv)
{
UNUSED(argc);
UNUSED(argv);
int rval = 0;
int *edges = 0;
double *capacities = 0;
double *flow = 0;
double flow_value;
struct Edge **cut_edges = 0;
FILE *f = fopen("tmp/flow.in", "r");
abort_if(!f, "could not open input file");
struct Graph graph;
graph_init(&graph);
int node_count, edge_count;
rval = fscanf(f, "%d %d ", &node_count, &edge_count);
abort_if(rval != 2, "invalid input format (node count, edge count)");
int sink, source;
rval = fscanf(f, "%d %d", &source, &sink);
abort_if(rval != 2, "invalid input format (source, sink)");
edges = (int *) malloc(4 * edge_count * sizeof(int));
abort_if(!edges, "could not allocate edges\n");
capacities = (double *) malloc(2 * edge_count * sizeof(double));
abort_if(!capacities, "could not allocate capacities");
for (int i = 0; i < edge_count; i++)
{
int from, to;
double cap;
rval = fscanf(f, "%d %d %lf ", &from, &to, &cap);
abort_if(rval != 3, "invalid input format (edge specification)");
edges[i * 4] = edges[i * 4 + 3] = from;
edges[i * 4 + 1] = edges[i * 4 + 2] = to;
capacities[2 * i] = cap;
capacities[2 * i + 1] = 0;
}
rval = graph_build(node_count, 2 * edge_count, edges, 1, &graph);
abort_if(rval, "graph_build failed");
for (int i = 0; i < edge_count; i++)
{
graph.edges[2 * i].reverse = &graph.edges[2 * i + 1];
graph.edges[2 * i + 1].reverse = &graph.edges[2 * i];
}
flow = (double *) malloc(graph.edge_count * sizeof(double));
abort_if(!flow, "could not allocate flow");
struct Node *from = &graph.nodes[source];
struct Node *to = &graph.nodes[sink];
rval = flow_find_max_flow(&graph, capacities, from, to, flow, &flow_value);
abort_if(rval, "flow_find_max_flow failed");
log_info("Optimal flow has value %f\n", flow_value);
for (int i = 0; i < graph.edge_count; i++)
{
struct Edge *e = &graph.edges[i];
if (flow[e->index] <= 0) continue;
log_info(" %d %d %6.2f / %6.2f\n", e->from->index, e->to->index,
flow[e->index], capacities[e->index]);
}
log_info("Nodes reachable from origin on residual graph:\n");
for (int i = 0; i < graph.node_count; i++)
{
struct Node *n = &graph.nodes[i];
if (n->mark)
log_info(" %d\n", n->index);
}
int cut_edges_count = 0;
cut_edges =
(struct Edge **) malloc(graph.edge_count * sizeof(struct Edge *));
abort_if(!cut_edges, "could not allocate cut_edges");
rval = get_cut_edges_from_marks(&graph, &cut_edges_count, cut_edges);
abort_if(rval, "get_cut_edges_from_marks failed");
log_info("Min cut edges:\n");
for (int i = 0; i < cut_edges_count; i++)
{
struct Edge *e = cut_edges[i];
if (capacities[e->index] <= 0) continue;
log_info(" %d %d\n", e->from->index, e->to->index);
}
CLEANUP:
if (cut_edges) free(cut_edges);
if (capacities) free(capacities);
if (edges) free(edges);
if (flow) free(flow);
return rval;
}

4
src/flow.h
View File

@@ -23,8 +23,4 @@ int flow_find_max_flow(
int flow_mark_reachable_nodes( int flow_mark_reachable_nodes(
const struct Graph *graph, double *residual_caps, struct Node *from); const struct Graph *graph, double *residual_caps, struct Node *from);
int flow_main(int argc, char **argv);
#include "graph.h"
#endif //_PROJECT_FLOW_H_ #endif //_PROJECT_FLOW_H_

4
src/geometry.h
View File

@@ -24,5 +24,7 @@ int get_euclidean_distance(
int generate_dist_matrix( int generate_dist_matrix(
int node_count, int node_count,
double *x_coordinates, double *x_coordinates,
double *y_coordinates, int** dist_matrix); double *y_coordinates,
int **dist_matrix);
#endif //_PROJECT_GEOMETRY_H_ #endif //_PROJECT_GEOMETRY_H_

45
src/graph.c
View File

@@ -1,7 +1,6 @@
#include <malloc.h> #include <malloc.h>
#include "graph.h" #include "graph.h"
#include "util.h" #include "util.h"
#include "lp.h"
void graph_init(struct Graph *graph) void graph_init(struct Graph *graph)
{ {
@@ -109,50 +108,6 @@ int graph_build(
return rval; return rval;
} }
void graph_dfs(
int n, struct Graph *G, double *x, int *island_size, int *island_nodes)
{
*(island_nodes + (*island_size)) = n;
(*island_size)++;
struct Node *pn = &G->nodes[n];
pn->mark = 1;
for (int i = 0; i < pn->degree; i++)
{
if (x[pn->adj[i].edge_index] > LP_EPSILON)
{
int neighbor = pn->adj[i].neighbor_index;
if (G->nodes[neighbor].mark == 0)
graph_dfs(neighbor, G, x, island_size, island_nodes);
}
}
}
void get_delta(
int island_node_count,
int *island_nodes,
int edge_count,
int *edges,
int *delta_count,
int *delta,
int *marks)
{
for (int i = 0; i < island_node_count; i++)
marks[island_nodes[i]] = 1;
int k = 0;
for (int i = 0; i < edge_count; i++)
if (marks[edges[2 * i]] + marks[edges[2 * i + 1]] == 1)
delta[k++] = i;
*delta_count = k;
for (int i = 0; i < island_node_count; i++)
marks[island_nodes[i]] = 0;
}
int get_cut_edges_from_marks( int get_cut_edges_from_marks(
struct Graph *graph, int *cut_edges_count, struct Edge **cut_edges) struct Graph *graph, int *cut_edges_count, struct Edge **cut_edges)
{ {

12
src/graph.h
View File

@@ -47,9 +47,6 @@ struct Graph
struct Adjacency *adj; struct Adjacency *adj;
}; };
void graph_dfs(
int n, struct Graph *G, double *x, int *icount, int *island);
void graph_init(struct Graph *graph); void graph_init(struct Graph *graph);
void graph_free(struct Graph *graph); void graph_free(struct Graph *graph);
@@ -61,15 +58,6 @@ int graph_build(
int is_directed, int is_directed,
struct Graph *graph); struct Graph *graph);
void get_delta(
int nsize,
int *nlist,
int ecount,
int *elist,
int *deltacount,
int *delta,
int *marks);
int get_cut_edges_from_marks( int get_cut_edges_from_marks(
struct Graph *graph, int *cut_edges_count, struct Edge **cut_edges); struct Graph *graph, int *cut_edges_count, struct Edge **cut_edges);

123
src/gtsp-comb.c
View File

@@ -2,7 +2,7 @@
#include "util.h" #include "util.h"
#include <assert.h> #include <assert.h>
int add_comb_cut( static int add_comb_cut(
struct LP *lp, struct LP *lp,
struct Graph *graph, struct Graph *graph,
int current_component, int current_component,
@@ -73,7 +73,7 @@ int add_comb_cut(
{ {
for (int i = 0; i < nz; i++) for (int i = 0; i < nz; i++)
{ {
if (OPTIMAL_X[rmatind[i]] < LP_EPSILON) continue; if (OPTIMAL_X[rmatind[i]] < EPSILON) continue;
if (rmatind[i] >= node_count) if (rmatind[i] >= node_count)
{ {
@@ -98,7 +98,7 @@ int add_comb_cut(
for (int i = 0; i < nz; i++) for (int i = 0; i < nz; i++)
sum += rmatval[i] * OPTIMAL_X[rmatind[i]]; sum += rmatval[i] * OPTIMAL_X[rmatind[i]];
log_verbose("%.2lf >= %.2lf\n", sum, rhs); log_verbose("%.2lf >= %.2lf\n", sum, rhs);
abort_if(sum <= rhs - LP_EPSILON, "cannot add invalid cut"); abort_if(sum <= rhs - EPSILON, "cannot add invalid cut");
} }
double lhs = 0.0; double lhs = 0.0;
@@ -107,7 +107,7 @@ int add_comb_cut(
log_verbose("Violation: %.4lf >= %.4lf\n", lhs, rhs); log_verbose("Violation: %.4lf >= %.4lf\n", lhs, rhs);
if (lhs + LP_EPSILON > rhs) if (lhs + EPSILON > rhs)
{ {
free(rmatind); free(rmatind);
free(rmatval); free(rmatval);
@@ -130,7 +130,7 @@ int add_comb_cut(
return rval; return rval;
} }
int find_components( static int find_components(
struct Graph *graph, double *x, int *components, int *component_sizes) struct Graph *graph, double *x, int *components, int *component_sizes)
{ {
int rval = 0; int rval = 0;
@@ -168,8 +168,8 @@ int find_components(
if (neighbor->mark) continue; if (neighbor->mark) continue;
double x_e = x[adj->edge->index]; double x_e = x[adj->edge->index];
if (x_e < LP_EPSILON) continue; if (x_e < EPSILON) continue;
if (x_e > 1 - LP_EPSILON) continue; if (x_e > 1 - EPSILON) continue;
stack[stack_top++] = neighbor; stack[stack_top++] = neighbor;
neighbor->mark = 1; neighbor->mark = 1;
@@ -196,7 +196,7 @@ int find_components(
return rval; return rval;
} }
int find_teeth( static int find_teeth(
struct Graph *graph, struct Graph *graph,
double *x, double *x,
int current_component, int current_component,
@@ -221,7 +221,7 @@ int find_teeth(
struct Node *from = e->from; struct Node *from = e->from;
struct Node *to = e->to; struct Node *to = e->to;
if (x[e->index] < 1 - LP_EPSILON) continue; if (x[e->index] < 1 - EPSILON) continue;
if (to->mark || from->mark) continue; if (to->mark || from->mark) continue;
@@ -242,11 +242,6 @@ int find_teeth(
return 0; return 0;
} }
int write_shrunken_graph(
double *shrunken_x,
struct Graph *shrunken_graph,
int const cluster_count);
static int shrink_clusters( static int shrink_clusters(
const struct GTSP *data, const struct GTSP *data,
double *x, double *x,
@@ -346,10 +341,53 @@ static int shrink_clusters(
return rval; return rval;
} }
int find_comb_cuts(struct LP *lp, struct GTSP *data) static int write_shrunken_graph(
double *shrunken_x,
struct Graph *shrunken_graph,
int const cluster_count)
{ {
int rval = 0; int rval = 0;
double initial_time = get_current_time();
FILE *file = 0;
file = fopen("gtsp-shrunken.in", "w");
abort_if(!file, "could not open file");
fprintf(file, "%d %d\n", (*shrunken_graph).node_count, cluster_count);
for (int i = 0; i < (*shrunken_graph).node_count; i++)
fprintf(file, "%.2lf %.2lf %d\n", (*shrunken_graph).x_coordinates[i],
(*shrunken_graph).y_coordinates[i], i);
fclose(file);
file = fopen("gtsp-shrunken.out", "w");
abort_if(!file, "could not open file");
int positive_edge_count = 0;
for (int i = 0; i < (*shrunken_graph).edge_count; i++)
if (shrunken_x[i] > EPSILON)
positive_edge_count++;
fprintf(file, "%d %d\n", (*shrunken_graph).node_count,
(*shrunken_graph).edge_count);
fprintf(file, "%d\n", positive_edge_count);
for (int i = 0; i < (*shrunken_graph).edge_count; i++)
if (shrunken_x[i] > EPSILON)
fprintf(file, "%d %d %.4lf\n",
(*shrunken_graph).edges[i].from->index,
(*shrunken_graph).edges[i].to->index, shrunken_x[i]);
fclose(file);
CLEANUP:
return rval;
}
int GTSP_find_comb_cuts(struct LP *lp, struct GTSP *data)
{
int rval = 0;
double initial_time = get_user_time();
double *x = 0; double *x = 0;
double *shrunken_x = 0; double *shrunken_x = 0;
@@ -380,6 +418,8 @@ int find_comb_cuts(struct LP *lp, struct GTSP *data)
#if LOG_LEVEL >= LOG_LEVEL_DEBUG #if LOG_LEVEL >= LOG_LEVEL_DEBUG
rval = write_shrunken_graph(shrunken_x, &shrunken_graph, cluster_count); rval = write_shrunken_graph(shrunken_x, &shrunken_graph, cluster_count);
abort_if(rval, "write_shrunken_graph failed"); abort_if(rval, "write_shrunken_graph failed");
#else
UNUSED(write_shrunken_graph);
#endif #endif
teeth = (int *) malloc(cluster_count * sizeof(int)); teeth = (int *) malloc(cluster_count * sizeof(int));
@@ -414,8 +454,8 @@ int find_comb_cuts(struct LP *lp, struct GTSP *data)
if (tooth_count % 2 == 0) if (tooth_count % 2 == 0)
{ {
for (int i = 0; i < cluster_count; i++) for (int k = 0; k < cluster_count; k++)
if (teeth[i] == tooth_count - 1) teeth[i] = -1; if (teeth[k] == tooth_count - 1) teeth[k] = -1;
tooth_count--; tooth_count--;
} }
@@ -428,7 +468,7 @@ int find_comb_cuts(struct LP *lp, struct GTSP *data)
int added_cuts_count = lp->cut_pool_size - original_cut_pool_size; int added_cuts_count = lp->cut_pool_size - original_cut_pool_size;
log_debug(" %d combs\n", added_cuts_count); log_debug(" %d combs\n", added_cuts_count);
COMBS_TIME += get_current_time() - initial_time; COMBS_TIME += get_user_time() - initial_time;
COMBS_COUNT += added_cuts_count; COMBS_COUNT += added_cuts_count;
CLEANUP: CLEANUP:
@@ -440,48 +480,3 @@ int find_comb_cuts(struct LP *lp, struct GTSP *data)
if (x) free(x); if (x) free(x);
return rval; return rval;
} }
int write_shrunken_graph(
double *shrunken_x,
struct Graph *shrunken_graph,
int const cluster_count)
{
int rval = 0;
FILE *file = 0;
file = fopen("gtsp-shrunken.in", "w");
abort_if(!file, "could not open file");
fprintf(file, "%d %d\n", (*shrunken_graph).node_count, cluster_count);
for (int i = 0; i < (*shrunken_graph).node_count; i++)
{
fprintf(file, "%.2lf %.2lf %d\n", (*shrunken_graph).x_coordinates[i],
(*shrunken_graph).y_coordinates[i], i);
}
fclose(file);
file = fopen("gtsp-shrunken.out", "w");
abort_if(!file, "could not open file");
int positive_edge_count = 0;
for (int i = 0; i < (*shrunken_graph).edge_count; i++)
if (shrunken_x[i] > LP_EPSILON)
positive_edge_count++;
fprintf(file, "%d %d\n", (*shrunken_graph).node_count,
(*shrunken_graph).edge_count);
fprintf(file, "%d\n", positive_edge_count);
for (int i = 0; i < (*shrunken_graph).edge_count; i++)
if (shrunken_x[i] > LP_EPSILON)
fprintf(file, "%d %d %.4lf\n",
(*shrunken_graph).edges[i].from->index,
(*shrunken_graph).edges[i].to->index, shrunken_x[i]);
fclose(file);
CLEANUP:
return rval;
}

2
src/gtsp-comb.h
View File

@@ -1,6 +1,6 @@
#ifndef PROJECT_GTSP_COMB_H #ifndef PROJECT_GTSP_COMB_H
#define PROJECT_GTSP_COMB_H #define PROJECT_GTSP_COMB_H
int find_comb_cuts(struct LP *lp, struct GTSP *data); int GTSP_find_comb_cuts(struct LP *lp, struct GTSP *data);
#endif //PROJECT_GTSP_COMB_H #endif //PROJECT_GTSP_COMB_H

170
src/gtsp-subtour.c
View File

@@ -1,12 +1,32 @@
#include "gtsp-subtour.h" #include "gtsp-subtour.h"
#include <assert.h> #include <assert.h>
#include <float.h>
#include "util.h" #include "util.h"
#include "flow.h" #include "flow.h"
extern double SUBTOUR_TIME; static void deactivate_cluster_node(double *capacities, struct Node *cluster_node)
{
for (int i = 0; i < cluster_node->degree; i++)
{
struct Adjacency *adj = &cluster_node->adj[i];
struct Edge *e = adj->edge;
int static build_flow_digraph( capacities[e->index] = 0;
}
}
static void activate_cluster_node(double *capacities, struct Node *cluster_node)
{
for (int i = 0; i < cluster_node->degree; i++)
{
struct Adjacency *adj = &cluster_node->adj[i];
struct Edge *e = adj->edge;
capacities[e->index] = 1e10;
capacities[e->reverse->index] = 1e10;
}
}
static int build_flow_digraph(
struct GTSP *data, double *x, struct Graph *digraph, double *capacities) struct GTSP *data, double *x, struct Graph *digraph, double *capacities)
{ {
int rval = 0; int rval = 0;
@@ -16,9 +36,8 @@ int static build_flow_digraph(
int node_count = data->graph->node_count; int node_count = data->graph->node_count;
struct Graph *graph = data->graph; struct Graph *graph = data->graph;
int digraph_node_count = node_count + data->cluster_count + 1; int digraph_node_count = node_count + data->cluster_count;
int digraph_edge_count = 4 * graph->edge_count + 2 * graph->node_count + int digraph_edge_count = 4 * graph->edge_count + 2 * graph->node_count;
2 * data->cluster_count;
digraph_edges = (int *) malloc(2 * digraph_edge_count * sizeof(int)); digraph_edges = (int *) malloc(2 * digraph_edge_count * sizeof(int));
abort_if(!digraph_edges, "could not allocate digraph_edges"); abort_if(!digraph_edges, "could not allocate digraph_edges");
@@ -28,7 +47,7 @@ int static build_flow_digraph(
int kc = 0; int kc = 0;
for (int i = 0; i < graph->edge_count; i++) for (int i = 0; i < graph->edge_count; i++)
{ {
if (x[i] < LP_EPSILON) continue; if (x[i] < EPSILON) continue;
struct Edge *e = &graph->edges[i]; struct Edge *e = &graph->edges[i];
int from = e->from->index; int from = e->from->index;
@@ -67,19 +86,6 @@ int static build_flow_digraph(
capacities[kc++] = 0; capacities[kc++] = 0;
} }
// Create an extra root node and connect it to each cluster node through
// some edge with zero capacity
for (int i = 0; i < data->cluster_count; i++)
{
digraph_edges[ke++] = node_count + i;
digraph_edges[ke++] = node_count + data->cluster_count;
capacities[kc++] = 0;
digraph_edges[ke++] = node_count + data->cluster_count;
digraph_edges[ke++] = node_count + i;
capacities[kc++] = 0;
}
assert(ke <= 2 * digraph_edge_count); assert(ke <= 2 * digraph_edge_count);
assert(kc <= digraph_edge_count); assert(kc <= digraph_edge_count);
@@ -99,10 +105,8 @@ int static build_flow_digraph(
return rval; return rval;
} }
int add_subtour_cut( static int add_subtour_cut(
struct LP *lp, struct LP *lp, struct Edge **cut_edges, int cut_edges_count)
struct Edge **cut_edges,
int cut_edges_count)
{ {
int rval = 0; int rval = 0;
@@ -135,7 +139,7 @@ int add_subtour_cut(
double sum = 0; double sum = 0;
for (int i = 0; i < newnz; i++) for (int i = 0; i < newnz; i++)
sum += rmatval[i] * OPTIMAL_X[rmatind[i]]; sum += rmatval[i] * OPTIMAL_X[rmatind[i]];
abort_if(sum <= rhs - LP_EPSILON, "cannot add invalid cut"); abort_if(sum <= rhs - EPSILON, "cannot add invalid cut");
} }
struct Row *cut = 0; struct Row *cut = 0;
@@ -155,19 +159,27 @@ int add_subtour_cut(
return rval; return rval;
} }
int find_exact_subtour_cuts( int GTSP_find_exact_subtour_cuts(struct LP *lp, struct GTSP *data)
struct LP *lp, struct GTSP *data, double min_cut_violation)
{ {
int rval = 0; int rval = 0;
double *x = 0; double *x = 0;
double *flow = 0;
double *capacities = 0; double *capacities = 0;
double initial_time = get_current_time(); struct Edge **cut_edges = 0;
double initial_time = get_user_time();
int added_cuts_count = 0; int added_cuts_count = 0;
struct Graph *graph = data->graph; struct Graph *graph = data->graph;
const int edge_count = graph->edge_count;
const int node_count = graph->node_count;
const int cluster_count = data->cluster_count;
int num_cols = LP_get_num_cols(lp); int num_cols = LP_get_num_cols(lp);
int digraph_edge_count = 4 * edge_count + 2 * node_count;
int original_cut_pool_size = lp->cut_pool_size;
x = (double *) malloc(num_cols * sizeof(double)); x = (double *) malloc(num_cols * sizeof(double));
abort_if(!x, "could not allocate x"); abort_if(!x, "could not allocate x");
@@ -181,78 +193,36 @@ int find_exact_subtour_cuts(
struct Graph digraph; struct Graph digraph;
graph_init(&digraph); graph_init(&digraph);
int digraph_edge_count = 4 * graph->edge_count + 2 * graph->node_count +
2 * data->cluster_count;
int original_cut_pool_size = lp->cut_pool_size;
flow = (double *) malloc(digraph_edge_count * sizeof(double));
capacities = (double *) malloc(digraph_edge_count * sizeof(double)); capacities = (double *) malloc(digraph_edge_count * sizeof(double));
cut_edges = (struct Edge **) malloc(edge_count * sizeof(struct Edge *));
abort_if(!flow, "could not allocate flow");
abort_if(!capacities, "could not allocate capacities"); abort_if(!capacities, "could not allocate capacities");
abort_if(!cut_edges, "could not allocate cut_edges");
rval = build_flow_digraph(data, x, &digraph, capacities); rval = build_flow_digraph(data, x, &digraph, capacities);
abort_if(rval, "build_flow_digraph failed"); abort_if(rval, "build_flow_digraph failed");
rval = find_exact_subtour_cuts_cluster_to_cluster(lp, data, &digraph,
capacities, min_cut_violation);
abort_if(rval, "find_exact_subtour_cuts_cluster_to_cluster failed");
added_cuts_count = lp->cut_pool_size - original_cut_pool_size;
log_debug(" %d cluster-to-cluster\n", added_cuts_count);
SUBTOUR_CLUSTER_CLUSTER_COUNT += added_cuts_count;
if (added_cuts_count > 0)
goto CLEANUP;
SUBTOUR_TIME += get_current_time() - initial_time;
CLEANUP:
graph_free(&digraph);
if (capacities) free(capacities);
if (x) free(x);
return rval;
}
int find_exact_subtour_cuts_cluster_to_cluster(
struct LP *lp,
struct GTSP *data,
struct Graph *digraph,
double *capacities,
double min_cut_violation)
{
int rval = 0;
double *flow = 0;
struct Edge **cut_edges = 0;
int cuts_count = 0; int cuts_count = 0;
struct Graph *graph = data->graph; for (int i = 0; i < cluster_count; i++)
cut_edges = (struct Edge **) malloc(
graph->edge_count * sizeof(struct Edge *));
flow = (double *) malloc(digraph->edge_count * sizeof(double));
abort_if(!cut_edges, "could not allocate cut_edges");
abort_if(!flow, "could not allocate flow");
struct Node *root_node = &digraph->nodes[graph->node_count +
data->cluster_count];
for (int i = 0; i < data->cluster_count; i++)
{ {
for (int j = i + 1; j < data->cluster_count; j++) for (int j = i + 1; j < cluster_count; j++)
{ {
struct Node *from = &digraph->nodes[graph->node_count + i]; struct Node *from = &digraph.nodes[node_count + i];
struct Node *to = &digraph->nodes[graph->node_count + j]; struct Node *to = &digraph.nodes[node_count + j];
double flow_value; double flow_value;
int cut_edges_count; int cut_edges_count;
activate_cluster_node(capacities, from); activate_cluster_node(capacities, from);
activate_cluster_node(capacities, to); activate_cluster_node(capacities, to);
deactivate_cluster_node(capacities, root_node);
log_verbose("Sending flow from cluster %d to cluster %d\n", i, j); log_verbose("Sending flow from cluster %d to cluster %d\n", i, j);
rval = flow_find_max_flow(digraph, capacities, from, to, flow, rval = flow_find_max_flow(&digraph, capacities, from, to, flow,
&flow_value); &flow_value);
abort_if(rval, "flow_find_max_flow failed"); abort_if(rval, "flow_find_max_flow failed");
@@ -262,13 +232,13 @@ int find_exact_subtour_cuts_cluster_to_cluster(
deactivate_cluster_node(capacities, from); deactivate_cluster_node(capacities, from);
deactivate_cluster_node(capacities, to); deactivate_cluster_node(capacities, to);
if (flow_value >= 2 - min_cut_violation) continue; if (flow_value >= 2 - EPSILON) continue;
log_verbose("Marked nodes:\n"); log_verbose("Marked nodes:\n");
for (int k = 0; k < graph->node_count; k++) for (int k = 0; k < node_count; k++)
{ {
graph->nodes[k].mark = digraph->nodes[k].mark; graph->nodes[k].mark = digraph.nodes[k].mark;
if (digraph->nodes[k].mark) log_verbose(" %d\n", k); if (digraph.nodes[k].mark) log_verbose(" %d\n", k);
} }
rval = get_cut_edges_from_marks(graph, &cut_edges_count, cut_edges); rval = get_cut_edges_from_marks(graph, &cut_edges_count, cut_edges);
@@ -286,31 +256,17 @@ int find_exact_subtour_cuts_cluster_to_cluster(
} }
} }
added_cuts_count = lp->cut_pool_size - original_cut_pool_size;
log_debug(" %d cluster-to-cluster\n", added_cuts_count);
SUBTOUR_COUNT += added_cuts_count;
SUBTOUR_TIME += get_user_time() - initial_time;
CLEANUP: CLEANUP:
graph_free(&digraph);
if (capacities) free(capacities);
if (cut_edges) free(cut_edges); if (cut_edges) free(cut_edges);
if (flow) free(flow); if (flow) free(flow);
if (x) free(x);
return rval; return rval;
} }
void deactivate_cluster_node(double *capacities, struct Node *cluster_node)
{
for (int i = 0; i < cluster_node->degree; i++)
{
struct Adjacency *adj = &cluster_node->adj[i];
struct Edge *e = adj->edge;
capacities[e->index] = 0;
}
}
void activate_cluster_node(double *capacities, struct Node *cluster_node)
{
for (int i = 0; i < cluster_node->degree; i++)
{
struct Adjacency *adj = &cluster_node->adj[i];
struct Edge *e = adj->edge;
capacities[e->index] = 1e10;
capacities[e->reverse->index] = 1e10;
}
}

32
src/gtsp-subtour.h
View File

@@ -3,36 +3,6 @@
#include "gtsp.h" #include "gtsp.h"
void deactivate_cluster_node(double *capacities, struct Node *cluster_node); int GTSP_find_exact_subtour_cuts(struct LP *lp, struct GTSP *data);
void activate_cluster_node(double *capacities, struct Node *cluster_node);
int find_exact_subtour_cuts_cluster_to_cluster(
struct LP *lp,
struct GTSP *data,
struct Graph *digraph,
double *capacities,
double min_cut_violation);
int find_exact_subtour_cuts_node_to_cluster(
struct LP *lp,
struct GTSP *data,
double *x,
struct Graph *digraph,
double *capacities,
double min_cut_violation);
int find_exact_subtour_cuts_node_to_node(
struct LP *lp,
struct GTSP *data,
double *x,
struct Graph *digraph,
double *capacities,
double min_cut_violation);
int find_exact_subtour_cuts(
struct LP *lp,
struct GTSP *data,
double min_cut_violation);
#endif //_PROJECT_GTSP_SUBTOUR_H_ #endif //_PROJECT_GTSP_SUBTOUR_H_

409
src/gtsp.c
View File

@@ -1,17 +1,16 @@
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <getopt.h>
#include <math.h>
#include <assert.h> #include <assert.h>
#include <sys/stat.h>
#include <float.h>
#include "gtsp.h" #include "gtsp.h"
#include "geometry.h" #include "geometry.h"
#include "util.h" #include "util.h"
#include "gtsp-subtour.h" #include "gtsp-subtour.h"
#include "gtsp-comb.h" #include "gtsp-comb.h"
int large_neighborhood_search(struct Tour *tour, struct GTSP *data, int *tour_cost); int large_neighborhood_search(
struct Tour *tour,
struct GTSP *data,
int *tour_cost);
int build_edge_map(struct GTSP *gtsp, int *edge_map); int build_edge_map(struct GTSP *gtsp, int *edge_map);
@@ -231,7 +230,8 @@ int GTSP_add_cutting_planes(struct LP *lp, struct GTSP *data)
} }
current_round++; current_round++;
abort_if(get_current_time() - INITIAL_TIME >= MAX_TOTAL_TIME, "time limit exceeded"); abort_if(get_user_time() - INITIAL_TIME >= MAX_TOTAL_TIME,
"time limit exceeded");
int original_cut_pool_size; int original_cut_pool_size;
int added_cuts_count; int added_cuts_count;
@@ -239,24 +239,22 @@ int GTSP_add_cutting_planes(struct LP *lp, struct GTSP *data)
original_cut_pool_size = lp->cut_pool_size; original_cut_pool_size = lp->cut_pool_size;
log_debug("Finding subtour cuts, round %d...\n", current_round); log_debug("Finding subtour cuts, round %d...\n", current_round);
rval = find_exact_subtour_cuts(lp, data, LP_EPSILON); rval = GTSP_find_exact_subtour_cuts(lp, data);
abort_if(rval, "find_exact_subtour_cuts failed"); abort_if(rval, "GTSP_find_exact_subtour_cuts failed");
added_cuts_count = lp->cut_pool_size - original_cut_pool_size; added_cuts_count = lp->cut_pool_size - original_cut_pool_size;
if (added_cuts_count > 0) if (added_cuts_count > 0)
continue; continue;
#ifdef ENABLE_COMB_INEQUALITIES
original_cut_pool_size = lp->cut_pool_size; original_cut_pool_size = lp->cut_pool_size;
log_debug("Finding comb cuts, round %d...\n", current_round); log_debug("Finding comb cuts, round %d...\n", current_round);
rval = find_comb_cuts(lp, data); rval = GTSP_find_comb_cuts(lp, data);
abort_if(rval, "find_comb_cuts failed"); abort_if(rval, "GTSP_find_comb_cuts failed");
added_cuts_count = lp->cut_pool_size - original_cut_pool_size; added_cuts_count = lp->cut_pool_size - original_cut_pool_size;
if (added_cuts_count > 0) if (added_cuts_count > 0)
continue; continue;
#endif
break; break;
} }
@@ -303,7 +301,7 @@ int GTSP_write_solution(struct GTSP *data, char *filename, double *x)
int positive_edge_count = 0; int positive_edge_count = 0;
for (int i = 0; i < edge_count; i++) for (int i = 0; i < edge_count; i++)
if (x[i] > LP_EPSILON) if (x[i] > EPSILON)
positive_edge_count++; positive_edge_count++;
fprintf(file, "%d %d\n", node_count, edge_count); fprintf(file, "%d %d\n", node_count, edge_count);
@@ -311,7 +309,7 @@ int GTSP_write_solution(struct GTSP *data, char *filename, double *x)
fprintf(file, "%d\n", positive_edge_count); fprintf(file, "%d\n", positive_edge_count);
for (int i = 0; i < edge_count; i++) for (int i = 0; i < edge_count; i++)
if (x[i] > LP_EPSILON) if (x[i] > EPSILON)
fprintf(file, "%d %d %.4lf\n", edges[i].from->index, fprintf(file, "%d %d %.4lf\n", edges[i].from->index,
edges[i].to->index, x[i]); edges[i].to->index, x[i]);
@@ -368,7 +366,7 @@ int GTSP_read_solution(struct GTSP *gtsp, char *filename, double **p_x)
for (int i = 0; i < edge_count; i++) for (int i = 0; i < edge_count; i++)
{ {
if (x[i] <= LP_EPSILON) continue; if (x[i] <= EPSILON) continue;
log_debug(" x%-5d = %.6f\n", i, x[i]); log_debug(" x%-5d = %.6f\n", i, x[i]);
} }
@@ -420,10 +418,10 @@ int GTSP_check_solution(struct GTSP *data, double *x)
for (int i = 0; i < edge_count; i++) for (int i = 0; i < edge_count; i++)
{ {
abort_iff(x[i] < 1.0 - LP_EPSILON && x[i] > LP_EPSILON, abort_iff(x[i] < 1.0 - EPSILON && x[i] > EPSILON,
"solution is not integral: x%d = %.4lf", i, x[i]); "solution is not integral: x%d = %.4lf", i, x[i]);
abort_iff(x[i] > 1.0 + LP_EPSILON || x[i] < 0.0 - LP_EPSILON, abort_iff(x[i] > 1.0 + EPSILON || x[i] < 0.0 - EPSILON,
"value out of bounds: x%d = %.4lf", i, x[i]); "value out of bounds: x%d = %.4lf", i, x[i]);
} }
@@ -435,7 +433,7 @@ int GTSP_check_solution(struct GTSP *data, double *x)
int initial; int initial;
for (initial = 0; initial < edge_count; initial++) for (initial = 0; initial < edge_count; initial++)
if (x[initial] > 1.0 - LP_EPSILON) break; if (x[initial] > 1.0 - EPSILON) break;
abort_if(initial == edge_count, "no initial node"); abort_if(initial == edge_count, "no initial node");
@@ -453,7 +451,7 @@ int GTSP_check_solution(struct GTSP *data, double *x)
struct Node *neighbor = adj->neighbor; struct Node *neighbor = adj->neighbor;
if (neighbor->mark) continue; if (neighbor->mark) continue;
if (x[adj->edge->index] < LP_EPSILON) continue; if (x[adj->edge->index] < EPSILON) continue;
stack[stack_top++] = neighbor; stack[stack_top++] = neighbor;
neighbor->mark = 1; neighbor->mark = 1;
@@ -499,7 +497,8 @@ int GTSP_solution_found(struct BNC *bnc, struct GTSP *data, double *x)
rval = large_neighborhood_search(tour, data, &tour_cost); rval = large_neighborhood_search(tour, data, &tour_cost);
abort_if(rval, "large_neighborhood_search failed"); abort_if(rval, "large_neighborhood_search failed");
if(tour_cost + LP_EPSILON < *best_val){ if (tour_cost + EPSILON < *best_val)
{
log_info("Local search improve the integral solution\n"); log_info("Local search improve the integral solution\n");
log_info(" obj val = %f\n", *best_val); log_info(" obj val = %f\n", *best_val);
log_info(" after LNS = %d\n", tour_cost); log_info(" after LNS = %d\n", tour_cost);
@@ -514,298 +513,6 @@ int GTSP_solution_found(struct BNC *bnc, struct GTSP *data, double *x)
return rval; return rval;
} }
static const struct option options_tab[] = {{"help", no_argument, 0, 'h'},
{"nodes", required_argument, 0, 'n'},
{"clusters", required_argument, 0, 'm'},
{"grid-size", required_argument, 0, 'g'},
{"optimal", required_argument, 0, 'x'},
{"seed", required_argument, 0, 's'},
{(char *) 0, (int) 0, (int *) 0, (int) 0}};
static int input_node_count = -1;
static int input_cluster_count = -1;
static int grid_size = 100;
static char input_x_filename[1000] = {0};
static void GTSP_print_usage()
{
printf("Parameters:\n");
printf("%4s %-13s %s\n", "-n", "--nodes", "number of nodes");
printf("%4s %-13s %s\n", "-m", "--clusters", "number of clusters");
printf("%4s %-13s %s\n", "-s", "--seed", "random seed");
printf("%4s %-13s %s\n", "-g", "--grid-size",
"size of the box used for generating random points");
printf("%4s %-13s %s\n", "-x", "--optimal",
"file containg valid solution (used to assert validity of cuts)");
}
static int GTSP_parse_args(int argc, char **argv)
{
int rval = 0;
opterr = 0;
while (1)
{
int c = 0;
int option_index = 0;
c = getopt_long(argc, argv, "n:m:g:x:s:", options_tab, &option_index);
if (c < 0) break;
switch (c)
{
case 'n':
input_node_count = atoi(optarg);
break;
case 'm':
input_cluster_count = atoi(optarg);
break;
case 'g':
grid_size = atoi(optarg);
break;
case 'x':
strcpy(input_x_filename, optarg);
break;
case 's':
SEED = (unsigned) atoi(optarg);
break;
case ':':
fprintf(stderr, "option '-%c' requires an argument\n", optopt);
rval = 1;
goto CLEANUP;
case '?':
default:
fprintf(stderr, "option '-%c' is invalid\n", optopt);
rval = 1;
goto CLEANUP;
}
}
if (input_cluster_count < 0)
{
input_cluster_count = (int) ceil(input_node_count / 5.0);
if (input_cluster_count < 3) input_cluster_count = 3;
}
if (input_node_count < 0)
{
printf("You must specify the number of nodes.\n");
rval = 1;
}
if (input_cluster_count > input_node_count)
{
printf("Number of clusters must be at most number of nodes.\n");
rval = 1;
}
if (rval)
{
GTSP_print_usage();
rval = 1;
}
CLEANUP:
return rval;
}
int GTSP_main(int argc, char **argv)
{
int rval = 0;
struct BNC bnc;
struct GTSP data;
double *initial_x = 0;
double initial_time = get_current_time();
SEED = (unsigned int) get_real_time() % 1000;
rval = GTSP_init_data(&data);
abort_if(rval, "GTSP_init_data failed");
rval = BNC_init(&bnc);
abort_if(rval, "BNC_init failed");
rval = GTSP_parse_args(argc, argv);
if (rval) return 1;
srand(SEED);
log_info("Generating random GTSP instance...\n");
log_info(" seed = %d\n", SEED);
log_info(" input_node_count = %d\n", input_node_count);
log_info(" input_cluster_count = %d\n", input_cluster_count);
log_info(" grid_size = %d\n", grid_size);
char instance_name[1000];
sprintf(instance_name, "gtsp-m%d-n%d-s%d", input_cluster_count,
input_node_count, SEED);
rval = GTSP_create_random_problem(input_node_count, input_cluster_count,
grid_size, &data);
abort_if(rval, "GTSP_create_random_problem failed");
char filename[100];
sprintf(filename, "input/%s.in", instance_name);
log_info("Writing random instance to file %s\n", filename);
rval = GTSP_write_problem(&data, filename);
abort_if(rval, "GTSP_write_problem failed");
#if LOG_LEVEL >= LOG_LEVEL_DEBUG
log_info("Writing random instance to file gtsp.in\n");
rval = GTSP_write_problem(&data, "gtsp.in");
#endif
int init_val = 0;
//
initial_x = (double *) malloc(
(data.graph->node_count + data.graph->edge_count) * sizeof(double));
abort_if(!initial_x, "could not allocate initial_x");
rval = inital_tour_value(&data, &init_val, initial_x);
abort_if(rval, "initial_tour_value failed");
rval = GTSP_solution_found(&bnc, &data, initial_x);
abort_if(rval, "check_sol failed");
bnc.best_x = initial_x;
bnc.best_obj_val = init_val;
bnc.problem_data = (void *) &data;
bnc.problem_init_lp = (int (*)(struct LP *, void *)) GTSP_init_lp;
bnc.problem_add_cutting_planes = (int (*)(
struct LP *, void *)) GTSP_add_cutting_planes;
bnc.problem_solution_found = (int (*)(
struct BNC *, void *, double *)) GTSP_solution_found;
double opt_val = 0.0;
if (strlen(input_x_filename) == 0)
{
sprintf(input_x_filename, "optimal/%s.out", instance_name);
FILE *file = fopen(input_x_filename, "r");
if (!file)
input_x_filename[0] = 0;
else
fclose(file);
}
if (strlen(input_x_filename) > 0)
{
rval = GTSP_read_solution(&data, input_x_filename, &OPTIMAL_X);
abort_if(rval, "GTSP_read_solution failed");
log_info("Optimal solution is available. Cuts will be checked.\n");
for (int i = 0; i < data.graph->edge_count; i++)
{
struct Edge *e = &data.graph->edges[i];
opt_val += OPTIMAL_X[i] * e->weight;
}
log_info(" opt = %.2lf\n", opt_val);
}
log_info("Initializing LP...\n");
rval = BNC_init_lp(&bnc);
abort_if(rval, "BNC_init_lp failed");
// log_info("Writing LP to file gtsp.lp...\n");
// rval = LP_write(bnc.lp, "gtsp.lp");
// abort_if(rval, "LP_write failed");
log_info("Starting branch-and-cut solver...\n");
rval = BNC_solve(&bnc);
abort_if(rval, "BNC_solve_node failed");
abort_if(!bnc.best_x, "problem has no feasible solution");
log_info("Optimal integral solution:\n");
log_info(" obj value = %.2lf **\n", bnc.best_obj_val);
if (OPTIMAL_X)
{
abort_iff(bnc.best_obj_val - LP_EPSILON > opt_val,
"Solution is not optimal: %.4lf > %.4lf", bnc.best_obj_val,
opt_val);
}
TOTAL_TIME = get_current_time() - initial_time;
log_info("Branch-and-bound nodes: %d\n", BNC_NODE_COUNT);
log_info("LP optimize calls: %d\n", LP_SOLVE_COUNT);
log_info("LP solving time: %.2lf\n", LP_SOLVE_TIME);
log_info("LP cut pool management time: %.2lf\n", CUT_POOL_TIME);
FILE *file = fopen("stats.tab", "a");
abort_if(!file, "could not open stats.tab");
struct stat st;
stat("stats.tab", &st);
if(st.st_size == 0)
{
fprintf(file, "%-20s ", "instance");
fprintf(file, "%-8s ", "time");
fprintf(file, "%-8s ", "subt-t");
fprintf(file, "%-8s ", "combs-t");
fprintf(file, "%-8s ", "pool-t");
fprintf(file, "%-8s ", "pool-m");
fprintf(file, "%-8s ", "lp-count");
fprintf(file, "%-8s ", "lp-time");
fprintf(file, "%-8s ", "lp-rows");
fprintf(file, "%-8s ", "lp-cols");
fprintf(file, "%-8s ", "init-v");
fprintf(file, "%-8s ", "opt-v");
fprintf(file, "%-8s ", "root-v");
fprintf(file, "%-8s ", "nodes");
fprintf(file, "%-8s ", "subt-cc");
fprintf(file, "%-8s ", "subt-nc");
fprintf(file, "%-8s ", "subt-nn");
fprintf(file, "%-8s ", "combs");
fprintf(file, "\n");
}
fprintf(file, "%-20s ", instance_name);
fprintf(file, "%-8.2lf ", TOTAL_TIME);
fprintf(file, "%-8.2lf ", SUBTOUR_TIME);
fprintf(file, "%-8.2lf ", COMBS_TIME);
fprintf(file, "%-8.2lf ", CUT_POOL_TIME);
fprintf(file, "%-8ld ", CUT_POOL_MAX_MEMORY/1024/1024);
fprintf(file, "%-8d ", LP_SOLVE_COUNT);
fprintf(file, "%-8.2lf ", LP_SOLVE_TIME);
fprintf(file, "%-8d ", LP_MAX_ROWS);
fprintf(file, "%-8d ", LP_MAX_COLS);
fprintf(file, "%-8d ", init_val);
fprintf(file, "%-8.0lf ", bnc.best_obj_val);
fprintf(file, "%-8.0lf ", ROOT_VALUE);
fprintf(file, "%-8d ", BNC_NODE_COUNT);
fprintf(file, "%-8d ", SUBTOUR_CLUSTER_CLUSTER_COUNT);
fprintf(file, "%-8d ", SUBTOUR_NODE_CLUSTER_COUNT);
fprintf(file, "%-8d ", SUBTOUR_NODE_NODE_COUNT);
fprintf(file, "%-8d ", COMBS_COUNT);
fprintf(file, "\n");
fclose(file);
CLEANUP:
if (OPTIMAL_X) free(OPTIMAL_X);
GTSP_free(&data);
BNC_free(&bnc);
return rval;
}
int build_x_from_tour(struct GTSP *data, struct Tour *tour, double *x) int build_x_from_tour(struct GTSP *data, struct Tour *tour, double *x)
{ {
int rval = 0; int rval = 0;
@@ -886,7 +593,8 @@ int inital_tour_value(struct GTSP *data, int *tour_cost, double *x)
{ {
cost = data->dist_matrix[i][tour[k].vertex]; cost = data->dist_matrix[i][tour[k].vertex];
if (cost < min_cost) if (cost < min_cost)
{ min_cost = cost; {
min_cost = cost;
min_dist_vertex = i; min_dist_vertex = i;
} }
} }
@@ -900,21 +608,30 @@ int inital_tour_value(struct GTSP *data, int *tour_cost, double *x)
int insertion_cost = -1, best_pose = -1, best_vertex = -1; int insertion_cost = -1, best_pose = -1, best_vertex = -1;
for (int k=0 ; k < data->clusters[cluster_to_insert].size; k++){ 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]; for (int j = 0; j < new_vertex; j++)
if (new_vertex ==1){ {
int vertex_to_insert = data->clusters[cluster_to_insert]
.nodes[k];
if (new_vertex == 1)
{
int vertex1 = tour[0].vertex; int vertex1 = tour[0].vertex;
insertion_cost = data->dist_matrix[vertex_to_insert][vertex1] insertion_cost =
+ data->dist_matrix[vertex1][vertex_to_insert]; data->dist_matrix[vertex_to_insert][vertex1] +
}else{ data->dist_matrix[vertex1][vertex_to_insert];
}
else
{
int vertex1 = tour[j].vertex; int vertex1 = tour[j].vertex;
int vertex2 = tour[tour[j].next].vertex; int vertex2 = tour[tour[j].next].vertex;
insertion_cost = data->dist_matrix[vertex1][vertex_to_insert] insertion_cost =
+ data->dist_matrix[vertex_to_insert][vertex2] - data->dist_matrix[vertex1][vertex_to_insert] +
data->dist_matrix[vertex_to_insert][vertex2] -
data->dist_matrix[vertex1][vertex2]; data->dist_matrix[vertex1][vertex2];
} }
if(insertion_cost < min_cost){ if (insertion_cost < min_cost)
{
min_cost = insertion_cost; min_cost = insertion_cost;
best_pose = j; best_pose = j;
best_vertex = vertex_to_insert; best_vertex = vertex_to_insert;
@@ -925,10 +642,13 @@ int inital_tour_value(struct GTSP *data, int *tour_cost, double *x)
tour[new_vertex].vertex = best_vertex; tour[new_vertex].vertex = best_vertex;
tour[new_vertex].prev = best_pose; tour[new_vertex].prev = best_pose;
if (new_vertex ==1){ if (new_vertex == 1)
{
tour[new_vertex].next = best_pose; tour[new_vertex].next = best_pose;
tour[best_pose].prev = new_vertex; tour[best_pose].prev = new_vertex;
}else{ }
else
{
int temp_vertex = tour[best_pose].next; int temp_vertex = tour[best_pose].next;
tour[new_vertex].next = temp_vertex; tour[new_vertex].next = temp_vertex;
tour[temp_vertex].prev = new_vertex; tour[temp_vertex].prev = new_vertex;
@@ -940,7 +660,6 @@ int inital_tour_value(struct GTSP *data, int *tour_cost, double *x)
tour[data->cluster_count].vertex = 0; tour[data->cluster_count].vertex = 0;
rval = large_neighborhood_search(tour, data, tour_cost); rval = large_neighborhood_search(tour, data, tour_cost);
abort_if(rval, "large_neighborhood_search failed"); abort_if(rval, "large_neighborhood_search failed");
@@ -1027,7 +746,10 @@ int two_opt(struct Tour *tour, struct GTSP *data)
return 0; return 0;
} }
int large_neighborhood_search(struct Tour *tour, struct GTSP *data, int *tour_cost) int large_neighborhood_search(
struct Tour *tour,
struct GTSP *data,
int *tour_cost)
{ {
int rval = 0; int rval = 0;
@@ -1043,7 +765,6 @@ int large_neighborhood_search(struct Tour *tour, struct GTSP *data, int *tour_co
int delete_vertex = rand() % (cluster_count - 1) + 1; int delete_vertex = rand() % (cluster_count - 1) + 1;
int prev_vertex = tour[delete_vertex].prev; int prev_vertex = tour[delete_vertex].prev;
int next_vertex = tour[delete_vertex].next; int next_vertex = tour[delete_vertex].next;
@@ -1164,7 +885,8 @@ int K_opt(struct Tour* tour, struct GTSP *data)
int cluster_count = data->cluster_count; int cluster_count = data->cluster_count;
int l; int l;
for(int i = 0 ; i < cluster_count - 3; i++){ for (int i = 0; i < cluster_count - 3; i++)
{
int location_in_path = 0; int location_in_path = 0;
for (int k = 0; k < i; k++) for (int k = 0; k < i; k++)
@@ -1172,12 +894,13 @@ int K_opt(struct Tour* tour, struct GTSP *data)
int current_vertex = tour[location_in_path].vertex; int current_vertex = tour[location_in_path].vertex;
for (int k = 3; k < cluster_count - i - 2; k++)
for(int k = 3; k < cluster_count - i - 2; k++){ {
int first_next_location = tour[location_in_path].next; int first_next_location = tour[location_in_path].next;
int first_next_vertex = tour[tour[location_in_path].next].vertex; int first_next_vertex = tour[tour[location_in_path].next].vertex;
int next_location = tour[location_in_path].next; int next_location = tour[location_in_path].next;
for(l = 0; l < k; l++){ for (l = 0; l < k; l++)
{
next_location = tour[next_location].next; next_location = tour[next_location].next;
} }
int next_vertex = tour[next_location].vertex; int next_vertex = tour[next_location].vertex;
@@ -1243,7 +966,7 @@ int build_tour_from_x(struct GTSP *data, struct Tour *tour, double *x)
int initial; int initial;
for (initial = 0; initial < edge_count; initial++) for (initial = 0; initial < edge_count; initial++)
if (x[initial] > 1.0 - LP_EPSILON) break; if (x[initial] > 1.0 - EPSILON) break;
initial = graph->edges[initial].from->index; initial = graph->edges[initial].from->index;
@@ -1266,7 +989,7 @@ int build_tour_from_x(struct GTSP *data, struct Tour *tour, double *x)
struct Node *neighbor = adj->neighbor; struct Node *neighbor = adj->neighbor;
if (neighbor->mark) continue; if (neighbor->mark) continue;
if (x[adj->edge->index] < LP_EPSILON) continue; if (x[adj->edge->index] < EPSILON) continue;
stack[stack_top++] = neighbor; stack[stack_top++] = neighbor;
tour[next_vertex].vertex = neighbor->index; tour[next_vertex].vertex = neighbor->index;
@@ -1276,16 +999,22 @@ int build_tour_from_x(struct GTSP *data, struct Tour *tour, double *x)
} }
} }
for (int i = 0; i < data->cluster_count; i++)
for(int i =0; i < data->cluster_count; i++){ {
if(i == 0){ if (i == 0)
{
tour[i].prev = data->cluster_count - 1; tour[i].prev = data->cluster_count - 1;
}else{ }
else
{
tour[i].prev = i - 1; tour[i].prev = i - 1;
} }
if(i == data->cluster_count - 1){ if (i == data->cluster_count - 1)
{
tour[i].next = 0; tour[i].next = 0;
}else{ }
else
{
tour[i].next = i + 1; tour[i].next = i + 1;
} }
} }

8
src/gtsp.h
View File

@@ -46,8 +46,6 @@ int GTSP_write_problem(struct GTSP *data, char *filename);
int GTSP_write_solution(struct GTSP *data, char *filename, double *x); int GTSP_write_solution(struct GTSP *data, char *filename, double *x);
int GTSP_main(int argc, char **argv);
int optimize_vertex_in_cluster(struct Tour * tour, struct GTSP *data); int optimize_vertex_in_cluster(struct Tour * tour, struct GTSP *data);
int two_opt(struct Tour * tour, struct GTSP *data); int two_opt(struct Tour * tour, struct GTSP *data);
@@ -64,4 +62,10 @@ void print_list(struct Tour * tour, struct GTSP* data);
int build_tour_from_x(struct GTSP *data, struct Tour *tour, double *x); int build_tour_from_x(struct GTSP *data, struct Tour *tour, double *x);
int GTSP_solution_found(struct BNC *bnc, struct GTSP *data, double *x);
int GTSP_check_solution(struct GTSP *data, double *x);
int GTSP_read_solution(struct GTSP *gtsp, char *filename, double **p_x);
#endif //_PROJECT_GTSP_H_ #endif //_PROJECT_GTSP_H_

397
src/lp.c
View File

@@ -7,6 +7,170 @@
#include "util.h" #include "util.h"
#include "main.h" #include "main.h"
static int compress_cut_pool(struct LP *lp)
{
int delete_count = 0;
for (int i = 0; i < lp->cut_pool_size; i++)
{
if (lp->cut_pool[i]->cplex_row_index < 0)
{
free(lp->cut_pool[i]->rmatind);
free(lp->cut_pool[i]->rmatval);
free(lp->cut_pool[i]);
delete_count++;
}
else
{
lp->cut_pool[i - delete_count] = lp->cut_pool[i];
}
}
lp->cut_pool_size -= delete_count;
return 0;
}
static int remove_old_cuts(struct LP *lp)
{
int rval = 0;
int *should_remove = 0;
int numrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp);
log_verbose(" numrows=%d\n", numrows);
should_remove = (int *) malloc((numrows + 1) * sizeof(int));
abort_if(!should_remove, "could not allocate should_remove");
for (int i = 0; i < numrows; i++)
should_remove[i] = 0;
should_remove[numrows] = 0;
log_verbose("Old cplex row index:\n");
for (int i = 0; i < lp->cut_pool_size; i++)
log_verbose(" %d\n", lp->cut_pool[i]->cplex_row_index);
log_verbose("Should remove:\n");
for (int i = 0; i < lp->cut_pool_size; i++)
{
struct Row *cut = lp->cut_pool[i];
if (cut->age <= MAX_CUT_AGE) continue;
if (cut->cplex_row_index < 0) continue;
should_remove[cut->cplex_row_index] = 1;
log_verbose(" %d\n", cut->cplex_row_index);
}
// Update cut->cplex_row_index
int count = 0;
for (int i = 0; i < numrows; i++)
{
if (!should_remove[i]) continue;
for (int j = 0; j < lp->cut_pool_size; j++)
{
struct Row *cut = lp->cut_pool[j];
if (cut->cplex_row_index == i - count) cut->cplex_row_index = -1;
else if (cut->cplex_row_index > i - count) cut->cplex_row_index--;
}
count++;
}
log_verbose("New cplex row index:\n");
for (int i = 0; i < lp->cut_pool_size; i++)
log_verbose(" %d\n", lp->cut_pool[i]->cplex_row_index);
// Remove from CPLEX, finding the largest intervals of contiguous rows
// that should be removed.
int start = 0;
int end = -1;
count = 0;
for (int i = 0; i < numrows + 1; i++)
{
if (should_remove[i])
end++;
else
{
if (end >= start)
{
rval = CPXdelrows(lp->cplex_env, lp->cplex_lp, start - count,
end - count);
abort_if(rval, "CPXdelrows failed");
log_verbose(" %d %d (%d)\n", start, end, end - start + 1);
count += end - start + 1;
}
start = i + 1;
end = i;
}
}
log_debug(" removed %d old cuts\n", count);
rval = CPXdualopt(lp->cplex_env, lp->cplex_lp);
abort_if(rval, "CPXoptimize failed");
log_debug("Compressing cut pool...\n");
compress_cut_pool(lp);
long nz = 0;
long size = 0;
for (int i = 0; i < lp->cut_pool_size; i++)
{
size += sizeof(struct Row);
nz += lp->cut_pool[i]->nz;
size += lp->cut_pool[i]->nz * sizeof(double);
size += lp->cut_pool[i]->nz * sizeof(int);
}
log_debug(" %ld cuts (%ld nz, %ld MiB)\n", lp->cut_pool_size, nz, size/1024/1024);
if(size > CUT_POOL_MAX_MEMORY)
CUT_POOL_MAX_MEMORY = size;
CLEANUP:
if (should_remove) free(should_remove);
return rval;
}
static int update_cut_ages(struct LP *lp)
{
int rval = 0;
double *slacks = 0;
int numrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp);
slacks = (double *) malloc(numrows * sizeof(double));
abort_if(!slacks, "could not allocate slacks");
rval = CPXgetslack(lp->cplex_env, lp->cplex_lp, slacks, 0, numrows - 1);
abort_if(rval, "CPXgetslack failed");
int count = 0;
for (int i = 0; i < lp->cut_pool_size; i++)
{
struct Row *cut = lp->cut_pool[i];
if (cut->cplex_row_index < 0) continue;
assert(cut->cplex_row_index < numrows);
if (slacks[cut->cplex_row_index] < -EPSILON)
cut->age++;
else
cut->age = 0;
if (cut->age > 5) count++;
}
CLEANUP:
if (slacks) free(slacks);
return rval;
}
int LP_open(struct LP *lp) int LP_open(struct LP *lp)
{ {
int rval = 0; int rval = 0;
@@ -50,32 +214,6 @@ void LP_free(struct LP *lp)
LP_free_cut_pool(lp); LP_free_cut_pool(lp);
} }
int LP_compress_cut_pool(struct LP *lp)
{
int rval = 0;
int delete_count = 0;
for (int i = 0; i < lp->cut_pool_size; i++)
{
if (lp->cut_pool[i]->cplex_row_index < 0)
{
free(lp->cut_pool[i]->rmatind);
free(lp->cut_pool[i]->rmatval);
free(lp->cut_pool[i]);
delete_count++;
}
else
{
lp->cut_pool[i - delete_count] = lp->cut_pool[i];
}
}
lp->cut_pool_size -= delete_count;
CLEANUP:
return rval;
}
int LP_create(struct LP *lp, const char *name) int LP_create(struct LP *lp, const char *name)
{ {
int rval = 0; int rval = 0;
@@ -98,7 +236,6 @@ int LP_new_row(struct LP *lp, char sense, double rhs)
int rval = 0; int rval = 0;
rval = CPXnewrows(lp->cplex_env, lp->cplex_lp, 1, &rhs, &sense, 0, 0); rval = CPXnewrows(lp->cplex_env, lp->cplex_lp, 1, &rhs, &sense, 0, 0);
abort_if(rval, "CPXnewrows failed"); abort_if(rval, "CPXnewrows failed");
CLEANUP: CLEANUP:
@@ -119,7 +256,6 @@ int LP_add_rows(
rval = CPXaddrows(lp->cplex_env, lp->cplex_lp, 0, newrows, newnz, rhs, rval = CPXaddrows(lp->cplex_env, lp->cplex_lp, 0, newrows, newnz, rhs,
sense, rmatbeg, rmatind, rmatval, 0, 0); sense, rmatbeg, rmatind, rmatval, 0, 0);
abort_if(rval, "CPXaddrows failed"); abort_if(rval, "CPXaddrows failed");
CLEANUP: CLEANUP:
@@ -159,46 +295,6 @@ int LP_change_bound(struct LP *lp, int col, char lower_or_upper, double bnd)
return rval; return rval;
} }
extern int LP_SOLVE_COUNT;
extern double LP_SOLVE_TIME;
extern double CUT_POOL_TIME;
int LP_update_cut_ages(struct LP *lp)
{
int rval = 0;
double *slacks = 0;
int numrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp);
slacks = (double *) malloc(numrows * sizeof(double));
abort_if(!slacks, "could not allocate slacks");
rval = CPXgetslack(lp->cplex_env, lp->cplex_lp, slacks, 0, numrows - 1);
abort_if(rval, "CPXgetslack failed");
int count = 0;
for (int i = 0; i < lp->cut_pool_size; i++)
{
struct Row *cut = lp->cut_pool[i];
if (cut->cplex_row_index < 0) continue;
assert(cut->cplex_row_index < numrows);
if (slacks[cut->cplex_row_index] < -LP_EPSILON)
cut->age++;
else
cut->age = 0;
if (cut->age > 5) count++;
}
CLEANUP:
if (slacks) free(slacks);
return rval;
}
int LP_optimize(struct LP *lp, int *infeasible) int LP_optimize(struct LP *lp, int *infeasible)
{ {
LP_SOLVE_COUNT++; LP_SOLVE_COUNT++;
@@ -215,12 +311,12 @@ int LP_optimize(struct LP *lp, int *infeasible)
log_debug("Optimizing LP (%d rows %d cols)...\n", numrows, numcols); log_debug("Optimizing LP (%d rows %d cols)...\n", numrows, numcols);
double initial_time = get_current_time(); double initial_time = get_user_time();
rval = CPXdualopt(lp->cplex_env, lp->cplex_lp); rval = CPXdualopt(lp->cplex_env, lp->cplex_lp);
abort_if(rval, "CPXdualopt failed"); abort_if(rval, "CPXdualopt failed");
LP_SOLVE_TIME += get_current_time() - initial_time; LP_SOLVE_TIME += get_user_time() - initial_time;
solstat = CPXgetstat(lp->cplex_env, lp->cplex_lp); solstat = CPXgetstat(lp->cplex_env, lp->cplex_lp);
if (solstat == CPX_STAT_INFEASIBLE) if (solstat == CPX_STAT_INFEASIBLE)
@@ -239,131 +335,22 @@ int LP_optimize(struct LP *lp, int *infeasible)
abort_if(rval, "LP_get_obj_val failed"); abort_if(rval, "LP_get_obj_val failed");
log_debug(" obj val = %.4lf\n", objval); log_debug(" obj val = %.4lf\n", objval);
log_debug(" time = %.4lf\n", get_current_time() - initial_time); log_debug(" time = %.4lf\n", get_user_time() - initial_time);
initial_time = get_current_time(); initial_time = get_user_time();
rval = LP_update_cut_ages(lp); rval = update_cut_ages(lp);
abort_if(rval, "LP_update_cut_ages failed"); abort_if(rval, "update_cut_ages failed");
log_debug("Removing old cuts...\n"); log_debug("Removing old cuts...\n");
rval = LP_remove_old_cuts(lp); rval = remove_old_cuts(lp);
abort_if(rval, "LP_remove_old_cuts failed"); abort_if(rval, "LP_remove_old_cuts failed");
CUT_POOL_TIME += get_current_time() - initial_time; CUT_POOL_TIME += get_user_time() - initial_time;
CLEANUP: CLEANUP:
return rval; return rval;
} }
int LP_remove_old_cuts(struct LP *lp)
{
int rval = 0;
int *should_remove = 0;
int numrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp);
log_verbose(" numrows=%d\n", numrows);
should_remove = (int *) malloc((numrows + 1) * sizeof(int));
abort_if(!should_remove, "could not allocate should_remove");
for (int i = 0; i < numrows; i++)
should_remove[i] = 0;
should_remove[numrows] = 0;
log_verbose("Old cplex row index:\n");
for (int i = 0; i < lp->cut_pool_size; i++)
log_verbose(" %d\n", lp->cut_pool[i]->cplex_row_index);
log_verbose("Should remove:\n");
for (int i = 0; i < lp->cut_pool_size; i++)
{
struct Row *cut = lp->cut_pool[i];
if (cut->age <= MAX_CUT_AGE) continue;
if (cut->cplex_row_index < 0) continue;
should_remove[cut->cplex_row_index] = 1;
log_verbose(" %d\n", cut->cplex_row_index);
}
int count = 0;
for (int i = 0; i < numrows; i++)
{
if (!should_remove[i]) continue;
for (int j = 0; j < lp->cut_pool_size; j++)
{
struct Row *cut = lp->cut_pool[j];
if (cut->cplex_row_index == i - count) cut->cplex_row_index = -1;
else if (cut->cplex_row_index > i - count) cut->cplex_row_index--;
}
count++;
}
log_verbose("New cplex row index:\n");
for (int i = 0; i < lp->cut_pool_size; i++)
log_verbose(" %d\n", lp->cut_pool[i]->cplex_row_index);
int start = 0;
int end = -1;
count = 0;
for (int i = 0; i < numrows + 1; i++)
{
if (should_remove[i])
{
end++;
}
else
{
if (end >= start)
{
rval = CPXdelrows(lp->cplex_env, lp->cplex_lp, start - count,
end - count);
abort_if(rval, "CPXdelrows failed");
log_verbose(" %d %d (%d)\n", start, end, end - start + 1);
count += end - start + 1;
}
start = i + 1;
end = i;
}
}
numrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp);
log_verbose(" new numrows=%d\n", numrows);
log_debug(" removed %d old cuts\n", count);
rval = CPXdualopt(lp->cplex_env, lp->cplex_lp);
abort_if(rval, "CPXoptimize failed");
log_debug("Compressing cut pool...\n");
LP_compress_cut_pool(lp);
long nz = 0;
long size = 0;
for (int i = 0; i < lp->cut_pool_size; i++)
{
size += sizeof(struct Row);
nz += lp->cut_pool[i]->nz;
size += lp->cut_pool[i]->nz * sizeof(double);
size += lp->cut_pool[i]->nz * sizeof(int);
}
log_debug(" %ld cuts (%ld nz, %ld MiB)\n", lp->cut_pool_size, nz, size/1024/1024);
if(size > CUT_POOL_MAX_MEMORY)
CUT_POOL_MAX_MEMORY = size;
CLEANUP:
if (should_remove) free(should_remove);
return rval;
}
int LP_get_obj_val(struct LP *lp, double *obj) int LP_get_obj_val(struct LP *lp, double *obj)
{ {
int rval = 0; int rval = 0;
@@ -414,8 +401,8 @@ int LP_write(struct LP *lp, const char *fname)
if((a)>(b)) return 1; if((a)>(b)) return 1;
#define return_if_neq_epsilon(a, b) \ #define return_if_neq_epsilon(a, b) \
if((a+LP_EPSILON)<(b)) return -1; \ if((a+EPSILON)<(b)) return -1; \
if((a-LP_EPSILON)>(b)) return 1; if((a-EPSILON)>(b)) return 1;
int compare_cuts(struct Row *cut1, struct Row *cut2) int compare_cuts(struct Row *cut1, struct Row *cut2)
{ {
@@ -430,12 +417,27 @@ int compare_cuts(struct Row *cut1, struct Row *cut2)
return 0; return 0;
} }
static int update_hash(struct Row *cut)
{
unsigned long hash = 0;
for (int i = 0; i < cut->nz; i++)
{
hash += cut->rmatind[i] * 65521;
hash %= 4294967291;
}
cut->hash = hash;
return 0;
}
int LP_add_cut(struct LP *lp, struct Row *cut) int LP_add_cut(struct LP *lp, struct Row *cut)
{ {
int rval = 0; int rval = 0;
double initial_time = get_current_time(); double initial_time = get_user_time();
rval = LP_update_hash(cut); rval = update_hash(cut);
abort_if(rval, "LP_update_hash failed"); abort_if(rval, "LP_update_hash failed");
for (int i = 0; i < lp->cut_pool_size; i++) for (int i = 0; i < lp->cut_pool_size; i++)
@@ -465,22 +467,7 @@ int LP_add_cut(struct LP *lp, struct Row *cut)
cut->cplex_row_index = CPXgetnumrows(lp->cplex_env, lp->cplex_lp) - 1; cut->cplex_row_index = CPXgetnumrows(lp->cplex_env, lp->cplex_lp) - 1;
cut->age = 0; cut->age = 0;
CUT_POOL_TIME += get_current_time() - initial_time; CUT_POOL_TIME += get_user_time() - initial_time;
CLEANUP: CLEANUP:
return rval; return rval;
} }
int LP_update_hash(struct Row *cut)
{
unsigned long hash = 0;
for (int i = 0; i < cut->nz; i++)
{
hash += cut->rmatind[i] * 65521;
hash %= 4294967291;
}
cut->hash = hash;
return 0;
}

4
src/lp.h
View File

@@ -69,10 +69,6 @@ int LP_get_x(struct LP *lp, double *x);
int LP_get_num_cols(struct LP *lp); int LP_get_num_cols(struct LP *lp);
int LP_remove_old_cuts(struct LP *lp);
int LP_update_hash(struct Row *cut);
int LP_add_cut(struct LP *lp, struct Row *cut); int LP_add_cut(struct LP *lp, struct Row *cut);
#endif #endif

319
src/main.c
View File

@@ -1,14 +1,12 @@
#include <getopt.h> #include <getopt.h>
#include <string.h>
#include <math.h>
#include <sys/stat.h>
#include "main.h" #include "main.h"
#include "gtsp.h" #include "gtsp.h"
#include "tsp.h" #include "util.h"
#include "flow.h"
char *INPUT_FILENAME = 0;
unsigned int SEED = 0; unsigned int SEED = 0;
int GEOMETRIC_DATA = 0;
int NODE_COUNT_RAND = 0;
int GRID_SIZE_RAND = 100;
double SUBTOUR_TIME = 0; double SUBTOUR_TIME = 0;
double COMBS_TIME = 0; double COMBS_TIME = 0;
@@ -25,53 +23,306 @@ int LP_SOLVE_COUNT = 0;
double TOTAL_TIME = 0; double TOTAL_TIME = 0;
double ROOT_VALUE = 0; double ROOT_VALUE = 0;
int SUBTOUR_CLUSTER_CLUSTER_COUNT = 0; int SUBTOUR_COUNT = 0;
int SUBTOUR_NODE_CLUSTER_COUNT = 0;
int SUBTOUR_NODE_NODE_COUNT = 0;
int COMBS_COUNT = 0; int COMBS_COUNT = 0;
double INITIAL_TIME = 0; double INITIAL_TIME = 0;
static const struct option options_tab[] = { int BNC_NODE_COUNT = 0;
{"help", no_argument, 0, 'h'}, {"tsp", no_argument, 0, 't'},
{"gtsp", no_argument, 0, 'g'}, {"flow", no_argument, 0, 'f'},
{(char *) 0, (int) 0, (int *) 0, (int) 0}
};
void GTSP_print_usage() static int input_node_count = -1;
static int input_cluster_count = -1;
static int grid_size = 100;
static const struct option options_tab[] = {{"help", no_argument, 0, 'h'},
{"nodes", required_argument, 0, 'n'},
{"clusters", required_argument, 0, 'm'},
{"grid-size", required_argument, 0, 'g'},
{"optimal", required_argument, 0, 'x'},
{"seed", required_argument, 0, 's'},
{(char *) 0, (int) 0, (int *) 0, (int) 0}};
static char input_x_filename[1000] = {0};
static void print_usage(int argc, char **argv)
{ {
printf("wrong usage\n"); printf("Usage: %s [OPTION]...\n", argv[0]);
printf("Solves the Generalized Traveling Salesman problem for the input graph.\n\n");
printf("Parameters:\n");
printf("%4s %-13s %s\n", "-n", "--nodes", "number of nodes");
printf("%4s %-13s %s\n", "-m", "--clusters", "number of clusters");
printf("%4s %-13s %s\n", "-s", "--seed", "random seed");
printf("%4s %-13s %s\n", "-g", "--grid-size",
"size of the box used for generating random points");
printf("%4s %-13s %s\n", "-x", "--optimal",
"file containg valid solution (used to assert validity of cuts)");
}
static int parse_args(int argc, char **argv)
{
int rval = 0;
opterr = 0;
while (1)
{
int c = 0;
int option_index = 0;
c = getopt_long(argc, argv, "n:m:g:x:s:h:", options_tab, &option_index);
if (c < 0) break;
switch (c)
{
case 'n':
input_node_count = atoi(optarg);
break;
case 'm':
input_cluster_count = atoi(optarg);
break;
case 'g':
grid_size = atoi(optarg);
break;
case 'x':
strcpy(input_x_filename, optarg);
break;
case 's':
SEED = (unsigned) atoi(optarg);
break;
case 'h':
print_usage(argc, argv);
exit(0);
case ':':
fprintf(stderr, "%s: option '-%c' requires an argument\n",
argv[0], optopt);
rval = 1;
goto CLEANUP;
case '?':
default:
fprintf(stderr, "%s: option '-%c' is invalid\n", argv[0],
optopt);
rval = 1;
goto CLEANUP;
}
}
if (input_cluster_count < 0)
{
input_cluster_count = (int) ceil(input_node_count / 5.0);
if (input_cluster_count < 3) input_cluster_count = 3;
}
if (input_node_count < 0)
{
fprintf(stderr, "You must specify the number of nodes.\n");
rval = 1;
}
else if (input_cluster_count > input_node_count)
{
printf("Number of clusters must be at most number of nodes.\n");
rval = 1;
}
CLEANUP:
if (rval)
fprintf(stderr, "Try '%s --help' for more information.\n", argv[0]);
return rval;
} }
int main(int argc, char **argv) int main(int argc, char **argv)
{ {
int c = 0; int rval = 0;
int option_index = 0;
c = getopt_long(argc, argv, "htgf", options_tab, &option_index); struct BNC bnc;
struct GTSP data;
if (c < 0) double *initial_x = 0;
double initial_time = get_user_time();
SEED = (unsigned int) get_real_time() % 1000;
rval = GTSP_init_data(&data);
abort_if(rval, "GTSP_init_data failed");
rval = BNC_init(&bnc);
abort_if(rval, "BNC_init failed");
rval = parse_args(argc, argv);
if (rval) return 1;
srand(SEED);
log_info("Generating random GTSP instance...\n");
log_info(" seed = %d\n", SEED);
log_info(" input_node_count = %d\n", input_node_count);
log_info(" input_cluster_count = %d\n", input_cluster_count);
log_info(" grid_size = %d\n", grid_size);
char instance_name[1000];
sprintf(instance_name, "gtsp-m%d-n%d-s%d", input_cluster_count,
input_node_count, SEED);
rval = GTSP_create_random_problem(input_node_count, input_cluster_count,
grid_size, &data);
abort_if(rval, "GTSP_create_random_problem failed");
char filename[100];
sprintf(filename, "input/%s.in", instance_name);
log_info("Writing random instance to file %s\n", filename);
rval = GTSP_write_problem(&data, filename);
abort_if(rval, "GTSP_write_problem failed");
#if LOG_LEVEL >= LOG_LEVEL_DEBUG
log_info("Writing random instance to file gtsp.in\n");
rval = GTSP_write_problem(&data, "gtsp.in");
#endif
int init_val = 0;
initial_x = (double *) malloc(
(data.graph->node_count + data.graph->edge_count) * sizeof(double));
abort_if(!initial_x, "could not allocate initial_x");
rval = inital_tour_value(&data, &init_val, initial_x);
abort_if(rval, "initial_tour_value failed");
rval = GTSP_solution_found(&bnc, &data, initial_x);
abort_if(rval, "check_sol failed");
bnc.best_x = initial_x;
bnc.best_obj_val = init_val;
bnc.problem_data = (void *) &data;
bnc.problem_init_lp = (int (*)(struct LP *, void *)) GTSP_init_lp;
bnc.problem_add_cutting_planes = (int (*)(
struct LP *, void *)) GTSP_add_cutting_planes;
bnc.problem_solution_found = (int (*)(
struct BNC *, void *, double *)) GTSP_solution_found;
double opt_val = 0.0;
if (strlen(input_x_filename) == 0)
{ {
GTSP_print_usage(); sprintf(input_x_filename, "optimal/%s.out", instance_name);
return 1;
FILE *file = fopen(input_x_filename, "r");
if (!file)
input_x_filename[0] = 0;
else
fclose(file);
} }
switch (c) if (strlen(input_x_filename) > 0)
{ {
case 'f': rval = GTSP_read_solution(&data, input_x_filename, &OPTIMAL_X);
return flow_main(argc, argv); abort_if(rval, "GTSP_read_solution failed");
case 't': log_info("Optimal solution is available. Cuts will be checked.\n");
return TSP_main(argc, argv);
case 'g': for (int i = 0; i < data.graph->edge_count; i++)
return GTSP_main(argc, argv); {
struct Edge *e = &data.graph->edges[i];
default: opt_val += OPTIMAL_X[i] * e->weight;
case 'h':
GTSP_print_usage();
return 1;
} }
log_info(" opt = %.2lf\n", opt_val);
} }
log_info("Initializing LP...\n");
rval = BNC_init_lp(&bnc);
abort_if(rval, "BNC_init_lp failed");
// log_info("Writing LP to file gtsp.lp...\n");
// rval = LP_write(bnc.lp, "gtsp.lp");
// abort_if(rval, "LP_write failed");
log_info("Starting branch-and-cut solver...\n");
rval = BNC_solve(&bnc);
abort_if(rval, "BNC_solve_node failed");
abort_if(!bnc.best_x, "problem has no feasible solution");
log_info("Optimal integral solution:\n");
log_info(" obj value = %.2lf **\n", bnc.best_obj_val);
if (OPTIMAL_X)
{
abort_iff(bnc.best_obj_val - EPSILON > opt_val,
"Solution is not optimal: %.4lf > %.4lf", bnc.best_obj_val,
opt_val);
}
TOTAL_TIME = get_user_time() - initial_time;
log_info("Branch-and-bound nodes: %d\n", BNC_NODE_COUNT);
log_info("LP optimize calls: %d\n", LP_SOLVE_COUNT);
log_info("LP solving time: %.2lf\n", LP_SOLVE_TIME);
log_info("LP cut pool management time: %.2lf\n", CUT_POOL_TIME);
// Write statistics to a file
FILE *file = fopen("stats.tab", "a");
abort_if(!file, "could not open stats.tab");
struct stat st;
stat("stats.tab", &st);
if (st.st_size == 0)
{
fprintf(file, "%-20s ", "instance");
fprintf(file, "%-8s ", "time");
fprintf(file, "%-8s ", "subt-t");
fprintf(file, "%-8s ", "combs-t");
fprintf(file, "%-8s ", "pool-t");
fprintf(file, "%-8s ", "pool-m");
fprintf(file, "%-8s ", "lp-count");
fprintf(file, "%-8s ", "lp-time");
fprintf(file, "%-8s ", "lp-rows");
fprintf(file, "%-8s ", "lp-cols");
fprintf(file, "%-8s ", "init-v");
fprintf(file, "%-8s ", "opt-v");
fprintf(file, "%-8s ", "root-v");
fprintf(file, "%-8s ", "nodes");
fprintf(file, "%-8s ", "subt-cc");
fprintf(file, "%-8s ", "subt-nc");
fprintf(file, "%-8s ", "subt-nn");
fprintf(file, "%-8s ", "combs");
fprintf(file, "\n");
}
fprintf(file, "%-20s ", instance_name);
fprintf(file, "%-8.2lf ", TOTAL_TIME);
fprintf(file, "%-8.2lf ", SUBTOUR_TIME);
fprintf(file, "%-8.2lf ", COMBS_TIME);
fprintf(file, "%-8.2lf ", CUT_POOL_TIME);
fprintf(file, "%-8ld ", CUT_POOL_MAX_MEMORY / 1024 / 1024);
fprintf(file, "%-8d ", LP_SOLVE_COUNT);
fprintf(file, "%-8.2lf ", LP_SOLVE_TIME);
fprintf(file, "%-8d ", LP_MAX_ROWS);
fprintf(file, "%-8d ", LP_MAX_COLS);
fprintf(file, "%-8d ", init_val);
fprintf(file, "%-8.0lf ", bnc.best_obj_val);
fprintf(file, "%-8.0lf ", ROOT_VALUE);
fprintf(file, "%-8d ", BNC_NODE_COUNT);
fprintf(file, "%-8d ", SUBTOUR_COUNT);
fprintf(file, "%-8d ", COMBS_COUNT);
fprintf(file, "\n");
fclose(file);
CLEANUP:
if (OPTIMAL_X) free(OPTIMAL_X);
GTSP_free(&data);
BNC_free(&bnc);
return rval;
}

9
src/main.h
View File

@@ -1,11 +1,7 @@
#ifndef _PROJECT_MAIN_H_ #ifndef _PROJECT_MAIN_H_
#define _PROJECT_MAIN_H_ #define _PROJECT_MAIN_H_
extern char *INPUT_FILENAME;
extern unsigned int SEED; extern unsigned int SEED;
extern int GEOMETRIC_DATA;
extern int NODE_COUNT_RAND;
extern int GRID_SIZE_RAND;
extern double *OPTIMAL_X; extern double *OPTIMAL_X;
extern double SUBTOUR_TIME; extern double SUBTOUR_TIME;
@@ -25,10 +21,7 @@ extern double TOTAL_TIME;
extern double INITIAL_TIME; extern double INITIAL_TIME;
extern double ROOT_VALUE; extern double ROOT_VALUE;
extern int SUBTOUR_CLUSTER_CLUSTER_COUNT; extern int SUBTOUR_COUNT;
extern int SUBTOUR_NODE_CLUSTER_COUNT;
extern int SUBTOUR_NODE_NODE_COUNT;
extern int COMBS_COUNT; extern int COMBS_COUNT;
#endif #endif

26
src/params.h
View File

@@ -1,17 +1,31 @@
#ifndef PROJECT_PARAMS_H #ifndef PROJECT_PARAMS_H
#define PROJECT_PARAMS_H #define PROJECT_PARAMS_H
#define LP_EPSILON 0.000001 /*
* Error margin for floating point comparisons.
*/
#define EPSILON 0.000001
/*
* Cut pool settings.
*/
#define MAX_CUT_AGE 5 #define MAX_CUT_AGE 5
#define MAX_CUT_POOL_SIZE 1000000 #define MAX_CUT_POOL_SIZE 1000000
/*
* Available log levels, in decreasing level of verboseness, are:
*
* LOG_LEVEL_VERBOSE
* LOG_LEVEL_DEBUG
* LOG_LEVEL_INFO
* LOG_LEVEL_WARNING
* LOG_LEVEL_ERROR
*/
#define LOG_LEVEL LOG_LEVEL_INFO #define LOG_LEVEL LOG_LEVEL_INFO
#define ENABLE_COMB_INEQUALITIES /*
* Time limit for the computation (user time, in seconds).
#define MAX_TOTAL_TIME 999999 */
#define MAX_TOTAL_TIME 3600
//#define ALLOW_FRACTIONAL_SOLUTIONS
#endif //PROJECT_PARAMS_H #endif //PROJECT_PARAMS_H

577
src/tsp.c
View File

@@ -1,577 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
#include "tsp.h"
#include "util.h"
#include "geometry.h"
#include "branch_and_cut.h"
static int TSP_parse_arguments(int argc, char **argv);
static void TSP_print_usage(char *f);
int TSP_init_data(struct TSPData *data)
{
data->node_count = 0;
data->edge_count = 0;
data->edge_weights = 0;
data->edge_list = 0;
return 0;
}
void TSP_free_data(struct TSPData *data)
{
if (!data) return;
if (data->edge_list) free(data->edge_list);
if (data->edge_weights) free(data->edge_weights);
}
int TSP_init_lp(struct LP *lp, struct TSPData *data)
{
int node_count = data->node_count;
int edge_count = data->edge_count;
int *edge_list = data->edge_list;
int *edge_weights = data->edge_weights;
int rval = 0;
/* Build a row for each degree equation */
for (int i = 0; i < node_count; i++)
{
rval = LP_new_row(lp, 'E', 2.0);
abort_if(rval, "LP_new_row failed");
}
/* Build a column for each edge_index of the graph */
double lb = 0.0;
double ub = 1.0;
int cmatbeg = 0;
double cmatval[] = {1.0, 1.0};
for (int j = 0; j < edge_count; j++)
{
double obj = (double) edge_weights[j];
int cmatind[] = {edge_list[2 * j], edge_list[2 * j + 1]};
rval = LP_add_cols(lp, 1, 2, &obj, &cmatbeg, cmatind, cmatval, &lb,
&ub);
abort_if(rval, "LP_add_cols failed");
}
CLEANUP:
return rval;
}
int TSP_find_violated_subtour_elimination_cut(
struct LP *lp, struct TSPData *data)
{
int ncount = data->node_count;
int edge_count = data->edge_count;
int *edges = data->edge_list;
int rval = 0;
int is_infeasible = 0;
double *x = 0;
int *delta = 0;
int *marks = 0;
int *island_nodes = 0;
int *island_start = 0;
int *island_sizes = 0;
struct Graph G;
graph_init(&G);
rval = LP_optimize(lp, &is_infeasible);
abort_if(rval, "LP_optimize failed");
abort_if(is_infeasible, "LP is infeasible");
rval = graph_build(ncount, edge_count, edges, 0, &G);
abort_if(rval, "graph_build failed");
x = (double *) malloc(edge_count * sizeof(double));
delta = (int *) malloc(edge_count * sizeof(int));
marks = (int *) malloc(ncount * sizeof(int));
abort_if(!x, "Could not allocate memory for x");
abort_if(!delta, "Could not allocate memory for delta");
abort_if(!marks, "Could not allocate memory for marks");
island_nodes = (int *) malloc(ncount * sizeof(int));
island_start = (int *) malloc(ncount * sizeof(int));
island_sizes = (int *) malloc(ncount * sizeof(int));
abort_if(!island_nodes, "Could not allocate memory for island_nodes");
abort_if(!island_start, "Could not allocate memory for island_start");
abort_if(!island_sizes, "Could not allocate memory for island_sizes");
for (int i = 0; i < ncount; i++)
marks[i] = 0;
rval = LP_get_x(lp, x);
abort_if(rval, "LP_get_x failed");
#if LOG_LEVEL >= LOG_LEVEL_VERBOSE
int current_round = 0;
#endif
int delta_count = 0;
int island_count = 0;
while (1)
{
rval = graph_find_connected_components(&G, x, &island_count,
island_sizes, island_start, island_nodes);
abort_if(rval, "graph_find_connected_components failed");
if(island_count == 1) break;
log_verbose("Adding %d subtour inequalities...\n", island_count);
for (int i = 0; i < island_count; i++)
{
get_delta(island_sizes[i], island_nodes + island_start[i],
edge_count, edges, &delta_count, delta, marks);
rval = TSP_add_subtour_elimination_cut(lp, delta_count, delta);
}
log_verbose("Reoptimizing (round %d)...\n", ++current_round);
abort_if(rval, "TSP_add_subtour_elimination_cut failed");
rval = LP_optimize(lp, &is_infeasible);
abort_if(rval, "LP_optimize failed");
if (is_infeasible) goto CLEANUP;
double objval = 0;
rval = LP_get_obj_val(lp, &objval);
abort_if(rval, "LP_get_obj_val failed");
rval = LP_get_x(lp, x);
abort_if(rval, "LP_get_x failed");
}
log_verbose(" graph is connected\n");
CLEANUP:
graph_free(&G);
if (x) free(x);
if (island_nodes) free(island_nodes);
if (delta) free(delta);
if (marks) free(marks);
return rval;
}
int TSP_add_subtour_elimination_cut(struct LP *lp, int delta_length, int *delta)
{
int rval = 0;
char sense = 'G';
double rhs = 2.0;
int rmatbeg = 0;
double *rmatval;
int *rmatind = delta;
rmatval = (double *) malloc(delta_length * sizeof(double));
abort_if(!rmatval, "out of memory for rmatval");
for (int i = 0; i < delta_length; i++)
rmatval[i] = 1.0;
rval = LP_add_rows(lp, 1, delta_length, &rhs, &sense, &rmatbeg, rmatind,
rmatval);
abort_if(rval, "LP_add_rows failed");
CLEANUP:
if (rmatval) free(rmatval);
return rval;
}
int graph_find_connected_components(
struct Graph *G,
double *x,
int *island_count,
int *island_sizes,
int *island_start,
int *island_nodes)
{
for (int i = 0; i < G->node_count; i++)
{
G->nodes[i].mark = 0;
island_nodes[i] = -1;
}
int k = 0, current_island = 0;
for (int i = 0; i < G->node_count; i++)
{
if (G->nodes[i].mark != 0) continue;
island_sizes[current_island] = 0;
graph_dfs(i, G, x, island_sizes + current_island, island_nodes + k);
island_start[current_island] = k;
k += island_sizes[current_island];
current_island++;
}
(*island_count) = current_island;
return 0;
}
int TSP_find_closest_neighbor_tour(
int start,
int node_count,
int edge_count,
int *edges,
int *elen,
int *path_length)
{
int rval;
int current_node = start;
struct Graph G;
graph_init(&G);
rval = graph_build(node_count, edge_count, edges, 0, &G);
abort_if(rval, "graph_build failed");
for (int j = 0; j < node_count; j++)
G.nodes[j].mark = 0;
for (int j = 0; j < node_count; j++)
{
if (j == node_count - 1)
G.nodes[start].mark = 0;
struct Node *pn = &G.nodes[current_node];
pn->mark = 1;
int closest_neighbor = -1;
int closest_edge_length = 10000000;
for (int i = 0; i < pn->degree; i++)
{
int edge = pn->adj[i].edge_index;
int neighbor = pn->adj[i].neighbor_index;
if (G.nodes[neighbor].mark == 1) continue;
if (elen[edge] > closest_edge_length) continue;
closest_neighbor = neighbor;
closest_edge_length = elen[edge];
}
*path_length += closest_edge_length;
current_node = closest_neighbor;
}
CLEANUP:
graph_free(&G);
return rval;
}
int TSP_read_problem(char *filename, struct TSPData *data)
{
int *p_node_count = &data->node_count;
int *p_edge_count = &data->edge_count;
int **p_edge_list = &data->edge_list;
int **p_edge_weights = &data->edge_weights;
struct _IO_FILE *f = (struct _IO_FILE *) NULL;
int i, j, end1, end2, w, rval = 0, node_count, edge_count;
int *edge_list = 0, *edge_weights = 0;
double *x = (double *) NULL, *y = (double *) NULL;
if (filename)
{
if ((f = fopen(filename, "r")) == NULL)
{
fprintf(stderr, "Unable to open %s for input\n", filename);
rval = 1;
goto CLEANUP;
}
}
if (filename && GEOMETRIC_DATA == 0)
{
if (fscanf(f, "%d %d", &node_count, &edge_count) != 2)
{
fprintf(stderr, "Input file %s has invalid format\n", filename);
rval = 1;
goto CLEANUP;
}
printf("Nodes: %d Edges: %d\n", node_count, edge_count);
fflush(stdout);
edge_list = (int *) malloc(2 * edge_count * sizeof(int));
if (!edge_list)
{
fprintf(stderr, "out of memory for edge_list\n");
rval = 1;
goto CLEANUP;
}
edge_weights = (int *) malloc(edge_count * sizeof(int));
if (!edge_weights)
{
fprintf(stderr, "out of memory for edge_weights\n");
rval = 1;
goto CLEANUP;
}
for (i = 0; i < edge_count; i++)
{
if (fscanf(f, "%d %d %d", &end1, &end2, &w) != 3)
{
fprintf(stderr, "%s has invalid input format\n", filename);
rval = 1;
goto CLEANUP;
}
edge_list[2 * i] = end1;
edge_list[2 * i + 1] = end2;
edge_weights[i] = w;
}
}
else
{
if (filename)
{
if (fscanf(f, "%d", &node_count) != 1)
{
fprintf(stderr, "Input file %s has invalid format\n", filename);
rval = 1;
goto CLEANUP;
}
}
else
{
node_count = NODE_COUNT_RAND;
}
x = (double *) malloc(node_count * sizeof(double));
y = (double *) malloc(node_count * sizeof(double));
if (!x || !y)
{
fprintf(stdout, "out of memory for x or y\n");
rval = 1;
goto CLEANUP;
}
if (filename)
{
for (i = 0; i < node_count; i++)
{
if (fscanf(f, "%lf %lf", &x[i], &y[i]) != 2)
{
fprintf(stderr, "%s has invalid input format\n", filename);
rval = 1;
goto CLEANUP;
}
}
}
else
{
rval = generate_random_points_2d(node_count, GRID_SIZE_RAND, x, y);
if (rval)
{
fprintf(stderr, "generate_random_points_2d failed\n");
goto CLEANUP;
}
printf("%d\n", node_count);
for (i = 0; i < node_count; i++)
{
printf("%.0f %.0f\n", x[i], y[i]);
}
printf("\n");
}
edge_count = (node_count * (node_count - 1)) / 2;
log_verbose("Complete graph: %d nodes, %d edges\n", node_count,
edge_count);
edge_list = (int *) malloc(2 * edge_count * sizeof(int));
if (!edge_list)
{
fprintf(stderr, "out of memory for edge_list\n");
rval = 1;
goto CLEANUP;
}
edge_weights = (int *) malloc(edge_count * sizeof(int));
if (!edge_weights)
{
fprintf(stderr, "out of memory for edge_weights\n");
rval = 1;
goto CLEANUP;
}
edge_count = 0;
for (i = 0; i < node_count; i++)
{
for (j = i + 1; j < node_count; j++)
{
edge_list[2 * edge_count] = i;
edge_list[2 * edge_count + 1] = j;
edge_weights[edge_count] = get_euclidean_distance(x, y, i, j);
edge_count++;
}
}
}
*p_node_count = node_count;
*p_edge_count = edge_count;
*p_edge_list = edge_list;
*p_edge_weights = edge_weights;
CLEANUP:
if (f) fclose(f);
if (x) free(x);
if (y) free(y);
return rval;
}
int TSP_add_cutting_planes(struct LP *lp, struct TSPData *data)
{
int rval = 0;
rval = TSP_find_violated_subtour_elimination_cut(lp, data);
abort_if (rval, "TSP_find_violated_subtour_elimination_cut failed\n");
CLEANUP:
return rval;
}
double TSP_find_initial_solution(struct TSPData *data)
{
double best_val = 1e99;
log_verbose("Finding closest neighbor_index tour\n");
for (int i = 0; i < data->node_count; i++)
{
int path_length = 0;
TSP_find_closest_neighbor_tour(i, data->node_count, data->edge_count,
data->edge_list, data->edge_weights, &path_length);
if (best_val > path_length) best_val = path_length;
}
log_verbose(" length = %lf\n", best_val);
return best_val;
}
int TSP_main(int argc, char **argv)
{
int rval = 0;
SEED = (unsigned int) get_real_time();
struct BNC bnc;
struct TSPData data;
rval = TSP_init_data(&data);
abort_if(rval, "TSP_init_data failed");
rval = BNC_init(&bnc);
abort_if(rval, "BNC_init failed");
rval = TSP_parse_arguments(argc, argv);
abort_if(rval, "Failed to parse arguments.");
printf("Seed = %d\n", SEED);
srand(SEED);
rval = TSP_read_problem(INPUT_FILENAME, &data);
abort_if(rval, "TSP_read_problem failed");
bnc.best_obj_val = TSP_find_initial_solution(&data);
bnc.problem_data = (void *) &data;
bnc.problem_init_lp = (int (*)(struct LP *, void *)) TSP_init_lp;
bnc.problem_add_cutting_planes = (int (*)(
struct LP *,
void *)) TSP_add_cutting_planes;
rval = BNC_init_lp(&bnc);
abort_if(rval, "BNC_init_lp failed");
rval = BNC_solve(&bnc);
abort_if(rval, "BNC_solve_node failed");
log_info("Optimal integral solution:\n");
log_info(" obj value = %.2lf **\n", bnc.best_obj_val);
CLEANUP:
BNC_free(&bnc);
TSP_free_data(&data);
return rval;
}
static int TSP_parse_arguments(int argc, char **argv)
{
int rval = 0;
int c;
if (argc == 1)
{
TSP_print_usage(argv[0]);
return 1;
}
while ((c = getopt(argc, argv, "ab:gk:s:")) != EOF)
{
switch (c)
{
case 'a':;
break;
case 'b':
GRID_SIZE_RAND = atoi(optarg);
break;
case 'g':
GEOMETRIC_DATA = 1;
break;
case 'k':
NODE_COUNT_RAND = atoi(optarg);
break;
case 's':
SEED = (unsigned) atoi(optarg);
break;
case '?':
default:
TSP_print_usage(argv[0]);
return 1;
}
}
if (optind < argc) INPUT_FILENAME = argv[optind++];
if (optind != argc)
{
TSP_print_usage(argv[0]);
return 1;
}
abort_if(!INPUT_FILENAME && !NODE_COUNT_RAND,
"Must specify an input file or use -k for random problem\n");
CLEANUP:
return rval;
}
static void TSP_print_usage(char *f)
{
fprintf(stderr, "Usage: %s [-see below-] [prob_file]\n"
" -a add all subtours cuts at once\n"
" -b d gridsize d for random problems\n"
" -g prob_file has x-y coordinates\n"
" -k d generate problem with d cities\n"
" -s d random SEED\n", f);
}

50
src/tsp.h
View File

@@ -1,50 +0,0 @@
#ifndef _PROJECT_TSP_H_
#define _PROJECT_TSP_H_
#include "lp.h"
#include "graph.h"
struct TSPData
{
int node_count;
int edge_count;
int *edge_list;
int *edge_weights;
};
int TSP_init_data(struct TSPData *data);
void TSP_free_data(struct TSPData *data);
int TSP_find_violated_subtour_elimination_cut
(struct LP *lp, struct TSPData *data);
int graph_find_connected_components(
struct Graph *G,
double *x,
int *island_count,
int *island_sizes,
int *island_start,
int *island_nodes);
int TSP_find_closest_neighbor_tour(
int start,
int node_count,
int edge_count,
int *edges,
int *elen,
int *path_length);
int TSP_add_subtour_elimination_cut(struct LP *lp, int delta_length, int *delta);
int TSP_read_problem(char *filename, struct TSPData *data);
int TSP_add_cutting_planes(struct LP *lp, struct TSPData *data);
int TSP_init_lp(struct LP *lp, struct TSPData *data);
double TSP_find_initial_solution(struct TSPData *data);
int TSP_main(int argc, char **argv);
#endif

6
src/util.c
View File

@@ -5,7 +5,7 @@
#include "util.h" #include "util.h"
#include "main.h" #include "main.h"
double get_current_time() double get_user_time()
{ {
struct rusage ru; struct rusage ru;
@@ -23,9 +23,9 @@ double get_real_time()
void time_printf(const char *fmt, ...) void time_printf(const char *fmt, ...)
{ {
if (INITIAL_TIME == 0) if (INITIAL_TIME == 0)
INITIAL_TIME = get_current_time(); INITIAL_TIME = get_user_time();
printf("[%10.2lf] ", get_current_time() - INITIAL_TIME); printf("[%10.2lf] ", get_user_time() - INITIAL_TIME);
va_list args; va_list args;
va_start(args, fmt); va_start(args, fmt);

9
src/util.h
View File

@@ -48,18 +48,11 @@
fprintf(stderr, "%28s:%d " msg "\n", __FILE__, __LINE__, __VA_ARGS__); \ fprintf(stderr, "%28s:%d " msg "\n", __FILE__, __LINE__, __VA_ARGS__); \
rval = 1; goto CLEANUP; } rval = 1; goto CLEANUP; }
#define swap(x, y) do \
{ unsigned char swap_temp[sizeof(x)]; \
memcpy(swap_temp,&y,sizeof(x)); \
memcpy(&y,&x, sizeof(x)); \
memcpy(&x,swap_temp,sizeof(x)); \
} while(0)
#define UNUSED(x) (void)(x) #define UNUSED(x) (void)(x)
void time_printf(const char *fmt, ...); void time_printf(const char *fmt, ...);
double get_current_time(void); double get_user_time(void);
double get_real_time(); double get_real_time();