Sudoku – Random IT Utensils https://blog.adamfurmanek.pl IT, operating systems, maths, and more. Fri, 19 Feb 2021 07:28:34 +0000 en-US hourly 1 https://wordpress.org/?v=6.6.2 ILP Part 65 – Fancy sudoku https://blog.adamfurmanek.pl/2021/05/08/ilp-part-65/ https://blog.adamfurmanek.pl/2021/05/08/ilp-part-65/#comments Sat, 08 May 2021 08:00:17 +0000 https://blog.adamfurmanek.pl/?p=3876 Continue reading ILP Part 65 – Fancy sudoku]]>

This is the sixtieth fifth part of the ILP series. For your convenience you can find other parts in the table of contents in Part 1 – Boolean algebra

Today we’re going to solve a modification of classical Sudoku riddle shown in Precz z cyframi.

We need to fill the board with digits as in Sudoku and meet some additional constraints:

  • Field with X indicates that number in that field is equal to the number of different values in corner neighbours of the field (so top-left field, top-right, bottom-left and bottom-right)
  • Field with + indicates that among four straight neighbours (field above, below, to the left, to the right) there are two which added values are equal to the value in the field
  • Field with circle indicates that field is equal to the number of different values in all neighbours

Additionally, the board is filled with all crosses, pluses and circles completely. This means that if given field is empty then no cross/field/circle can be placed in it.

Solution:

var board = new int[][]{
	new int []{0,0,1,2,3,0,2,0,0},
	new int []{0,2,0,0,4,3,0,2,2},
	new int []{0,2,0,1,2,0,0,0,0},
	new int []{1,2,0,0,2,1,0,0,0},
	new int []{2,0,0,2,2,0,2,1,3},
	new int []{0,1,1,2,0,2,0,2,0},
	new int []{4,0,2,1,0,0,2,4,0},
	new int []{0,2,0,0,1,2,1,2,1},
	new int []{4,1,0,0,2,0,3,0,0}
};

var width = board[0].Length;
var height = board.Length;

var solver = new CplexMilpSolver(new CplexMilpSolverSettings
{
	IntegerWidth = 15,
	Epsilon = 0.1,
	StoreDebugExpressions = false,
	CacheConstants = false,
	StoreDebugConstraints = false
});

Console.WriteLine("Create variables");
var fields = Enumerable.Range(0, height).Select(r => Enumerable.Range(0, width).Select(c => solver.CreateAnonymous(Domain.PositiveOrZeroInteger)).ToArray()).ToArray();

Console.WriteLine("Set ranges");
for(int row = 0; row < height;++row){
	for(int column = 0; column < width; ++ column){
		solver.Set<LessOrEqual>(fields[row][column], solver.FromConstant(9));
		solver.Set<GreaterOrEqual>(fields[row][column], solver.FromConstant(1));
	}
}

Console.WriteLine("Rows");
for(int row = 0; row < height;++row){
	solver.Set<AllDifferent>(fields[row][0], fields[row].Skip(1).ToArray());
}

Console.WriteLine("Columns");
for(int column = 0; column < width; ++ column){
	solver.Set<AllDifferent>(fields[0][column], Enumerable.Range(1, height-1).Select(row => fields[row][column]).ToArray());
}

Console.WriteLine("Quadrants");
for(int row = 0; row < height / 3;++row){
	for(int column = 0; column < width / 3; ++ column){
		var quandrant = Enumerable.Range(0, 3).SelectMany(i => Enumerable.Range(0, 3).Select(j => fields[row * 3 + i][column * 3 + j])).ToArray();
		solver.Set<AllDifferent>(quandrant[0], quandrant.Skip(1).ToArray());
	}
}

Console.WriteLine("Shapes");
for(int row = 0; row < height;++row){
	for(int column = 0; column < width; ++ column){
		var corners = Enumerable.Range(-1, 3).SelectMany(i => Enumerable.Range(-1, 3)
		.Where(j => i % 2 != 0 && j % 2 != 0)
		.Where(j => row + i >= 0 && row + i < height)
		.Where(j => column + j >= 0 && column + j < width)
		.Select(j => 
			fields[row+i][column+j]
		)).ToArray();
		
		var straights = Enumerable.Range(-1, 3).SelectMany(i => Enumerable.Range(-1, 3)
		.Where(j => (i % 2 == 0 && j % 2 != 0) || (i % 2 != 0 && j % 2 == 0))
		.Where(j => row + i >= 0 && row + i < height)
		.Where(j => column + j >= 0 && column + j < width)
		.Select(j => 
			fields[row+i][column+j]
		)).ToArray();
		
		var neighbours = corners.Concat(straights).ToArray();
		
		var allSums = Enumerable.Range(0, straights.Length).SelectMany(i => Enumerable.Range(0, straights.Length).Where(j => j > i).Select(j => solver.Operation<Addition>(straights[i], straights[j]))).ToArray();
		var isEqual = solver.Operation<Addition>(allSums.Select(s => solver.Operation<IsEqual>(fields[row][column], s)).ToArray());
		
		switch(board[row][column]){
			case 0:
				solver.Set<NotEqual>(solver.Operation<DifferentValuesCount>(corners), fields[row][column]);
				solver.Set<Equal>(isEqual, solver.FromConstant(0));
				solver.Set<NotEqual>(solver.Operation<DifferentValuesCount>(neighbours), fields[row][column]);
			break;
			case 1:
				solver.Set<Equal>(solver.Operation<DifferentValuesCount>(corners), fields[row][column]);
				solver.Set<Equal>(isEqual, solver.FromConstant(0));
				solver.Set<NotEqual>(solver.Operation<DifferentValuesCount>(neighbours), fields[row][column]);
			break;
			case 2:
				solver.Set<NotEqual>(solver.Operation<DifferentValuesCount>(corners), fields[row][column]);
				solver.Set<GreaterOrEqual>(isEqual, solver.FromConstant(1));
				solver.Set<NotEqual>(solver.Operation<DifferentValuesCount>(neighbours), fields[row][column]);
			break;
			case 3:
				solver.Set<NotEqual>(solver.Operation<DifferentValuesCount>(corners), fields[row][column]);
				solver.Set<Equal>(isEqual, solver.FromConstant(0));
				solver.Set<Equal>(solver.Operation<DifferentValuesCount>(neighbours), fields[row][column]);
			break;
			case 4:
				solver.Set<NotEqual>(solver.Operation<DifferentValuesCount>(corners), fields[row][column]);
				solver.Set<GreaterOrEqual>(isEqual, solver.FromConstant(1));
				solver.Set<Equal>(solver.Operation<DifferentValuesCount>(neighbours), fields[row][column]);
			break;
		}
	}
}


solver.AddGoal("goal", solver.FromConstant(0));
solver.Solve();

for(int row = 0; row < height;++row){
	if(row % 3 == 0)Console.WriteLine("-------------");
	for(int column = 0; column < width; ++ column){
		if(column % 3 == 0)Console.Write("|");
		Console.Write(Math.Round(solver.GetValue(fields[row][column])));
	}
	Console.Write("|");
	Console.WriteLine();
}
Console.WriteLine("-------------");

Lines 25-52 are regular Sudoku rules.

Next, in lines 54-106 we encode additional rules. First, we extract corner neighbours (lines 57-63), straight neighbours (65-71), and all neighbours together (line 73). Next, we prepare sums of pairs of straight neighbours (line 75) and then compare them with current field value.

Finally, in lines 78-103 we encode the rules. If field is of type 0 then it means field is empty so number of different values for corner neighbours must not be equal to the current value, number of straight neighbours summing to the current field value is zero (so no sum matches), and similarly for all neighbours – number of different values among them does not match current value.

Type 1 indicates X, type 2 is +, type 3 is circle, type 4 is + with circle. Notice that X with circle is not possible, nor is X with + with circle.

Output:

Tried aggregator 2 times.
MIP Presolve eliminated 5032 rows and 2922 columns.
MIP Presolve modified 32939 coefficients.
Aggregator did 915 substitutions.
Reduced MIP has 22057 rows, 10412 columns, and 68119 nonzeros.
Reduced MIP has 10331 binaries, 81 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.06 sec. (63.34 ticks)
Probing fixed 2358 vars, tightened 75 bounds.
Probing changed sense of 858 constraints.
Probing time = 0.41 sec. (93.68 ticks)
Cover probing fixed 0 vars, tightened 803 bounds.
Tried aggregator 2 times.
MIP Presolve eliminated 6449 rows and 2938 columns.
MIP Presolve modified 8576 coefficients.
Aggregator did 854 substitutions.
Reduced MIP has 14754 rows, 6620 columns, and 46066 nonzeros.
Reduced MIP has 6553 binaries, 67 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.13 sec. (100.36 ticks)
Probing fixed 148 vars, tightened 7 bounds.
Probing changed sense of 23 constraints.
Probing time = 0.16 sec. (42.71 ticks)
Cover probing fixed 0 vars, tightened 41 bounds.
Tried aggregator 2 times.
MIP Presolve eliminated 532 rows and 257 columns.
MIP Presolve modified 1534 coefficients.
Aggregator did 64 substitutions.
Reduced MIP has 14158 rows, 6299 columns, and 44200 nonzeros.
Reduced MIP has 6232 binaries, 67 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.09 sec. (57.17 ticks)
Probing fixed 6232 vars, tightened 112 bounds.
Probing changed sense of 39 constraints.
Probing time = 0.75 sec. (187.24 ticks)
Clique table members: 18414.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 8 threads.
Root relaxation solution time = 0.02 sec. (5.76 ticks)

        Nodes                                         Cuts/
   Node  Left     Objective  IInf  Best Integer    Best Bound    ItCnt     Gap

*     0     0      integral     0        0.0000        0.0000        0    0.00%
Elapsed time = 1.88 sec. (690.31 ticks, tree = 0.00 MB, solutions = 1)

Root node processing (before b&c):
  Real time             =    1.88 sec. (690.79 ticks)
Parallel b&c, 8 threads:
  Real time             =    0.00 sec. (0.00 ticks)
  Sync time (average)   =    0.00 sec.
  Wait time (average)   =    0.00 sec.
                          ------------
Total (root+branch&cut) =    1.88 sec. (690.79 ticks)
-------------
|962|753|814|
|451|286|739|
|783|491|526|
-------------
|295|674|183|
|817|932|645|
|634|518|297|
-------------
|548|327|961|
|176|849|352|
|329|165|478|
-------------

]]> https://blog.adamfurmanek.pl/2021/05/08/ilp-part-65/feed/ 1 ILP Part 30 — Sudoku https://blog.adamfurmanek.pl/2016/03/12/ilp-part-30/ https://blog.adamfurmanek.pl/2016/03/12/ilp-part-30/#comments Sat, 12 Mar 2016 09:00:17 +0000 https://blog.adamfurmanek.pl/?p=1580 Continue reading ILP Part 30 — Sudoku]]>

This is the thirtieth part of the ILP series. For your convenience you can find other parts in the table of contents in Part 1 – Boolean algebra

Hello! Today we are going to solve Sudoku using ILP.

Description

In Sudoku we need to fill 9 \times 9 grid with numbers. In every column, in every row, and in nine smaller subgrids we can use every number from 1 to 9 exactly once. We usually start with few fields filled and we need to find the rest. See the picture below (from wikipedia):
Sudoku

In this post we will use basic set operations a lot so make sure that you understand them.

Variables

Let’s start with variables for the problem. We usually use binary variables to indicate decision, however, this time we will use integer variables in range [1; 9] to represent the solutions. It is quite intuitive: after solving the problem we obtain the result by checking the value of variable for given board field. So we define the following variables:

    \begin{gather*} \displaystyle\mathop{\forall}_{1 \le i \le 9}\displaystyle\mathop{\forall}_{1 \le j \le 9} v_{i,j} \end{gather*}

We also make sure that they are in the correct range:

    \begin{gather*} \displaystyle\mathop{\forall}_{1 \le i \le 9}\displaystyle\mathop{\forall}_{1 \le j \le 9} 1 \le v_{i,j} \le 9 \end{gather*}

Let’s now define constraints.

Constraints

This part is very easy. We already know how to define “all different” constraint, so all we need to do is to add constraints for every row, every column, and every subgrid:

    \begin{gather*} \displaystyle\mathop{\forall}_{1 \le i \le 9} \text{ AllDifferent}\left( v_{i, 1}, v_{i, 2}, \ldots, v_{i, 9} \right) \\ \displaystyle\mathop{\forall}_{1 \le j \le 9} \text{ AllDifferent}\left( v_{1, j}, v_{2, j}, \ldots, v_{9, j} \right) \\ \displaystyle\mathop{\forall}_{0 \le x \le 2} \displaystyle\mathop{\forall}_{0 \le y \le 2} \text{ AllDifferent}\left( v_{3x+1, 3y+1}, \ldots, v_{3x+3, 3y+3} \right) \end{gather*}

Code

var board = new int[][]{
	new int []{5,3,0,0,7,0,0,0,0},
	new int []{6,0,0,1,9,5,0,0,0},
	new int []{0,9,8,0,0,0,0,6,0},
	new int []{8,0,0,0,6,0,0,0,3},
	new int []{4,0,0,8,0,3,0,0,1},
	new int []{7,0,0,0,2,0,0,0,6},
	new int []{0,6,0,0,0,0,2,8,0},
	new int []{0,0,0,4,1,9,0,0,5},
	new int []{0,0,0,0,8,0,0,7,9}
};

var width = board[0].Length;
var height = board.Length;

var solver = new CplexMilpSolver(new CplexMilpSolverSettings
{
	IntegerWidth = 15,
	Epsilon = 0.1,
	StoreDebugExpressions = false,
	CacheConstants = false,
	StoreDebugConstraints = false
});

Console.WriteLine("Create variables");
var fields = Enumerable.Range(0, height).Select(r => Enumerable.Range(0, width).Select(c => solver.CreateAnonymous(Domain.PositiveOrZeroInteger)).ToArray()).ToArray();

Console.WriteLine("Set ranges");
for(int row = 0; row < height;++row){
	for(int column = 0; column < width; ++ column){
		solver.Set<LessOrEqual>(fields[row][column], solver.FromConstant(9));
		solver.Set<GreaterOrEqual>(fields[row][column], solver.FromConstant(1));
	}
}

Console.WriteLine("Set known");
for(int row = 0; row < height;++row){
	for(int column = 0; column < width; ++ column){
		if(board[row][column] != 0)solver.Set<Equal>(fields[row][column], solver.FromConstant(board[row][column]));
	}
}

Console.WriteLine("Rows");
for(int row = 0; row < height;++row){
	solver.Set<AllDifferent>(fields[row][0], fields[row].Skip(1).ToArray());
}

Console.WriteLine("Columns");
for(int column = 0; column < width; ++ column){
	solver.Set<AllDifferent>(fields[0][column], Enumerable.Range(1, height-1).Select(row => fields[row][column]).ToArray());
}

Console.WriteLine("Quadrants");
for(int row = 0; row < height / 3;++row){
	for(int column = 0; column < width / 3; ++ column){
		var quandrant = Enumerable.Range(0, 3).SelectMany(i => Enumerable.Range(0, 3).Select(j => fields[row * 3 + i][column * 3 + j])).ToArray();
		solver.Set<AllDifferent>(quandrant[0], quandrant.Skip(1).ToArray());
	}
}

solver.AddGoal("goal", solver.FromConstant(0));
solver.Solve();

for(int row = 0; row < height;++row){
	if(row % 3 == 0)Console.WriteLine("-------------");
	for(int column = 0; column < width; ++ column){
		if(column % 3 == 0)Console.Write("|");
		Console.Write(solver.GetValue(fields[row][column]));
	}
	Console.Write("|");
	Console.WriteLine();
}
Console.WriteLine("-------------");

Output:

Tried aggregator 1 time.
MIP Presolve eliminated 7995 rows and 3969 columns.
MIP Presolve modified 23909 coefficients.
All rows and columns eliminated.
Presolve time = 0.02 sec. (10.42 ticks)

Root node processing (before b&c):
  Real time             =    0.02 sec. (10.97 ticks)
Parallel b&c, 8 threads:
  Real time             =    0.00 sec. (0.00 ticks)
  Sync time (average)   =    0.00 sec.
  Wait time (average)   =    0.00 sec.
                          ------------
Total (root+branch&cut) =    0.02 sec. (10.97 ticks)
-------------
|534|678|912|
|672|195|348|
|198|342|567|
-------------
|859|761|423|
|426|853|791|
|713|924|856|
-------------
|961|537|284|
|287|419|635|
|345|286|179|
-------------

Summary

As we can see, our toolkit for building ILP formulas is quite powerfull and we are able to solve Sudoku almost instantly just by translating rules of the game (nine different numbers in each column, nine different numbers in each row, nine different numbers in each subgrid) to mathematical operators.

]]> https://blog.adamfurmanek.pl/2016/03/12/ilp-part-30/feed/ 1