032: Maximum density still life

/*

  Maximum density still life in Comet.

  CSPLib 032: http://www.csplib.org/Problems/prob032


  This model was inspired by the OPL model from
  Toni Mancini, Davide Micaletto, Fabio Patrizi, Marco Cadoli: 
  "Evaluating ASP and commercial solvers on the CSPLib"
  http://www.dis.uniroma1.it/~tmancini/index.php?problemid=032&solver=OPL&spec=BASE&currItem=research.publications.webappendices.csplib2x.problemDetails#listing


  This Comet model was created by Hakan Kjellerstrand (hakank@gmail.com)
  Also, see my Comet page: http://www.hakank.org/comet 

 */

// Licenced under CC-BY-4.0 : http://creativecommons.org/licenses/by/4.0/

import cotfd;

int t0 = System.getCPUTime();

int size = 6; // to change

range objFunctionBoardCoord      = 2..size+1;
range checkConstraintsBoardCoord = 1..size+2;
range augmentedBoardCoord        = 0..size+3;


Solver<CP> m();

// Search space: The set of all possible assignments of 0s (dead) and 1s (live) 
// to the cells of the board section. However, to be able to easily express 
// constraints on "boundary" cells, we take as search space the set of 0/1 
// boards of size n+4 by n+4: the actual stable pattern appears in the sub-board 
// defined by ignoring the first/last two rows/columns.
var<CP>{int} grid[augmentedBoardCoord,augmentedBoardCoord](m, 0..1);

// Objective function: Maximize the number of live cells in the sub-board defined 
// by ignoring the first/last two/ rows/columns.

Integer num_solutions(0);

// exploreall<m> {
maximize<m> 
   sum(r in objFunctionBoardCoord, c in objFunctionBoardCoord) (grid[r,c])
subject to {

  // C1: Cells in the first/last two rows/columns are all 0 (dead)
  forall(x in augmentedBoardCoord) {
    m.post(grid[0,x] == 0);
    m.post(grid[1,x] == 0);
    m.post(grid[size+2,x] == 0);  
    m.post(grid[size+3,x] == 0);
    m.post(grid[x,0] == 0);       
    m.post(grid[x,1] == 0);
    m.post(grid[x,size+2] == 0);  
    m.post(grid[x,size+3] == 0);
  }
  
  forall(r in checkConstraintsBoardCoord,c in checkConstraintsBoardCoord) { 
    // C2: Each cell of the board (except those of the first/last row/column) 
    //     that has exactly three live neighbors is alive. 
    //     Together with constraint C1, this implies that cells in the
    //     second/last-but-one row/column cannot have three live neighbors.
    m.post(( ( grid[r-1,c-1] + grid[r-1,c] + grid[r-1,c+1] + 
        grid[r,c-1] + grid[r,c+1] + 
        grid[r+1,c-1] + grid[r+1,c] + grid[r+1,c+1]
        ) == 3 
             ) => (grid[r,c] == 1));
    
    // C3: Each live cell must have 2 or 3 live neighbors (cells of the first/last 
    // row/column may be ignored by this constraint)
    m.post((grid[r,c] == 1) => (
                        2 <= 
                        ( grid[r-1,c-1] + grid[r-1,c] + grid[r-1,c+1] +
                          grid[r,c-1] + grid[r,c+1] +
                          grid[r+1,c-1] + grid[r+1,c] + grid[r+1,c+1] 
                          )
                        && 
                        ( grid[r-1,c-1] + grid[r-1,c] + grid[r-1,c+1] +
                          grid[r,c-1] + grid[r,c+1] +
                          grid[r+1,c-1] + grid[r+1,c] + grid[r+1,c+1] 
                          ) <= 3
                                ));
  }
  

  // SBSO: Symmetry-breaking by selective ordering
  // The assignment is forced to respect an ordering on the values that occur in corner entries
  // of the board. In particular:  
  // - if the NW-corner cell is dead, the SE-corner cell
  // must be dead too 
  // - if the NE-corner cell is dead, the SW-corner cell must be dead too
  // 
  m.post(grid[2,2] >= grid[size+1,size+1]);
  m.post(grid[2,size+1] >= grid[size+1,2]);

} using {
      
  // labelFF(m);

  forall(i in augmentedBoardCoord,  j in augmentedBoardCoord: !grid[i,j].bound()) {
   label(grid[i,j]);
  }

  num_solutions := num_solutions + 1;

  forall(i in augmentedBoardCoord) {
    forall(j in augmentedBoardCoord) {
      cout << grid[i,j] << " ";
    }
    cout << endl;
  }
  cout << endl;


}

cout << "\nnum_solutions: " << num_solutions << endl;
    
int t1 = System.getCPUTime();
cout << "time:      " << (t1-t0) << endl;
cout << "#choices = " << m.getNChoice() << endl;
cout << "#fail    = " << m.getNFail() << endl;
cout << "#propag  = " << m.getNPropag() << endl;