construct.c

Lin-Kernighan heuristic for the TSP and minimum weight perfect matching

#define init_max_per_city_nn_cache 30
#define max(X,Y) ((X) <(Y) ?(Y) :(X) ) 
#define min(X,Y) ((X) >(Y) ?(Y) :(X) )  \

/*1:*/
#line 213 "./construct.w"

#include <config.h>
#include "lkconfig.h"
/*25:*/
#line 548 "./construct.w"

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#define FIXINCLUDES_NEED_NRAND48 
#include "fixincludes.h"
#undef FIXINCLUDES_NEED_NRAND48

/*:25*/
#line 216 "./construct.w"

/*3:*/
#line 237 "./construct.w"

#include "construct.h"

/*:3*//*7:*/
#line 281 "./construct.w"

#include "error.h"

/*:7*//*12:*/
#line 346 "./construct.w"

#include "prng.h"
#include "read.h"


/*:12*//*24:*/
#line 540 "./construct.w"

#include "memory.h"
#include "pq.h"
#include "pool.h"

/*:24*//*28:*/
#line 612 "./construct.w"

#include "lk.h"
#include "kdtree.h"

/*:28*//*42:*/
#line 801 "./construct.w"

#include "dsort.h"


/*:42*/
#line 217 "./construct.w"

/*19:*/
#line 483 "./construct.w"

typedef struct{
length_t len;
int this_end;
int other_end;
}pq_edge_t;

/*:19*//*45:*/
#line 822 "./construct.w"

typedef struct{
int count;
int city[2];
}adj_entry_t;

/*:45*/
#line 218 "./construct.w"

/*20:*/
#line 491 "./construct.w"

static void pq_edge_print(FILE*out,void*edgep);
static void
pq_edge_print(FILE*out,void*edgep)
{
pq_edge_t*edge= edgep;
if(edge==NULL){fprintf(out,"(null)");}
else{
fprintf(out,"{%p,%d,%d,"length_t_spec"}",
edge,edge->this_end,edge->other_end,length_t_pcast(edge->len));
}
}

/*:20*//*26:*/
#line 574 "./construct.w"

static int
cmp_pq_edge(const void*a,const void*b){
length_t len_diff= ((const pq_edge_t*)a)->len-((const pq_edge_t*)b)->len;
return len_diff<0?-1:(len_diff>0?1:
(int)(((const pq_edge_t*)a)-((const pq_edge_t*)b)));
}

/*:26*/
#line 219 "./construct.w"

/*4:*/
#line 256 "./construct.w"

length_t
construct(const int n,int*tour,const int heuristic,const long heur_param,
const long random_seed)
{
switch(heuristic){
/*8:*/
#line 288 "./construct.w"

case CONSTRUCT_CANONICAL:{
int i;length_t len;
for(i= 0;i<n;i++)tour[i]= i;

len= cost(tour[0],tour[n-1]);
for(i= 1;i<n;i++)len+= cost(tour[i-1],tour[i]);
return len;
break;
}

/*:8*//*10:*/
#line 312 "./construct.w"

case CONSTRUCT_RANDOM:{
int i;length_t len;
prng_t*random_stream;

for(i= 0;i<n;i++){
tour[i]= i;
}

random_stream= prng_new(PRNG_DEFAULT,heur_param);
for(i= 0;i<n;i++){
const int next= prng_unif_int(random_stream,n-i);
const int t= tour[next];
tour[next]= tour[n-1-i];
tour[n-1-i]= t;
}
prng_free(random_stream);

len= cost(tour[0],tour[n-1]);
for(i= 1;i<n;i++)len+= cost(tour[i-1],tour[i]);

return len;
break;
}


/*:10*//*16:*/
#line 424 "./construct.w"

case CONSTRUCT_GREEDY:
{
length_t len= 0;
const int E2_case= E2_supports(tsp_instance);
#define DO_TOUR
/*21:*/
#line 523 "./construct.w"

pool_t*nn_link_pool= pool_create(sizeof(pq_edge_t),n);
pq_t*pq_edge= pq_create_size(cmp_pq_edge,n);

/*:21*//*29:*/
#line 623 "./construct.w"

int*farthest_in_queue= NULL;

/*:29*//*32:*/
#line 644 "./construct.w"

int*unsaturated= NULL,num_unsaturated= 0,*inv_unsaturated= NULL;

/*:32*//*37:*/
#line 699 "./construct.w"

int sqrt_n= (int)sqrt((double)n);

/*:37*//*43:*/
#line 806 "./construct.w"

pq_edge_t*nn_work= NULL;


/*:43*//*46:*/
#line 830 "./construct.w"

#if defined(DO_TOUR)
adj_entry_t*adj= new_arr_of(adj_entry_t,n);
#endif

/*:46*//*49:*/
#line 859 "./construct.w"

#if defined(DO_TOUR)
int*tail= new_arr_of(int,n);
#endif

/*:49*//*55:*/
#line 957 "./construct.w"

#if defined(DO_RANDOM)
prng_t*random_stream= NULL;
#endif

/*:55*/
#line 430 "./construct.w"

/*22:*/
#line 528 "./construct.w"

errorif(pq_edge==NULL,"Couldn't create the priority queue!");
pq_set_print_func(pq_edge,pq_edge_print);

/*:22*//*27:*/
#line 596 "./construct.w"

if(E2_case){
int i;
for(i= 0;i<n;i++){
pq_edge_t*e= pool_alloc(nn_link_pool);
e->this_end= i;
e->other_end= E2_nn(i);
e->len= cost(i,e->other_end);
pq_insert(pq_edge,e);
/*48:*/
#line 842 "./construct.w"

#if defined(DO_TOUR)
adj[i].count= 0;
#endif

/*:48*//*51:*/
#line 871 "./construct.w"

#if defined(DO_TOUR)
tail[i]= i;
#endif

/*:51*//*68:*/
#line 1197 "./construct.w"

#if defined(DO_MATCHING)
mate[i]= -1;
#endif


/*:68*/
#line 605 "./construct.w"

}
}else{
/*30:*/
#line 627 "./construct.w"

farthest_in_queue= new_arr_of(int,n);

/*:30*//*33:*/
#line 648 "./construct.w"

unsaturated= new_arr_of(int,n);
inv_unsaturated= new_arr_of(int,n);
num_unsaturated= n;
{int i;
for(i= 0;i<n;i++)inv_unsaturated[i]= unsaturated[i]= i;
}

/*:33*//*36:*/
#line 687 "./construct.w"

{
int i;
nn_work= new_arr_of(pq_edge_t,n);
for(i= 0;i<n;i++){
const int x= i,not_me= -1;
/*38:*/
#line 732 "./construct.w"

if(num_unsaturated>sqrt_n){
/*41:*/
#line 767 "./construct.w"

{
int i,num_chosen,farthest_city;
size_t w;
length_t farthest_len;
/*81:*/
#line 1407 "./construct.w"

#if defined(CONSTRUCT_MAX_VERBOSE)
#if CONSTRUCT_MAX_VERBOSE >= 750
fprintf(stderr,"construct: Select fresh neighbours for %d, excepting %d\n",
x,not_me);fflush(stderr);
#endif
#endif

/*:81*/
#line 772 "./construct.w"

for(i= w= 0;i<num_unsaturated;i++){
const int c= unsaturated[i];
if(c==x||c==not_me)continue;
nn_work[w].this_end= x;
nn_work[w].other_end= c;
nn_work[w].len= cost(x,c);
w++;
}
num_chosen= min(w,init_max_per_city_nn_cache);
errorif(num_chosen==0,"Bug!");
(void)select_range(nn_work,w,sizeof(pq_edge_t),cmp_pq_edge,
0,num_chosen,0);
farthest_len= nn_work[num_chosen-1].len;
farthest_city= nn_work[num_chosen-1].other_end;
for(i= 0;i<num_chosen;i++){
pq_edge_t*e;
e= pool_alloc(nn_link_pool);
*e= nn_work[i];
pq_insert(pq_edge,e);
if(farthest_len<e->len){
farthest_city= e->other_end;
farthest_len= e->len;
}
}
farthest_in_queue[x]= farthest_city;
}

/*:41*/
#line 734 "./construct.w"
;
}else{
/*39:*/
#line 740 "./construct.w"

{int i;
/*78:*/
#line 1378 "./construct.w"

#if defined(CONSTRUCT_MAX_VERBOSE)
#if CONSTRUCT_MAX_VERBOSE >= 500
fprintf(stderr,"construct: All neighbours for %d, excepting %d\n",
x,not_me);fflush(stderr);
#endif
#endif


/*:78*/
#line 742 "./construct.w"

for(i= 0;i<num_unsaturated;i++){
const int y= unsaturated[i];
length_t len;
if(y==x||y==not_me)continue;
len= cost(x,y);
/*40:*/
#line 754 "./construct.w"

{
pq_edge_t*e;
e= pool_alloc(nn_link_pool);
e->len= len;
e->this_end= x;
e->other_end= y;
pq_insert(pq_edge,e);
/*80:*/
#line 1397 "./construct.w"

#if defined(CONSTRUCT_MAX_VERBOSE)
#if CONSTRUCT_MAX_VERBOSE >= 750
fprintf(stderr,"construct: alloc new link (%d,%d) "
length_t_native_spec"\n",
e->this_end,e->other_end,length_t_native_pcast(e->len));fflush(stderr);
#endif
#endif

/*:80*/
#line 762 "./construct.w"

}


/*:40*/
#line 748 "./construct.w"

}
}


/*:39*/
#line 736 "./construct.w"

}

/*:38*/
#line 693 "./construct.w"

/*48:*/
#line 842 "./construct.w"

#if defined(DO_TOUR)
adj[i].count= 0;
#endif

/*:48*//*51:*/
#line 871 "./construct.w"

#if defined(DO_TOUR)
tail[i]= i;
#endif

/*:51*//*68:*/
#line 1197 "./construct.w"

#if defined(DO_MATCHING)
mate[i]= -1;
#endif


/*:68*/
#line 694 "./construct.w"

}
}

/*:36*/
#line 608 "./construct.w"

}

/*:27*//*56:*/
#line 969 "./construct.w"

#if defined(DO_RANDOM)
random_stream= prng_new(PRNG_DEFAULT,random_seed^(6502*6510));
#endif

/*:56*/
#line 431 "./construct.w"

/*52:*/
#line 887 "./construct.w"

{int i,x,y,tx,ty;
errorif(n<3,"Only %d cities.  Can't build a tour.",n);
tx= ty= -1;
for(i= 0;i<n-1;i++){
pq_edge_t*e;
/*53:*/
#line 925 "./construct.w"

#if defined(DO_RANDOM)
{
pq_edge_t*candidate[2];
/*58:*/
#line 993 "./construct.w"

/*75:*/
#line 1338 "./construct.w"

#if defined(CONSTRUCT_MAX_VERBOSE)
#if CONSTRUCT_MAX_VERBOSE >= 2000
if(verbose>=2000){
printf("Priority queue before Get one short valid edge or NULL:\n");
pq_print(pq_edge,stdout);
printf("\n");
}
#endif
#endif

/*:75*/
#line 994 "./construct.w"

while(1){
e= pq_delete_min(pq_edge);
if(e==NULL)break;
/*76:*/
#line 1350 "./construct.w"

#if defined(CONSTRUCT_MAX_VERBOSE)
#if CONSTRUCT_MAX_VERBOSE >= 2000
if(verbose>=2000){
printf("Extract ");
pq_edge_print(stdout,e);

printf("\nPriority queue just after extract:\n");
pq_print(pq_edge,stdout);
printf("\n");
}
#endif
#endif

/*:76*/
#line 998 "./construct.w"

x= e->this_end;y= e->other_end;
if(adj[x].count==2){
continue;
}
if(adj[y].count<2&&y!=tail[x])break;
/*60:*/
#line 1032 "./construct.w"

if(E2_case){
E2_hide(tail[x]);
e->other_end= E2_nn(x);
e->len= cost(x,e->other_end);
E2_unhide(tail[x]);
pq_insert(pq_edge,e);
}else{
pool_free(nn_link_pool,e);
if(farthest_in_queue[x]==y){
const int not_me= tail[x];
/*82:*/
#line 1416 "./construct.w"

#if defined(CONSTRUCT_MAX_VERBOSE)
#if CONSTRUCT_MAX_VERBOSE >= 750
fprintf(stderr,"Need to refresh because I got x=%d %d y=%d %d "
length_t_native_spec"\n",
x,e->this_end,y,e->other_end,length_t_native_pcast(e->len));
#endif
#endif

/*:82*/
#line 1043 "./construct.w"

/*38:*/
#line 732 "./construct.w"

if(num_unsaturated>sqrt_n){
/*41:*/
#line 767 "./construct.w"

{
int i,num_chosen,farthest_city;
size_t w;
length_t farthest_len;
/*81:*/
#line 1407 "./construct.w"

#if defined(CONSTRUCT_MAX_VERBOSE)
#if CONSTRUCT_MAX_VERBOSE >= 750
fprintf(stderr,"construct: Select fresh neighbours for %d, excepting %d\n",
x,not_me);fflush(stderr);
#endif
#endif

/*:81*/
#line 772 "./construct.w"

for(i= w= 0;i<num_unsaturated;i++){
const int c= unsaturated[i];
if(c==x||c==not_me)continue;
nn_work[w].this_end= x;
nn_work[w].other_end= c;
nn_work[w].len= cost(x,c);
w++;
}
num_chosen= min(w,init_max_per_city_nn_cache);
errorif(num_chosen==0,"Bug!");
(void)select_range(nn_work,w,sizeof(pq_edge_t),cmp_pq_edge,
0,num_chosen,0);
farthest_len= nn_work[num_chosen-1].len;
farthest_city= nn_work[num_chosen-1].other_end;
for(i= 0;i<num_chosen;i++){
pq_edge_t*e;
e= pool_alloc(nn_link_pool);
*e= nn_work[i];
pq_insert(pq_edge,e);
if(farthest_len<e->len){
farthest_city= e->other_end;
farthest_len= e->len;
}
}
farthest_in_queue[x]= farthest_city;
}

/*:41*/
#line 734 "./construct.w"
;
}else{
/*39:*/
#line 740 "./construct.w"

{int i;
/*78:*/
#line 1378 "./construct.w"

#if defined(CONSTRUCT_MAX_VERBOSE)
#if CONSTRUCT_MAX_VERBOSE >= 500
fprintf(stderr,"construct: All neighbours for %d, excepting %d\n",
x,not_me);fflush(stderr);
#endif
#endif


/*:78*/
#line 742 "./construct.w"

for(i= 0;i<num_unsaturated;i++){
const int y= unsaturated[i];
length_t len;
if(y==x||y==not_me)continue;
len= cost(x,y);
/*40:*/
#line 754 "./construct.w"

{
pq_edge_t*e;
e= pool_alloc(nn_link_pool);
e->len= len;
e->this_end= x;
e->other_end= y;
pq_insert(pq_edge,e);
/*80:*/
#line 1397 "./construct.w"

#if defined(CONSTRUCT_MAX_VERBOSE)
#if CONSTRUCT_MAX_VERBOSE >= 750
fprintf(stderr,"construct: alloc new link (%d,%d) "
length_t_native_spec"\n",
e->this_end,e->other_end,length_t_native_pcast(e->len));fflush(stderr);
#endif
#endif

/*:80*/
#line 762 "./construct.w"

}


/*:40*/
#line 748 "./construct.w"

}
}


/*:39*/
#line 736 "./construct.w"

}

/*:38*/
#line 1044 "./construct.w"

}
}

/*:60*/
#line 1004 "./construct.w"

}
/*77:*/
#line 1365 "./construct.w"

#if defined(CONSTRUCT_MAX_VERBOSE)
#if CONSTRUCT_MAX_VERBOSE >= 1000
if(verbose>=2000){
printf("Priority queue AFTER Get one short valid edge or NULL:\n");
pq_print(pq_edge,stdout);
printf("\n");
}
if(verbose>=1000)printf(" valid");
#endif
#endif

/*:77*/
#line 1006 "./construct.w"


/*:58*/
#line 929 "./construct.w"
candidate[0]= e;
/*58:*/
#line 993 "./construct.w"

/*75:*/
#line 1338 "./construct.w"

#if defined(CONSTRUCT_MAX_VERBOSE)
#if CONSTRUCT_MAX_VERBOSE >= 2000
if(verbose>=2000){
printf("Priority queue before Get one short valid edge or NULL:\n");
pq_print(pq_edge,stdout);
printf("\n");
}
#endif
#endif

/*:75*/
#line 994 "./construct.w"

while(1){
e= pq_delete_min(pq_edge);
if(e==NULL)break;
/*76:*/
#line 1350 "./construct.w"

#if defined(CONSTRUCT_MAX_VERBOSE)
#if CONSTRUCT_MAX_VERBOSE >= 2000
if(verbose>=2000){
printf("Extract ");
pq_edge_print(stdout,e);

printf("\nPriority queue just after extract:\n");
pq_print(pq_edge,stdout);
printf("\n");
}
#endif
#endif

/*:76*/
#line 998 "./construct.w"

x= e->this_end;y= e->other_end;
if(adj[x].count==2){
continue;
}
if(adj[y].count<2&&y!=tail[x])break;
/*60:*/
#line 1032 "./construct.w"

if(E2_case){
E2_hide(tail[x]);
e->other_end= E2_nn(x);
e->len= cost(x,e->other_end);
E2_unhide(tail[x]);
pq_insert(pq_edge,e);
}else{
pool_free(nn_link_pool,e);
if(farthest_in_queue[x]==y){
const int not_me= tail[x];
/*82:*/
#line 1416 "./construct.w"

#if defined(CONSTRUCT_MAX_VERBOSE)
#if CONSTRUCT_MAX_VERBOSE >= 750
fprintf(stderr,"Need to refresh because I got x=%d %d y=%d %d "
length_t_native_spec"\n",
x,e->this_end,y,e->other_end,length_t_native_pcast(e->len));
#endif
#endif

/*:82*/
#line 1043 "./construct.w"

/*38:*/
#line 732 "./construct.w"

if(num_unsaturated>sqrt_n){
/*41:*/
#line 767 "./construct.w"

{
int i,num_chosen,farthest_city;
size_t w;