Merge branch 'charm' of charmgit:charm into charm
[charm.git] / src / ck-ldb / tm_bucket.c
1 #ifdef _WIN32
2 #include <windows.h>
3 #include <winbase.h>
4 #define random() rand()
5 #define srandom(x)  srand(x)
6 #endif
7 #include <stdio.h>
8 #include <float.h>
9 #include <math.h>
10 #include <assert.h>
11 #include "treematch.h"
12 #include "treebucket.h"
13 #include "treetimings.h" 
14
15 #undef DEBUG
16 bucket_list_t global_bl; 
17
18 int tab_cmp(const void* x1,const void* x2){ 
19   int *e1,*e2,i1,i2,j1,j2;
20   double **tab;
21   bucket_list_t bl;
22
23   bl=global_bl;
24
25   e1=((int *)x1);
26   e2=((int *)x2);
27
28   
29   tab=bl->tab;
30
31   i1=e1[0];
32   j1=e1[1];
33   i2=e2[0];
34   j2=e2[1];
35
36   return tab[i1][j1]>tab[i2][j2]?-1:1;
37 }
38
39
40 int old_bucket_id(int i,int j,bucket_list_t bucket_list){
41   double *pivot,val;
42   int n,sup,inf,p;
43   pivot=bucket_list->pivot;
44   n=bucket_list->nb_buckets;  
45   val=bucket_list->tab[i][j];
46   
47   inf=-1;
48   sup=n;
49
50   while(sup-inf>1){
51     p=(sup+inf)/2;
52     //printf("%f [%d,%d,%d]=%f\n",val,inf,p,sup,pivot[p]);
53     if(val<pivot[p]){
54       inf=p;
55       if(inf==sup)
56         inf--;
57     }else{
58       sup=p;
59       if(sup==inf)
60         sup++;
61     }
62   }
63   //exit(-1);
64   return sup;
65 }
66
67
68 int bucket_id(int i,int j,bucket_list_t bucket_list){
69   double *pivot_tree,val;
70   int p,k;
71   pivot_tree=bucket_list->pivot_tree;
72   val=bucket_list->tab[i][j];
73
74
75
76   p=1;
77   for(k=0;k<bucket_list->max_depth;k++){
78     if(val>pivot_tree[p])
79       p=p*2;
80     else
81       p=p*2+1;
82   }
83
84   return  (int)pivot_tree[p];
85 }
86   
87
88
89 void  display_bucket(bucket_t *b){
90   printf("\tb.bucket=%p\n",b->bucket);
91   printf("\tb.bucket_len=%d\n",(int)b->bucket_len);
92   printf("\tb.nb_elem=%d\n",(int)b->nb_elem);
93
94 }
95
96 void check_bucket(bucket_t *b,double **tab,double inf, double sup,int N){
97   int k,i,j;
98   for(k=0;k<b->nb_elem;k++){
99     i=b->bucket[k].i;
100     j=b->bucket[k].j;
101     if((tab[i][j]<inf)||(tab[i][j]>sup)){
102       printf("[%d] (%d,%d):%f not in [%f,%f]\n",k,i,j,tab[i][j],inf,sup);
103       exit(-1);
104     }
105   }
106 }
107
108 void  display_bucket_list(bucket_list_t bucket_list){
109   int i;
110   double inf,sup;
111
112   for(i=0;i<bucket_list->nb_buckets-1;i++){
113     printf("pivot[%d]=%f\n",i,bucket_list->pivot[i]);
114   }
115   printf("\n");
116
117   for(i=0;i<bucket_list->nb_buckets;i++){
118     inf=bucket_list->pivot[i];
119     sup=bucket_list->pivot[i-1];
120     if(i==0)
121       sup=DBL_MAX;
122     if(i==bucket_list->nb_buckets-1)
123       inf=0;
124     printf("Bucket %d:\n",i);
125     display_bucket(bucket_list->bucket_tab[i]);
126     printf("\n");
127     check_bucket(bucket_list->bucket_tab[i],bucket_list->tab,inf,sup,bucket_list->N);
128   }
129   
130 }
131
132 void add_to_bucket(int id,int i,int j,bucket_list_t bucket_list){
133   bucket_t *bucket;
134   int N,n,size;
135
136   bucket=bucket_list->bucket_tab[id];
137   //display_bucket(bucket);
138   
139   if(bucket->bucket_len==bucket->nb_elem){
140     N=bucket_list->N;
141     n=bucket_list->nb_buckets;  
142     size=N*N/n;
143     //display_bucket(bucket);
144     bucket->bucket=(coord*)realloc(bucket->bucket,sizeof(coord)*(size+bucket->bucket_len));
145     bucket->bucket_len+=size;
146 #ifdef DEBUG
147     printf("malloc/realloc: %d\n",id);
148     printf("(%d,%d)\n",i,j);
149     display_bucket(bucket);
150     printf("\n");
151 #endif
152   }
153   
154  bucket->bucket[bucket->nb_elem].i=i;
155  bucket->bucket[bucket->nb_elem].j=j;
156  bucket->nb_elem++;
157
158   //printf("\n");
159   //exit(-1);
160 }
161
162 void dfs(int i,int inf,int sup,double *pivot,double *pivot_tree,int depth,int max_depth){
163   int p;
164   if(depth==max_depth)
165     return;
166
167   p=(inf+sup)/2;
168   pivot_tree[i]=pivot[p-1];
169
170   dfs(2*i,inf,p-1,pivot,pivot_tree,depth+1,max_depth);
171   dfs(2*i+1,p+1,sup,pivot,pivot_tree,depth+1,max_depth);
172 }
173
174 void  built_pivot_tree(bucket_list_t bucket_list){
175   double *pivot_tree,*pivot;
176   int n,i,k;
177   pivot=bucket_list->pivot;
178   n=bucket_list->nb_buckets;
179   pivot_tree=(double*)malloc(sizeof(double)*2*n);
180   bucket_list->max_depth=(int)log2(n);
181
182   dfs(1,1,n-1,pivot,pivot_tree,0,bucket_list->max_depth);
183
184   k=0;
185   for(i=n;i<2*n;i++)
186     pivot_tree[i]=k++;
187
188   bucket_list->pivot_tree=pivot_tree;  
189   /*
190   for(i=0;i<2*n;i++)
191     printf("%d:%f\t",i,pivot_tree[i]);
192   printf("\n");
193   */
194 }
195
196 void fill_buckets(bucket_list_t bucket_list){
197   int N,i,j,id;
198
199   N=bucket_list->N;
200
201   for(i=0;i<N;i++){
202     for(j=i+1;j<N;j++){
203       id=bucket_id(i,j,bucket_list);
204       add_to_bucket(id,i,j,bucket_list);
205     }
206   }
207   
208   
209 }
210
211 int is_power_of_2(int val){
212   int n=1;
213   do{
214     if(n==val)
215       return 1;
216     n<<=1;
217   }while(n>0);
218   return 0;
219 }
220
221
222 void partial_sort(bucket_list_t *bl,double **tab,int N,int nb_buckets){
223   int *sample;
224   int i,j,k,n;
225   int id;
226   double *pivot;
227   bucket_list_t bucket_list;
228
229
230   if(!is_power_of_2(nb_buckets)){
231     fprintf(stderr,"Error! Paramater nb_buckets is: %d and should be a power of 2\n",nb_buckets);
232     exit(-1);
233   }
234
235
236   bucket_list=(bucket_list_t)malloc(sizeof(_bucket_list_t));
237
238   bucket_list->tab=tab;
239   bucket_list->N=N;
240
241
242   n=pow(nb_buckets,2);
243   n=N;
244   sample=(int*)malloc(2*sizeof(int)*n);
245   
246   for(k=0;k<n;k++){
247     i=random()%(N-2)+1;
248     if(i==N-2)
249       j=N-1;
250     else
251       j=random()%(N-i-2)+i+1;
252     assert(i!=j);
253     assert(i<j);
254     assert(i<N);
255     assert(j<N);
256     sample[2*k]=i;
257     sample[2*k+1]=j;
258   }
259   
260   global_bl=bucket_list;
261   qsort(sample,n,2*sizeof(int),tab_cmp);
262
263   /*  for(k=0;k<N;k++){
264     i=sample[k]/N;
265     j=sample[k]%N;
266     printf("%f\n",tab[i][j]);
267   }
268   */
269   pivot=(double*)malloc(sizeof(double)*nb_buckets-1);
270   id=2;
271   for(k=1;k<nb_buckets;k++){
272     i=sample[2*(id-1)];
273     j=sample[2*(id-1)+1];
274     id*=2;
275     
276
277     /*    i=sample[k*N/nb_buckets]/N;
278           j=sample[k*N/nb_buckets]%N;*/
279     pivot[k-1]=tab[i][j];
280     //printf("pivot[%d]=%f\n",k-1,tab[i][j]);
281   }
282
283   bucket_list->pivot=pivot;
284   bucket_list->nb_buckets=nb_buckets;
285   built_pivot_tree(bucket_list);
286   
287   bucket_list->bucket_tab=(bucket_t**)malloc(nb_buckets*sizeof(bucket_t*));
288   for(i=0;i<nb_buckets;i++){
289     bucket_list->bucket_tab[i]=(bucket_t*)calloc(1,sizeof(bucket_t));
290   }
291
292   fill_buckets(bucket_list);
293   
294   //display_bucket_list(bucket_list);
295
296   bucket_list->cur_bucket=0;
297   bucket_list->bucket_indice=0;
298   
299   free(sample);
300
301   *bl=bucket_list;
302 }
303
304 void next_bucket_elem(bucket_list_t bucket_list,int *i,int *j){
305   int N;
306   bucket_t *bucket=bucket_list->bucket_tab[bucket_list->cur_bucket];
307
308
309   while(bucket->nb_elem<=bucket_list->bucket_indice){
310     bucket_list->bucket_indice=0;
311     bucket=bucket_list->bucket_tab[bucket_list->cur_bucket++];
312
313     printf("### From bucket %d to bucket %d\n",bucket_list->cur_bucket-1,bucket_list->cur_bucket);
314   }
315
316   if(!bucket->sorted){
317     global_bl=bucket_list;
318     qsort(bucket->bucket,bucket->nb_elem,2*sizeof(int),tab_cmp);
319     bucket->sorted=1;
320   }
321
322
323   
324   N=bucket_list->N;
325
326   *i=bucket->bucket[bucket_list->bucket_indice].i;
327   *j=bucket->bucket[bucket_list->bucket_indice].j;
328   bucket_list->bucket_indice++;
329 }
330
331
332 int add_edge_3(double **tab,tree_t *tab_node, tree_t *parent,int i,int j,int N,int *nb_groups){
333   //printf("%d <-> %d ?\n",tab_node[i].id,tab_node[j].id);
334
335   if((!tab_node[i].parent) && (!tab_node[j].parent)){
336     if(parent){
337       parent->child[0]=&tab_node[i];
338       parent->child[1]=&tab_node[j];
339       tab_node[i].parent=parent;
340       tab_node[j].parent=parent;
341 #ifdef DEBUG
342       printf("%d: %d-%d\n",*nb_groups,parent->child[0]->id,parent->child[1]->id);
343 #endif
344       return 1;
345     } 
346     return 0;
347   }
348   
349   if(tab_node[i].parent && (!tab_node[j].parent)){
350     parent=tab_node[i].parent;
351     if(!parent->child[2]){
352       parent->child[2]=&tab_node[j];
353       tab_node[j].parent=parent;
354 #ifdef DEBUG
355       printf("%d: %d-%d-%d\n",*nb_groups,parent->child[0]->id,parent->child[1]->id,parent->child[2]->id);
356 #endif
357       (*nb_groups)++;
358     }
359     return 0;
360   }
361
362   if(tab_node[j].parent && (!tab_node[i].parent)){
363     parent=tab_node[j].parent;
364     if(!parent->child[2]){
365       parent->child[2]=&tab_node[i];
366       tab_node[i].parent=parent;
367 #ifdef DEBUG
368       printf("%d: %d-%d-%d\n",*nb_groups,parent->child[0]->id,parent->child[1]->id,parent->child[2]->id);
369 #endif      
370       (*nb_groups)++;
371     }
372     return 0;
373   }
374
375   return 0;
376 }
377
378 int try_add_edge(double **tab,tree_t *tab_node, tree_t *parent,int arity,int i,int j,int N,int *nb_groups){
379
380   assert(i!=j);
381
382   
383   switch(arity){
384   case 2:
385     if(tab_node[i].parent)
386       return 0;
387     if(tab_node[j].parent)
388       return 0;
389
390     parent->child[0]=&tab_node[i];
391     parent->child[1]=&tab_node[j];
392     tab_node[i].parent=parent;
393     tab_node[j].parent=parent;
394     
395     (*nb_groups)++;
396
397     return 1;
398   case 3:
399     return add_edge_3(tab,tab_node,parent,i,j,N,nb_groups);
400   default:
401     fprintf(stderr,"Cannot handle arity %d\n",parent->arity);
402     exit(-1);
403   }
404 }
405
406
407 void free_bucket(bucket_t *bucket){
408   free(bucket->bucket);
409   free(bucket);
410
411 }
412
413 void free_tab_bucket(bucket_t **bucket_tab,int N){
414   int i;
415   for(i=0;i<N;i++){
416     free_bucket(bucket_tab[i]);
417   }
418   free(bucket_tab);
419 }
420
421
422 void free_bucket_list(bucket_list_t bucket_list){
423
424   // Do not free the tab field it is used elsewhere
425
426   free_tab_bucket(bucket_list->bucket_tab,bucket_list->nb_buckets);
427   free(bucket_list->pivot);
428   free(bucket_list->pivot_tree);
429   free(bucket_list);
430 }
431
432 void bucket_grouping(double **tab,tree_t *tab_node, tree_t *new_tab_node, int arity,int N, int M,long int k){
433   bucket_list_t bucket_list;
434   double duration;
435   int l,i,j,nb_groups;
436   double val=0;
437  
438   TIC;
439   partial_sort(&bucket_list,tab,N,8);
440   duration=TOC;
441   printf("Partial sorting=%fs\n",duration);  
442  
443   TIC;
444   l=0;
445   i=0;
446   nb_groups=0;
447   while(l<M){
448     next_bucket_elem(bucket_list,&i,&j);
449     if(try_add_edge(tab,tab_node,&new_tab_node[l],arity,i,j,N,&nb_groups)){
450       l++;
451     }
452   }
453
454 #ifdef DEBUG
455   printf("l=%d,nb_groups=%d\n",l,nb_groups);
456 #endif
457
458   while(nb_groups<M){
459     next_bucket_elem(bucket_list,&i,&j);
460     try_add_edge(tab,tab_node,NULL,arity,i,j,N,&nb_groups);
461   }
462
463 #ifdef DEBUG
464   printf("l=%d,nb_groups=%d\n",l,nb_groups);
465 #endif
466
467   for(l=0;l<M;l++){
468     update_val(tab,&new_tab_node[l],N);      
469     val+=new_tab_node[l].val;
470   }
471
472
473       
474
475
476   duration=TOC;
477   printf("Grouping =%fs\n",duration);  
478
479   printf("Bucket: %d, indice:%d\n",bucket_list->cur_bucket,bucket_list->bucket_indice);
480
481   printf("val=%f\n",val);
482   free_bucket_list(bucket_list);
483
484   //  exit(-1);
485
486   //  display_grouping(new_tab_node,M,arity,val);
487   
488
489 }
490