[TECH] The shopping problem is as hard as TSP.

The Shopping Problem : Let I={i₁,i₂...i_k} be the items on your shopping list. Let S = {s₁,s₂,....,s_m} be the list of stores each store s_i sells only a subset of the items on your shopping list, and also the cost the items sold by the shops may vary i.e item i_j may be sold with different prices in different shops, so let C[i,j] be the cost of item i in store j. Since the shops are located in different places let D[p,q] be the cost of reaching store q from store p. So our aim is to start at home (H) and buy all the items we want and return to home with minimum possible cost. As I mentioned in the previous post this problem is was mentioned in Google Code Jam practice problems, however I did not see many submissions to this problem.

Firstly the problem is NP-hard because if items sold by each of the stores are disjoint and they are singleton sets then we should visit all the stores in the right order to minimize the total distance traveled, this is nothing by TSP. The following is the dynamic programming formulation for this problem. Let S(V,j) be a optimal solution to buy all the items in set V ending at store j. Then the following is the DP formulation.
S(V,j) = min { S(V-{k},p) + D[p,j] + C[k,j] } ; k \in V, 1<=p<=m Final Answer = min { S(I,j) + D[j,Home] } Initialization of the DP is simple you may want to look at the code below.

Further complexity to the problem can be added by saying that the items are classified as PerishableP and Non-PerishableN. So when ever you buy a perishable item you need to go back home to put in the refrigerator. The above DP formulation can be extended to this easily as follows. Let S(V,j,P) be the optimal sub problem of shopping all the items in the set V ending at the store j and having a Perishable item in the shopping cart, on the same lines S(V,j,N) optimal sub problem of shopping all the items in the set V ending at store j with NO perishable item in the shopping cart. Then the following formulation will give the optimal solution.

S(V,j,P) = min {S(V-{k},p,P) + D[p,Home] + D[Home,j] + C[k,j] } or
              min {S(V-{k},p,N) + D[p,j] + C[k,j] } 
If k is a perishable item

We can write similar formulations for Non-Perishable items also, get the code here implementing the shopping problem.

  1 /*Vamsi Kundeti (vamsi(dot)krishnak(at)gmail*/
  2 #include<stdio.h>
  3 #include<stdlib.h>
  4 #include<string.h>
  5 #include<assert.h>
  6 #include<math.h>
  7 #include "List.h"
  8 unsigned int **GRAPH;
  9 //#define VERBOSE_MODE 1
 10 unsigned int item_bits;
 11 char **item_names;
 12 unsigned int item_subset;
 13 /*has a unsigned int indicating
 14 *if this shop sells the item or not
 15 */
 16 unsigned int *shop_sells;
 17 /*dist[i][j] indicates the distance 
 18  *between i and j*/
 19 double **dist;
 20 /*cost[i][j] indicates the cost of
 21 *item i in shop j
 22 */
 23 double **cost;
 24 /*The number of the items*/
 25 unsigned int icount;
 26 /*The number of test cases*/
 27 unsigned int tcount;
 28 unsigned int caseno;
 29 /*The number of shops**/
 30 unsigned int scount;
 31 typedef struct _pair_ {
 32     int x;
 33     int y;
 34 }Coor;
 35 Coor *shop_coor;
 36 void EnumSubSetsBFS(unsigned int n);
 37 unsigned int LookUpItem(const char *key){
 38     unsigned int i=0;
 39     for(i=0;i<icount;i++){
 40         if(!strcmp(key,item_names[i])){
 41             return i+1;
 42         }
 43     }
 44     /*Ignore this item*/
 45     return 0;
 46 }
 47 unsigned int gas_price;
 48 unsigned int perish;
 49 inline char CheckPerish(unsigned int i){
 50     return (perish & (1<<(i)))?1:0;
 51 }
 52 void ReadInput(void){
 53     unsigned int i,j;
 54     unsigned int len;
 55     unsigned int icost;
 56     char item_buffer[64];
 57     char delim; 
 58     scanf("%u",&tcount);
 59     while(tcount--){
 60         perish=0;
 61         scanf("%u %u %u",&icount,&scount,&gas_price);
 62         assert(icount && scount);
 63         /*Read the items to be bought*/
 64         item_names = (char **) malloc(sizeof(char *)*icount);
 65         for(i=0;i<icount;i++){
 66             scanf("%s",&item_buffer);
 67             len = strlen(item_buffer);
 68             len = (len > 64)?64:len;
 69             item_names[i] = (char *) malloc(sizeof(char)*(len+1));
 70             assert(memcpy(item_names[i],item_buffer,(len+1))==item_names[i]);
 71             if(len && (item_names[i])[len-1]=='!'){
 72                 perish = perish | (1 << (i));
 73                 (item_names[i])[len-1]='\0';
 74             }
 75         }
 76         shop_coor = (Coor *) malloc(sizeof(Coor)*(scount+1));
 77         cost = (double **) malloc(sizeof(double *)*(scount+1));
 78         shop_sells = (unsigned int *)malloc(sizeof(unsigned int)*(scount+1));
 79         shop_coor[0].x = 0.0;
 80         shop_coor[0].y = 0.0;
 81         for(i=1;i<scount+1;i++){
 82             cost[i] = (double *) malloc(sizeof(double)*icount);
 83             scanf("%d %d",&(shop_coor[i].x),&(shop_coor[i].y));
 84             shop_sells[i] = 0;
 85             for(j=0;j<icount;j++){
 86                 cost[i][j] = 0;
 87             }
 88             while(1){
 89                 scanf(" %[^:]:%u",item_buffer,&icost);
 90                 if(j=LookUpItem(item_buffer)){
 91                     cost[i][j-1] = icost;
 92                     /*Also set the sells bit_set*/
 93                     shop_sells[i] = shop_sells[i] | (1<<(j-1));
 94                 }
 95                 scanf("%c",&delim);
 96                 if(delim == '\n'){
 97                     break;
 98                 }
 99             }
100         }
101         /*Now Compute the distance matrix*/
102         dist = (double **)malloc(sizeof(double *)*(scount+1));
103         for(i=0;i<scount+1;i++){
104             dist[i] = (double *)malloc(sizeof(double)*(scount+1));
105             for(j=0;j<scount+1;j++){
106                 dist[i][j] = sqrtf((shop_coor[i].x - shop_coor[j].x)*
107                     (shop_coor[i].x - shop_coor[j].x) +(shop_coor[i].y - shop_coor[j].y)*
108                         (shop_coor[i].y - shop_coor[j].y));
109             }
110         }
111 #if 0
112         printf("Items..\n");
113         for(i=0;i<icount;i++){
114             printf("%s ",item_names[i]);
115         }
116         printf("\n");
117         printf("Cost Matrix\n");
118         for(i=0;i<scount+1;i++){
119             printf("(%d,%d)",shop_coor[i].x,shop_coor[i].y);
120             for(j=0;j<icount;j++){
121                 if(cost[i][j]){
122                     printf(" %s[%c]:%lf",item_names[j],((CheckPerish(j))?'P':'U'),cost[i][j]);
123                 }
124             }
125             printf("\n");
126         }
127         printf("\n");
128 #endif
129         /*Make the call to the dynamic programming routine*/
130         EnumSubSetsBFS(icount);
131 
132         /*Freeup the stuff*/
133         for(i=0;i<icount;i++){
134             free(item_names[i]); 
135         }
136         free(item_names);
137         free(shop_coor);
138         free(shop_sells);
139         for(i=1;i<scount+1;i++){
140             free(cost[i]);
141         }
142         free(cost);
143         for(i=0;i<scount+1;i++){
144             free(dist[i]);
145         }
146         free(dist);
147     }
148 }
149 
150 void CreateGraph(unsigned int *n){
151     unsigned int i,j;
152     printf("Enter the Graph Size\n");
153     scanf("%u",n);
154     printf("Enter the weighted graph\n");
155     GRAPH = (unsigned int **) malloc(sizeof(unsigned int *)*(*n));
156     for(i=0;i<*n;i++){
157         GRAPH[i] = (unsigned int *) malloc(sizeof(unsigned int)*(*n));
158         for(j=0;j<*n;j++){
159             scanf("%u",&(GRAPH[i][j]));
160         }
161     }
162 }
163 
164 
165 unsigned int subset_count=0;
166 int main(){
167     /*Create the Graph*/
168     unsigned int n;
169     caseno=0;
170     ReadInput();
171     //EnumSubSetsBFS(n);
172 
173 }
174 /*Returns the new tail of the list*/
175 inline List* ListAppend(List **tail,void *_data){
176     List *new_node = (List *)malloc(sizeof(new_node)*1);
177     assert(new_node);
178     if(!(*tail)){
179         (*tail) = new_node;
180         (*tail)->_data = _data;
181         (*tail)->next = NULL;
182         return *tail;
183     }
184     new_node->_data = _data;
185     new_node->next = NULL;
186     (*tail)->next = new_node;
187     return new_node;
188 }
189 typedef struct _subset_{
190     /*if i^th bit is set then subset has i^th
191      *element*/
192     unsigned int ebits;
193     unsigned int max; /*Maximum element in the subset*/
194     unsigned int size; /*Size of subset*/
195     /*Pointer to the sub-problem S({},i,0) S({},i,1)*/
196     void *sub_problem;
197 }SubSet;
198 
199 /*Enumerate the Subset using BFS, n is 
200  *the set size [1,2...n] */
201 void EnumSubSetsBFS(unsigned int n){
202     unsigned int check_sub_set,k,i,j;
203     unsigned int current_level = 0;
204     List *bfs_list = NULL;
205     List *tail = NULL; List *free_node,*ptr;
206     SubSet *sub_set,*new_sub_set;
207     SubSet *one_sets = NULL;
208 
209     one_sets = (SubSet *)malloc(sizeof(SubSet)*icount); 
210     assert(one_sets);
211     /*Add the 1-subsets i=items*/
212     for(i=0;i<icount;i++){
213         one_sets[i].ebits = 0; one_sets[i].ebits |= (1<<i);
214         one_sets[i].max = i; one_sets[i].size = 1;
215         /****SUBPROBLEM SPECIFIC*****/
216         one_sets[i].sub_problem = (double **)malloc(sizeof(double *)*2);
217         /* Perishable S({i},X,0) */
218         ((double **)one_sets[i].sub_problem)[0] = (double *) malloc(sizeof(double)*(scount+1));
219         /* Unperishable S({i},X,1) */
220         ((double **)one_sets[i].sub_problem)[1] = (double *) malloc(sizeof(double)*(scount+1));
221         double **sub_problem_i = (double **)one_sets[i].sub_problem;
222 
223         /*j is the shop*/
224         for(j=1;j<scount+1;j++){
225             /*if item i can be bought at shop j
226              *depending on perishable or not initialize sub_problem[0] 
227              *sub_problem[1]
228              */
229              sub_problem_i[0][j] = HUGE_VAL; sub_problem_i[1][j] = HUGE_VAL;
230              if(shop_sells[j]&(1<<i)){
231                  if(CheckPerish(i)){
232                      sub_problem_i[0][j] = dist[0][j]*(gas_price) + cost[j][i];
233                  }else{
234                      sub_problem_i[1][j] = dist[0][j]*(gas_price) + cost[j][i];
235                  }
236              }
237 #ifdef VERBOSE_MODE
238             printf("S({%u},%u,0) = %lf\n",i+1,j+1,sub_problem_i[0][j]);
239             printf("S({%u},%u,1) = %lf\n",i+1,j+1,sub_problem_i[1][j]);
240 #endif
241         }
242         /*****************************/
243         if(!tail){
244             ListAppend(&bfs_list,(void *)&(one_sets[i]));
245             tail = bfs_list;
246             continue;
247         }
248         tail = ListAppend(&tail,(void *)&(one_sets[i]));
249     }
250 
251 
252 
253     /*Now do the BFS*/
254     while(bfs_list){
255         sub_set = (SubSet *)bfs_list->_data; 
256         free_node = bfs_list;
257         bfs_list=bfs_list->next;
258         subset_count++;
259         assert(sub_set->max < icount);
260         for(i=((sub_set->max)+1);i<icount;i++){
261             new_sub_set = (SubSet *) malloc(sizeof(SubSet)*1);
262             new_sub_set->size = (sub_set->size)+1;
263             new_sub_set->max = i;
264             new_sub_set->ebits = (sub_set->ebits | (1<<i));
265             /*This subset enumeration is critical to Subset based D.P*/
266 
267             /*****<begin>**SUBPROBLEM SPECIFIC******/
268             /* S(X,j) = min_k {S(X-k,k) + d[k][j]}*/
269             ptr = free_node;
270             double current_min = HUGE_VAL;
271             double **new_sub_problem;
272             double **sub_problem; double local_min;
273             double local_min_perish,local_min_unperish;
274             unsigned int k1,perish,shop1;
275             List *ptr_prev;
276             
277             new_sub_set->sub_problem = (double **) malloc(sizeof(double *)*2);
278             ((double **)new_sub_set->sub_problem)[0] = (double *) malloc(sizeof(double)*(scount+1));
279             ((double **)new_sub_set->sub_problem)[1] = (double *) malloc(sizeof(double)*(scount+1));
280             new_sub_problem = (double **)new_sub_set->sub_problem;
281 
282             for(perish=0;perish<=1;perish++){
283                 for(j=1;j<scount+1;j++){
284                     /*Objective here is to fill up S({new_sub_set},j,0) 
285                      *and S({new_sub_set},j,1)*/
286                     ptr = free_node;
287                     current_min = HUGE_VAL;
288                     for(k=1;k<=new_sub_set->size;k++){
289                         do{
290                             check_sub_set = ((SubSet *)ptr->_data)->ebits; 
291                             check_sub_set ^= new_sub_set->ebits; 
292                             k1 = (int)log2(check_sub_set);
293                             ptr_prev = ptr;
294                             ptr=ptr->next;
295                         }while(check_sub_set & (check_sub_set-1));
296 
297                         /*j is the shop and k1 is the item which you want
298                          *to buy at this shop, depending on
299                          *perishablity create the value of current min.
300                         */
301                         if(!(shop_sells[j] & check_sub_set)){
302                              continue; /*The shop don't sell this*/
303                         }
304                         sub_problem = (double **)
305                             ((SubSet *)ptr_prev->_data)->sub_problem;
306                         /*item k1 is sold by shop j*/
307                         local_min = HUGE_VAL;
308                         if((perish && !CheckPerish(k1)) 
309                                 || (!perish && CheckPerish(k1))){
310                             for(shop1=1;shop1<scount+1;shop1++){
311                                 /*S{{new_sub_set}-{k1},shop1} + 
312                                 d(shop1,h) + d(h,j) + c[k1][j]*/
313                                 if(j != shop1 || perish){
314                                     local_min_perish = 
315                                         sub_problem[0][shop1] + 
316                                         (dist[shop1][0] + dist[0][j])*
317                                             (gas_price) + cost[j][k1];
318                                 }else if(!perish){
319                                     local_min_perish = sub_problem[0][shop1] + 
320                                         cost[j][k1];
321                                 }
322                                 local_min_unperish = sub_problem[1][shop1] + 
323                                     (dist[shop1][j])*(gas_price) + cost[j][k1];
324                                 local_min_perish = 
325                                  (local_min_perish < local_min_unperish)? 
326                                     local_min_perish:local_min_unperish;
327 
328                                 if(local_min_perish < local_min){
329                                     local_min = local_min_perish;
330                                 }
331                             }
332                         }
333                         if(local_min < current_min){
334                             current_min = local_min; 
335                         }
336                     }
337                     new_sub_problem[perish][j] = current_min;
338                 }
339             }
340 
341 #ifdef VERBOSE_MODE
342             printf("Printing S(V,i) size=%u \n",sub_set->size);
343             /*Print SubSet*/
344             unsigned int i1;
345             for(i1=1;i1<scount+1;i1++){
346                 printf("S({");
347                 unsigned int j1;
348                 /*Get the elements back*/
349                 check_sub_set = new_sub_set->ebits;
350                 for(j1=0;j1<icount;j1++){
351                     if(check_sub_set & (1<<j1)){
352                         printf(" %u ",j1+1);
353                     }
354                 }
355                 printf("},%u,0)=%lf\n",i1+1,((double **)(new_sub_set->sub_problem))[0][i1]);
356                 printf("S({ABOVE},%u,1)=%lf\n",i1+1,((double **)(new_sub_set->sub_problem))[1][i1]);
357             }
358 #endif
359 
360             /**********<END>****************/
361             tail = ListAppend(&tail,(void *)new_sub_set);
362         }
363 
364         if(sub_set->size == icount){
365             double final_answer=HUGE_VAL;
366             double local_min; double **sub_problem;
367             sub_problem = (double **) sub_set->sub_problem;
368             for(j=1;j<scount+1;j++){
369                 local_min = (sub_problem[0][j]<sub_problem[1][j])?
370                     sub_problem[0][j]:sub_problem[1][j];
371                 local_min += (dist[j][0])*(gas_price);
372                 if(local_min < final_answer){
373                     final_answer = local_min;
374                 }
375             }
376             printf("Case #%u: %0.7lf\n",++caseno,final_answer);
377         }
378 
379         if(!(sub_set >= one_sets && sub_set <= &(one_sets[icount-1]))){
380             free(((double **)sub_set->sub_problem)[0]);
381             free(((double **)sub_set->sub_problem)[1]);
382             free(sub_set->sub_problem);
383             free(sub_set); 
384         }
385         free(free_node);
386     }
387     for(i=0;i<icount;i++){
388         free(((double **)one_sets[i].sub_problem)[0]);
389         free(((double **)one_sets[i].sub_problem)[1]);
390         free(one_sets[i].sub_problem);
391     }
392     free(one_sets);
393 //    printf("SubsetCount=%u\n",subset_count);
394 }
395 

[TECH] A O(log(n)) time and O(1) space algorithm to find out Nth Fibonacci number.

We often come across Fibonacci numbers in several contexts, if you see Knuth's "Concrete Mathematics" in the chapter on "Generating functions" one interesting example illustrates the possible brick patterns to build a wall whose width is 'n' using bricks of dimension 2x1. Also see the UVA problem here

We know a simple dynamic programming based algorithm can lead us to a O(n) time algorithm with O(1). But we can do much better than that in fact we can do it in O(log(n)) time and O(1) space. We can rewrite the dynamic programming as matrix formulation as follows.

Fn Fn-1

=

1 1 1 0

Fn-1 Fn-2

This is really great we can get the closed form for F_n by repeated substitution. In fact the following is the closed form for F_n is as follows.

Fn Fn-1

=

1 1 1 0

n-1

1 0

So the problem of finding the n^th Fibonacci number is reduces to evaluating X^n-1 where X is a 2x2 matrix as shown above. Now the question is how we can evaluate X^n-1 as efficiently as possible. In fact to compute any a^b we need exactly |_{_}log(n)_{_}| multiplications by making use of the bit representation of b = (2^i₁ + 2^i₂ ...) where i_j are the positions of bits which have 1 in the binary representation of b. I used this trick to get a O(log(n)) time and O(1) time algorithm to compute the n^th Fibonacci number. Please see the code below get this here

  1 /*Find the n^th fibonacci number in O(log(n)) 
  2  *time, using the following idea of matrix formulation 
  3  * 
  4  * [   F_n  ]    [ 1      1  ] [F_{n-1}]
  5  * [        ] = [         ] [       ]
  6  * [ F_{n-1}]    [ 1      0  ] [F_{n-2}]
  7  *
  8  * vamsik(at)engr(dot)uconn July 11, 2008
  9  */
 10 #include<stdio.h>
 11 #include<stdlib.h>
 12 /*Fib returns the n^th Fibonacci number*/
 13 unsigned int Fib(unsigned int n){
 14     unsigned int deg=1;
 15     unsigned int fib[2][2] = {{1,1},{1,0}};
 16     unsigned int fibAns[2][2];
 17     unsigned char initAns=0;
 18     unsigned int mul[4];
 19     if(n<=2){ return (n<2)?n:2;}
 20     while(deg < n-1){
 21          mul[0] = fib[0][0]*fib[0][0] + fib[0][1]*fib[1][0];
 22          mul[1] = fib[0][0]*fib[0][1] + fib[0][1]*fib[1][1];
 23          mul[2] = fib[1][0]*fib[0][0] + fib[1][1]*fib[1][0];
 24          mul[3] = fib[1][0]*fib[0][1] + fib[1][1]*fib[1][1];
 25          fib[0][0] = mul[0]; fib[0][1] = mul[1]; 
 26          fib[1][0] = mul[2]; fib[1][1] = mul[3];
 27          deg <<=1;
 28          if(deg & n-1){ /*i^th bit set*/
 29              if(!initAns){
 30                  fibAns[0][0] = fib[0][0]; fibAns[0][1] = fib[0][1]; 
 31                  fibAns[1][0] = fib[1][0]; fibAns[1][1] = fib[1][1];
 32                 initAns=1;
 33              }else{
 34                  mul[0] = fibAns[0][0]*fib[0][0] + fibAns[0][1]*fib[1][0];
 35                  mul[1] = fibAns[0][0]*fib[0][1] + fibAns[0][1]*fib[1][1];
 36                  mul[2] = fibAns[1][0]*fib[0][0] + fibAns[1][1]*fib[1][0];
 37                  mul[3] = fibAns[1][0]*fib[0][1] + fibAns[1][1]*fib[1][1];
 38                  fibAns[0][0] = mul[0]; fibAns[0][1] = mul[1]; 
 39                  fibAns[1][0] = mul[2]; fibAns[1][1] = mul[3];
 40              }
 41          }
 42     }
 43     if(n-1 & 1){
 44         return fibAns[0][0] + fibAns[0][1];
 45     }
 46     return fibAns[0][0];
 47 }
 48 int main(int argc,char **argv){
 49     unsigned int n;
 50     while(1){
 51         scanf("%u",&n);
 52         if(!n){
 53             exit(0);
 54         }
 55         printf("%u\n",n,Fib(n));
 56     }
 57 }
 58 

[TECH] Illustration of how to use the BFS based subset enumeration algorithmic framework to solve TSP.

In the previous post I posted on the idea of using BFS to enumerate the subsets incrementally based on the size (see the previous post). The power of this core framework is so immense let me show how we can implement the TSP based on the framework. The extra code is between the lines 82-90 (TSP Initialization) and 113-175. The same core can be easily adapted to solve several sub-set based Dynamic programming optimizations like the job-scheduling, vehicle routing problem (VRP). If you have not looked into one of the Google Code Jam 2008 (practice) problems the problem#4 which is a shopping problem can be easily solved by this core framework. Check this code TSP.C (use -std=c99 to compile).

  1 /*******************************
  2  * Subset Enumeration Algorithmic 
  3  * Framework and its illustration in 
  4  * finding TSP. 
  5  ******************************/
  6 #include<stdio.h>
  7 #include<stdlib.h>
  8 #include<string.h>
  9 #include<assert.h>
 10 #include<math.h>
 11 #include "List.h"
 12 unsigned int **GRAPH;
 13 unsigned int subset_count=0;
 14 //#define VERBOSE_MODE 1
 15 void EnumSubSetsBFS(unsigned int n);
 16 
 17 void CreateGraph(unsigned int *n){
 18     unsigned int i,j;
 19     printf("Enter the Graph Size\n");
 20     scanf("%u",n);
 21     if(*n>32){
 22         fprintf(stderr,"TODO: The Bitset currently supports only 32 bits\n");
 23         exit(1);
 24     }
 25     printf("Enter the weighted graph\n");
 26     GRAPH = (unsigned int **) malloc(sizeof(unsigned int *)*(*n));
 27     for(i=0;i<*n;i++){
 28         GRAPH[i] = (unsigned int *) malloc(sizeof(unsigned int)*(*n));
 29         for(j=0;j<*n;j++){
 30             scanf("%u",&(GRAPH[i][j]));
 31         }
 32     }
 33 }
 34 
 35 int main(){
 36     /*Create the Graph*/
 37     unsigned int n;
 38     CreateGraph(&n);
 39     EnumSubSetsBFS(n);
 40 
 41 }
 42 /*Returns the new tail of the list*/
 43 inline List* ListAppend(List **tail,void *_data){
 44     List *new_node = (List *)malloc(sizeof(new_node)*1);
 45     assert(new_node);
 46     if(!(*tail)){
 47         (*tail) = new_node;
 48         (*tail)->_data = _data;
 49         (*tail)->next = NULL;
 50         return *tail;
 51     }
 52     new_node->_data = _data;
 53     new_node->next = NULL;
 54     (*tail)->next = new_node;
 55     return new_node;
 56 }
 57 typedef struct _subset_{
 58     /*if i^th bit is set then subset has i^th
 59      *element*/
 60     unsigned int ebits;
 61     unsigned int max; /*Maximum element in the subset*/
 62     unsigned int size; /*Size of subset*/
 63     /*Pointer to the sub-problem */
 64     void *sub_problem;
 65 }SubSet;
 66 /*Enumerate the Subset using BFS, n is 
 67  *the set size [1,2...n] */
 68 void EnumSubSetsBFS(unsigned int n){
 69     unsigned int check_sub_set,k,i,j;
 70     unsigned int current_level = 0;
 71     List *bfs_list = NULL;
 72     List *tail = NULL; List *free_node,*ptr;
 73     SubSet *sub_set,*new_sub_set;
 74     SubSet *one_sets = NULL;
 75 
 76     one_sets = (SubSet *)malloc(sizeof(SubSet)*n); 
 77     assert(one_sets);
 78     /*Add the 1-subsets*/
 79     for(i=1;i<n;i++){
 80         one_sets[i].ebits = 0; one_sets[i].ebits |= (1<<i);
 81         one_sets[i].max = i; one_sets[i].size = 1;
 82         /****SUBPROBLEM SPECIFIC*****/
 83         one_sets[i].sub_problem = (float *)malloc(sizeof(float)*n);
 84         for(j=1;j<n;j++){
 85             if(i==j){ continue; }
 86             /* TSP({i},j) */
 87             ((float *)one_sets[i].sub_problem)[j] = GRAPH[0][i]+GRAPH[i][j];
 88             printf("S({%u},%u) = %f\n",i+1,j+1,((float *)one_sets[i].sub_problem)[j]);
 89         }
 90         /*****************************/
 91         if(!tail){
 92             ListAppend(&bfs_list,(void *)&(one_sets[i]));
 93             tail = bfs_list;
 94             continue;
 95         }
 96         tail = ListAppend(&tail,(void *)&(one_sets[i]));
 97     }
 98     /*Now do the BFS*/
 99     while(bfs_list){
100         sub_set = (SubSet *)bfs_list->_data; 
101         free_node = bfs_list;
102         bfs_list=bfs_list->next;
103         subset_count++;
104         assert(sub_set->max < n);
105         for(i=((sub_set->max)+1);i<n;i++){
106             new_sub_set = (SubSet *) malloc(sizeof(SubSet)*1);
107             new_sub_set->size = (sub_set->size)+1;
108             new_sub_set->max = i;
109             new_sub_set->ebits = (sub_set->ebits | (1<<i));
110             
111             
112             
113             /*****<begin>**SUBPROBLEM SPECIFIC******/
114             /* S(X,j) = min_k {S(X-k,k) + d[k][j]}*/
115             ptr = free_node;
116             float current_min = HUGE_VAL;
117             float *new_sub_problem;
118             float *sub_problem; float local_min;
119             unsigned int k1,jstart,jend;
120             List *ptr_prev;
121             
122             new_sub_set->sub_problem = (float *) malloc(sizeof(float)*n);
123             new_sub_problem = (float *)new_sub_set->sub_problem;
124 
125             if(new_sub_set->size == n-1){
126                 jstart = 0; jend = 1;
127             }else{
128                 jstart = 1; jend = n;
129             }
130 
131             for(j=jstart;j<jend;j++){
132                 ptr = free_node;
133                 current_min = HUGE_VAL;
134                 if(new_sub_set->ebits & (1<<j)){
135                     /*Avoid computation of redundant subproblems like
136                      *S({....,x_i,x_j,x_k,....},x_i) */
137                     new_sub_problem[j]=HUGE_VAL;
138                     continue;
139                 }
140 
141                 for(k=1;k<=new_sub_set->size;k++){
142                     do{
143                         check_sub_set = ((SubSet *)ptr->_data)->ebits; 
144                         check_sub_set ^= new_sub_set->ebits; 
145                         k1 = (int)log2(check_sub_set);
146                         ptr_prev = ptr;
147                         ptr=ptr->next;
148                     }while(check_sub_set & (check_sub_set-1));
149                     sub_problem = (float *)((SubSet *)ptr_prev->_data)->sub_problem;
150                     local_min = sub_problem[k1] + GRAPH[k1][j];
151 
152                     if(local_min < current_min){
153                         current_min = local_min; 
154                         new_sub_problem[j] = current_min;
155                     }
156 
157                 }
158             }
159             /*Print Solution to the Subproblem */
160             unsigned int i1,istart,iend;
161             unsigned int j1;
162             
163             for(i1=jstart;i1<jend;i1++){
164                 check_sub_set = new_sub_set->ebits;
165                 if(check_sub_set & (1<<i1)){ continue; }
166                 printf("S(%u,{ ",i1+1);
167                 /*Get the elements back*/
168                 for(j1=0;j1<n;j1++){
169                     if(check_sub_set & (1<<j1)){
170                         printf(" %u ",j1+1);
171                     }
172                 }
173                 printf("})=%f\n",((float *)(new_sub_set->sub_problem))[i1]);
174             }
175             /**********<END>***SUBPROBLEM SPECIFIC*************/
176 
177 
178             tail = ListAppend(&tail,(void *)new_sub_set);
179         }
180         if(!(sub_set >= one_sets && sub_set <= &(one_sets[n-1]))){
181             free(sub_set->sub_problem);
182             free(sub_set); 
183         }
184         free(free_node);
185     }
186     for(i=0;i<n;i++){
187         free(one_sets[i].sub_problem);
188     }
189     free(one_sets);
190     printf("SubsetCount=%u\n",subset_count);
191 }
192 

Sample Input.
---------test5.txt----------
5
0 3 1 5 4 
1 0 5 4 3 
5 4 0 2 1 
3 1 3 0 3 
5 2 4 1 0
---------------------
[vamsik@abadon ~]$ ./a.out < test5.txt 
Enter the Graph Size
Enter the weighted graph
S({2},3) = 8.000000
S({2},4) = 7.000000
S({2},5) = 6.000000
S({3},2) = 5.000000
S({3},4) = 3.000000
S({3},5) = 2.000000
S({4},2) = 6.000000
S({4},3) = 8.000000
S({4},5) = 8.000000
S({5},2) = 6.000000
S({5},3) = 8.000000
S({5},4) = 5.000000
S(4,{  2  3 })=9.000000
S(5,{  2  3 })=8.000000
S(3,{  2  4 })=10.000000
S(5,{  2  4 })=9.000000
S(3,{  2  5 })=10.000000
S(4,{  2  5 })=7.000000
S(2,{  3  4 })=4.000000
S(5,{  3  4 })=6.000000
S(2,{  3  5 })=4.000000
S(4,{  3  5 })=3.000000
S(2,{  4  5 })=6.000000
S(3,{  4  5 })=8.000000
S(5,{  2  3  4 })=7.000000
S(4,{  2  3  5 })=8.000000
S(3,{  2  4  5 })=10.000000
S(2,{  3  4  5 })=4.000000
S(1,{  2  3  4  5 })=5.000000
SubsetCount=15

[TECH] A note on enumerating subsets using BFS.

I have several things which I want to make a note on my blog, but unfortunately I'm not getting enough time to write the blog posts. Any way.

Recently I came across several Dynamic Programming formulations which involve subsets in the Subproblem for example take the TSP problem, the Subproblem T(S-{j},j) is often used in the D.P formulation which indicates the optimal path starting from node '1' and covering all the nodes in the subset 'S' and ending in node 'j'. Here 'S' is a subset of V (set of all the nodes/vertices's of the given graph.). The final answer for the TSP problem is found by finding a 'k' such that T(V-{k},k) + d(k,1) is minimum.

If we observe we need a small algorithm framework which moves the solution from all the subsets of size '1' to '|V|'. In simple terms we need to enumerate the subsets in the order of their sizes. I found that a BFS can really help in enumerating the Subsets. The following code should be very handy in providing an enumeration framework and will be greatly useful to solve optimization problems which involve subsets in their D.P formulation.

---

#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<assert.h>
#include<math.h>
#include "List.h"
#define VERBOSE_MODE 1

/*Returns the new tail of the list*/
inline List* ListAppend(List **tail,void *_data){
    List *new_node = (List *)malloc(sizeof(new_node)*1);
    assert(new_node);
    if(!(*tail)){
        (*tail) = new_node;
        (*tail)->_data = _data;
        (*tail)->next = NULL;
        return *tail;
    }
    new_node->_data = _data;
    new_node->next = NULL;
    (*tail)->next = new_node;
    return new_node;
}
typedef struct _subset_{
    /*if i^th bit is set then subset has i^th
     *element*/
    unsigned int ebits;
    unsigned int max; /*Maximum element in the subset*/
    unsigned int size; /*Size of subset*/
}SubSet;
/*Enumerate the Subset using BFS, n is 
 *the set size [1,2...n] */
void EnumSubSetsBFS(unsigned int n){
    unsigned int i;
    List *bfs_list = NULL;
    List *tail = NULL; List *free_node;
    SubSet *sub_set,*new_sub_set;
    SubSet *one_sets = (SubSet *)malloc(sizeof(SubSet)*n);
    unsigned int current_level = 0;
    /*Add the 1-subsets*/
    for(i=0;i<n;i++){
        one_sets[i].ebits = 0; one_sets[i].ebits |= (1<<i);
        one_sets[i].max = i; one_sets[i].size = 1;
        if(!tail){
            ListAppend(&bfs_list,(void *)&(one_sets[i]));
            tail = bfs_list;
            continue;
        }
        tail = ListAppend(&tail,(void *)&(one_sets[i]));
    }
    /*Now do the BFS*/
    while(bfs_list){
        sub_set = (SubSet *)bfs_list->_data; 
        free_node = bfs_list;
        bfs_list=bfs_list->next;

#ifdef VERBOSE_MODE 
        if(sub_set->size != current_level){
            current_level = sub_set->size;
            printf("Printing SubSets of Size [%u] \n",current_level);
            printf("===============================\n");
        }
        /*print this set*/
        printf("[ ");
        for(i=0;i<n;i++){
            if(sub_set->ebits & (1<<i)){
                printf("%u ",i);
            }
        }
        printf("]\n");
#endif

        assert(sub_set->max < n);
        for(i=((sub_set->max)+1);i<n;i++){
            new_sub_set = (SubSet *) malloc(sizeof(SubSet)*1);
            new_sub_set->size = (sub_set->size)+1;
            new_sub_set->max = i;
            new_sub_set->ebits = (sub_set->ebits | (1<<i));
            tail = ListAppend(&tail,(void *)new_sub_set);
        }
        if(!(sub_set >= one_sets && sub_set <= &(one_sets[n-1]))){
            free(sub_set); 
        }
        free(free_node);
    }
    free(one_sets);
}
int main(int argc,char **argv){
    EnumSubSetsBFS(5);
    return 0;
}

Thanks to this http://puzzleware.net/codehtmler/default.aspx for converting my code.

Cheers! Vamsi.