Maximum Density Still Life

From CometPublic

Contents 1 Description of the problem 2 Definition of the model 2.1 The board 2.2 The constraints 2.2.1 Violations 3 The search 3.1 Objective function 3.2 Initial State 3.3 Meta-heuristic 3.4 Heuristic 4 Results 4.1 Example of solution 4.2 Some statistics

[edit] Description of the problem

This problem was proposed by Barbara Smith on CSPLib (prob032). From CSPLib:

"This problem arises from the Game of Life, invented by John Horton Conway in the 1960s and popularized by Martin Gardner in his Scientific American columns."

The purpose of this problem is to find the densest stable pattern, or still-life, that can fit on a section of the board (in the Game of Life, the board is supposed to be infinite in all directions, so we consider the rest of the board to be dead). A still-life is an arrangement of live cells that does not change over time.

[edit] Definition of the model

[edit] The board

The model we propose here is based on an (or more generally an ) array of integers taking values in {0,1} (we did not choose booleans as we did not know their behaviour as array indices or in arithmetic expressions), representing the state of each cell: dead (0) or alive (1). Every border cell (type a) and b) in the below representation) is set to 0 and won't be taken into account in the search, they are only used to simplify the model description and represent the rest of the board: this allows to avoid considering the border cells in a different way.

Model representation
	Legend: a) Border cell, dead and unconstrained b) Border cell, dead but constrained c) Modifiable cell, constrained

Our Comet implementation is the LifeGrid class. It allows to specify the height and width of the board section to consider, and if the user wants to start with an empty (without live cells) or full (all cells initialized alive) board through initNull. This might have an effect on the search and the results as we will see later.

LifeGrid provides several methods for general use (getNumAlive(), getNumDead(), getNumCells()...) or to manipulate the model at a lower level (getLocalSolver(), getGrid() etc.). This allows the user to define additional constraints and invariants before closing the model (through the close() method) and of course to manipulate the model efficiently during the search.

[edit] The constraints

There is only one constraint in this problem, which is defined in this way:

• if a cell has 3 live neighbours, then it must be alive
• if a cell has strictly more than 3 live neighbours or strictly less than 2 live neighbours, then it must be dead

In this model, we have chosen to implement the Constraint interface with a simple constraint StillCell, and apply this constraint on each cell. This simplifies the implementation as we only have to consider one cell and it's neighbours at a time instead of the whole board, without taking into account the board disposition.

This constraint is applied to all cells of type b) and c) (see the Model representation above). Cells of type a) do not need this constraint as they always keep dead as well as their neighbours, and thus always check the constraint.

[edit] Violations

We have chosen to define the violations count as the number of neighbours that should be changed if we wanted to keep _cell in it's current state. For better efficiency and easier implementation, we build once and for all an array that gives the number of violations in function of the state of _cell and the number of live neighbours (_numAliveNeighb):

			_violations = new int[j in 0..1, i in _numNeighb]
					= (j == 1 ? (i <= 2 ? 2 - i : i - 3)
						  : (i == 3 ? 1 : 0));

This allows us to easily define an invariant on the number of violations:

			_violation <- _violations[_cell, _numAliveNeighb];

The variable violations are defined in a similar way, but we have considered that we could count only one violation for each variable. However, decrease(var{int}), getAssignDelta(var{int}, int) and getSwapDelta(var{int}, var{int}) give the true (exact) variation in the number of violations. We have not implemented the other methods as we didn't need them.

[edit] The search

[edit] Objective function

This is an optimization problem: the objective that must be minimized depends on the total number of violations (which counts positively) and the number of live cells (which counts negatively). Thus, the objective function looks like with . Generally, a is greater than b because we want the final solution to have no violations.

We have chosen a = max(2,w * h / 200) and b = 1 because it seemed to give better results from our own experience. However the value for a should perhaps be tweaked for great board sections.

[edit] Initial State

To find the optimal solution, the search considers each cell individually and can toggle its value. If the search starts from the empty configuration then it can be trapped in poor local optima : the initial state has no violations and similar states (with only a few live cells) have also no violations. These states have great objective values and some of them are local optima. To avoid this, we have choosen to start the search with the initial board filled with live cells.

[edit] Meta-heuristic

We use here the variable-depth neighborhood search (VDNS, see Constraint-Based Local Search p. 169). This search explores a sequence of moves and then returns the best solution encountered during the exploration (with respect to the objective function). For this purpose we use the selectBest function (see Constraint-Based Local Search p. 172). Notice that the version presented here differs from the original one presented in the Constraint-based Local Search book in that it uses Solutions instead of Checkpoints. We use Solutions because using events and Checkpoints together don't seem to work.

This VDNS is made several times and starts at each time with the best solution found so far. This is a kind of intensification that works better than doing just one VDNS (with the same computation time).

forall (i in 1..maxit) {
	selectBest(m, 100, obj, step);
}

In addition, a tabu component is added to avoid being trapped in poor local optima.

int tabu[Height, Width] = 0;
Counter it(m, 0);
forall (i in Height, j in Width)
	whenever alive[i, j]@changes(int o, int n) {
		tabu[i, j] = it + 8;
		it++;
	}

Tweaking the tabu length did not seem to give better results.

[edit] Heuristic

A move consists in taking the most violated variable of the board and toggle its value. Once there are no more violations, we can try to toggle the dead cell that least increases the number of violations in order to find a better solution. We do this explicitly because we want to guaranty an increase in the number of live cells, even if it is worse for the number of violations. Else it could lead to a live cell becoming dead.

Notice that we do not implement the heuristic by simply looking at the objective function, because it would never guaranty increases in the number of live cells while making the search very dependant on the chosen weights.

Closure step = closure {
	if (S.violations() == 0)
		selectMin(i in Height, j in Width : tabu[i, j] <= it && alive[i, j] == 0)
			(S.getAssignDelta(alive[i, j], 1))
			alive[i, j] := 1;
 	else
		selectMax(i in Height, j in Width : tabu[i, j] <= it)
			(S.violations(alive[i, j]))
			alive[i, j] := 1 - alive[i, j];
};

The entire implementation of the search can be found here.

[edit] Results

[edit] Example of solution

Here is an example of solution for a section of the board: (a "[]" represents one live cell)

violations : 0; num alive : 48
[][]  []  []    [][]
  []  [][]  []    []
[]          [][][]  
[][][][][][]        
            []      
[][]  [][][]    [][]
[]  [][]    [][]  []
          []        
  [][][][]  [][]  []
  []    []  []  [][]

[edit] Some statistics

The following table presents the number of live cells and violations over samples of five solutions.

Board size	live cells			violations
	min	max	mean	(mean)
7x7	24	28	26.2	0
10x10	48	51	49.4	0.2
13x13	81	88	83.4	0.2
15x15	107	114	110.6	0