Implement column generation

master
Alinson S. Xavier 11 years ago
parent 274a10777e
commit 42360c79a8

@ -100,9 +100,9 @@ else:
c = white.blend(red, 0.1 + 0.9 * edges[k][2]) c = white.blend(red, 0.1 + 0.9 * edges[k][2])
plot = plot + line([all_points[edges[k][0]], all_points[edges[k][1]]], color=c) plot = plot + line([all_points[edges[k][0]], all_points[edges[k][1]]], color=c)
print ('Drawing labels...') #print ('Drawing labels...')
for i in range(node_count): #for i in range(node_count):
plot = plot + text(str(i), all_points[i] + text_offset, color='gray') # plot = plot + text(str(i), all_points[i] + text_offset, color='gray')
print ('Drawing clusters...') print ('Drawing clusters...')
for i in range(cluster_count): for i in range(cluster_count):

@ -22,13 +22,13 @@
#include "util.h" #include "util.h"
#include "gtsp.h" #include "gtsp.h"
static int BNC_solve_node(struct BNC *bnc, int depth); static int solve_node(struct BNC *bnc, int depth);
static int BNC_branch_node(struct BNC *bnc, double *x, int depth); static int branch_node(struct BNC *bnc, double *x, int depth);
static int BNC_is_integral(double *x, int num_cols); static int is_integral(double *x, int num_cols);
static int BNC_find_best_branching_var(double *x, int num_cols); static int find_best_branching_var(double *x, int num_cols);
int BNC_init(struct BNC *bnc) int BNC_init(struct BNC *bnc)
{ {
@ -86,10 +86,10 @@ int BNC_init_lp(struct BNC *bnc)
int BNC_solve(struct BNC *bnc) int BNC_solve(struct BNC *bnc)
{ {
return BNC_solve_node(bnc, 1); return solve_node(bnc, 1);
} }
static int BNC_solve_node(struct BNC *bnc, int depth) static int solve_node(struct BNC *bnc, int depth)
{ {
struct LP *lp = bnc->lp; struct LP *lp = bnc->lp;
double *best_val = &bnc->best_obj_val; double *best_val = &bnc->best_obj_val;
@ -151,11 +151,17 @@ static int BNC_solve_node(struct BNC *bnc, int depth)
goto CLEANUP; goto CLEANUP;
} }
free(x);
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); rval = LP_get_x(lp, x);
abort_if(rval, "LP_get_x failed"); abort_if(rval, "LP_get_x failed");
} }
if (BNC_is_integral(x, num_cols)) if (is_integral(x, num_cols))
{ {
log_debug("Solution is integral\n"); log_debug("Solution is integral\n");
@ -179,8 +185,8 @@ static int BNC_solve_node(struct BNC *bnc, int depth)
else else
{ {
log_debug("Solution is fractional\n"); log_debug("Solution is fractional\n");
rval = BNC_branch_node(bnc, x, depth); rval = branch_node(bnc, x, depth);
abort_if(rval, "BNC_branch_node failed"); abort_if(rval, "branch_node failed");
} }
CLEANUP: CLEANUP:
@ -188,14 +194,14 @@ static int BNC_solve_node(struct BNC *bnc, int depth)
return rval; return rval;
} }
static int BNC_branch_node(struct BNC *bnc, double *x, int depth) static int branch_node(struct BNC *bnc, double *x, int depth)
{ {
int rval = 0; int rval = 0;
struct LP *lp = bnc->lp; struct LP *lp = bnc->lp;
int num_cols = LP_get_num_cols(lp); int num_cols = LP_get_num_cols(lp);
int best_branch_var = BNC_find_best_branching_var(x, num_cols); int best_branch_var = find_best_branching_var(x, num_cols);
log_debug("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); best_branch_var, x[best_branch_var], depth);
@ -204,8 +210,8 @@ static int BNC_branch_node(struct BNC *bnc, double *x, int depth)
rval = LP_change_bound(lp, best_branch_var, 'L', 1.0); rval = LP_change_bound(lp, best_branch_var, 'L', 1.0);
abort_if(rval, "LP_change_bound failed"); abort_if(rval, "LP_change_bound failed");
rval = BNC_solve_node(bnc, depth + 1); rval = solve_node(bnc, depth + 1);
abort_if(rval, "BNC_solve_node failed"); abort_if(rval, "solve_node failed");
rval = LP_change_bound(lp, best_branch_var, 'L', 0.0); rval = LP_change_bound(lp, best_branch_var, 'L', 0.0);
abort_if(rval, "LP_change_bound failed"); abort_if(rval, "LP_change_bound failed");
@ -214,8 +220,8 @@ static int BNC_branch_node(struct BNC *bnc, double *x, int depth)
rval = LP_change_bound(lp, best_branch_var, 'U', 0.0); rval = LP_change_bound(lp, best_branch_var, 'U', 0.0);
abort_if(rval, "LP_change_bound failed"); abort_if(rval, "LP_change_bound failed");
rval = BNC_solve_node(bnc, depth + 1); rval = solve_node(bnc, depth + 1);
abort_if(rval, "BNC_solve_node failed"); abort_if(rval, "solve_node failed");
rval = LP_change_bound(lp, best_branch_var, 'U', 1.0); rval = LP_change_bound(lp, best_branch_var, 'U', 1.0);
abort_if(rval, "LP_change_bound failed"); abort_if(rval, "LP_change_bound failed");
@ -226,7 +232,7 @@ static int BNC_branch_node(struct BNC *bnc, double *x, int depth)
return rval; return rval;
} }
static int BNC_is_integral(double *x, int num_cols) static int is_integral(double *x, int num_cols)
{ {
for (int i = 0; i < num_cols; i++) for (int i = 0; i < num_cols; i++)
if (x[i] > EPSILON && x[i] < 1.0 - EPSILON) if (x[i] > EPSILON && x[i] < 1.0 - EPSILON)
@ -235,7 +241,7 @@ static int BNC_is_integral(double *x, int num_cols)
return 1; return 1;
} }
static int BNC_find_best_branching_var(double *x, int num_cols) static int find_best_branching_var(double *x, int num_cols)
{ {
int best_index = 0; int best_index = 0;
double best_index_frac = 1.0; double best_index_frac = 1.0;

@ -21,18 +21,44 @@
struct BNC struct BNC
{ {
struct LP *lp; struct LP *lp;
/*
* Best integral solution found so far.
*/
double *best_x; double *best_x;
/*
* Value of the best integral solution found so far.
*/
double best_obj_val; double best_obj_val;
int *problem_data; /*
* Pointer to problem-specific data. This pointer is passed to other
* problem-specific functions, such as problem_init_lp and
* problem_add_cutting_planes.
*/
void *problem_data;
int (*problem_init_lp)(struct LP *, void *); /*
* This callback is called at the beginning of the branch-and-cut search.
* It should initialize the LP, create the initial rows and columns.
*/
int (*problem_init_lp)(struct LP *lp, void *data);
int (*problem_add_cutting_planes)(struct LP *, void *); /*
* This callback is called after an LP relaxation is solved and before
* checking whether the solution is integral. It should find and add
* any problem-specific cutting planes to the LP.
*/
int (*problem_add_cutting_planes)(struct LP *lp, void *data);
int (*problem_solution_found)(struct BNC*, void *data, double *x); /*
* This callback is called after a better integral solution has been found.
* This solution satisfies all the problem-specific cuts.
*/
int (*problem_solution_found)(struct BNC *, void *data, double *x);
}; };
int BNC_init(struct BNC *bnc); int BNC_init(struct BNC *bnc);
@ -43,6 +69,4 @@ int BNC_init_lp(struct BNC *bnc);
void BNC_free(struct BNC *bnc); void BNC_free(struct BNC *bnc);
extern int BNC_NODE_COUNT;
#endif //_PROJECT_BRANCH_AND_CUT_H_ #endif //_PROJECT_BRANCH_AND_CUT_H_

@ -22,53 +22,17 @@
int FLOW_MAX_FLOW_COUNT = 0; int FLOW_MAX_FLOW_COUNT = 0;
int flow_mark_reachable_nodes( static int mark_reachable_nodes(
const struct Graph *graph, double *residual_caps, struct Node *from) const struct Graph *graph, double *residual_caps, struct Node *from);
{
int rval = 0;
struct Node **stack;
int stack_top = 0;
int *parents = 0;
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: static int find_augmenting_path(
if(parents) free(parents); const struct Graph *graph,
if (stack) free(stack); const double *residual_caps,
return rval; struct Node *from,
} struct Node *to,
int *path_length,
struct Edge **path_edges,
double *path_capacity);
int flow_find_max_flow( int flow_find_max_flow(
const struct Graph *digraph, const struct Graph *digraph,
@ -85,27 +49,18 @@ int flow_find_max_flow(
for (int i = 0; i < digraph->node_count; i++) for (int i = 0; i < digraph->node_count; i++)
digraph->nodes[i].mark = 0; digraph->nodes[i].mark = 0;
log_verbose("Input graph:\n");
#if LOG_LEVEL >= LOG_LEVEL_VERBOSE
graph_dump(digraph);
#endif
log_verbose("Solving flow problem:\n"); log_verbose("Solving flow problem:\n");
log_verbose("%d %d\n", digraph->node_count, digraph->edge_count); log_verbose("%d %d\n", digraph->node_count, digraph->edge_count);
log_verbose("%d %d\n", from->index, to->index); log_verbose("%d %d\n", from->index, to->index);
for (int i = 0; i < digraph->edge_count; i++) for (int i = 0; i < digraph->edge_count; i++)
{ log_verbose("%d %d %.4lf\n", digraph->edges[i].from->index,
log_verbose("%d %d %.4lf\n", digraph->edges[i].from->index, digraph->edges[i].to->index, capacities[i]);
digraph->edges[i].to->index, capacities[i]);
}
int path_length;
struct Edge **path_edges = 0;
double path_capacity;
int path_length = 0;
double path_capacity = 0;
double *residual_caps = 0; double *residual_caps = 0;
struct Edge **path_edges = 0;
residual_caps = (double *) malloc(digraph->edge_count * sizeof(double)); residual_caps = (double *) malloc(digraph->edge_count * sizeof(double));
abort_if(!residual_caps, "could not allocate residual_caps"); abort_if(!residual_caps, "could not allocate residual_caps");
@ -128,8 +83,8 @@ int flow_find_max_flow(
while (1) while (1)
{ {
flow_find_augmenting_path(digraph, residual_caps, from, to, find_augmenting_path(digraph, residual_caps, from, to, &path_length,
&path_length, path_edges, &path_capacity); path_edges, &path_capacity);
if (path_length == 0) break; if (path_length == 0) break;
@ -142,7 +97,8 @@ int flow_find_max_flow(
{ {
struct Edge *e = &digraph->edges[path_edges[i]->index]; struct Edge *e = &digraph->edges[path_edges[i]->index];
log_verbose(" %d %d (%d)\n", e->from->index, e->to->index, e->index); log_verbose(" %d %d (%d)\n", e->from->index, e->to->index,
e->index);
residual_caps[e->index] -= path_capacity; residual_caps[e->index] -= path_capacity;
residual_caps[e->reverse->index] += path_capacity; residual_caps[e->reverse->index] += path_capacity;
@ -151,24 +107,24 @@ int flow_find_max_flow(
flow[e->reverse->index] -= path_capacity; flow[e->reverse->index] -= path_capacity;
} }
#if LOG_LEVEL >= LOG_LEVEL_VERBOSE
log_verbose("New residual capacities:\n"); log_verbose("New residual capacities:\n");
for (int i = 0; i < digraph->edge_count; i++) for (int i = 0; i < digraph->edge_count; i++)
{ {
#if LOG_LEVEL >= LOG_LEVEL_VERBOSE
struct Edge *e = &digraph->edges[i]; struct Edge *e = &digraph->edges[i];
#endif
if (residual_caps[i] < 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,
residual_caps[e->index]); e->index, residual_caps[e->index]);
} }
#endif
} }
log_verbose("No more paths found.\n"); log_verbose("No more paths found.\n");
rval = flow_mark_reachable_nodes(digraph, residual_caps, from); rval = mark_reachable_nodes(digraph, residual_caps, from);
abort_if(rval, "flow_mark_reachable_nodes failed"); abort_if(rval, "mark_reachable_nodes failed");
CLEANUP: CLEANUP:
if (path_edges) free(path_edges); if (path_edges) free(path_edges);
@ -176,7 +132,56 @@ int flow_find_max_flow(
return rval; return rval;
} }
int flow_find_augmenting_path( int static 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 = 0;
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 static find_augmenting_path(
const struct Graph *graph, const struct Graph *graph,
const double *residual_caps, const double *residual_caps,
struct Node *from, struct Node *from,

@ -19,15 +19,6 @@
#include "graph.h" #include "graph.h"
int flow_find_augmenting_path(
const struct Graph *graph,
const double *residual_caps,
struct Node *from,
struct Node *to,
int *path_length,
struct Edge **path_edges,
double *path_capacity);
int flow_find_max_flow( int flow_find_max_flow(
const struct Graph *digraph, const struct Graph *digraph,
const double *capacities, const double *capacities,
@ -36,7 +27,4 @@ int flow_find_max_flow(
double *flow, double *flow,
double *value); double *value);
int flow_mark_reachable_nodes(
const struct Graph *graph, double *residual_caps, struct Node *from);
#endif //_PROJECT_FLOW_H_ #endif //_PROJECT_FLOW_H_

@ -21,8 +21,6 @@
#include "geometry.h" #include "geometry.h"
#include "util.h" #include "util.h"
/* function for creating a random set of points in unit square */
int generate_random_points_2d( int generate_random_points_2d(
int node_count, int node_count,
int grid_size, int grid_size,
@ -110,8 +108,8 @@ int generate_random_clusters_2d(
for (int j = 0; j < cluster_count; j++) for (int j = 0; j < cluster_count; j++)
{ {
int distance = int distance = get_euclidean_distance(x_coordinates, y_coordinates,
get_euclidean_distance(x_coordinates, y_coordinates, i, j); i, j);
if (distance < closest_distance) if (distance < closest_distance)
{ {
@ -128,18 +126,21 @@ int generate_random_clusters_2d(
} }
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)
{ {
int i,j; int i, j;
for (i = 0; i < node_count; i++){ for (i = 0; i < node_count; i++)
for (j = 0; j < node_count; j++){ {
dist_matrix[i][j] = for (j = 0; j < node_count; j++)
get_euclidean_distance(x_coordinates, y_coordinates, i, j); {
} dist_matrix[i][j] = get_euclidean_distance(x_coordinates,
} y_coordinates, i, j);
return 0; }
}
return 0;
} }
int get_euclidean_distance( int get_euclidean_distance(

@ -78,6 +78,7 @@ int graph_build(
graph->edges[i].index = i; graph->edges[i].index = i;
graph->edges[i].from = &graph->nodes[a]; graph->edges[i].from = &graph->nodes[a];
graph->edges[i].to = &graph->nodes[b]; graph->edges[i].to = &graph->nodes[b];
graph->edges[i].column = -1;
} }
p = graph->adj; p = graph->adj;
@ -136,29 +137,3 @@ int get_cut_edges_from_marks(
return 0; return 0;
} }
int graph_dump(const struct Graph *graph)
{
(void) graph;
#if LOG_LEVEL >= LOG_LEVEL_DEBUG
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 (e->reverse) printf("reverse: %d ", e->reverse->index);
printf("\n");
}
#endif
return 0;
}

@ -37,12 +37,30 @@ struct Node
struct Edge struct Edge
{ {
/*
* Index of the edge. Each edge is numbered from 0 to graph->edge_count.
*/
int index; int index;
/*
* Weight of the edge.
*/
int weight; int weight;
/*
* If this edge corresponds to a column of the LP, this field contains
* the index of that column. Otherwise, this field contains a negative
* value.
*/
int column;
struct Node *from; struct Node *from;
struct Node *to; struct Node *to;
/*
* Pointer to an edge that points in the opposite direction of this one.
* Used by flow algorithms.
*/
struct Edge *reverse; struct Edge *reverse;
}; };
@ -71,9 +89,11 @@ int graph_build(
int is_directed, int is_directed,
struct Graph *graph); struct Graph *graph);
/*
* Returns the list of edges e=uv such that either u or v (but not both) are
* marked nodes.
*/
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);
int graph_dump(const struct Graph *graph);
#endif #endif

@ -0,0 +1,111 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "gtsp-cols.h"
#include "util.h"
int GTSP_find_columns(struct LP *lp, struct GTSP *data)
{
int rval = 0;
double *y = 0;
double initial_time = get_user_time();
int num_rows = LP_get_num_rows(lp);
y = (double *) malloc(num_rows * sizeof(double));
abort_if(!y, "could not allocate y");
rval = LP_get_y(lp, y);
abort_if(rval, "LP_get_y failed");
int violated_count = 0;
int edge_count = data->graph->edge_count;
log_debug("Finding new columns...\n");
for (int i = 0; i < edge_count; i++)
{
struct Edge *e = &data->graph->edges[i];
if (e->column >= 0) continue;
double sum = 0;
for (int j = 0; j < lp->cut_pool_size; j++)
{
struct Row *cut = lp->cut_pool[j];
if (!cut->edges[e->index]) continue;
sum += y[cut->cplex_row_index];
}
if (sum - EPSILON < e->weight) continue;
rval = GTSP_add_column(lp, e);
abort_if(rval, "GTSP_add_column failed");
violated_count++;
}
log_debug(" %d of %d edges are violated\n", violated_count, edge_count);
COLUMNS_TIME += get_user_time() - initial_time;
CLEANUP:
if (y) free(y);
return rval;
}
int GTSP_add_column(struct LP *lp, struct Edge *e)
{
int rval = 0;
int *cmatind = 0;
double *cmatval = 0;
int num_rows = LP_get_num_rows(lp);
int num_cols = LP_get_num_cols(lp);
cmatind = (int *) malloc(num_rows * sizeof(int));
cmatval = (double *) malloc(num_rows * sizeof(double));
abort_if(!cmatind, "could not allocate cmatbeg");
abort_if(!cmatval, "could not allocate cmatval");
int nz = 0;
int cmatbeg = 0;
double lb = 0.0;
double ub = 1.0;
double obj = e->weight;
for (int j = 0; j < lp->cut_pool_size; j++)
{
struct Row *cut = lp->cut_pool[j];
if (!cut->edges[e->index]) continue;
cmatind[nz] = cut->cplex_row_index;
cmatval[nz] = cut->edge_val;
nz++;
}
e->column = num_cols;
LP_add_cols(lp, 1, nz, &obj, &cmatbeg, cmatind, cmatval, &lb, &ub);
CLEANUP:
if (cmatind) free(cmatind);
if (cmatval) free(cmatval);
return rval;
}

@ -0,0 +1,27 @@
/* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef PROJECT_GTSP_COLS_H
#define PROJECT_GTSP_COLS_H
#include "lp.h"
#include "gtsp.h"
int GTSP_find_columns(struct LP *lp, struct GTSP *data);
int GTSP_add_column(struct LP *lp, struct Edge *e);
#endif //PROJECT_GTSP_COLS_H

@ -45,6 +45,15 @@ static int add_comb_cut(
abort_if(!rmatind, "could not allocate rmatind"); abort_if(!rmatind, "could not allocate rmatind");
abort_if(!rmatval, "could not allocate rmatval"); abort_if(!rmatval, "could not allocate rmatval");
cut = (struct Row *) malloc(sizeof(struct Row));
cut->edges = (char *) malloc(graph->edge_count * sizeof(char));
abort_if(!cut, "could not allocate cut");
abort_if(!cut->edges, "could not allocate cut->edges");
for (int i = 0; i < graph->edge_count; i++)
cut->edges[i] = 0;
double rhs = -component_sizes[current_component] - tooth_count + double rhs = -component_sizes[current_component] - tooth_count +
(tooth_count + 1) / 2; (tooth_count + 1) / 2;
@ -57,7 +66,10 @@ static int add_comb_cut(
if (components[clusters[e->from->index]] != current_component) continue; if (components[clusters[e->from->index]] != current_component) continue;
if (components[clusters[e->to->index]] != current_component) continue; if (components[clusters[e->to->index]] != current_component) continue;
rmatind[nz] = e->index; cut->edges[e->index] = 1;
if (e->column < 0) continue;
rmatind[nz] = e->column;
rmatval[nz] = -1.0; rmatval[nz] = -1.0;
nz++; nz++;
@ -76,11 +88,16 @@ static int add_comb_cut(
if (teeth[clusters[from->index]] != teeth[clusters[to->index]]) if (teeth[clusters[from->index]] != teeth[clusters[to->index]])
continue; continue;
cut->edges[e->index] = 1;
if (e->column < 0) continue;
log_verbose(" tooth (%d %d)\n", e->from->index, e->to->index); log_verbose(" tooth (%d %d)\n", e->from->index, e->to->index);
rmatind[nz] = e->index; rmatind[nz] = e->column;
rmatval[nz] = -1.0; rmatval[nz] = -1.0;
nz++; nz++;
} }
#if LOG_LEVEL >= LOG_LEVEL_VERBOSE #if LOG_LEVEL >= LOG_LEVEL_VERBOSE
@ -123,21 +140,20 @@ static int add_comb_cut(
log_verbose("Violation: %.4lf >= %.4lf\n", lhs, rhs); log_verbose("Violation: %.4lf >= %.4lf\n", lhs, rhs);
if (lhs + EPSILON > rhs) // if (lhs + EPSILON > rhs)
{ // {
free(rmatind); // free(rmatind);
free(rmatval); // free(rmatval);
goto CLEANUP; // goto CLEANUP;
} // }
cut = (struct Row *) malloc(sizeof(struct Row));
abort_if(!cut, "could not allocate cut");
cut->nz = nz; cut->nz = nz;
cut->sense = sense; cut->sense = sense;
cut->rhs = rhs; cut->rhs = rhs;
cut->rmatval = rmatval; cut->rmatval = rmatval;
cut->rmatind = rmatind; cut->rmatind = rmatind;
cut->edge_val = -1.0;
cut->edge_count = graph->edge_count;
rval = LP_add_cut(lp, cut); rval = LP_add_cut(lp, cut);
abort_if(rval, "LP_add_cut failed"); abort_if(rval, "LP_add_cut failed");
@ -184,6 +200,7 @@ static 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 < EPSILON) continue; if (x_e < EPSILON) continue;
if (x_e > 1 - EPSILON) continue; if (x_e > 1 - EPSILON) continue;
@ -238,7 +255,6 @@ static int find_teeth(
struct Node *to = e->to; struct Node *to = e->to;
if (x[e->index] < 1 - EPSILON) continue; if (x[e->index] < 1 - EPSILON) continue;
if (to->mark || from->mark) continue; if (to->mark || from->mark) continue;
int z = 0; int z = 0;
@ -342,12 +358,13 @@ static int shrink_clusters(
for (int i = 0; i < graph->edge_count; i++) for (int i = 0; i < graph->edge_count; i++)
{ {
struct Edge *e = &graph->edges[i]; struct Edge *e = &graph->edges[i];
if(e->column < 0) continue;
int from = clusters[e->from->index]; int from = clusters[e->from->index];
int to = clusters[e->to->index]; int to = clusters[e->to->index];
int shunk_e_index = edge_map[from * cluster_count + to]; int shunk_e_index = edge_map[from * cluster_count + to];
shrunken_x[shunk_e_index] += x[e->index]; shrunken_x[shunk_e_index] += x[e->column];
} }
CLEANUP: CLEANUP:

@ -0,0 +1,565 @@
/* Copyright (c) 2015 Armin Sadeghi
* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include <assert.h>
#include "util.h"
#include "gtsp.h"
#include "gtsp-heur.h"
static int large_neighborhood_search(
struct Tour *tour, struct GTSP *data, int *tour_cost);
static int build_tour_from_x(struct GTSP *data, struct Tour *tour, double *x);
static int build_x_from_tour(struct GTSP *data, struct Tour *tour, double *x);
static int optimize_vertex_in_cluster(struct Tour *tour, struct GTSP *data);
static int two_opt(struct Tour *tour, struct GTSP *data);
//static int tour_length(int *tour, struct GTSP *data);
static int get_tour_length(struct Tour *tour, struct GTSP *data);
static int K_opt(struct Tour *tour, struct GTSP *data);
int GTSP_find_initial_tour(struct GTSP *data, int *tour_cost, double *x)
{
int rval = 0;
int cluster_count = data->cluster_count;
struct Tour *tour = 0;
int *uncovered_sets = 0;
int *cluster_in_tour = 0;
tour = (struct Tour *) malloc((cluster_count + 1) * sizeof(struct Tour));
uncovered_sets = (int *) malloc((cluster_count - 1) * sizeof(int));
cluster_in_tour = (int *) malloc(cluster_count * sizeof(int));
abort_if(!tour, "could not allocate tour");
abort_if(!uncovered_sets, "could not allocate uncovered_sets");
abort_if(!cluster_in_tour, "could not allocate cluster_in_tour");
int cluster_num = 0;
for (int i = 0; i < cluster_count; i++)
{
cluster_in_tour[i] = 0;
if (data->node_to_cluster[0] != i)
{
uncovered_sets[cluster_num] = i;
cluster_num++;
}
}
int new_vertex = 1, cost;
tour[0].vertex = 0;
cluster_in_tour[0] = 1;
while (new_vertex < data->cluster_count)
{
int min_dist_vertex = -1;
int min_cost = INT_MAX;
for (int i = 1; i < data->graph->node_count; i++)
{
if (cluster_in_tour[data->node_to_cluster[i]]) continue;
for (int k = 0; k < new_vertex; k++)
{
cost = data->dist_matrix[i][tour[k].vertex];
if (cost < min_cost)
{
min_cost = cost;
min_dist_vertex = i;
}
}
}
assert(min_dist_vertex >= 0);
assert(min_dist_vertex < data->graph->node_count);
int cluster_to_insert = data->node_to_cluster[min_dist_vertex];
cluster_in_tour[cluster_to_insert] = 1;
min_cost = INT_MAX;
int insertion_cost = -1, best_pose = -1, best_vertex = -1;
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];
if (new_vertex == 1)
{
int vertex1 = tour[0].vertex;
insertion_cost =
data->dist_matrix[vertex_to_insert][vertex1] +
data->dist_matrix[vertex1][vertex_to_insert];
}
else
{
int vertex1 = tour[j].vertex;
int vertex2 = tour[tour[j].next].vertex;
insertion_cost =
data->dist_matrix[vertex1][vertex_to_insert] +
data->dist_matrix[vertex_to_insert][vertex2] -
data->dist_matrix[vertex1][vertex2];
}
if (insertion_cost < min_cost)
{
min_cost = insertion_cost;
best_pose = j;
best_vertex = vertex_to_insert;
}
}
}
tour[new_vertex].vertex = best_vertex;
tour[new_vertex].prev = best_pose;
if (new_vertex == 1)
{
tour[new_vertex].next = best_pose;
tour[best_pose].prev = new_vertex;
}
else
{
int temp_vertex = tour[best_pose].next;
tour[new_vertex].next = temp_vertex;
tour[temp_vertex].prev = new_vertex;
}
tour[best_pose].next = new_vertex;
new_vertex += 1;
}
tour[data->cluster_count].vertex = 0;
rval = large_neighborhood_search(tour, data, tour_cost);
abort_if(rval, "large_neighborhood_search failed");
log_info("Initial upper-bound: %d \n", *tour_cost);
rval = build_x_from_tour(data, tour, x);
abort_if(rval, "build_x_from_tour failed");
CLEANUP:
if (tour) free(tour);
if (cluster_in_tour) free(cluster_in_tour);
if (uncovered_sets) free(uncovered_sets);
return rval;
}
int GTSP_improve_solution(struct BNC *bnc, struct GTSP *data, double *x)
{
int rval = 0;
int tour_cost;
double *best_val = &bnc->best_obj_val;
double **best_x = &bnc->best_x;
struct Tour *tour;
tour = (struct Tour *) malloc(
(data->cluster_count + 1) * sizeof(struct Tour));
rval = build_tour_from_x(data, tour, x);
abort_if(rval, "build_tour_from_x failed");
for (int i = 0; i < data->cluster_count; i++)
{
log_debug(" %d\n", tour[i]);
}
rval = large_neighborhood_search(tour, data, &tour_cost);
abort_if(rval, "large_neighborhood_search failed");
if (tour_cost + EPSILON < *best_val)
{
log_info("Local search improved the integral solution\n");
log_info(" before = %f\n", *best_val);
log_info(" after = %f\n", tour_cost);
build_x_from_tour(data, tour, x);
*best_val = tour_cost;
*best_x = x;
}
CLEANUP:
return rval;
}
int static optimize_vertex_in_cluster(struct Tour *tour, struct GTSP *data)
{
int current_cluster;
int insertion_cost;
int **dist_matrix = data->dist_matrix;
int cluster_count = data->cluster_count;
struct Cluster *vertex_set = data->clusters;
for (int i = 0; i < cluster_count; i++)
{
int vertex = tour[i].vertex;
int prev_vertex = tour[tour[i].prev].vertex;
int next_vertex = tour[tour[i].next].vertex;
current_cluster = data->node_to_cluster[vertex];
insertion_cost = dist_matrix[prev_vertex][vertex] +
dist_matrix[vertex][next_vertex];
for (int j = 0; j < vertex_set[current_cluster].size; j++)
{
int vertex_in_cluster = vertex_set[current_cluster].nodes[j];
int cost = dist_matrix[vertex_in_cluster][prev_vertex] +
dist_matrix[vertex_in_cluster][next_vertex];
if (insertion_cost > cost)
{
insertion_cost = cost;
tour[i].vertex = vertex_in_cluster;
}
}
}
return 0;
}
int static two_opt(struct Tour *tour, struct GTSP *data)
{
int **dist_matrix = data->dist_matrix;
for (int i = 0; i < data->cluster_count; i++)
{
int v1 = tour[i].vertex;
int v2 = tour[tour[i].prev].vertex;
int v3 = tour[tour[i].next].vertex;
int v4 = tour[tour[tour[i].next].next].vertex;
int current_cost = dist_matrix[v2][v1] + dist_matrix[v3][v4];
int temp_cost = dist_matrix[v2][v3] + dist_matrix[v1][v4];
if (current_cost > temp_cost)
{
int temp_next = tour[i].next;
int temp_prev = tour[i].prev;
tour[i].next = tour[temp_next].next;
tour[i].prev = temp_next;
tour[tour[temp_next].next].prev = i;
tour[temp_next].next = i;
tour[temp_next].prev = temp_prev;
tour[temp_prev].next = temp_next;
}
}
return 0;
}
static int large_neighborhood_search(
struct Tour *tour, struct GTSP *data, int *tour_cost)
{
int rval = 0;
int cluster_count = data->cluster_count;
int *node_to_cluster = data->node_to_cluster;
int **dist_matrix = data->dist_matrix;
struct Cluster *vertex_set = data->clusters;
//LNS starts
for (int iter = 0; iter < 2000; iter++)
{
//Delete a random vertex
int delete_vertex = rand() % (cluster_count - 1) + 1;
int prev_vertex = tour[delete_vertex].prev;
int next_vertex = tour[delete_vertex].next;
tour[prev_vertex].next = next_vertex;
tour[next_vertex].prev = prev_vertex;
int cluster_to_insert = node_to_cluster[tour[delete_vertex].vertex];
int best_pose = -1;
int best_vertex = -1;
int min_cost = INT_MAX;
for (int i = 0; i < vertex_set[cluster_to_insert].size; i++)
{
int vertex_to_insert = vertex_set[cluster_to_insert].nodes[i];
int next_edge = tour[0].next;
for (int j = 1; j < cluster_count; j++)
{
int vertex1 = tour[next_edge].vertex;
int vertex2 = tour[tour[next_edge].next].vertex;
int insert_cost = dist_matrix[vertex1][vertex_to_insert] +
dist_matrix[vertex_to_insert][vertex2] -
dist_matrix[vertex1][vertex2];
if (insert_cost < min_cost)
{
min_cost = insert_cost;
best_pose = next_edge;
best_vertex = vertex_to_insert;
}
next_edge = tour[next_edge].next;
}
}
assert(best_pose >= 0);
assert(best_vertex >= 0);
next_vertex = tour[best_pose].next;
tour[delete_vertex].prev = best_pose;
tour[delete_vertex].vertex = best_vertex;
tour[delete_vertex].next = next_vertex;
tour[best_pose].next = delete_vertex;
tour[next_vertex].prev = delete_vertex;
rval = optimize_vertex_in_cluster(tour, data);
abort_if(rval, "optimize_vertex_in_cluster failed");
}
rval = K_opt(tour, data);
abort_if(rval, "two_opt failed");
rval = two_opt(tour, data);
abort_if(rval, "two_opt failed");
*tour_cost = get_tour_length(tour, data);
CLEANUP:
//if (vertex_seq) free(vertex_seq);
return rval;
}
int static get_tour_length(struct Tour *tour, struct GTSP *data)
{
int tour_cost = 0;
for (int i = 0; i < data->cluster_count; i++)
{
int vertex1 = tour[i].vertex;
int vertex2 = tour[tour[i].next].vertex;
tour_cost += data->dist_matrix[vertex1][vertex2];
}
return tour_cost;
}
int static K_opt(struct Tour *tour, struct GTSP *data)
{
int cluster_count = data->cluster_count;
int l;
for (int i = 0; i < cluster_count - 3; i++)
{
int location_in_path = 0;
for (int k = 0; k < i; k++)
location_in_path = tour[location_in_path].next;
int current_vertex = tour[location_in_path].vertex;
for (int k = 3; k < cluster_count - i - 2; k++)
{
int first_next_location = tour[location_in_path].next;
int first_next_vertex = tour[tour[location_in_path].next].vertex;
int next_location = tour[location_in_path].next;
for (l = 0; l < k; l++)
{
next_location = tour[next_location].next;
}
int next_vertex = tour[next_location].vertex;
int next_next_location = tour[next_location].next;
int next_next_vertex = tour[next_next_location].vertex;
if (next_next_vertex == current_vertex) break;
int cost1 = data->dist_matrix[current_vertex][first_next_vertex] +
data->dist_matrix[next_vertex][next_next_vertex];
int cost2 = data->dist_matrix[current_vertex][next_vertex] +
data->dist_matrix[first_next_vertex][next_next_vertex];
if (cost2 < cost1)
{
int tmp_location = next_location;
for (int m = 0; m < l + 1; m++)
{
int tmp_vertex = tour[tmp_location].next;
tour[tmp_location].next = tour[tmp_location].prev;
tour[tmp_location].prev = tmp_vertex;
tmp_location = tour[tmp_location].next;
}
tour[location_in_path].next = next_location;
tour[next_location].prev = location_in_path;
tour[first_next_location].next = next_next_location;
tour[next_next_location].prev = first_next_location;
}
}
}
return 0;
}
int static build_tour_from_x(struct GTSP *data, struct Tour *tour, double *x)
{
int rval = 0;
int *cluster_mark = 0;
struct Node **stack = 0;
int stack_top = 0;
struct Graph *graph = data->graph;
const int node_count = graph->node_count;
const int edge_count = graph->edge_count;
cluster_mark = (int *) malloc(data->cluster_count * sizeof(int));
abort_if(!cluster_mark, "could not allocate cluster_mark");
stack = (struct Node **) malloc(graph->node_count * sizeof(struct Node *));
abort_if(!stack, "could not allocate stack");
for (int i = 0; i < node_count; i++)
graph->nodes[i].mark = 0;
for (int i = 0; i < data->cluster_count; i++)
cluster_mark[i] = 0;
int initial = -1;
for (int i = 0; i < edge_count; i++)
{
int col = graph->edges[i].column;
if (col < 0) continue;
if (x[col] > 1.0 - EPSILON)
{
initial = i;
break;
}
}
initial = graph->edges[initial].from->index;
abort_if(initial < 0, "no initial node");
stack[stack_top++] = &graph->nodes[initial];
graph->nodes[initial].mark = 1;
tour[0].vertex = graph->nodes[initial].index;
int next_vertex = 1;
while (stack_top > 0)
{
struct Node *n = stack[--stack_top];
cluster_mark[data->node_to_cluster[n->index]]++;
for (int i = 0; i < n->degree; i++)
{
struct Adjacency *adj = &n->adj[i];
struct Node *neighbor = adj->neighbor;
if (neighbor->mark) continue;
if (adj->edge->column < 0) continue;
if (x[adj->edge->column] < EPSILON) continue;
stack[stack_top++] = neighbor;
tour[next_vertex].vertex = neighbor->index;
next_vertex++;
neighbor->mark = 1;
}
}
for (int i = 0; i < data->cluster_count; i++)
{
if (i == 0)
{
tour[i].prev = data->cluster_count - 1;
}
else
{
tour[i].prev = i - 1;
}
if (i == data->cluster_count - 1)
{
tour[i].next = 0;
}
else
{
tour[i].next = i + 1;
}
}
CLEANUP:
if (stack) free(stack);
if (cluster_mark) free(cluster_mark);
return rval;
}
static int build_x_from_tour(struct GTSP *data, struct Tour *tour, double *x)
{
int rval = 0;
int *edge_map = 0;
int node_count = data->graph->node_count;
int edge_count = data->graph->edge_count;
edge_map = (int *) malloc(node_count * node_count * sizeof(int));
abort_if(!edge_map, "could not allocate edge_map");
rval = GTSP_build_edge_map(data, edge_map);
abort_if(rval, "GTSP_build_edge_map failed");
for (int i = 0; i < edge_count; i++)
x[i] = 0.0;
int k = 0;
int next_vertex = tour[0].next;
int current_vertex = tour[0].vertex;
for (int i = 0; i < data->cluster_count; i++)
{
int from = current_vertex;
int to = tour[next_vertex].vertex;
current_vertex = tour[next_vertex].vertex;
next_vertex = tour[next_vertex].next;
struct Edge *e = &data->graph->edges[edge_map[from * node_count + to]];
if (e->column < 0) e->column = k++;
x[e->column] = 1.0;
}
CLEANUP:
if (edge_map) free(edge_map);
return rval;
}

@ -0,0 +1,27 @@
/* Copyright (c) 2015 Armin Sadeghi
* Copyright (c) 2015 Alinson Xavier
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef PROJECT_GTSP_HEUR_H
#define PROJECT_GTSP_HEUR_H
#include "gtsp.h"
int GTSP_find_initial_tour(struct GTSP *data, int *value, double *x);
int GTSP_improve_solution(struct BNC *bnc, struct GTSP *data, double *x);
#endif //PROJECT_GTSP_HEUR_H

@ -19,7 +19,8 @@
#include "util.h" #include "util.h"
#include "flow.h" #include "flow.h"
static void deactivate_cluster_node(double *capacities, struct Node *cluster_node) static void deactivate_cluster_node(
double *capacities, struct Node *cluster_node)
{ {
for (int i = 0; i < cluster_node->degree; i++) for (int i = 0; i < cluster_node->degree; i++)
{ {
@ -63,7 +64,10 @@ static int 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] < EPSILON) continue; int col = graph->edges[i].column;
if (col < 0) continue;
if (x[col] < EPSILON) continue;
struct Edge *e = &graph->edges[i]; struct Edge *e = &graph->edges[i];
int from = e->from->index; int from = e->from->index;
@ -71,7 +75,7 @@ static int build_flow_digraph(
digraph_edges[ke++] = from; digraph_edges[ke++] = from;
digraph_edges[ke++] = to; digraph_edges[ke++] = to;
capacities[kc++] = x[i]; capacities[kc++] = x[col];
digraph_edges[ke++] = to; digraph_edges[ke++] = to;
digraph_edges[ke++] = from; digraph_edges[ke++] = from;
@ -79,7 +83,7 @@ static int build_flow_digraph(
digraph_edges[ke++] = to; digraph_edges[ke++] = to;
digraph_edges[ke++] = from; digraph_edges[ke++] = from;
capacities[kc++] = x[i]; capacities[kc++] = x[col];
digraph_edges[ke++] = from; digraph_edges[ke++] = from;
digraph_edges[ke++] = to; digraph_edges[ke++] = to;
@ -122,51 +126,68 @@ static int build_flow_digraph(
} }
static int add_subtour_cut( static int add_subtour_cut(
struct LP *lp, struct Edge **cut_edges, int cut_edges_count) struct LP *lp,
struct Edge **cut_edges,
int cut_edges_count,
struct Graph *graph)
{ {
int rval = 0; int rval = 0;
char sense = 'G'; char sense = 'G';
double rhs = 2.0; double rhs = 2.0;
int newnz = cut_edges_count; int nz = cut_edges_count;
int *rmatind = 0; int *rmatind = 0;
double *rmatval = 0; double *rmatval = 0;
struct Row *cut = 0;
rmatind = (int *) malloc(newnz * sizeof(int)); rmatind = (int *) malloc(nz * sizeof(int));
abort_if(!rmatind, "could not allocate rmatind"); abort_if(!rmatind, "could not allocate rmatind");
rmatval = (double *) malloc(newnz * sizeof(double)); rmatval = (double *) malloc(nz * sizeof(double));
abort_if(!rmatval, "could not allocate rmatval"); abort_if(!rmatval, "could not allocate rmatval");
nz = 0;
for (int i = 0; i < cut_edges_count; i++) for (int i = 0; i < cut_edges_count; i++)
{ {
rmatind[i] = cut_edges[i]->index; if (cut_edges[i]->column < 0) continue;
rmatval[i] = 1.0; rmatind[nz] = cut_edges[i]->column;
rmatval[nz] = 1.0;
nz++;
} }
log_verbose("Generated cut:\n"); log_verbose("Generated cut:\n");
for (int i = 0; i < newnz; i++) for (int i = 0; i < nz; i++)
log_verbose(" %8.2f x%d\n", rmatval[i], rmatind[i]); log_verbose(" %8.2f x%d\n", rmatval[i], rmatind[i]);
log_verbose(" %c %.2lf\n", sense, rhs); log_verbose(" %c %.2lf\n", sense, rhs);
if (OPTIMAL_X) if (OPTIMAL_X)
{ {
double sum = 0; double sum = 0;
for (int i = 0; i < newnz; i++) for (int i = 0; i < nz; i++)
sum += rmatval[i] * OPTIMAL_X[rmatind[i]]; sum += rmatval[i] * OPTIMAL_X[rmatind[i]];
abort_if(sum <= rhs - EPSILON, "cannot add invalid cut"); abort_if(sum <= rhs - EPSILON, "cannot add invalid cut");
} }
struct Row *cut = 0;
cut = (struct Row *) malloc(sizeof(struct Row)); cut = (struct Row *) malloc(sizeof(struct Row));
abort_if(!cut, "could not allocate cut"); abort_if(!cut, "could not allocate cut");
cut->nz = newnz; cut->edges = (char *) malloc(graph->edge_count * sizeof(char));
abort_if(!cut->edges, "could not allocate cut->edges");
for (int i = 0; i < graph->edge_count; i++)
cut->edges[i] = 0;
for (int i = 0; i < cut_edges_count; i++)
cut->edges[cut_edges[i]->index] = 1;
cut->nz = nz;
cut->sense = sense; cut->sense = sense;
cut->rhs = rhs; cut->rhs = rhs;
cut->rmatval = rmatval; cut->rmatval = rmatval;
cut->rmatind = rmatind; cut->rmatind = rmatind;
cut->edge_count = graph->edge_count;
cut->edge_val = 1.0;
rval = LP_add_cut(lp, cut); rval = LP_add_cut(lp, cut);
abort_if(rval, "LP_add_cut failed"); abort_if(rval, "LP_add_cut failed");
@ -203,7 +224,7 @@ int GTSP_find_exact_subtour_cuts(struct LP *lp, struct GTSP *data)
abort_if(rval, "LP_get_x failed"); abort_if(rval, "LP_get_x failed");
#if LOG_LEVEL >= LOG_LEVEL_DEBUG #if LOG_LEVEL >= LOG_LEVEL_DEBUG
if(strlen(FRAC_SOLUTION_FILENAME) > 0) if (strlen(FRAC_SOLUTION_FILENAME) > 0)
{ {
rval = GTSP_write_solution(data, FRAC_SOLUTION_FILENAME, x); rval = GTSP_write_solution(data, FRAC_SOLUTION_FILENAME, x);
abort_if(rval, "GTSP_write_solution failed"); abort_if(rval, "GTSP_write_solution failed");
@ -268,7 +289,7 @@ int GTSP_find_exact_subtour_cuts(struct LP *lp, struct GTSP *data)
log_verbose(" %d %d (%d)\n", cut_edges[k]->from->index, log_verbose(" %d %d (%d)\n", cut_edges[k]->from->index,
cut_edges[k]->to->index, cut_edges[k]->index); cut_edges[k]->to->index, cut_edges[k]->index);
rval = add_subtour_cut(lp, cut_edges, cut_edges_count); rval = add_subtour_cut(lp, cut_edges, cut_edges_count, graph);
abort_if(rval, "add_subtour_cut failed"); abort_if(rval, "add_subtour_cut failed");
cuts_count++; cuts_count++;

@ -17,19 +17,12 @@
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <assert.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"
#include "gtsp-cols.h"
int large_neighborhood_search(
struct Tour *tour, struct GTSP *data, int *tour_cost);
int build_edge_map(struct GTSP *gtsp, int *edge_map);
double *OPTIMAL_X = 0;
int GTSP_init_data(struct GTSP *data) int GTSP_init_data(struct GTSP *data)
{ {
@ -76,6 +69,7 @@ int GTSP_create_random_problem(
int *weights = 0; int *weights = 0;
int *node_to_cluster = 0; int *node_to_cluster = 0;
struct Cluster *clusters = 0; struct Cluster *clusters = 0;
int *edge_map = 0;
int **dist_matrix = 0; int **dist_matrix = 0;
@ -180,6 +174,7 @@ int GTSP_create_random_problem(
data->clusters = clusters; data->clusters = clusters;
CLEANUP: CLEANUP:
if (edge_map) free(edge_map);
if (weights) free(weights); if (weights) free(weights);
if (edges) free(edges); if (edges) free(edges);
if (rval) if (rval)
@ -196,31 +191,22 @@ int GTSP_init_lp(struct LP *lp, struct GTSP *data)
int rval = 0; int rval = 0;
int edge_count = data->graph->edge_count; int edge_count = data->graph->edge_count;
int cluster_count = data->cluster_count;
int *clusters = data->node_to_cluster;
struct Edge *edges = data->graph->edges; struct Edge *edges = data->graph->edges;
for (int i = 0; i < cluster_count; i++)
{
rval = LP_new_row(lp, 'E', 2.0);
abort_if(rval, "LP_new_row failed");
}
double lb = 0.0; double lb = 0.0;
double ub = 1.0; double ub = 1.0;
int cmatbeg = 0; int empty = 0;
double emptyd = 0.0;
int k = 0;
for (int i = 0; i < edge_count; i++) for (int i = 0; i < edge_count; i++)
{ {
struct Node *from = edges[i].from; if (edges[i].column < 0) continue;
struct Node *to = edges[i].to; edges[i].column = k++;
double obj = (double) edges[i].weight; double obj = (double) edges[i].weight;
double cmatval[] = {1.0, 1.0};
int cmatind[] = {clusters[from->index], clusters[to->index]};
rval = LP_add_cols(lp, 1, 2, &obj, &cmatbeg, cmatind, cmatval, &lb, rval = LP_add_cols(lp, 1, 0, &obj, &empty, &empty, &emptyd, &lb, &ub);
&ub);
abort_if(rval, "LP_add_cols failed"); abort_if(rval, "LP_add_cols failed");
} }
@ -251,6 +237,7 @@ int GTSP_add_cutting_planes(struct LP *lp, struct GTSP *data)
int original_cut_pool_size; int original_cut_pool_size;
int added_cuts_count; int added_cuts_count;
// Generalized subtour cuts
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);
@ -261,6 +248,7 @@ int GTSP_add_cutting_planes(struct LP *lp, struct GTSP *data)
if (added_cuts_count > 0) if (added_cuts_count > 0)
continue; continue;
// Generalized comb cuts
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);
@ -271,6 +259,17 @@ int GTSP_add_cutting_planes(struct LP *lp, struct GTSP *data)
if (added_cuts_count > 0) if (added_cuts_count > 0)
continue; continue;
// Column generation
int original_cols_count = LP_get_num_cols(lp);
rval = GTSP_find_columns(lp, data);
abort_if(rval, "GTSP_find_columns failed");
int added_cols_count = LP_get_num_cols(lp) - original_cols_count;
if (added_cols_count > 0)
continue;
break; break;
} }
@ -309,9 +308,13 @@ int GTSP_print_solution(struct GTSP *data, double *x)
int edge_count = data->graph->edge_count; int edge_count = data->graph->edge_count;
for (int i = 0; i < edge_count; i++) for (int i = 0; i < edge_count; i++)
if (x[i] > EPSILON) {
int col = edges[i].column;
if (col < 0) continue;
if (x[col] > EPSILON)
log_info(" %-3d %-3d %8.4lf\n", edges[i].from->index, log_info(" %-3d %-3d %8.4lf\n", edges[i].from->index,
edges[i].to->index, x[i]); edges[i].to->index, x[col]);
}
return 0; return 0;
} }
@ -330,15 +333,27 @@ 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] > EPSILON) {
positive_edge_count++; struct Edge *e = &edges[i];
if (e->column < 0) continue;
if (x[e->column] < EPSILON) continue;
positive_edge_count++;
}
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] > EPSILON) {
fprintf(file, "%d %d %.4lf\n", edges[i].from->index, struct Edge *e = &edges[i];
edges[i].to->index, x[i]);
if (e->column < 0) continue;
if (x[e->column] < EPSILON) continue;
fprintf(file, "%d %d %.4lf\n", edges[i].from->index, edges[i].to->index,
x[e->column]);
}
CLEANUP: CLEANUP:
if (file) fclose(file); if (file) fclose(file);
@ -376,8 +391,8 @@ int GTSP_read_solution(struct GTSP *gtsp, char *filename, double **p_x)
edge_map = (int *) malloc(node_count * node_count * sizeof(int)); edge_map = (int *) malloc(node_count * node_count * sizeof(int));
abort_if(!edge_map, "could not allocate edge_map"); abort_if(!edge_map, "could not allocate edge_map");
rval = build_edge_map(gtsp, edge_map); rval = GTSP_build_edge_map(gtsp, edge_map);
abort_if(rval, "build_edge_map failed"); abort_if(rval, "GTSP_build_edge_map failed");
for (int i = 0; i < edge_count; i++) for (int i = 0; i < edge_count; i++)
{ {
@ -406,25 +421,6 @@ int GTSP_read_solution(struct GTSP *gtsp, char *filename, double **p_x)
return rval; return rval;
} }
int build_edge_map(struct GTSP *gtsp, int *edge_map)
{
int node_count = gtsp->graph->node_count;
int k = 0;
for (int i = 0; i < node_count; i++)
{
for (int j = i + 1; j < node_count; j++)
{
if (gtsp->node_to_cluster[i] == gtsp->node_to_cluster[j]) continue;
edge_map[i * node_count + j] = k;
edge_map[j * node_count + i] = k;
k++;
}
}
return 0;
}
int GTSP_check_solution(struct GTSP *data, double *x) int GTSP_check_solution(struct GTSP *data, double *x)
{ {
int rval = 0; int rval = 0;
@ -445,11 +441,14 @@ 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 - EPSILON && x[i] > EPSILON, int col = graph->edges[i].column;
"solution is not integral: x%d = %.4lf", i, x[i]); if (col < 0) continue;
abort_iff(x[i] > 1.0 + EPSILON || x[i] < 0.0 - EPSILON, abort_iff(x[col] < 1.0 - EPSILON && x[col] > EPSILON,
"value out of bounds: x%d = %.4lf", i, x[i]); "solution is not integral: x%d = %.4lf", i, x[col]);
abort_iff(x[col] > 1.0 + EPSILON || x[col] < 0.0 - EPSILON,
"value out of bounds: x%d = %.4lf", i, x[col]);
} }
for (int i = 0; i < node_count; i++) for (int i = 0; i < node_count; i++)
@ -464,8 +463,15 @@ int GTSP_check_solution(struct GTSP *data, double *x)
abort_if(initial == edge_count, "no initial node"); abort_if(initial == edge_count, "no initial node");
stack[stack_top++] = graph->edges[initial].from; for (int i = 0; i < edge_count; i++)
graph->edges[initial].from->mark = 1; {
if (graph->edges[i].column == initial)
{
stack[stack_top++] = graph->edges[i].from;
graph->edges[i].from->mark = 1;
break;
}
}
while (stack_top > 0) while (stack_top > 0)
{ {
@ -478,7 +484,9 @@ 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] < EPSILON) continue;
if (adj->edge->column < 0) continue;
if (x[adj->edge->column] < EPSILON) continue;
stack[stack_top++] = neighbor; stack[stack_top++] = neighbor;
neighbor->mark = 1; neighbor->mark = 1;
@ -486,10 +494,7 @@ int GTSP_check_solution(struct GTSP *data, double *x)
} }
for (int i = 0; i < data->cluster_count; i++) for (int i = 0; i < data->cluster_count; i++)
{
abort_iff(cluster_mark[i] > 1, "cluster %d visited multiple times", i);
abort_iff(cluster_mark[i] < 1, "cluster %d not visited", i); abort_iff(cluster_mark[i] < 1, "cluster %d not visited", i);
}
log_info(" solution is valid\n"); log_info(" solution is valid\n");
@ -501,14 +506,7 @@ int GTSP_check_solution(struct GTSP *data, double *x)
int GTSP_solution_found(struct BNC *bnc, struct GTSP *data, double *x) int GTSP_solution_found(struct BNC *bnc, struct GTSP *data, double *x)
{ {
UNUSED(bnc);
int rval = 0; int rval = 0;
int tour_cost;
double *best_val = &bnc->best_obj_val;
struct Tour *tour;
tour = (struct Tour *) malloc(data->cluster_count * sizeof(struct Tour));
if (strlen(SOLUTION_FILENAME) > 0) if (strlen(SOLUTION_FILENAME) > 0)
{ {
@ -517,20 +515,6 @@ int GTSP_solution_found(struct BNC *bnc, struct GTSP *data, double *x)
abort_if(rval, "GTSP_write_solution failed"); abort_if(rval, "GTSP_write_solution failed");
} }
rval = build_tour_from_x(data, tour, x);
abort_if(rval, "build_tour_from_x failed");
rval = large_neighborhood_search(tour, data, &tour_cost);
abort_if(rval, "large_neighborhood_search failed");
if (tour_cost + EPSILON < *best_val)
{
log_info("Local search improve the integral solution\n");
log_info(" obj val = %f\n", *best_val);
log_info(" after LNS = %d\n", tour_cost);
*best_val = tour_cost;
}
log_info("Checking solution...\n"); log_info("Checking solution...\n");
rval = GTSP_check_solution(data, x); rval = GTSP_check_solution(data, x);
abort_if(rval, "GTSP_check_solution failed"); abort_if(rval, "GTSP_check_solution failed");
@ -539,513 +523,21 @@ int GTSP_solution_found(struct BNC *bnc, struct GTSP *data, double *x)
return rval; return rval;
} }
int build_x_from_tour(struct GTSP *data, struct Tour *tour, double *x) int GTSP_build_edge_map(struct GTSP *gtsp, int *edge_map)
{ {
int rval = 0; int node_count = gtsp->graph->node_count;
int *edge_map = 0;
int node_count = data->graph->node_count;
int edge_count = data->graph->edge_count;
edge_map = (int *) malloc(node_count * node_count * sizeof(int));
abort_if(!edge_map, "could not allocate edge_map");
rval = build_edge_map(data, edge_map);
abort_if(rval, "build_edge_map failed");
for (int i = 0; i < edge_count; i++)
x[i] = 0.0;
int next_vertex = tour[0].next;
int current_vertex = tour[0].vertex;
for (int i = 0; i < data->cluster_count; i++)
{
int from = current_vertex;
int to = tour[next_vertex].vertex;
current_vertex = tour[next_vertex].vertex;
next_vertex = tour[next_vertex].next;
x[edge_map[from * node_count + to]] = 1.0;
}
CLEANUP:
if (edge_map) free(edge_map);
return rval;
}
int inital_tour_value(struct GTSP *data, int *tour_cost, double *x)
{
int rval = 0;
int cluster_count = data->cluster_count;
struct Tour *tour = 0;
int *uncovered_sets = 0;
int *cluster_in_tour = 0;
tour = (struct Tour *) malloc((cluster_count + 1) * sizeof(struct Tour));
uncovered_sets = (int *) malloc((cluster_count - 1) * sizeof(int));
cluster_in_tour = (int *) malloc(cluster_count * sizeof(int));
abort_if(!tour, "could not allocate tour");
abort_if(!uncovered_sets, "could not allocate uncovered_sets");
abort_if(!cluster_in_tour, "could not allocate cluster_in_tour");
int cluster_num = 0;
for (int i = 0; i < cluster_count; i++)
{
cluster_in_tour[i] = 0;
if (data->node_to_cluster[0] != i)
{
uncovered_sets[cluster_num] = i;
cluster_num++;
}
}
int new_vertex = 1, cost;
tour[0].vertex = 0;
cluster_in_tour[0] = 1;
while (new_vertex < data->cluster_count)
{
int min_dist_vertex = -1;
int min_cost = INT_MAX;
for (int i = 1; i < data->graph->node_count; i++)
{
if (cluster_in_tour[data->node_to_cluster[i]]) continue;
for (int k = 0; k < new_vertex; k++)
{
cost = data->dist_matrix[i][tour[k].vertex];
if (cost < min_cost)
{
min_cost = cost;
min_dist_vertex = i;
}
}
}
assert(min_dist_vertex >= 0);
int cluster_to_insert = data->node_to_cluster[min_dist_vertex];
cluster_in_tour[cluster_to_insert] = 1;
min_cost = INT_MAX;
int insertion_cost = -1, best_pose = -1, best_vertex = -1;
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];
if (new_vertex == 1)
{
int vertex1 = tour[0].vertex;
insertion_cost =
data->dist_matrix[vertex_to_insert][vertex1] +
data->dist_matrix[vertex1][vertex_to_insert];
}
else
{
int vertex1 = tour[j].vertex;
int vertex2 = tour[tour[j].next].vertex;
insertion_cost =
data->dist_matrix[vertex1][vertex_to_insert] +
data->dist_matrix[vertex_to_insert][vertex2] -
data->dist_matrix[vertex1][vertex2];
}
if (insertion_cost < min_cost)
{
min_cost = insertion_cost;
best_pose = j;
best_vertex = vertex_to_insert;
}
}
}
tour[new_vertex].vertex = best_vertex;
tour[new_vertex].prev = best_pose;
if (new_vertex == 1)
{
tour[new_vertex].next = best_pose;
tour[best_pose].prev = new_vertex;
}
else
{
int temp_vertex = tour[best_pose].next;
tour[new_vertex].next = temp_vertex;
tour[temp_vertex].prev = new_vertex;
}
tour[best_pose].next = new_vertex;
new_vertex += 1;
}
tour[data->cluster_count].vertex = 0;
rval = large_neighborhood_search(tour, data, tour_cost);
abort_if(rval, "large_neighborhood_search failed");
log_info("Initial upper-bound: %d \n", *tour_cost);
rval = build_x_from_tour(data, tour, x);
abort_if(rval, "build_x_from_tour failed");
CLEANUP:
if (tour) free(tour);
if (cluster_in_tour) free(cluster_in_tour);
if (uncovered_sets) free(uncovered_sets);
return rval;
}
int optimize_vertex_in_cluster(struct Tour *tour, struct GTSP *data)
{
int current_cluster;
int insertion_cost;
int **dist_matrix = data->dist_matrix;
int cluster_count = data->cluster_count;
struct Cluster *vertex_set = data->clusters;
for (int i = 0; i < cluster_count; i++)
{
int vertex = tour[i].vertex;
int prev_vertex = tour[tour[i].prev].vertex;
int next_vertex = tour[tour[i].next].vertex;
current_cluster = data->node_to_cluster[vertex];
insertion_cost = dist_matrix[prev_vertex][vertex] +
dist_matrix[vertex][next_vertex];
for (int j = 0; j < vertex_set[current_cluster].size; j++)
{
int vertex_in_cluster = vertex_set[current_cluster].nodes[j];
int cost = dist_matrix[vertex_in_cluster][prev_vertex] +
dist_matrix[vertex_in_cluster][next_vertex];
if (insertion_cost > cost)
{
insertion_cost = cost;
tour[i].vertex = vertex_in_cluster;
}
}
}
return 0;
}
int two_opt(struct Tour *tour, struct GTSP *data)
{
int **dist_matrix = data->dist_matrix;
for (int i = 0; i < data->cluster_count; i++)
{
int v1 = tour[i].vertex;
int v2 = tour[tour[i].prev].vertex;
int v3 = tour[tour[i].next].vertex;
int v4 = tour[tour[tour[i].next].next].vertex;
int current_cost = dist_matrix[v2][v1] + dist_matrix[v3][v4];
int temp_cost = dist_matrix[v2][v3] + dist_matrix[v1][v4];
if (current_cost > temp_cost)
{
int temp_next = tour[i].next;
int temp_prev = tour[i].prev;
tour[i].next = tour[temp_next].next;
tour[i].prev = temp_next;
tour[tour[temp_next].next].prev = i;
tour[temp_next].next = i;
tour[temp_next].prev = temp_prev;
tour[temp_prev].next = temp_next;
}
}
return 0;
}
int large_neighborhood_search(
struct Tour *tour, struct GTSP *data, int *tour_cost)
{
int rval = 0;
int cluster_count = data->cluster_count;
int *clusters = data->node_to_cluster;
int **dist_matrix = data->dist_matrix;
struct Cluster *vertex_set = data->clusters;
//LNS starts
for (int iter = 0; iter < 2000; iter++)
{
//Delete a random vertex
int delete_vertex = rand() % (cluster_count - 1) + 1;
int prev_vertex = tour[delete_vertex].prev;
int next_vertex = tour[delete_vertex].next;
tour[prev_vertex].next = next_vertex;
tour[next_vertex].prev = prev_vertex;
int cluster_to_insert = clusters[tour[delete_vertex].vertex];
int best_pose = -1;
int best_vertex = -1;
int min_cost = INT_MAX;
for (int i = 0; i < vertex_set[cluster_to_insert].size; i++)
{
int vertex_to_insert = vertex_set[cluster_to_insert].nodes[i];
int next_edge = tour[0].next;
for (int j = 1; j < cluster_count; j++)
{
int vertex1 = tour[next_edge].vertex;
int vertex2 = tour[tour[next_edge].next].vertex;
int insert_cost = dist_matrix[vertex1][vertex_to_insert] +
dist_matrix[vertex_to_insert][vertex2] -
dist_matrix[vertex1][vertex2];
if (insert_cost < min_cost)
{
min_cost = insert_cost;
best_pose = next_edge;
best_vertex = vertex_to_insert;
}
next_edge = tour[next_edge].next;
}
}
assert(best_pose >= 0);
assert(best_vertex >= 0);
next_vertex = tour[best_pose].next;
tour[delete_vertex].prev = best_pose;
tour[delete_vertex].vertex = best_vertex;
tour[delete_vertex].next = next_vertex;
tour[best_pose].next = delete_vertex;
tour[next_vertex].prev = delete_vertex;
rval = optimize_vertex_in_cluster(tour, data);
abort_if(rval, "optimize_vertex_in_cluster failed");
}
rval = K_opt(tour, data);
abort_if(rval, "two_opt failed");
rval = two_opt(tour, data);
abort_if(rval, "two_opt failed");
*tour_cost = list_length(tour, data);
CLEANUP:
//if (vertex_seq) free(vertex_seq);
return rval;
}
int tour_length(int *tour, struct GTSP *data)
{
int tour_cost = 0;
for (int i = 0; i < data->cluster_count; i++)
{
if (i == data->cluster_count - 1)
tour_cost += data->dist_matrix[tour[i]][tour[0]];
else
tour_cost += data->dist_matrix[tour[i]][tour[i + 1]];
}
return tour_cost;
}
int list_length(struct Tour *tour, struct GTSP *data)
{
int tour_cost = 0;
for (int i = 0; i < data->cluster_count; i++)
{
int vertex1 = tour[i].vertex;
int vertex2 = tour[tour[i].next].vertex;
tour_cost += data->dist_matrix[vertex1][vertex2];
}
return tour_cost;
}
void print_tour(int *tour, struct GTSP *data)
{
for (int i = 0; i < data->cluster_count; i++)
{
printf("%d\t", tour[i]);
}
printf("\n");
}
void print_list(struct Tour *tour, struct GTSP *data)
{
printf("%d\t", tour[0].vertex);
int vertex_next = tour[0].next;
for (int i = 1; i < data->cluster_count; i++)
{
printf("%d\t", tour[vertex_next].vertex);
vertex_next = tour[vertex_next].next;
}
printf("\n");
}
int K_opt(struct Tour *tour, struct GTSP *data)
{
int cluster_count = data->cluster_count;
int l;
for (int i = 0; i < cluster_count - 3; i++)
{
int location_in_path = 0;
for (int k = 0; k < i; k++)
location_in_path = tour[location_in_path].next;
int current_vertex = tour[location_in_path].vertex;
for (int k = 3; k < cluster_count - i - 2; k++)
{
int first_next_location = tour[location_in_path].next;
int first_next_vertex = tour[tour[location_in_path].next].vertex;
int next_location = tour[location_in_path].next;
for (l = 0; l < k; l++)
{
next_location = tour[next_location].next;
}
int next_vertex = tour[next_location].vertex;
int next_next_location = tour[next_location].next;
int next_next_vertex = tour[next_next_location].vertex;
if (next_next_vertex == current_vertex) break;
int cost1 = data->dist_matrix[current_vertex][first_next_vertex] +
data->dist_matrix[next_vertex][next_next_vertex];
int cost2 = data->dist_matrix[current_vertex][next_vertex] +
data->dist_matrix[first_next_vertex][next_next_vertex];
if (cost2 < cost1)
{
int tmp_location = next_location;
for (int m = 0; m < l + 1; m++)
{
int tmp_vertex = tour[tmp_location].next;
tour[tmp_location].next = tour[tmp_location].prev;
tour[tmp_location].prev = tmp_vertex;
tmp_location = tour[tmp_location].next;
}
tour[location_in_path].next = next_location;
tour[next_location].prev = location_in_path;
tour[first_next_location].next = next_next_location;
tour[next_next_location].prev = first_next_location;
}
}
}
return 0;
}
int build_tour_from_x(struct GTSP *data, struct Tour *tour, double *x)
{
int rval = 0;
int *cluster_mark = 0;
struct Node **stack = 0;
int stack_top = 0;
struct Graph *graph = data->graph;
const int node_count = graph->node_count;
const int edge_count = graph->edge_count;
cluster_mark = (int *) malloc(data->cluster_count * sizeof(int));
abort_if(!cluster_mark, "could not allocate cluster_mark");
stack = (struct Node **) malloc(graph->node_count * sizeof(struct Node *));
abort_if(!stack, "could not allocate stack");
int k = 0;
for (int i = 0; i < node_count; i++) for (int i = 0; i < node_count; i++)
graph->nodes[i].mark = 0;
for (int i = 0; i < data->cluster_count; i++)
cluster_mark[i] = 0;
int initial;
for (initial = 0; initial < edge_count; initial++)
if (x[initial] > 1.0 - EPSILON) break;
initial = graph->edges[initial].from->index;
abort_if(initial == edge_count, "no initial node");
stack[stack_top++] = &graph->nodes[initial];
graph->nodes[initial].mark = 1;
tour[0].vertex = graph->nodes[initial].index;
int next_vertex = 1;
while (stack_top > 0)
{ {
struct Node *n = stack[--stack_top]; for (int j = i + 1; j < node_count; j++)
cluster_mark[data->node_to_cluster[n->index]]++;
for (int i = 0; i < n->degree; i++)
{
struct Adjacency *adj = &n->adj[i];
struct Node *neighbor = adj->neighbor;
if (neighbor->mark) continue;
if (x[adj->edge->index] < EPSILON) continue;
stack[stack_top++] = neighbor;
tour[next_vertex].vertex = neighbor->index;
next_vertex++;
neighbor->mark = 1;
}
}
for (int i = 0; i < data->cluster_count; i++)
{
if (i == 0)
{
tour[i].prev = data->cluster_count - 1;
}
else
{
tour[i].prev = i - 1;
}
if (i == data->cluster_count - 1)
{
tour[i].next = 0;
}
else
{ {
tour[i].next = i + 1; if (gtsp->node_to_cluster[i] == gtsp->node_to_cluster[j]) continue;
edge_map[i * node_count + j] = k;
edge_map[j * node_count + i] = k;
k++;
} }
} }
CLEANUP: return 0;
if (stack) free(stack);
if (cluster_mark) free(cluster_mark);
return rval;
} }

@ -31,26 +31,36 @@ struct Tour
struct Cluster struct Cluster
{ {
int size; /*
int* nodes; * Number of nodes inside the cluster.
*/
int size;
/*
* List of nodes inside the cluster.
*/
int *nodes;
}; };
struct GTSP struct GTSP
{ {
struct Graph *graph; struct Graph *graph;
int** dist_matrix; int **dist_matrix;
int cluster_count; /*
* Mapping between a node and the cluster which contains it. If a node
* has index i and is contained in cluster k, then node_to_cluster[i] = k.
*/
int *node_to_cluster; int *node_to_cluster;
int cluster_count;
struct Cluster *clusters; struct Cluster *clusters;
}; };
int GTSP_create_random_problem( int GTSP_create_random_problem(
int node_count, int cluster_count, int grid_size, struct GTSP *data); int node_count, int cluster_count, int grid_size, struct GTSP *data);
int inital_tour_value(struct GTSP *data, int *value, double *x);
void GTSP_free(struct GTSP *data); void GTSP_free(struct GTSP *data);
int GTSP_init_data(struct GTSP *data); int GTSP_init_data(struct GTSP *data);
@ -63,30 +73,15 @@ 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 optimize_vertex_in_cluster(struct Tour * tour, struct GTSP *data);
int two_opt(struct Tour * tour, struct GTSP *data);
int K_opt(struct Tour* tour, struct GTSP *data);
int tour_length(int* tour, struct GTSP* data);
void print_tour(int* tour, struct GTSP* data);
int list_length(struct Tour * tour, struct GTSP* data);
void print_list(struct Tour * tour, struct GTSP* data);
int GTSP_print_solution(struct GTSP *data, double *x); int GTSP_print_solution(struct GTSP *data, 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_solution_found(struct BNC *bnc, struct GTSP *data, double *x);
int GTSP_check_solution(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); int GTSP_read_solution(struct GTSP *gtsp, char *filename, double **p_x);
int GTSP_build_edge_map(struct GTSP *gtsp, int *edge_map);
#endif //_PROJECT_GTSP_H_ #endif //_PROJECT_GTSP_H_

@ -21,18 +21,19 @@
#include <assert.h> #include <assert.h>
#include "lp.h" #include "lp.h"
#include "util.h" #include "util.h"
#include "main.h" #include "gtsp-subtour.h"
static int compress_cut_pool(struct LP *lp) static int compress_cut_pool(struct LP *lp)
{ {
int delete_count = 0; int delete_count = 0;
for (int i = 0; i < lp->cut_pool_size; i++) for (int i = 0; i < lp->cut_pool_size; i++)
{ {
if (lp->cut_pool[i]->cplex_row_index < 0) struct Row *cut = lp->cut_pool[i];
if (cut->cplex_row_index < 0)
{ {
free(lp->cut_pool[i]->rmatind); free(cut->edges);
free(lp->cut_pool[i]->rmatval); free(cut);
free(lp->cut_pool[i]);
delete_count++; delete_count++;
} }
else else
@ -87,8 +88,10 @@ static int remove_old_cuts(struct LP *lp)
{ {
struct Row *cut = lp->cut_pool[j]; struct Row *cut = lp->cut_pool[j];
if (cut->cplex_row_index == i - count) cut->cplex_row_index = -1; if (cut->cplex_row_index == i - count)
else if (cut->cplex_row_index > i - count) cut->cplex_row_index--; cut->cplex_row_index = -1;
else if (cut->cplex_row_index > i - count)
cut->cplex_row_index--;
} }
count++; count++;
@ -132,20 +135,23 @@ static int remove_old_cuts(struct LP *lp)
log_debug("Compressing cut pool...\n"); log_debug("Compressing cut pool...\n");
compress_cut_pool(lp); compress_cut_pool(lp);
long nz = 0; long nz = 0;
long size = 0; long size = 0;
for (int i = 0; i < lp->cut_pool_size; i++) for (int i = 0; i < lp->cut_pool_size; i++)
{ {
size += sizeof(struct Row); size += sizeof(struct Row);
nz += lp->cut_pool[i]->nz; struct Row *cut = lp->cut_pool[i];
size += lp->cut_pool[i]->nz * sizeof(double);
size += lp->cut_pool[i]->nz * sizeof(int); nz += cut->nz;
size += cut->nz * sizeof(double);
size += cut->nz * sizeof(int);
size += cut->edge_count * sizeof(char);
} }
log_debug(" %ld cuts (%ld nz, %ld MiB)\n", lp->cut_pool_size, nz, size/1024/1024); log_debug(" %ld cuts (%ld nz, %ld MiB)\n", lp->cut_pool_size, nz,
size / 1024 / 1024);
if(size > CUT_POOL_MAX_MEMORY) if (size > CUT_POOL_MAX_MEMORY)
CUT_POOL_MAX_MEMORY = size; CUT_POOL_MAX_MEMORY = size;
CLEANUP: CLEANUP:
@ -167,7 +173,6 @@ static int update_cut_ages(struct LP *lp)
rval = CPXgetslack(lp->cplex_env, lp->cplex_lp, slacks, 0, numrows - 1); rval = CPXgetslack(lp->cplex_env, lp->cplex_lp, slacks, 0, numrows - 1);
abort_if(rval, "CPXgetslack failed"); abort_if(rval, "CPXgetslack failed");
int count = 0;
for (int i = 0; i < lp->cut_pool_size; i++) for (int i = 0; i < lp->cut_pool_size; i++)
{ {
struct Row *cut = lp->cut_pool[i]; struct Row *cut = lp->cut_pool[i];
@ -178,8 +183,6 @@ static int update_cut_ages(struct LP *lp)
cut->age++; cut->age++;
else else
cut->age = 0; cut->age = 0;
if (cut->age > 5) count++;
} }
CLEANUP: CLEANUP:
@ -208,9 +211,9 @@ void LP_free_cut_pool(struct LP *lp)
{ {
for (int i = 0; i < lp->cut_pool_size; i++) for (int i = 0; i < lp->cut_pool_size; i++)
{ {
free(lp->cut_pool[i]->rmatind); struct Row *cut = lp->cut_pool[i];
free(lp->cut_pool[i]->rmatval); free(cut->edges);
free(lp->cut_pool[i]); free(cut);
} }
if (lp->cut_pool) free(lp->cut_pool); if (lp->cut_pool) free(lp->cut_pool);
@ -322,8 +325,8 @@ int LP_optimize(struct LP *lp, int *infeasible)
int numrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp); int numrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp);
int numcols = CPXgetnumcols(lp->cplex_env, lp->cplex_lp); int numcols = CPXgetnumcols(lp->cplex_env, lp->cplex_lp);
if(numrows > LP_MAX_ROWS) LP_MAX_ROWS = numrows; if (numrows > LP_MAX_ROWS) LP_MAX_ROWS = numrows;
if(numrows > LP_MAX_COLS) LP_MAX_COLS = numcols; if (numrows > LP_MAX_COLS) LP_MAX_COLS = numcols;
log_debug("Optimizing LP (%d rows %d cols)...\n", numrows, numcols); log_debug("Optimizing LP (%d rows %d cols)...\n", numrows, numcols);
@ -337,6 +340,8 @@ int LP_optimize(struct LP *lp, int *infeasible)
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)
{ {
log_debug(" infeasible\n");
*infeasible = 1; *infeasible = 1;
goto CLEANUP; goto CLEANUP;
} }
@ -357,9 +362,12 @@ int LP_optimize(struct LP *lp, int *infeasible)
rval = update_cut_ages(lp); rval = update_cut_ages(lp);
abort_if(rval, "update_cut_ages failed"); abort_if(rval, "update_cut_ages failed");
log_debug("Removing old cuts...\n"); if (LP_SOLVE_COUNT % MAX_CUT_AGE == 0)
rval = remove_old_cuts(lp); {
abort_if(rval, "LP_remove_old_cuts failed"); log_debug("Removing old cuts...\n");
rval = remove_old_cuts(lp);
abort_if(rval, "LP_remove_old_cuts failed");
}
CUT_POOL_TIME += get_user_time() - initial_time; CUT_POOL_TIME += get_user_time() - initial_time;
@ -392,11 +400,30 @@ int LP_get_x(struct LP *lp, double *x)
return rval; return rval;
} }
int LP_get_y(struct LP *lp, double *y)
{
int rval = 0;
int nrows = CPXgetnumrows(lp->cplex_env, lp->cplex_lp);
abort_if(!nrows, "No rows in LP");
rval = CPXgetpi(lp->cplex_env, lp->cplex_lp, y, 0, nrows - 1);
abort_iff(rval, "CPXgetpi failed (errno = %d)", rval);
CLEANUP:
return rval;
}
int LP_get_num_cols(struct LP *lp) int LP_get_num_cols(struct LP *lp)
{ {
return CPXgetnumcols(lp->cplex_env, lp->cplex_lp); return CPXgetnumcols(lp->cplex_env, lp->cplex_lp);
} }
int LP_get_num_rows(struct LP *lp)
{
return CPXgetnumrows(lp->cplex_env, lp->cplex_lp);
}
int LP_write(struct LP *lp, const char *fname) int LP_write(struct LP *lp, const char *fname)
{ {
int rval = 0; int rval = 0;
@ -427,12 +454,10 @@ int LP_write(struct LP *lp, const char *fname)
int compare_cuts(struct Row *cut1, struct Row *cut2) int compare_cuts(struct Row *cut1, struct Row *cut2)
{ {
return_if_neq(cut1->nz, cut2->nz);
for (int i = 0; i < cut1->nz; i++) for (int i = 0; i < cut1->edge_count; i++)
{ {
return_if_neq(cut1->rmatind[i], cut2->rmatind[i]); return_if_neq(cut1->edges[i], cut2->edges[i]);
return_if_neq_epsilon(cut1->rmatval[i], cut2->rmatval[i]);
} }
return 0; return 0;
@ -442,9 +467,9 @@ static int update_hash(struct Row *cut)
{ {
unsigned long hash = 0; unsigned long hash = 0;
for (int i = 0; i < cut->nz; i++) for (int i = 0; i < cut->edge_count; i++)
{ {
hash += cut->rmatind[i] * 65521; hash += cut->edges[i] * 65521;
hash %= 4294967291; hash %= 4294967291;
} }
@ -472,6 +497,7 @@ int LP_add_cut(struct LP *lp, struct Row *cut)
free(cut->rmatval); free(cut->rmatval);
free(cut->rmatind); free(cut->rmatind);
free(cut->edges);
free(cut); free(cut);
return 0; return 0;
} }
@ -488,6 +514,11 @@ 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;
free(cut->rmatval);
free(cut->rmatind);
cut->rmatind = 0;
cut->rmatval = 0;
CUT_POOL_TIME += get_user_time() - initial_time; CUT_POOL_TIME += get_user_time() - initial_time;
CLEANUP: CLEANUP:
return rval; return rval;

@ -40,6 +40,10 @@ struct Row
double rhs; double rhs;
double *rmatval; double *rmatval;
int *rmatind; int *rmatind;
int edge_count;
double edge_val;
char *edges;
}; };
static const int MAX_NAME_LENGTH = 100; static const int MAX_NAME_LENGTH = 100;
@ -83,8 +87,12 @@ int LP_get_obj_val(struct LP *lp, double *obj);
int LP_get_x(struct LP *lp, double *x); int LP_get_x(struct LP *lp, double *x);
int LP_get_y(struct LP *lp, double *y);
int LP_get_num_cols(struct LP *lp); int LP_get_num_cols(struct LP *lp);
int LP_get_num_rows(struct LP *lp);
int LP_add_cut(struct LP *lp, struct Row *cut); int LP_add_cut(struct LP *lp, struct Row *cut);
#endif #endif

@ -21,11 +21,11 @@
#include "main.h" #include "main.h"
#include "gtsp.h" #include "gtsp.h"
#include "util.h" #include "util.h"
#include "gtsp-heur.h"
unsigned int SEED = 0;
double SUBTOUR_TIME = 0; double SUBTOUR_TIME = 0;
double COMBS_TIME = 0; double COMBS_TIME = 0;
double COLUMNS_TIME = 0;
double CUT_POOL_TIME = 0; double CUT_POOL_TIME = 0;
long CUT_POOL_MAX_MEMORY = 0; long CUT_POOL_MAX_MEMORY = 0;
@ -52,19 +52,24 @@ char FRAC_SOLUTION_FILENAME[100] = {0};
char WRITE_GRAPH_FILENAME[100] = {0}; char WRITE_GRAPH_FILENAME[100] = {0};
char STATS_FILENAME[100] = {0}; char STATS_FILENAME[100] = {0};
double *OPTIMAL_X = 0;
unsigned int SEED = 0;
static int input_node_count = -1; static int input_node_count = -1;
static int input_cluster_count = -1; static int input_cluster_count = -1;
static int grid_size = 100; static int grid_size = 100;
static const struct option options_tab[] = {{"help", no_argument, 0, 'h'}, static const struct option options_tab[] = {{"help", no_argument, 0, 'h'},
{"nodes", required_argument, 0, 'n'}, {"nodes", required_argument, 0, 'n'},
{"clusters", required_argument, 0, 'm'}, {"clusters", required_argument, 0, 'm'},
{"grid-size", required_argument, 0, 'g'}, {"grid-size", required_argument, 0, 'g'},
{"optimal", required_argument, 0, 'x'}, {"optimal", required_argument, 0, 'x'},
{"seed", required_argument, 0, 's'}, {"out", required_argument, 0, 'o'}, {"seed", required_argument, 0, 's'},
{"stats", required_argument, 0, 't'}, {"lp", required_argument, 0, 'l'}, {"out", required_argument, 0, 'o'},
{"write-graph", required_argument, 0, 'w'}, {"stats", required_argument, 0, 't'},
{(char *) 0, (int) 0, (int *) 0, (int) 0}}; {"lp", required_argument, 0, 'l'},
{"write-graph", required_argument, 0, 'w'},
{(char *) 0, (int) 0, (int *) 0, (int) 0}};
static char input_x_filename[1000] = {0}; static char input_x_filename[1000] = {0};
@ -103,7 +108,7 @@ static int parse_args(int argc, char **argv)
{ {
int c = 0; int c = 0;
int option_index = 0; int option_index = 0;
c = getopt_long(argc, argv, "n:m:g:x:s:h:o:t:l:w:", options_tab, c = getopt_long(argc, argv, "n:m:g:x:s:h:o:t:l:w:f:", options_tab,
&option_index); &option_index);
if (c < 0) break; if (c < 0) break;
@ -146,6 +151,10 @@ static int parse_args(int argc, char **argv)
strcpy(WRITE_GRAPH_FILENAME, optarg); strcpy(WRITE_GRAPH_FILENAME, optarg);
break; break;
case 'f':
strcpy(FRAC_SOLUTION_FILENAME, optarg);
break;
case 'h': case 'h':
print_usage(argv); print_usage(argv);
exit(0); exit(0);
@ -235,11 +244,10 @@ int main(int argc, char **argv)
int init_val = 0; int init_val = 0;
initial_x = (double *) malloc( initial_x = (double *) malloc((data.graph->edge_count) * sizeof(double));
(data.graph->node_count + data.graph->edge_count) * sizeof(double));
abort_if(!initial_x, "could not allocate initial_x"); abort_if(!initial_x, "could not allocate initial_x");
rval = inital_tour_value(&data, &init_val, initial_x); rval = GTSP_find_initial_tour(&data, &init_val, initial_x);
abort_if(rval, "initial_tour_value failed"); abort_if(rval, "initial_tour_value failed");
rval = GTSP_solution_found(&bnc, &data, initial_x); rval = GTSP_solution_found(&bnc, &data, initial_x);
@ -289,7 +297,7 @@ int main(int argc, char **argv)
log_info("Starting branch-and-cut solver...\n"); log_info("Starting branch-and-cut solver...\n");
rval = BNC_solve(&bnc); rval = BNC_solve(&bnc);
abort_if(rval, "BNC_solve_node failed"); abort_if(rval, "solve_node failed");
abort_if(!bnc.best_x, "problem has no feasible solution"); abort_if(!bnc.best_x, "problem has no feasible solution");
@ -317,6 +325,7 @@ int main(int argc, char **argv)
fprintf(file, "%-8s ", "time"); fprintf(file, "%-8s ", "time");
fprintf(file, "%-8s ", "subt-t"); fprintf(file, "%-8s ", "subt-t");
fprintf(file, "%-8s ", "combs-t"); fprintf(file, "%-8s ", "combs-t");
fprintf(file, "%-8s ", "cols-t");
fprintf(file, "%-8s ", "pool-t"); fprintf(file, "%-8s ", "pool-t");
fprintf(file, "%-8s ", "pool-m"); fprintf(file, "%-8s ", "pool-m");
fprintf(file, "%-8s ", "lp-count"); fprintf(file, "%-8s ", "lp-count");
@ -337,6 +346,7 @@ int main(int argc, char **argv)
fprintf(file, "%-8.2lf ", TOTAL_TIME); fprintf(file, "%-8.2lf ", TOTAL_TIME);
fprintf(file, "%-8.2lf ", SUBTOUR_TIME); fprintf(file, "%-8.2lf ", SUBTOUR_TIME);
fprintf(file, "%-8.2lf ", COMBS_TIME); fprintf(file, "%-8.2lf ", COMBS_TIME);
fprintf(file, "%-8.2lf ", COLUMNS_TIME);
fprintf(file, "%-8.2lf ", CUT_POOL_TIME); fprintf(file, "%-8.2lf ", CUT_POOL_TIME);
fprintf(file, "%-8ld ", CUT_POOL_MAX_MEMORY / 1024 / 1024); fprintf(file, "%-8ld ", CUT_POOL_MAX_MEMORY / 1024 / 1024);
fprintf(file, "%-8d ", LP_SOLVE_COUNT); fprintf(file, "%-8d ", LP_SOLVE_COUNT);
@ -353,7 +363,7 @@ int main(int argc, char **argv)
fclose(file); fclose(file);
} }
if(strlen(SOLUTION_FILENAME) == 0) if (strlen(SOLUTION_FILENAME) == 0)
{ {
log_info("Optimal solution:\n"); log_info("Optimal solution:\n");
rval = GTSP_print_solution(&data, bnc.best_x); rval = GTSP_print_solution(&data, bnc.best_x);
@ -363,7 +373,6 @@ int main(int argc, char **argv)
log_info("Optimal solution value:\n"); log_info("Optimal solution value:\n");
log_info(" %.4lf\n", bnc.best_obj_val); log_info(" %.4lf\n", bnc.best_obj_val);
CLEANUP: CLEANUP:
if (OPTIMAL_X) free(OPTIMAL_X); if (OPTIMAL_X) free(OPTIMAL_X);
GTSP_free(&data); GTSP_free(&data);

@ -22,6 +22,7 @@ extern unsigned int SEED;
extern double *OPTIMAL_X; extern double *OPTIMAL_X;
extern double SUBTOUR_TIME; extern double SUBTOUR_TIME;
extern double COMBS_TIME; extern double COMBS_TIME;
extern double COLUMNS_TIME;
extern double LP_SOLVE_TIME; extern double LP_SOLVE_TIME;
extern int LP_SOLVE_COUNT; extern int LP_SOLVE_COUNT;
@ -44,4 +45,6 @@ extern char LP_FILENAME[100];
extern char SOLUTION_FILENAME[100]; extern char SOLUTION_FILENAME[100];
extern char FRAC_SOLUTION_FILENAME[100]; extern char FRAC_SOLUTION_FILENAME[100];
extern int BNC_NODE_COUNT;
#endif #endif

@ -42,6 +42,6 @@
/* /*
* Time limit for the computation (user time, in seconds). * Time limit for the computation (user time, in seconds).
*/ */
#define MAX_TOTAL_TIME 3600 #define MAX_TOTAL_TIME (24*3600)
#endif //PROJECT_PARAMS_H #endif //PROJECT_PARAMS_H