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

Fix subtour cuts, create separate main functions

Browse Source
This commit is contained in:
Alinson S. Xavier
2015-03-24 10:00:45 -04:00
parent 75a016c8c2
commit 5fe42908b4
13 changed files with 886 additions and 310 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
scripts/draw_solution.sage
View File

@@ -56,9 +56,10 @@ for i in range(cluster_count):
for i in range(node_count):
(x,y,cluster) = [ int(float(x)) for x in raw_input().split(' ') ]
points[cluster].append((x,y))
all_points.append((x,y))
all_points.append(vector([x,y]))
# solutions file
(node_count, edge_count) = [ int(x) for x in raw_input().split(' ') ]
edges_count = int(raw_input())
edges = []
for i in range(edges_count):
@@ -67,6 +68,11 @@ for i in range(edges_count):
plot = list_plot([], xmax=100, xmin=0, ymax=100, ymin=0)
#plot = plot + sum([text(str(i), points[i]) for i in range(len(points))])
for i in range(node_count):
plot = plot + text(str(i), all_points[i] + vector([0,-2]), color='gray')
for i in range(cluster_count):
plot = plot + list_plot(points[i], color='gray', figsize=FIGURE_SIZE,
pointsize=POINT_SIZE)

21
src/branch_and_cut.c
View File

@@ -45,7 +45,6 @@ void BNC_free(struct BNC *bnc)
int BNC_init_lp(struct BNC *bnc)
{
int rval = 0;
log_verbose("Initializing LP...\n");
rval = LP_open(bnc->lp);
abort_if(rval, "LP_open failed");
@@ -76,7 +75,7 @@ static int BNC_solve_node(struct BNC *bnc, int depth)
int rval = 0;
double *x = (double *) NULL;
log_verbose("Optimizing...\n");
log_debug("Optimizing...\n");
int is_infeasible;
rval = LP_optimize(lp, &is_infeasible);
@@ -84,7 +83,7 @@ static int BNC_solve_node(struct BNC *bnc, int depth)
if (is_infeasible)
{
log_verbose("Branch pruned by infeasibility.\n");
log_debug("Branch pruned by infeasibility.\n");
goto CLEANUP;
}
@@ -92,11 +91,11 @@ static int BNC_solve_node(struct BNC *bnc, int depth)
rval = LP_get_obj_val(lp, &objval);
abort_if(rval, "LP_get_obj_val failed\n");
log_verbose(" objective value = %.2f\n", objval);
log_debug(" objective value = %.2f\n", objval);
if (objval > *best_val)
{
log_verbose("Branch pruned by bound (%.2lf > %.2lf).\n", objval,
log_debug("Branch pruned by bound (%.2lf > %.2lf).\n", objval,
*best_val);
rval = 0;
goto CLEANUP;
@@ -127,7 +126,7 @@ static int BNC_solve_node(struct BNC *bnc, int depth)
if (BNC_is_integral(x, num_cols))
{
log_verbose(" solution is integral\n");
log_debug(" solution is integral\n");
if (objval < *best_val)
{
@@ -141,7 +140,7 @@ static int BNC_solve_node(struct BNC *bnc, int depth)
}
else
{
log_verbose(" solution is fractional\n");
log_debug(" solution is fractional\n");
rval = BNC_branch_node(bnc, x, depth);
abort_if(rval, "BNC_branch_node failed");
}
@@ -160,10 +159,10 @@ static int BNC_branch_node(struct BNC *bnc, double *x, int depth)
int num_cols = LP_get_num_cols(lp);
int best_branch_var = BNC_find_best_branching_var(x, num_cols);
log_verbose("Branching on variable x%d = %.6lf (depth %d)...\n",
log_debug("Branching on variable x%d = %.6lf (depth %d)...\n",
best_branch_var, x[best_branch_var], depth);
log_verbose("Fixing variable x%d to one...\n", best_branch_var);
log_debug("Fixing variable x%d to one...\n", best_branch_var);
rval = LP_change_bound(lp, best_branch_var, 'L', 1.0);
abort_if(rval, "LP_change_bound failed");
@@ -173,7 +172,7 @@ static int BNC_branch_node(struct BNC *bnc, double *x, int depth)
rval = LP_change_bound(lp, best_branch_var, 'L', 0.0);
abort_if(rval, "LP_change_bound failed");
log_verbose("Fixing variable x%d to zero...\n", best_branch_var);
log_debug("Fixing variable x%d to zero...\n", best_branch_var);
rval = LP_change_bound(lp, best_branch_var, 'U', 0.0);
abort_if(rval, "LP_change_bound failed");
@@ -183,7 +182,7 @@ static int BNC_branch_node(struct BNC *bnc, double *x, int depth)
rval = LP_change_bound(lp, best_branch_var, 'U', 1.0);
abort_if(rval, "LP_change_bound failed");
log_verbose("Finished branching on variable %d\n", best_branch_var);
log_debug("Finished branching on variable %d\n", best_branch_var);
CLEANUP:
return rval;

2
src/branch_and_cut.h
View File

@@ -10,6 +10,8 @@ struct BNC
double *best_x;
double best_obj_val;
double *optimal_x;
int *problem_data;
int (*problem_init_lp)(struct LP *, void *);

194
src/flow.c
View File

@@ -1,8 +1,57 @@
#include <malloc.h>
#include <float.h>
#include "flow.h"
#include "gtsp.h"
#include "util.h"
int flow_mark_reachable_nodes(
const struct Graph *graph, double *residual_caps, struct Node *from)
{
int rval = 0;
struct Node **stack;
int stack_top = 0;
int *parents;
stack = (struct Node **) malloc(graph->node_count * sizeof(struct Node *));
abort_if(!stack, "could not allocate stack");
parents = (int *) malloc(graph->node_count * sizeof(int ));
abort_if(!parents, "could not allocate parents");
stack[stack_top++] = from;
from->mark = 1;
while (stack_top > 0)
{
struct Node *n = stack[--stack_top];
for (int j = 0; j < n->degree; j++)
{
struct Edge *e = n->adj[j].edge;
struct Node *neighbor = n->adj[j].neighbor;
if (neighbor->mark) continue;
if (residual_caps[e->index] <= 0) continue;
stack[stack_top++] = neighbor;
neighbor->mark = 1;
parents[neighbor->index] = n->index;
}
}
log_verbose("Reachable nodes:\n");
for (int i = 0; i < graph->node_count; i++)
if (graph->nodes[i].mark)
log_verbose(" %d from %d\n", graph->nodes[i].index, parents[i]);
CLEANUP:
if(parents) free(parents);
if (stack) free(stack);
return rval;
}
int flow_find_max_flow(
const struct Graph *digraph,
const double *capacities,
@@ -13,6 +62,22 @@ int flow_find_max_flow(
{
int rval = 0;
for (int i = 0; i < digraph->node_count; i++)
digraph->nodes[i].mark = 0;
log_verbose("Input graph:\n");
graph_dump(digraph);
log_verbose("Solving flow problem:\n");
log_verbose("%d %d\n", digraph->node_count, digraph->edge_count);
log_verbose("%d %d\n", from->index, to->index);
for (int i = 0; i < digraph->edge_count; i++)
{
log_verbose("%d %d %.4lf\n", digraph->edges[i].from->index,
digraph->edges[i].to->index, capacities[i]);
}
int path_length;
struct Edge **path_edges = 0;
double path_capacity;
@@ -32,6 +97,8 @@ int flow_find_max_flow(
residual_caps[i] = capacities[i];
abort_if(!digraph->edges[i].reverse,
"digraph must have reverse edge information");
abort_if(digraph->edges[i].reverse->reverse != &digraph->edges[i],
"invalid reverse edge");
}
*value = 0;
@@ -43,19 +110,36 @@ int flow_find_max_flow(
if (path_length == 0) break;
log_verbose("Found augmenting path of capacity %.4lf:\n",
path_capacity);
(*value) += path_capacity;
for (int i = 0; i < path_length; i++)
{
struct Edge *e = &digraph->edges[path_edges[i]->index];
log_verbose(" %d %d (%d)\n", e->from->index, e->to->index, e->index);
residual_caps[e->index] -= path_capacity;
residual_caps[e->reverse->index] += path_capacity;
flow[e->index] += path_capacity;
flow[e->reverse->index] -= path_capacity;
}
log_verbose("New residual capacities:\n");
for (int i = 0; i < digraph->edge_count; i++)
{
struct Edge *e = &digraph->edges[i];
if (residual_caps[i] < LP_EPSILON) continue;
log_verbose("%d %d %.4lf (%d)\n", e->from->index, e->to->index, e->index,
residual_caps[e->index]);
}
}
log_verbose("No more paths found.\n");
rval = flow_mark_reachable_nodes(digraph, residual_caps, from);
abort_if(rval, "flow_mark_reachable_nodes failed");
@@ -155,42 +239,108 @@ int flow_find_augmenting_path(
return rval;
}
int flow_mark_reachable_nodes(
struct Graph *graph, double *residual_caps, struct Node *from)
int flow_main(int argc, char **argv)
{
int rval = 0;
struct Node **stack;
int stack_top = 0;
int *edges = 0;
double *capacities = 0;
double *flow = 0;
double flow_value;
stack = (struct Node**) malloc(graph->node_count * sizeof(struct Node*));
abort_if(!stack, "could not allocate stack");
struct Edge **cut_edges = 0;
for (int i = 0; i < graph->node_count; i++)
graph->nodes[i].mark = 0;
FILE *f = fopen("tmp/flow.in", "r");
abort_if(!f, "could not open input file");
from->mark = 1;
stack[stack_top++] = from;
struct Graph graph;
graph_init(&graph);
while(stack_top > 0)
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++)
{
struct Node *n = stack[--stack_top];
int from, to;
double cap;
for (int j = 0; j < n->degree; j++)
{
struct Edge *e = n->adj[j].edge;
struct Node *neighbor = n->adj[j].neighbor;
rval = fscanf(f, "%d %d %lf ", &from, &to, &cap);
abort_if(rval != 3, "invalid input format (edge specification)");
if(neighbor->mark) continue;
if(residual_caps[e->index] <= 0) continue;
stack[stack_top++] = neighbor;
neighbor->mark = 1;
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(stack) free(stack);
if (cut_edges) free(cut_edges);
if (capacities) free(capacities);
if (edges) free(edges);
if (flow) free(flow);
return rval;
}

11
src/flow.h
View File

@@ -1,10 +1,8 @@
//
// Created by isoron on 19/03/15.
//
#ifndef _PROJECT_FLOW_H_
#define _PROJECT_FLOW_H_
#include "graph.h"
int flow_find_augmenting_path(
const struct Graph *graph,
const double *residual_caps,
@@ -22,6 +20,11 @@ int flow_find_max_flow(
double *flow,
double *value);
int flow_mark_reachable_nodes(
const struct Graph *graph, double *residual_caps, struct Node *from);
int flow_main(int argc, char **argv);
#include "graph.h"
#endif //_PROJECT_FLOW_H_

26
src/graph.c
View File

@@ -196,3 +196,29 @@ int get_cut_edges_from_marks(
return 0;
}
int graph_dump(struct Graph *graph)
{
int rval = 0;
log_debug("node_count: %d edge_count: %d\n", graph->node_count, graph->edge_count);
for (int i = 0; i < graph->node_count; i++)
{
struct Node *n = &graph->nodes[i];
log_debug("%3d degree: %d mark: %d\n", n->index, n->degree, n->mark);
}
for (int i = 0; i < graph->edge_count; i++)
{
struct Edge *e = &graph->edges[i];
log_debug("%3d (%d, %d) weight: %d ", e->index, e->from->index, e->to->index, e->weight);
#if LOG_LEVEL >= LOG_LEVEL_VERBOSE
if (e->reverse) printf("reverse: %d ", e->reverse->index);
printf("\n");
#endif
}
CLEANUP:
return rval;
}

2
src/graph.h
View File

@@ -73,4 +73,6 @@ int graph_build_directed_from_undirected
int get_cut_edges_from_marks(
struct Graph *graph, int *cut_edges_count, struct Edge **cut_edges);
int graph_dump(struct Graph *graph);
#endif

503
src/gtsp.c
View File

@@ -1,24 +1,42 @@
#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <assert.h>
#include <getopt.h>
#include "gtsp.h"
#include "geometry.h"
#include "util.h"
#include "flow.h"
#include "branch_and_cut.h"
static double *OPTIMAL_X = 0;
int GTSP_init_data(struct GTSP *data)
{
int rval = 0;
data->clusters = 0;
data->cluster_count = 0;
data->x_coordinates = 0;
data->y_coordinates = 0;
data->graph = (struct Graph *) malloc(sizeof(struct Graph));
abort_if(!data->graph, "could not allocate data->graph");
graph_init(data->graph);
return 0;
CLEANUP:
return rval;
}
void GTSP_free(struct GTSP *data)
{
if (!data) return;
if (data->graph) graph_free(data->graph);
if (data->graph)
{
graph_free(data->graph);
free(data->graph);
}
if (data->clusters) free(data->clusters);
if (data->x_coordinates) free(data->x_coordinates);
if (data->y_coordinates) free(data->y_coordinates);
@@ -78,6 +96,9 @@ int GTSP_create_random_problem(
rval = graph_build(node_count, edge_count, edges, 0, graph);
abort_if(rval, "graph_build failed");
for (int i = 0; i < edge_count; i++)
graph->edges[i].weight = weights[i];
data->graph = graph;
data->clusters = clusters;
data->cluster_count = cluster_count;
@@ -145,7 +166,271 @@ int GTSP_init_lp(struct LP *lp, struct GTSP *data)
return rval;
}
int GTSP_write_data(struct GTSP *data, char *filename)
int GTSP_add_subtour_elimination_cut(
struct LP *lp,
struct Graph *graph,
struct Node *from,
struct Node *to,
struct Edge **cut_edges,
int cut_edges_count)
{
int rval = 0;
char sense = 'G';
double rhs = -2.0;
int newnz = cut_edges_count + 2;
int rmatbeg = 0;
int *rmatind = 0;
double *rmatval = 0;
rmatind = (int *) malloc(newnz * sizeof(int));
abort_if(!rmatind, "could not allocate rmatind");
rmatval = (double *) malloc(newnz * sizeof(double));
abort_if(!rmatval, "could not allocate rmatval");
for (int i = 0; i < cut_edges_count; i++)
{
rmatind[i] = cut_edges[i]->index + graph->node_count;
rmatval[i] = 1.0;
}
rmatind[cut_edges_count] = from->index;
rmatval[cut_edges_count] = -2.0;
rmatind[cut_edges_count + 1] = to->index;
rmatval[cut_edges_count + 1] = -2.0;
log_debug("Generated cut:\n");
for (int i = 0; i < newnz; i++)
log_debug("%8.2f x%d\n", rmatval[i], rmatind[i]);
log_debug(" %c %.2lf\n", sense, rhs);
if (OPTIMAL_X)
{
double sum = 0;
for (int i = 0; i < newnz; i++)
sum += rmatval[i] * OPTIMAL_X[rmatind[i]];
abort_if(sum <= rhs - LP_EPSILON, "cannot add invalid cut");
}
rval = LP_add_rows(lp, 1, newnz, &rhs, &sense, &rmatbeg, rmatind, rmatval);
abort_if(rval, "LP_add_rows failed");
CLEANUP:
if (rmatval) free(rmatval);
if (rmatind) free(rmatind);
return rval;
}
int GTSP_find_exact_subtour_elimination_cuts(
struct LP *lp, struct GTSP *data, int *added_cuts_count)
{
int rval = 0;
int *clusters = data->clusters;
double *x = 0;
double *capacities = 0;
double *flow = 0;
struct Edge **cut_edges = 0;
int *digraph_edges = 0;
struct Graph *graph = data->graph;
int node_count = graph->node_count;
int num_cols = LP_get_num_cols(lp);
x = (double *) malloc(num_cols * sizeof(double));
abort_if(!x, "could not allocate x");
rval = LP_get_x(lp, x);
abort_if(rval, "LP_get_x failed");
struct Graph digraph;
graph_init(&digraph);
digraph_edges = (int *) malloc(8 * graph->edge_count * sizeof(int));
flow = (double *) malloc(4 * graph->edge_count * sizeof(double));
capacities = (double *) malloc(4 * graph->edge_count * sizeof(double));
cut_edges =
(struct Edge **) malloc(
4 * graph->edge_count * sizeof(struct Edge *));
abort_if(!digraph_edges, "could not allocate digraph_edges");
abort_if(!flow, "could not allocate flow");
abort_if(!capacities, "could not allocate capacities");
abort_if(!cut_edges, "could not allocate cut_edges");
// Create four directed edges for each edge of the original graph.
for (int i = 0; i < graph->edge_count; i++)
{
assert(node_count + i < num_cols);
struct Edge *e = &graph->edges[i];
int from = e->from->index;
int to = e->to->index;
digraph_edges[8 * i] = from;
digraph_edges[8 * i + 1] = to;
capacities[4 * i] = x[node_count + i];
digraph_edges[8 * i + 2] = to;
digraph_edges[8 * i + 3] = from;
capacities[4 * i + 1] = 0;
digraph_edges[8 * i + 4] = to;
digraph_edges[8 * i + 5] = from;
capacities[4 * i + 2] = x[node_count + i];
digraph_edges[8 * i + 6] = from;
digraph_edges[8 * i + 7] = to;
capacities[4 * i + 3] = 0;
}
rval = graph_build(node_count, 4 * graph->edge_count, digraph_edges, 1, &digraph);
abort_if(rval, "graph_build failed");
for (int i = 0; i < graph->edge_count; i++)
{
digraph.edges[4 * i].reverse = &digraph.edges[4 * i + 1];
digraph.edges[4 * i + 1].reverse = &digraph.edges[4 * i];
digraph.edges[4 * i + 2].reverse = &digraph.edges[4 * i + 3];
digraph.edges[4 * i + 3].reverse = &digraph.edges[4 * i + 2];
}
int max_x_index = 0;
double max_x = DBL_MIN;
for (int i = 0; i < graph->node_count; i++)
{
struct Node *n = &graph->nodes[i];
if (x[n->index] > max_x)
{
max_x = x[n->index];
max_x_index = i;
}
}
int i = max_x_index;
for (int j = 0; j < digraph.node_count; j++)
{
if (i == j) continue;
if (clusters[i] == clusters[j]) continue;
if (x[i] + x[j] - 1 <= LP_EPSILON) continue;
struct Node *from = &digraph.nodes[i];
struct Node *to = &digraph.nodes[j];
log_verbose("Calculating max flow from %d to %to\n", from->index,
to->index);
double flow_value;
rval = flow_find_max_flow(&digraph, capacities, from, to, flow,
&flow_value);
abort_if(rval, "flow_find_max_flow failed");
log_verbose(" %.2lf\n", flow_value);
if (flow_value >= 2 * (x[i] + x[j] - 1) - LP_EPSILON) continue;
log_verbose("violation: %.2lf >= %.2lf\n", flow_value,
2 * (x[i] + x[j] - 1));
int cut_edges_count;
rval = get_cut_edges_from_marks(&digraph, &cut_edges_count, cut_edges);
abort_if(rval, "get_cut_edges_from_marks failed");
log_verbose("Adding cut for i=%d j=%d, cut edges:\n", i, j);
for (int k = 0; k < cut_edges_count/2; k++)
{
cut_edges[k] = &graph->edges[cut_edges[k*2]->index / 4];
log_verbose(" %d %d\n", cut_edges[k*2]->from->index,
cut_edges[k*2]->to->index);
}
rval = GTSP_add_subtour_elimination_cut(lp, graph, from, to, cut_edges,
cut_edges_count/2);
abort_if(rval, "GTSP_add_subtour_elimination_cut failed");
(*added_cuts_count)++;
}
CLEANUP:
if (digraph_edges) free(digraph_edges);
if (flow) free(flow);
if (cut_edges) free(cut_edges);
if (capacities) free(capacities);
if (x) free(x);
return rval;
}
int GTSP_add_cutting_planes(struct LP *lp, struct GTSP *data)
{
int rval = 0;
double *x = 0;
int num_cols = LP_get_num_cols(lp);
x = (double *) malloc(num_cols * sizeof(double));
abort_if(!x, "could not allocate x");
while (1)
{
int added_cuts_count = 0;
rval = GTSP_find_exact_subtour_elimination_cuts(lp, data,
&added_cuts_count);
abort_if(rval, "GTSP_find_exact_subtour_elimination_cuts failed");
log_verbose("Reoptimizing...\n");
int is_infeasible;
rval = LP_optimize(lp, &is_infeasible);
abort_if(rval, "LP_optimize failed");
if (is_infeasible) break;
double objval;
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("Current solution:\n");
for (int i = 0; i < data->graph->node_count; i++)
if (x[i] > LP_EPSILON) log_verbose(" node %d = %.2f\n", i, x[i]);
for (int i = 0; i < data->graph->edge_count; i++)
{
struct Edge *e = &data->graph->edges[i];
int idx = e->index + data->graph->node_count;
if (x[idx] > LP_EPSILON)
{
log_verbose(" edge (%d, %d) = %.2f\n", e->from->index,
e->to->index, x[idx]);
}
}
log_debug(" obj val = %f\n", objval);
if (added_cuts_count > 0)
{
log_debug("Found %d subtour elimination cuts using exact "
"separation\n", added_cuts_count);
}
else break;
}
CLEANUP:
if (x) free(x);
return rval;
}
int GTSP_write_input_data(struct GTSP *data, char *filename)
{
int rval = 0;
@@ -184,6 +469,8 @@ int GTSP_write_solution(struct GTSP *data, char *filename, double *x)
if (x[i + node_count] > 0.5)
positive_edge_count++;
fprintf(file, "%d %d\n", node_count, edge_count);
fprintf(file, "%d\n", positive_edge_count);
for (int i = 0; i < edge_count; i++)
@@ -195,3 +482,213 @@ int GTSP_write_solution(struct GTSP *data, char *filename, double *x)
return rval;
}
int get_edge_num(int node_count, int from, int to)
{
int idx = node_count;
for (int k = 0; k < from; k++)
idx += node_count - k - 1;
idx += to - from - 1;
return idx;
}
int GTSP_read_x(char *filename, double **p_x)
{
int rval = 0;
int node_count;
int edge_count;
double *x;
FILE *file;
log_info("Reading optimal solution from file %s\n", filename);
file = fopen(filename, "r");
abort_if(!file, "could not open file");
rval = fscanf(file, "%d %d", &node_count, &edge_count);
abort_if(rval != 2, "invalid input format (node and edge count)");
int num_cols = node_count + edge_count;
x = (double *) malloc(num_cols * sizeof(double));
abort_if(!x, "could not allocate x");
for (int i = 0; i < node_count + edge_count; i++) x[i] = 0.0;
rval = fscanf(file, "%d", &edge_count);
abort_if(rval != 1, "invalid input format (positive edge count)");
for (int i = 0; i < edge_count; i++)
{
int from, to, edge;
rval = fscanf(file, "%d %d", &from, &to);
abort_if(rval != 2, "invalid input format (edge endpoints)");
if (from > to) swap(from, to);
edge = get_edge_num(node_count, from, to);
abort_if(edge > num_cols, "invalid edge");
x[from] = x[
to] = 1.0;
x[edge] = 1;
}
for (int i = 0; i < num_cols; i++)
{
if (x[i] <= LP_EPSILON) continue;
log_verbose(" x%-3d = %.2f\n", i, x[i]);
}
*p_x = x;
rval = 0;
CLEANUP:
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 = 20;
static int input_cluster_count = 5;
static int grid_size = 100;
static void GTSP_print_usage(char **argv)
{
printf("wrong usage\n");
}
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':
rval = GTSP_read_x(optarg, &OPTIMAL_X);
abort_if(rval, "GTSP_read_x failed");
break;
case 's':
SEED = atoi(optarg);
break;
case ':':
fprintf(stderr, "option '-%c' requires an argument\n", optopt);
return 1;
case '?':
default:
fprintf(stderr, "option '-%c' is invalid\n", optopt);
return 1;
}
}
CLEANUP:
return rval;
}
int GTSP_main(int argc, char **argv)
{
int rval = 0;
struct BNC bnc;
struct GTSP data;
SEED = (unsigned int) get_real_time() % 10000;
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);
rval = GTSP_create_random_problem(input_node_count, input_cluster_count,
grid_size,
&data);
abort_if(rval, "GTSP_create_random_problem failed");
log_info("Writing random instance to file gtsp.in\n");
rval = GTSP_write_input_data(&data, "gtsp.in");
abort_if(rval, "GTSP_write_problem failed");
bnc.best_obj_val = DBL_MAX;
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;
if (OPTIMAL_X)
log_info("Optimal solution is available. Cuts will be checked.\n");
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);
rval = GTSP_write_solution(&data, "gtsp.out", bnc.best_x);
abort_if(rval, "GTSP_write_solution failed");
CLEANUP:
GTSP_free(&data);
BNC_free(&bnc);
return rval;
}

6
src/gtsp.h
View File

@@ -28,8 +28,12 @@ int GTSP_init_data(struct GTSP *data);
int GTSP_init_lp(struct LP *lp, struct GTSP *data);
int GTSP_write_data(struct GTSP *data, char *filename);
int GTSP_write_input_data(struct GTSP *data, char *filename);
int GTSP_write_solution(struct GTSP *data, char *filename, double *x);
int GTSP_add_cutting_planes(struct LP *lp, struct GTSP *data);
int GTSP_main(int argc, char **argv);
#endif //_PROJECT_GTSP_H_

287
src/main.c
View File

@@ -1,13 +1,7 @@
#include <stdio.h>
#include <stdlib.h>
#include <getopt.h>
#include <float.h>
#include "lp.h"
#include "util.h"
#include "main.h"
#include "tsp.h"
#include "branch_and_cut.h"
#include "gtsp.h"
#include "tsp.h"
#include "flow.h"
char *INPUT_FILENAME = 0;
@@ -16,276 +10,45 @@ int GEOMETRIC_DATA = 0;
int NODE_COUNT_RAND = 0;
int GRID_SIZE_RAND = 100;
static int parse_arguments_tsp(int ac, char **av);
static const struct option options_tab[] = {
{"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}
};
static void print_usage_tsp(char *f);
int test_max_flow()
void GTSP_print_usage(char **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");
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, cap;
rval = fscanf(f, "%d %d %d ", &from, &to, &cap);
abort_if(rval != 3, "invalid input format");
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;
printf("wrong usage\n");
}
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++)
int main(int argc, char **argv)
{
graph.edges[2 * i].reverse = &graph.edges[2 * i + 1];
graph.edges[2 * i + 1].reverse = &graph.edges[2 * i];
}
int c = 0;
int option_index = 0;
flow = (double *) malloc(graph.edge_count * sizeof(double));
abort_if(!flow, "could not allocate flow");
c = getopt_long(argc, argv, "htgf", options_tab, &option_index);
struct Node *from = &graph.nodes[0];
struct Node *to = &graph.nodes[graph.node_count - 1];
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++)
if (c < 0)
{
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;
}
int main_tsp(int ac, char **av)
{
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 = parse_arguments_tsp(ac, av);
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;
}
int main_gtsp(int ac, char **av)
{
int rval = 0;
SEED = (unsigned int) get_real_time();
srand(SEED);
int node_count = 50;
int cluster_count = node_count / 5;
int grid_size = 100;
struct BNC bnc;
struct GTSP data;
rval = GTSP_init_data(&data);
abort_if(rval, "GTSP_init_data failed");
rval = GTSP_create_random_problem(node_count, cluster_count, grid_size,
&data);
abort_if(rval, "GTSP_create_random_problem failed");
rval = GTSP_write_data(&data, "gtsp.in");
abort_if(rval, "GTSP_write_problem failed");
rval = BNC_init(&bnc);
abort_if(rval, "BNC_init failed");
log_info("Setting seed = %d\n", SEED);
srand(SEED);
bnc.best_obj_val = DBL_MAX;
bnc.problem_data = (void *) &data;
bnc.problem_init_lp = (int (*)(struct LP *, void *)) GTSP_init_lp;
rval = BNC_init_lp(&bnc);
abort_if(rval, "BNC_init_lp failed");
log_info("Starting branch-and-cut solver...\n");
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);
rval = GTSP_write_solution(&data, "gtsp.out", bnc.best_x);
abort_if(rval, "GTSP_write_solution failed");
CLEANUP:
GTSP_free(&data);
BNC_free(&bnc);
return rval;
}
static int parse_arguments_tsp(int ac, char **av)
{
int rval = 0;
int c;
if (ac == 1)
{
print_usage_tsp(av[0]);
GTSP_print_usage(argv);
return 1;
}
while ((c = getopt(ac, av, "ab:gk:s:")) != EOF)
{
switch (c)
{
case 'a':;
break;
case 'b':
GRID_SIZE_RAND = atoi(optarg);
break;
case 'h':
GTSP_print_usage(argv);
return 1;
case 'f':
return flow_main(argc, argv);
case 't':
return TSP_main(argc, argv);
case 'g':
GEOMETRIC_DATA = 1;
break;
case 'k':
NODE_COUNT_RAND = atoi(optarg);
break;
case 's':
SEED = (unsigned) atoi(optarg);
break;
case '?':
default:
print_usage_tsp(av[0]);
return 1;
}
return GTSP_main(argc, argv);
}
if (optind < ac) INPUT_FILENAME = av[optind++];
if (optind != ac)
{
print_usage_tsp(av[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;
}
int main(int ac, char **av)
{
return test_max_flow();
// return main_gtsp(ac, av);
// return main_tsp(ac, av);
}
static void print_usage_tsp(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);
}

117
src/tsp.c
View File

@@ -1,9 +1,14 @@
#include <stdio.h>
#include <stdlib.h>
#include "main.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)
{
@@ -113,7 +118,7 @@ int TSP_find_violated_subtour_elimination_cut(
while (!TSP_is_graph_connected(&G, x, &island_count, island_sizes,
island_start, island_nodes))
{
log_verbose("Adding %d BNC_solve_node inequalities...\n", island_count);
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],
@@ -152,6 +157,7 @@ int TSP_find_violated_subtour_elimination_cut(
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;
@@ -454,3 +460,110 @@ double TSP_find_initial_solution(struct TSPData *data)
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);
}

2
src/tsp.h
View File

@@ -45,4 +45,6 @@ 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

11
src/util.h
View File

@@ -1,13 +1,15 @@
#ifndef _PROJECT_UTIL_H_
#define _PROJECT_UTIL_H_
#include <string.h>
#define LOG_LEVEL_ERROR 10
#define LOG_LEVEL_WARNING 20
#define LOG_LEVEL_INFO 30
#define LOG_LEVEL_DEBUG 40
#define LOG_LEVEL_VERBOSE 50
#define LOG_LEVEL LOG_LEVEL_DEBUG
#define LOG_LEVEL LOG_LEVEL_INFO
#if LOG_LEVEL < LOG_LEVEL_DEBUG
#define log_debug(...)
@@ -43,6 +45,13 @@
fprintf(stderr, "%20s:%d " msg "\n", __FILE__, __LINE__); \
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)
void time_printf(const char *fmt, ...);
double get_current_time(void);