From: Ian Jackson <ijackson@chiark.greenend.org.uk>
Date: Fri, 7 Mar 2014 18:54:57 +0000 (+0000)
Subject: commentary
X-Git-Tag: v1~3
X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ianmdlvl/git?p=matchsticks-search.git;a=commitdiff_plain;h=55f78b65884b33c8e0d991a52a5335dcc8a00134

commentary
---

diff --git a/main.c b/main.c
index 5380737..43cfc3e 100644
--- a/main.c
+++ b/main.c
@@ -1,3 +1,26 @@
+/*
+ * Searches for "good" ways to divide n matchsticks up and reassemble them
+ * into m matchsticks.  "Good" means the smallest fragment is as big
+ * as possible.
+ *
+ * Invoke as   ./main n m
+ *
+ * The algorithm is faster if the arguments are ordered so that n > m.
+ */
+
+/*
+ * matchsticks/main.c  Copyright 2014 Ian Jackson
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ */
 
 #include <stdio.h>
 #include <stdint.h>
@@ -10,6 +33,48 @@
 #include <publib.h>
 #include <glpk.h>
 
+/*
+ * Algorithm.
+ *
+ * Each input match contributes, or does not contribute, to each
+ * output match; we do not need to consider multiple fragments
+ * relating to the same input/output pair this gives an n*m adjacency
+ * matrix (bitmap).  Given such an adjacency matrix, the problem of
+ * finding the best sizes for the fragments can be expressed as a
+ * linear programming problem.
+ *
+ * We search all possible adjacency matrices, and for each one we run
+ * GLPK's simplex solver.  We represent the adjacency matrix as an
+ * array of bitmaps.
+ *
+ * However, there are a couple of wrinkles:
+ *
+ * To best represent the problem as a standard LP problem, we separate
+ * out the size of each fragment into a common minimum size variable,
+ * plus a fragment-specific extra size variable.  This reduces the LP
+ * problem size at the cost of making the problem construction, and
+ * interpretation of the results, a bit fiddly.
+ *
+ * Many of the adjacency matrices are equivalent.  In particular,
+ * permutations of the columns, or of the rows, do not change the
+ * meaning.  It is only necessasry to consider any one permutation.
+ * We make use of this by considering only adjacency matrices whose
+ * bitmap array contains bitmap words whose numerical values are
+ * nondecreasing in array order.
+ *
+ * Once we have a solution, we also avoid considering any candidate
+ * which involves dividing one of the output sticks into so many
+ * fragment that the smallest fragment would necessarily be no bigger
+ * than our best solution.  That is, we reject candidates where any of
+ * the hamming weights of the adjacency bitmap words are too large.
+ *
+ * And, we want to do the search in order of increasing maximum
+ * hamming weight.  This is because in practice optimal solutions tend
+ * to have low hamming weight, and having found a reasonable solution
+ * early allows us to eliminate a lot of candidates without doing the
+ * full LP.
+ */
+
 typedef uint32_t AdjWord;
 #define PRADJ "08"PRIx32
 
@@ -45,9 +110,17 @@ static int count_set_adj_bits(AdjWord w) {
 }
 
 static void optimise(int doprint) {
+  /* Consider the best answer (if any) for a given adjacency matrix */
   glp_prob *prob = 0;
   int i, j, totalfrags;
 
+  /*
+   * Up to a certain point, optimise() can be restarted.  We use this
+   * to go back and print the debugging output if it turns out that we
+   * have an interesting case.  The HAVE_PRINTED macro does this: its
+   * semantics are to go back in time and make sure that we have
+   * printed the description of the search case.
+   */
 #define HAVE_PRINTED ({						\
       if (!doprint) { doprint = 1; goto retry_with_print; }	\
     })
@@ -74,6 +147,10 @@ static void optimise(int doprint) {
     }
   }
   if (!had_max) {
+    /* Skip this candidate as its max hamming weight is lower than
+     * we're currently looking for (which means we must have done it
+     * already).  (The recursive iteration ensures that none of the
+     * words have more than the max hamming weight.) */
     PRINTF(" nomaxham");
     goto out;
   }
@@ -128,8 +205,8 @@ static void optimise(int doprint) {
   int ME_totals_j__minimum = next_matrix_entry;
   for (j=0; j<m; j++) ADD_MATRIX_ENTRY(Y_totals_j+j, X_minimum);
 
-  /* \forall_i x_totals_i = m */
-  /* \forall_i x_totals_j = n */
+  /* \forall_i x_total_i = m */
+  /* \forall_i x_total_j = n */
   for (i=0; i<n; i++) glp_set_row_bnds(prob, Y_totals_i+i, GLP_FX, m,m);
   for (j=0; j<m; j++) glp_set_row_bnds(prob, Y_totals_j+j, GLP_FX, n,n);