decluster.c

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

#define child city
#define NO_CHILD (-1)  \
 \

#define NO_PARENT (-2)  \

#define lo_mask(X) (pow_2[(X) ]-1) 
#define hi_mask(X) (~lo_mask(i) )  \

#define min(A,B) ((A) <(B) ?(A) :(B) ) 
#define max(A,B) ((A) >(B) ?(A) :(B) )  \

/*9:*/
#line 530 "./decluster.w"

#include <config.h>
#include "lkconfig.h"
/*12:*/
#line 624 "./decluster.w"

#include <stdio.h>
#include <stdlib.h>
#include <stddef.h>

/*:12*/
#line 533 "./decluster.w"

/*15:*/
#line 651 "./decluster.w"

#include "error.h"
#include "memory.h"
#include "length.h"
#include "read.h"


/*:15*/
#line 534 "./decluster.w"

/*14:*/
#line 640 "./decluster.w"

#include "decluster.h"


/*:14*//*27:*/
#line 956 "./decluster.w"

#include "pq.h"
#include "kdtree.h"

/*:27*//*39:*/
#line 1128 "./decluster.w"

#include "lk.h"

/*:39*/
#line 535 "./decluster.w"


/*36:*/
#line 1088 "./decluster.w"

int decluster_discard_topology_tree= 1;

/*:36*/
#line 537 "./decluster.w"

/*41:*/
#line 1168 "./decluster.w"

#if SIZEOF_INT==8
# define copy_1_down(X) \
 (((X) |=  (X)>>1), \
  ((X) |=  (X)>>2), \
  ((X) |=  (X)>>4), \
  ((X) |=  (X)>>8), \
  ((X) |=  (X)>>16), \
  ((X) |=  (X)>>32) )
#else 
# define copy_1_down(X) \
 (((X) |=  (X)>>1), \
  ((X) |=  (X)>>2), \
  ((X) |=  (X)>>4), \
  ((X) |=  (X)>>8), \
  ((X) |=  (X)>>16) )
#endif


/*:41*//*75:*/
#line 1943 "./decluster.w"

#if DECLUSTER_DEBUG
#define print_tree(A,B) decluster_print_tree(stdout,A,B)
#else
#define print_tree(A,B)
#endif


/*:75*/
#line 538 "./decluster.w"

/*48:*/
#line 1416 "./decluster.w"

typedef struct{
int level,inlabel,ascendant;
length_t cost;
}digest_t;

/*:48*/
#line 539 "./decluster.w"

/*17:*/
#line 676 "./decluster.w"

static int n;

/*:17*//*34:*/
#line 1049 "./decluster.w"

static decluster_tree_t*T_prime= NULL;

/*:34*//*42:*/
#line 1190 "./decluster.w"


static const int pow_2[]= {
0x1,
0x2,
0x4,
0x8,
0x10,
0x20,
0x40,
0x80,
0x100,
0x200,
0x400,
0x800,
0x1000,
0x2000,
0x4000,
0x8000,
0x10000,
0x20000,
0x40000,
0x80000,
0x100000,
0x200000,
0x400000,
0x800000,
0x1000000,
0x2000000,
0x4000000,
0x8000000,
0x10000000,
0x20000000,
0x40000000,
0x80000000
#if SIZEOF_INT==8
,
0x100000000,
0x200000000,
0x400000000,
0x800000000,
0x1000000000,
0x2000000000,
0x4000000000,
0x8000000000,
0x10000000000,
0x20000000000,
0x40000000000,
0x80000000000,
0x100000000000,
0x200000000000,
0x400000000000,
0x800000000000,
0x1000000000000,
0x2000000000000,
0x4000000000000,
0x8000000000000,
0x10000000000000,
0x20000000000000,
0x40000000000000,
0x80000000000000,
0x100000000000000,
0x200000000000000,
0x400000000000000,
0x800000000000000,
0x1000000000000000,
0x2000000000000000,
0x4000000000000000,
0x8000000000000000
#endif
};


/*:42*//*43:*/
#line 1265 "./decluster.w"

static const int floor_log_2_small[16]= {
0,
0,
1,1,
2,2,2,2,
3,3,3,3,3,3,3,3};


/*:43*//*49:*/
#line 1423 "./decluster.w"

static digest_t*digest;

/*:49*//*61:*/
#line 1657 "./decluster.w"

int*parent_of_head;

/*:61*//*91:*/
#line 2250 "./decluster.w"

#if DECLUSTER_DEBUG
static char pp[]= "          ";
#endif

/*:91*/
#line 540 "./decluster.w"

/*44:*/
#line 1313 "./decluster.w"

#if SIZEOF_INT>8
#error "The bit twiddling handles integers of at most 64 bits."
#endif
static inline int
floor_log_2(unsigned int x)
{
int ans= 0;

#if SIZEOF_UNSIGNED_INT==8
if(x&(0xffffffff<<32))ans+= 32,x>>= 32;
#endif
if(x&(0xffff<<16))ans+= 16,x>>= 16;
if(x&(0xff<<8))ans+= 8,x>>= 8;
if(x&(0xf<<4))ans+= 4,x>>= 4;
ans+= floor_log_2_small[x];
return ans;
}


/*:44*//*64:*/
#line 1672 "./decluster.w"

static inline int
decluster_lca(int x,int y)
{
const digest_t xd= digest[x],yd= digest[y];
if(xd.inlabel==yd.inlabel)return(xd.level<=yd.level)?x:y;
else{
const int xil= xd.inlabel,yil= yd.inlabel;
int b,inlabel_z,common,jmask,xpp,ypp;
/*65:*/
#line 1735 "./decluster.w"

{
const unsigned int xfuzz= xil^(xil-1),yfuzz= yil^(yil-1);
const unsigned int lomask= max(xfuzz,yfuzz)>>1;
int hi_diff= (xil^yil)&(~lomask);
if(hi_diff){
int t= hi_diff;
copy_1_down(t);


b= (~t)&xil;
b|= t^(t>>1);
}else{
b= (xfuzz>=yfuzz)?xil:yil;
}
/*92:*/
#line 2256 "./decluster.w"

#if DECLUSTER_DEBUG
if(verbose>=2000){
printf("\n%s x=%d {lev=%d, inl=%d, asc=%d}\n",pp,
x,xd.level,xd.inlabel,xd.ascendant);
printf("%s y=%d {lev=%d, inl=%d, asc=%d}\n",pp,
y,yd.level,yd.inlabel,yd.ascendant);
printf("%s b=%d\n",pp,b);
}
#endif

/*:92*/
#line 1750 "./decluster.w"

}

/*:65*/
#line 1681 "./decluster.w"

/*66:*/
#line 1768 "./decluster.w"

{const int imask= (b^(b-1))>>1;
const int u= (common= xd.ascendant&yd.ascendant)&~imask;
jmask= (u^(u-1))>>1;

inlabel_z= (xd.inlabel&~jmask)|(jmask+1);
/*93:*/
#line 2268 "./decluster.w"

#if DECLUSTER_DEBUG
if(verbose>=2000){
printf("%s inlabel_z=%d, common=%d, jmask=%d\n",pp,inlabel_z,common,jmask);
}
#endif


/*:93*/
#line 1774 "./decluster.w"

}




/*:66*/
#line 1682 "./decluster.w"

/*67:*/
#line 1787 "./decluster.w"

if(xd.inlabel==inlabel_z)xpp= x;
else{
int inlabel_w_x;
int kmask= xd.ascendant&jmask;
copy_1_down(kmask);
inlabel_w_x= (~kmask)&xil;
inlabel_w_x|= kmask^(kmask>>1);
xpp= parent_of_head[inlabel_w_x];
/*94:*/
#line 2277 "./decluster.w"

#if DECLUSTER_DEBUG
if(verbose>=2000){
printf("%s inlabel_w_x=%d, xpp=%d\n",pp,inlabel_w_x,xpp);
}
#endif

/*:94*/
#line 1796 "./decluster.w"

}
if(yd.inlabel==inlabel_z)ypp= y;
else{
int inlabel_w_y;
int kmask= yd.ascendant&jmask;
copy_1_down(kmask);
inlabel_w_y= (~kmask)&yil;
inlabel_w_y|= kmask^(kmask>>1);
ypp= parent_of_head[inlabel_w_y];
/*95:*/
#line 2285 "./decluster.w"

#if DECLUSTER_DEBUG
if(verbose>=2000){
printf("%s inlabel_w_y=%d, ypp=%d\n",pp,inlabel_w_y,ypp);
}
#endif

/*:95*/
#line 1806 "./decluster.w"

}


/*:67*/
#line 1683 "./decluster.w"

return(digest[xpp].level<=digest[ypp].level)?xpp:ypp;
}
}


/*:64*/
#line 541 "./decluster.w"

/*16:*/
#line 664 "./decluster.w"

decluster_tree_t*
decluster_setup(int the_n)
{
decluster_tree_t*mst;
n= the_n;
/*20:*/
#line 708 "./decluster.w"

mst= new_of(decluster_tree_t);
mst->n= n-1;
mst->edge= new_arr_of(decluster_edge_t,n-1);


/*:20*/
#line 670 "./decluster.w"

/*31:*/
#line 1032 "./decluster.w"

/*30:*/
#line 1017 "./decluster.w"

if(T_prime==NULL||T_prime->n!=n+n-1){
int i;
/*32:*/
#line 1036 "./decluster.w"

if(T_prime){
const int n= T_prime->n;
free_mem(T_prime->edge);mem_deduct(n*sizeof(decluster_edge_t));
free_mem(T_prime);mem_deduct(sizeof(decluster_tree_t));
T_prime= NULL;
}

/*:32*/
#line 1020 "./decluster.w"

T_prime= new_of(decluster_tree_t);
T_prime->n= n+n-1;
T_prime->edge= new_arr_of(decluster_edge_t,n+n-1);
for(i= 0;i<n;i++){
T_prime->edge[i].child[0]= NO_CHILD;
T_prime->edge[i].child[1]= NO_CHILD;
T_prime->edge[i].cost= 0;
}
}

/*:30*/
#line 1033 "./decluster.w"


/*:31*//*50:*/
#line 1431 "./decluster.w"

digest= new_arr_of(digest_t,n+n-1);
{
int i;
for(i= 0;i<n;i++){
digest[i].cost= 0;
}
}

/*:50*//*62:*/
#line 1661 "./decluster.w"

parent_of_head= new_arr_of(int,n+n);

/*:62*/
#line 671 "./decluster.w"

return mst;
}

/*:16*//*18:*/
#line 680 "./decluster.w"

void
decluster_cleanup(void)
{
/*33:*/
#line 1045 "./decluster.w"

/*32:*/
#line 1036 "./decluster.w"

if(T_prime){
const int n= T_prime->n;
free_mem(T_prime->edge);mem_deduct(n*sizeof(decluster_edge_t));
free_mem(T_prime);mem_deduct(sizeof(decluster_tree_t));
T_prime= NULL;
}

/*:32*/
#line 1046 "./decluster.w"


/*:33*//*51:*/
#line 1441 "./decluster.w"

free_mem(digest);
mem_deduct((n+n-1)*sizeof(digest_t));

/*:51*//*63:*/
#line 1665 "./decluster.w"

free_mem(parent_of_head);mem_deduct((n+n)*sizeof(int));


/*:63*/
#line 684 "./decluster.w"

n= 0;
}

/*:18*//*21:*/
#line 716 "./decluster.w"

void
decluster_cleanup_tree(decluster_tree_t*T)
{
if(T){
size_t r= (T->n)*sizeof(decluster_edge_t)+sizeof(decluster_tree_t);
T->n= 0;
free_mem(T->edge);
free_mem(T);
mem_deduct(r);
}
}

/*:21*//*22:*/
#line 764 "./decluster.w"

length_t
decluster_mst(tsp_instance_t*tsp_instance,decluster_tree_t*T)
{
extern int verbose;
length_t total_len= 0;
errorif(T->n!=n-1,"decluster_mst: passed storage for a tree with %d"
" vertices instead of %d vertices",T->n,n-1);
if(E2_supports(tsp_instance)){
/*26:*/
#line 908 "./decluster.w"

{
decluster_edge_t*bridge= new_arr_of(decluster_edge_t,n);
pq_t*bridge_pq= pq_create_size(decluster_edge_cmp,n);
int next_edge;char*is_in_component= new_arr_of_zero(char,n);

errorif(bridge_pq==NULL,"Couldn't allocate a priority queue!");

E2_hide(0);
bridge[0].city[0]= 0;
bridge[0].city[1]= E2_nn(0);
bridge[0].cost= cost(0,bridge[0].city[1]);
is_in_component[0]= 1;
pq_insert(bridge_pq,&bridge[0]);

for(next_edge= 0;next_edge<n-1;next_edge++){
int in,out;
decluster_edge_t*short_bridge;
while(1){
short_bridge= pq_delete_min(bridge_pq);
in= short_bridge->city[0];

out= short_bridge->city[1];
if(!is_in_component[out])break;
bridge[in].city[1]= E2_nn(in);
bridge[in].cost= cost(in,bridge[in].city[1]);
pq_insert(bridge_pq,bridge+in);