main.c

   1
   2 #include <stdio.h>
   3 #include <stdint.h>
   4 #include <stdlib.h>
   5 #include <inttypes.h>
   6
   7 #include <publib.h>
   8
   9 typedef uint32_t AdjWord;
  10 #define PRADJ "08"PRIx32
  11
  12 static int n, m;
  13 static AdjWord *adjmatrix;
  14 static AdjWord adjall;
  15
  16 static double best;
  17
  18 static void prep(void) {
  19   adjall = ~((~(AdjWord)0) << m);
  20   adjmatrix = xmalloc(sizeof(*adjmatrix)*n);
  21 }
  22
  23 static AdjWord one_adj_bit(int bitnum) {
  24   return (AdjWord)1 << j;
  25 }
  26
  27 static int count_set_adj_bits(AdjWord w) {
  28   int j, total;
  29   for (j=0, total=0; j<m; j++)
  30     total += !!(w & one_adj_bit(j));
  31   return total;
  32 }
  33
  34 static void optimise(void) {
  35   int i, totalfrags;
  36   for (i=0, totalfrags=0; i<n; i++) {
  37     int frags = count_set_adj_bits(adjmatrix[i]);
  38     totalfrags += frags;
  39     printf("%"PRADJ" ", adjmatrix[i]);
  40     double maxminsize = (double)m / frags;
  41     if (maxminsize < best) {
  42       printf(" too fine\n");
  43       return;
  44     }
  45   }
  46
  47   /*
  48    * We formulate our problem as an LP problem as follows.
  49    * In this file "n" and "m" are the matchstick numbers.
  50    *
  51    * Each set bit in the adjacency matrix corresponds to taking a
  52    * fragment from old match i and making it part of new match j.
  53    *
  54    * The structural variables (columns) are:
  55    *   x_minimum        minimum size of any fragment (bounded below by 0)
  56    *   x_morefrag_i_j   the amount by which the size of the fragment
  57    *                     i,j exceeds the minimum size (bounded below by 0)
  58    *
  59    * The auxiliary variables (rows) are:
  60    *   x_total_i       total length for each input match (fixed variable)
  61    *   x_total_j       total length for each output match (fixed variable)
  62    *
  63    * The objective function is simply
  64    *   maximise x_minimum
  65    *
  66    * We use X_ and Y_ to refer to GLPK's (1-based) column and row indices.
  67    * ME_ refers to entries in the list of constraint matrix elements
  68    * which we build up as we go.
  69    */
  70
  71   glp_prob *prob = glp_create_prob();
  72
  73   int Y_totals_i = glp_add_rows(prob, i);
  74   int Y_totals_j = glp_add_rows(prob, j);
  75   int X_minimum = glp_add_cols(prob, 1);
  76   int rows = glp_get_num_rows(prob);
  77   int cols = glp_get_num_rows(cols);
  78
  79   int next_matrix_entry = 1; /* wtf GLPK! */
  80   int matrix_entries_size = next_matrix_entry + i + j + totalfrags*2;
  81   double matrix_entries[matrix_entries_size];
  82   int matrix_entries_XY[2][matrix_entries_size];
  83
  84 #define ADD_MATRIX_ENTRY(Y,X) ({                        \
  85       assert(matrix_entries_size < next_matrix_entry);  \
  86       matrix_entries_XY[0][next_matrix_entry] = X;      \
  87       matrix_entries_XY[1][next_matrix_entry] = Y;      \
  88       matrix_entries[next_matrix_entry] = 0;            \
  89       next_matrix_entry++;                              \
  90     })
  91
  92   int ME_totals_i__minimum = next_matrix_entry;
  93   for (i=0; i<n; i++) ADD_MATRIX_ENTRY(Y_totals_i+i, X_minimum);
  94
  95   int ME_totals_j__minimum = next_matrix_entry;
  96   for (j=0; j<m; j++) ADD_MATRIX_ENTRY(Y_totals_j+j, X_minimum);
  97
  98   /* \forall_i x_totals_i = m */
  99   /* \forall_i x_totals_j = n */
 100   for (i=0; i<n; i++) glp_set_row_bounds(prob, Y_totals_i+i, GLP_FX, m,m);
 101   for (j=0; j<m; j++) glp_set_row_bounds(prob, Y_totals_j+j, GLP_FX, n,n);
 102
 103   /* x_minimum >= 0 */
 104   glp_set_col_bounds(prob, X_minimum, GLP_LB, 0, 0);
 105
 106   /* objective is maximising x_minimum */
 107   glp_set_obj_dir(prob, GLP_MAX);
 108   glp_set_obj_coef(prob, X_minimum, 1);
 109
 110   for (i=0; i<n; j++) {
 111     for (j=0; j<m; j++) {
 112       if (!(adjmatrix[i] & one_adj_bit(j)))
 113         continue;
 114       /* x_total_i += x_minimum */
 115       /* x_total_j += x_minimum */
 116       matrix_entries[ ME_totals_i__minimum + i ] ++;
 117       matrix_entries[ ME_totals_j__minimum + j ] ++;
 118
 119       /* x_morefrag_i_j >= 0 */
 120       int X_morefrag_i_j = glp_add_cols(prob, 1);
 121       glp_set_col_bnds(prob, X_morefrag_i_j, GLP_LO, 0, 0);
 122
 123       /* x_total_i += x_morefrag_i_j */
 124       /* x_total_j += x_morefrag_i_j */
 125       int ME_totals_i__mf_i_j = ADD_MATRIX_ENTRY(Y_totals_i+i, X_morefrag_i_j);
 126       int ME_totals_j__mf_i_j = ADD_MATRIX_ENTRY(Y_totals_j+j, X_morefrag_i_j);
 127       matrix_entries[ME_totals_i__mf_i_j] = 1;
 128       matrix_entries[ME_totals_j__mf_i_j] = 1;
 129     }
 130   }
 131
 132   assert(next_matrix_entry == matrix_entries_size);
 133
 134   for (row=1; row<=rows; row++) {
 135     glp_load_matrix(prob, next_matrix_entry,
 136                     matrix_entries_XY[1], matrix_entries_XY[0],
 137                     matrix_entries);
 138
 139   printf("nyi\n");
 140 }
 141
 142 static void iterate_recurse(int i, AdjWord min) {
 143   if (i >= n) {
 144     optimise();
 145     return;
 146   }
 147   for (adjmatrix[i] = min;
 148        ;
 149        adjmatrix[i]++) {
 150     iterate_recurse(i+1, adjmatrix[i]);
 151     if (adjmatrix[i] == adjall)
 152       return;
 153   }
 154 }
 155
 156 static void iterate(void) {
 157   iterate_recurse(0, 1);
 158 }
 159
 160 int main(int argc, char **argv) {
 161   n = atoi(argv[1]);
 162   m = atoi(argv[2]);
 163   prep();
 164   iterate();
 165   if (ferror(stdout) || fclose(stdout)) { perror("stdout"); exit(-1); }
 166   return 0;
 167 }