I finished a first pass at my lattice regression library over the weekend. The idea with that is pretty straightforward. Essentially, there’s some function we want to model, and it’s unknown, but we have a bunch of observations of inputs and corresponding outputs. So, we throw down a lattice (i.e. a regularly-spaced grid) of points over the space of inputs, and we use the data to “learn” some values of the function at the lattice points. Then, we discard the training data but can predict new values of the function by interpolating between the values at the lattice points. For more details, see, e.g. this paper.

Code-wise, one challenge of the project was in representing and dealing with the lattice. For example, suppose the function we want to model has 4 inputs. Then, our learned values on a grid over the space of inputs might naturally be stored in something like a 4-D array,

double b[p][p][p][p];

where p is the number of grid points in any one direction. Then, if we want to do something to every lattice point, we can do something like

for (int i = 0; i < p; i++) {
  for (int j = 0; j < p; j++) {
    for (int k = 0; k < p; k++) {
      for (int l = 0; l < p; l++) {
        do_something_to(b[i][j][k][l]);
      }
    }
  }
}

Unfortunately, we want our code to be able to model functions with different numbers of inputs. Effectively, the number of dimensions needed in b array and the number of nested for loops necessary to act on it are unknown until runtime. To handle the general case, we need something a bit smarter than the above. Essentially, we need some way to simulate arbitrarily-deeply nested loops. There are good recursive solutions to this problem, but I didn’t go that route. Here are a couple approaches I used in the Lattice Regression library.

1) Put all the lattice points in a long 1-D array

double b[pow(t,d)]

where d is the number of dimensions. Then, loop over it like normal

for (int i = 0; i < pow(t,d), i++) {
  // do cool things that matter
}

If we need to know the “real” coordinates of the current point, we can extract them from the loop counter i by mapping i to the “real” coordinates . If we only care about the real coordinates of points, then any one-to-one mapping from i to real coordinates will work. I happen to like

Then, we can convert from i to the coordinates with something like

int indx = i;
int coords[d];
for (int k = d-1; k >= 0; k--) {
  coords[k] = indx % (p);
  indx /= p;
}

and we can go back again, by just performing the above sum.

Now, if we’re iterating over all points, then we’re essentially counting in base p from 0 to (p-1)(p-1)...(p-1), (where each (p-1) is one digit of a base-p number). We can make that more explicit and perhaps make a couple things easier by storing each “digit” separately…

2) Put artifical loop counters into a dynamically-sized array and “count” upwards.

int counters[d];
memset(counters, 0, d*sizeof(int));

We start out with each loop counter 0. Then we increment the first one until we reach its maximum value (p-1) in our case). Then we set it back to zero, increment the second loop counter once, and again increment the first counter until it maxes out. Then we increment the second loop counter once, and pass over the first counter again. Eventually, the second counter maxes out also, and we’ve covered all possible combinations of the first two counters once. We then increment the third counter once, reset the first and second, and repeat all combinations of the first two. And so on.

In pseudocode,

while (! all_the_counters_are_maxed_out()) {
  // do some stuff with the "loop indices" counters[0],counters[1],...,counters[d-1]

  // increment counters
  int i = 0;
  while (counters[i] == max_value_of_counter_i)
    counters[i] = 0;
  counters[i]++;
}

The lines here under //increment counters set all maxed-out counters back to 0 and then increment the first not-maxed-out counter found. Of course, the most common case is that the least-significant counter is not maxed out, in which case the body of the while loop doesn’t execute at all.