Genetically Seeking Sparse Rulers

Mathematica provides powerful computational and visualization tools for the study of a genetic algorithm.

Using an n-inch ruler with markings at every inch, one can measure any integral number of inches from 1 to n. But we can quickly see that not all of those marks are necessary. Rulers with minimal or nearly minimal numbers of markings are called “sparse rulers.” This notebook contains the code for finding sparse rulers using a genetic algorithm.

The general procedure in a genetic algorithm is to work iteratively toward a solution to a problem. At any given iteration, a family of possible solutions is evaluated for their “fitness,” that is, how close they come to satifying the conditions of the problem. The algorithm selects the fittest solutions and then recombines these solutions to produce a new generation of solutions. In this way a solution to the problem “evolves.”

References to genetic algorithms [Holland 1992], genetic programming [Koza 1992], and further general conepts of evolution and order [Kauffman 1993] can be found at the conclusion of this notebook.

Seeking sparse rulers to measure all the numbers from 0 to 36

The sparse ruler problem is ideal for illustrating the process of genetic solution. In the following notebook, a ruler for the number p is simply represented as a list beginning with 0 and ending with p. The case for p = 36 is interesting, because there is exactly one minimal sparse ruler for this value.

Defining a fitness function

By fixing the seed for the pseudorandom number generator, we can guarantee a specific outcome. This value can be changed, or the function left unevaluated to get new results.

SeedRandom[1]

As an example, consider the ruler for p = 36 for the following arbitrarily selected values:

p = 36;

r = {0, 1, 3, 4, 13, 23, 29, 36};

The fitness of this ruler is measured by comparing the number of marks in the ruler with the number of integers the ruler will measure. The number of marks in the ruler (counting 0) is just the Mathematica Length function. The set of all numbers that the ruler will measure is found with the defined function span. One way to find this set is to map the absolute value function across all differences generated by a given ruler list.

absDiff[{s_, t_}] := Abs[s - t];
span[s_List] := 
      Union[absDiff /@ Flatten[Outer[List, s, s], 1]]

For the example r,

span[r]

The function fitness can now be defined as a function that will increase when either the number of marks in the ruler decreases or the number of integers measured by the ruler increases.

fitness[
        s_List] := ((p + 1) - Length[s])/((p + 2) - 
            Length[span[s]])

For the example given here, the fitness value is about 2.

fitness[r]

A ruler containing marks at every point from 0 to 36 would have a fitness equal to 0. We know that a sparse ruler with ten marks exists for p = 36. The length of the span of this ruler would be 37, and the resulting fitness value would be 27.

Producing a new generation of solutions.

Through a rather loose usage of biological terminology, our rulers can be thought of as “chromosomes,” and each number in the ruler set as a “gene.” The function donate is a “crossing-over” function that pastes together parts of two rulers. This function occasionally leaves out some integer common to both rulers, resulting in a “mutation by omission.” The function newGene will result in “mutation by insertion” of a new integer. The rate is fixed at a rather high value, 0.7 in this case. The function mate combines both types of mutation. This function also guarantees that the values 0 and p will be retained in any new combination. Finally, the offspring function generates a new set of n rulers by “mating” a given pair of parent rulers.

donate[s_List, t_List] := 
      Block[{ri1 = Random[Integer, Min[Length[s], Length[t]]], 
          ri2 = Min[ri1 + Round[Random[Real, 0.7]], Length[t]]}, 
        Union[Take[s, ri1], Drop[t, ri2]]]

newGene[] := Random[Integer, p]*Round[Random[Real, 0.7]]

mate[{s_, t_}] := Union[{newGene[]}, {0, p}, donate[s, t]]

offspring[{s_, t_}, n_] := Table[mate[{s, t}], {n}]

The example below shows a generation of five new rulers resulting from the combination of two parent rulers. The fitness function tells us which two to select.

r1 = {0, 1, 2, 15, 29, 36};
r2 = {0, 3, 7, 9, 31, 33, 36};

offspring[{r1, r2}, 5]

Map[fitness, %] // N

Evolving a solution through multiple generations of rulers.

The scheme is to begin with two original rulers, labeled alphaRuler and betaRuler, then recursively generate new rulers in population sizes equal to 5. (This value, like the mutation rate of 0.7 above, was arrived at by experimenting. We will return to a discussion of these values later.) The two fittest rulers from each generation are indexed as newAlpha[i], and newBeta[i], where fitness[newAlpha[i]] ≤ fitness[newBeta[i]]. Note the use of delayed assignment (:=) followed by immediate assignment (=) in setting up the recursion.

alphaRuler = {0, p};
betaRuler = {0, p};

rulerGeneration[0] = {alphaRuler, betaRuler};

rulerGeneration[n_] := rulerGeneration[n] = 
          offspring[
            Take[Sort[rulerGeneration[n - 1], 
                fitness[#1] < fitness[#2] &], -2], 5];

newAlphaRuler[n_] := newAlphaRuler[n] = 
        First[Take[
            Sort[rulerGeneration[n], 
              fitness[#1] < fitness[#2] &], -2]]

newBetaRuler[n_] := newBetaRuler[n] = 
        Last[Take[
            Sort[rulerGeneration[n], 
              fitness[#1] < fitness[#2] &], -2]]

Since the ruler labeled newBetaRuler[i] will have the best fitness value, a plot of these values will indicate how well the generations of rulers are evolving toward a solution.

ListPlot[Table[fitness[newBetaRuler[i]], {i, 1, 20}], 
        PlotJoined -> True];

The graph shows that the rulers are evolving in terms of increasing fitness. The evolution is obviously not monotonic, however. The best ruler in this set of generations is newBetaRuler[19].

newBetaRuler[19]

The ruler is reasonably short.

Length[%]

The ruler will measure all but seven of the numbers between 0 and 36.

Length[span[newBetaRuler[19]]]

We can see from the fitness value that there is still considerable room for improvement. (Recall that we should be able to find fitness values up to 27.)

fitness[newBetaRuler[19]]

This fitness function seems to select quickly the rulers that will measure all or most of the numbers between 0 and a given value p. Here is a typical plot of 100 generations.

ListPlot[Table[fitness[newBetaRuler[i]], {i, 1, 100}], 
        PlotJoined -> True];

Here is a selection of eleven rulers from this set of generations. These will measure the complete set of numbers between 0 and 36.

testSet  =  Table[newBetaRuler[i], {i, 50, 60}]

Map[Length, Map[span, testSet]]

Map[fitness, testSet]

Union[testSet]

The above list, with sequences of identical rulers, is typical. The evolutionary process goes through periods of stability and periods of fluctuation, as we can see by examining a plot of 1,000 generations.

fitnessTable  =  
        Table[fitness[newBetaRuler[i]], {i, 1000}];

fitnessPlot  =  
        ListPlot[fitnessTable,  PlotJoined -> True,  
          PlotRange -> {0, 29}];

An alternative visualization

The next two functions give us a different "picture" of a given ruler.

bitPlane  =  Table[0, {i, 37}];

\(\(convertToBin[x_List] := 
          ReplacePart[bitPlane, 1, \ 
            Partition[Plus[1, x], 1]];\)\(\ \)\(
\)\)

To create an image of the rulers themselves, we first convert them to binary strings. Thus the ruler {0, 29, 36} would have the respective binary representation:

convertToBin[{0,  22,  36}]

We can use the Mathematica ListDensityPlot[15] function to transform numeric lists into “actual” rulers.

For example, if we view the first 15 best rulers from a particular series using the Mathematica ListDensityPlot[15] function, we obtain this picture:

This plot shows us at a glance that the first ruler is {0, 22, 36}; also, that the mark 24 is found on each of the second through the 15th rulers. The marks on any other ruler in the set, and the persistance of marks across several generations is easily visualized.

If we turn this plot sideways, we can match each ruler to its fitness value:

lpf15  = 
        ListPlot[Table[fitness[newAlphaRuler[i]], {i, 15}], 
          PlotJoined -> True];

ldp15  =  
        ListDensityPlot[-Transpose[ 
              Map[convertToBin, 
                Table[newBetaRuler[i],  {i, 15}]]]];

Show[
        GraphicsArray[{{ldp15}, {lpf15}}, 
          GraphicsSpacing -> 0.0]];

For an entire set of 1000 generations, we have the following picture:

ldp1000 = 
        ListDensityPlot[-Transpose[
              convertToBin /@ Table[newBetaRuler[i], {i, 1000}]], 
          Mesh -> False];

Show[
        GraphicsArray[{{ldp1000}, {fitnessPlot}}, 
          GraphicsSpacing -> 0.0]];

Exploring the fitness landscape for small values of p

For small values of p, we can look at an exhaustive list of rulers for the number p. A convenient way to generate these rulers is by using the GrayCode function from the DiscreteMath`Combinatorica` package. The GrayCode function generates a collection of sets such that adjacent sets in the list differ by only one element. From this list of sets the function rulerList constructs a set of rulers for the number p, such that two adjacent rulers in the list differ by only one integer.

Needs["DiscreteMath`Combinatorica`"]

setList[i_] := setList[i] = GrayCode[Range[i - 1]];

rulerList[i_] := rulerList[i] = 
          Map[Union[#, {0, i}] &, setList[i]];

rulerList[4]

The fitness values of the set of all rulers for a given p can be plotted as a “landscape” that shows how the fitness varies across a set of rulers. For example, when p = 10, the maximum fitness is 5. The plot below shows the fitness values for the (2^9) rulers arranged so that each ruler differs from its neighbor by one element.

p = 10;

fitnessList = Map[fitness, rulerList[10]];

landScape = 
        ListPlot[fitnessList, PlotJoined -> True, 
          PlotRange -> All];

The table also shows that for p = 10 there are multiple rulers with the maximum fitness value. The Position function shows the location of these rulers. The list of locations can be used as an index to show all the sparse rulers for p = 10.

bestRulers  =  Position[fitnessList, Max[fitnessList]]

sparserulers = rulerList[
            10][[Flatten[bestRulers]]];

Running a new experiment for p = 10 and plotting the fitness values of the best ruler at each generation shows that the evolution of rulers for this value of p can reach a steady state at about the 20th generation. We need to clear results of previous experiments. A simple way to do this is to clear and then reinitialize the following two functions:

Clear[rulerGeneration,  newBetaRuler]

[Reinitialize the functions rulerGeneration and newBetaRuler]

betaRulerList = Table[newBetaRuler[i], {i, 0, 100}];

ListPlot[Map[fitness, betaRulerList], PlotJoined -> True, 
        PlotRange -> All];

One can show how the evolving rulers migrate into the fitness landscape. betaRulerListReduced produces an abbreviated list of the best rulers, removing successive identical rulers while maintaining the order in which the rulers have evolved.

\(betaRulerListReduced = 
	
      Block[{bLR = {newBetaRuler[0]}}, 
		Do[
          If[betaRulerList[[i]] != 
              betaRulerList[[i - 1]], 
            bLR = Append[bLR, 
                betaRulerList[[
                  i]]]], {i, 1, Length[betaRulerList]}]; 
        bLR]\)

The set sites list the horizontal values of the positions of the rulers, while sitePoints pairs each horizontal value with its corresponding fitness. siteLines traces the evolution of the rulers through the landscape.

sites = 
      Flatten[Table[
          Position[rulerList[10], 
            betaRulerListReduced[[
              i]]], {i, 1, 
            Length[betaRulerListReduced]}]]

sitePoints = 
        Graphics[{PointSize[0.02], 
            Point /@ Table[{sites[[i]], 
                  fitness[
                    betaRulerListReduced[[
                      i]]]}, {i, 1, 
                  Length[betaRulerListReduced]}]}];

siteLines = 
        ListPlot[
          Table[{sites[[i]], 
              fitness[betaRulerListReduced[[
                  i]]]}, {i, 1, 
              Length[betaRulerListReduced]}], 
          PlotJoined -> True];

Lightening the color of the landscape makes it easier to distinguish the evolutionary path.

landscapeGray = 
        ListPlot[fitnessList, PlotJoined -> True, 
          PlotStyle -> GrayLevel[0.5], 
          PlotRange -> All];

Show[{landscapeGray}, {siteLines, sitePoints}];

Animation of the genetic algorithm

Finally, we have an animation showing the fitness values of the evolving rulers moving through the landscape.

Table[
        Show[{landscapeGray, siteLines, 
            Graphics[{PointSize[0.02], 
                Point[{sites[[i]], 
                    fitness[
                      betaRulerListReduced[[
                        i]]]}]}]}], {i, 1, 
          Length[betaRulerListReduced]}];

Conclusion

A more general way of thinking about sparse rulers and the fitness landscape is to define the rulers as binary strings. Each mark i on a ruler is represented by a 1 in the ith position of the binary string. For example, the ruler {0, 1, 3, 7, 8, 10} would be represented by the string {1, 1, 0, 1, 0, 0, 0, 1, 1, 0, 1}. Since the rulers to measure all numbers up to the value p always include the numbers 0 and p, the binary representations of the rulers always include a 1 in the first and pth positions. The fitness function can therefore be thought of as being defined on a hypercube of dimension p - 1. The genetic evolution of rulers corresponds to a tour through this hypercube.

Since the values of (2^(p - 1)) rapidly become too large, the entire fitness landscape for large values of p cannot be pictured. However, “sketches” of the landscape can provide useful information. For example, every fitness landscape will include the fitness value of the ruler containing all possible marks. This is the “zero valley” located in the graph above at the point 341 on the horizontal axis. There are also zones that are relatively rich in high fitness values, such as the region from 442 to 502 on the horizontal axis. Instead of always beginning with the zero generation rulers as {0, p}, collections of rulers could be seeded in the most promising zones of the landscape.

The numbers of offspring and the mutation rate could vary with the fitness of the evolving rulers. If generation size and mutation rate were to vary inversely with fitness, then the population should tend to evolve within an increasingly smaller zone. The appropriate metric here is the Hamming distance, which defines the distance between two binary strings of length n as the number of positions in which the two strings are different. Two strings could be a small Hamming distance apart and still not look close together in a graphical representation of the fitness landscape. For example, in the representation of the landscape above, the first fitness measure is for the ruler {0, 10} and the last measure is for the ruler {0, 9, 10}—the rulers at opposite sides of the landscape diagram are only one unit apart in Hamming Distance. This is a result of the way the Gray code generates the list of rulers.

An alternative algorithm could contain a timing device that would periodically raise and lower the mutation rate. Such an algorithm could produce occasional bursts of mutant rulers that would wander into unexplored areas of the landscape. Finally, the fitness function itself might change, in which case the landscape would be different at different stages of evolution.

(2^36)

There are 68,719,476,736 (about 70 billion) possible rulers of length 36, and only one of minimal length 9. Where in hyperspace does this unique solution live? How long would it take to find it? How good are the 36 adjacent rulers? Can the genetic algorithm find this absolute maximum? How good are the 94,143,280 other rulers of length 9? How are the best rulers of length 10, 11, and 12 distributed through hyperspace?

Binomial[36, 9]

References and links

Holland, John H. Adaptation in Natural and Artificial Systems. MIT Press, 1992.

Kauffman, Stuart A. The Origins of Order. Oxford Press, 1993.

Koza, John R. Genetic Programming. MIT Press, 1992.

Skiena, Steven S. The Algorithm Design Manual. TELOS-Springer Verlag, 1998.