# Sudoku¶

This tutorial explains how to solve a Sudoku puzzle using an Ising machine with Amplify.

## About Sudoku¶

Sudoku is a puzzle that fills a $9\times9$ grid with the numbers from $1$ to $9$ accordingly to the following rules:

- Place one number from $1$ to $9$ to each cell.
- No identical numbers may be placed in each row and column, and in each of the nine $3\times3$ subgrids.

First, we are given initial clues, where 17 or more numbers are placed (it has been proven that no solution can be found with 16 or less clues). The game is played by determining the numbers that go into the unfilled cells accordingly to the rules above. When the game's difficulty level is low, it is relatively easy to determine which numbers to place. As the difficulty of the game increases, and it requires some experience to solve the puzzles.

Various algorithms have been devised to solve Sudoku puzzles using computers, such as depth-first search, probabilistic methods, constraint satisfaction problems, exact cover problems, etc. By following these algorithms, Sudoku puzzles can be solved mechanically.

In this tutorial, we will show how to solve Sudoku puzzles using an Ising machine. Formulating the problem from the above Sudoku rules, we define the corresponding function, where finding its minimum value is equivalent to solving the Sudoku puzzle. In the case of Sudoku, only constraint conditions are needed, and no cost function appears.

Now let us take a look at how the code for solving a Sudoku puzzle is implemented using Amplify.

## Formulation of constraints¶

### How to express Sudoku constraints on annealing machines¶

First, we consider how to express the constraint conditions that correspond to the Sudoku rules, using binary variables whose values are $\{0, 1\}$.

The rules to be expressed are the followings:

- Each row contains the digits from $1$ to $9$ without overlap.
- Each column contains the digits $1$ to $9$ without overlap.
- Each $3\times3$ subgrid contains the digits $1$ to $9$ without overlap.

Without imposing the above rules 1-3, $9$ possible numbers can be placed in each of $9\times9=81$ cells. We thus consider assigning $9$ binary variables to each of the cells, so that $9\times9\times9=729$ variables are prepared in total. One way to describe this situation is that we stack $9$ layers of $9\times9=81$ grids, so that each layer corresponds to the numbers from $1$ to $9$, respectively. For the $9\times9=81$ grid in each layer, we let $i$ and $j$ be the indices expressing the rows and columns, respectively, with $i,j=0,\cdots,8$, corresponding to $1, \cdots, 9$ th row and column. We also assign the layer index $k=0,\dots,8$, so that the layer specified by the index corresponds to numbers $1,\cdots,9$, respectively.

From this observation, we can express the binary variables by $q_{i,j,k}$, where $i,j,k=0,\cdots,8$ correspond to the indices for rows, columns, and layers, respectively. For example, $q_{2,4,6}=1$ means that number $7$ is placed in the cell at the 3rd row and 5th column, and $q_{2,4,6}=0$ means that $7$ is not to be placed here.

Using these variables, the constraints 1-3 can be written as the following one-hot constraints, respectively.

$$ \left\{ \begin{align} &\begin{split} (a) \quad &\sum_{j=0}^8 q_{i,j,k}=1 \end{split}\\ &\begin{split} (b) \quad &\sum_{i=0}^8 q_{i,j,k}=1 \end{split}\\ &\begin{split} (c) \quad &\sum_{i,j\in 3\times3\,\text{subgrid}}q_{i,j,k}=1 \end{split}\\ &\begin{split} (d) \quad &\sum_{k=0}^8 q_{i,j,k}=1 \end{split} \end{align} \right. $$

The constraints $(a)$, $(b)$ and $(c)$ correspond to rules 1, 2 and 3, respectively. The $(d)$ constraint corresponds to the basic condition that each cell can only contain one number.

## Initialization¶

In Sudoku, clues are given as an initial placement with 17 or more numbers already placed. Here we use the following initial placement which is considered a difficult problem. In the following notation, the cells with no filled numbers are marked with $0$.

```
import numpy as np
# Indicate initial placement as a list
# Reference: http://www.sudokugame.org/archive/printable.php?nd=4&y=2018&m=03&d=1
initial = np.array(
[
[2, 0, 5, 1, 3, 0, 0, 0, 4],
[0, 0, 0, 0, 4, 8, 0, 0, 0],
[0, 0, 0, 0, 0, 7, 0, 2, 0],
[0, 3, 8, 5, 0, 0, 0, 9, 2],
[0, 0, 0, 0, 9, 0, 7, 0, 0],
[0, 0, 0, 0, 0, 0, 4, 5, 0],
[8, 6, 0, 9, 7, 0, 0, 0, 0],
[9, 5, 0, 0, 0, 0, 0, 3, 1],
[0, 0, 4, 0, 0, 0, 0, 0, 0],
]
)
# Display Sudoku grid with numbers
def print_sudoku(sudoku):
print("\n\n")
for i in range(len(sudoku)):
line = ""
if i == 3 or i == 6:
print("---------------------")
for j in range(len(sudoku[i])):
if j == 3 or j == 6:
line += "| "
line += str(sudoku[i][j]) + " "
print(line)
print_sudoku(initial)
```

```
from amplify import VariableGenerator
gen = VariableGenerator()
q = gen.array("Binary", 9, 9, 9)
```

This gives a three-dimensional array with $9^3=729$ variables. Each of `9, 9, 9`

represents
the number of elements in the rows, columns, and number layers, and their indices are `i`

,
`j`

, and `k`

, respectively. Each element can be accessed with
`q[i][j][k]`

. For example, the 9 variables for row number $0$ and column number $0$ can be
displayed as follows:

```
# variables at column 0 and row 0
print(q[0, 0])
# 9 variables at column 2 and number layer 5
print(q[2, :, 5])
```

From the initial placement stored in `initial`

, we can narrow down the unknown variables
that
can be determined accordingly to the rules of Sudoku. For example, the cell indicated by
`i=1`

and `j=5`

already has a number $8$ (`initial[1][5]=8`

), so the variable in the
layer
`k=7`

corresponding to the number $8$ can be specified by `q[1][5][7]=1`

. In
this
specific case, we can place some numbers accordingly to the rules in the following way.

Rules 1 and 2 allow us to further determine the values of the variables `q[i][5][7]=0`

($i\neq1$), `q[1][j][7]=0`

($j\neq5$), and `q[i][5][7]=0`

($j\neq5$), since the
row
and column containing the cell indicated by `i=1`

and `j=5`

in the layer
`k=7`

cannot contain the same number as $8$. This corresponds to imposing the constraints
$(a)$
and $(b)$.

Also, by rule 3, the $3\times3$ subgrid, which this number $8$ above belongs to, cannot contain $8$
any
more, so in $(i,j)\in\{(0,3),(0,4),(0,5),(1,3),(1,4),(2,3),(2,4),(2,5)\}$, we can let
`q[i ][j][7]=0`

. We can imposed the constraint $(c)$ in this manner.

Furthermore, we impose the constraint $(d)$ to place only one of the numbers in a cell where the
number
is fixed. In the above example, we have `q[1][5][k]=0`

($k\neq7$). By performing similar
operations for all the cells given from the initial layout, we can narrow down the required variables
so
that we can perform the calculations with fewer number of variables.

```
for i, j in zip(
*np.where(initial != 0)
): # Obtain the indices of non-zero rows, columns, and layers
k = initial[i, j] - 1 # Note that -1 is used to convert from value to index
q[i, :, k] = 0 # Constraint (a)
q[:, j, k] = 0 # Constraint (b)
q[i, j, :] = 0 # Constraint (d)
for m in range(9):
q[(3 * (i // 3) + m // 3), (3 * (j // 3) + m % 3), k] = 0 # Constraint (c)
q[i, j, k] = 1 # # Set the value of the variable to 1
```

We have now prepared the initial settings. As an example, if you display the 9 variables with the row number $0$ and the column number $0$, you will see that the second element is fixed as $1$, that is, the number $2$ is assigned to this cell.

```
q[0][0]
```

Similarly, if we display the 9 variables with the row number $0$ and column number $1$, we can see that the numbers $1,2,3,4,5,6$ are not candidates.

```
q[0][1]
```

## Setting constraints¶

Next, we show how the constraint parts are implemented.
The one-hot constraint in $(a)$-$(d)$ can be expressed using Amplify's `one_hot`

function.

We start by defining the constraint condition corresponding to $(a)$ that each row of cannot contain
the
same number.
Since the sum of the variables for all columns specified by row `i`

and layer
`k`

is
given by $\sum_{j=0}^{8}q_{i,j,k}$, the constraint can be obtained by setting this sum to $1$ as
follows.

```
from amplify import one_hot
# (a): Constraint that each row cannot contain the same number
# Imposing the one_hot constraint for the direction j (i.e. axis=1).
row_constraints = one_hot(q, axis=1)
```

Similarly, the constraint $(b)$ that each column has cannot contain the same number and $(d)$ that only one number can be placed in a given cell can be expressed as follows:

```
# (b): Constraint that each column cannot contain the same number
# Imposing the one_hot constraint for the direction i (i.e. axis=0).
col_constraints = one_hot(q, axis=0)
# (d): Constraint that only one number can be placed in a cell
# Imposing the one_hot constraint for the direction k (i.e. axis=2).
num_constraints = one_hot(q, axis=2)
```

Finally, we express the constraint $(c)$ that each $3\times3$ subgrid cannot contain the same number. We take the sum of the variables in each $3\times3$ subgrid for all layers and impose the one-hot constraint as follows:

```
from amplify import sum
# (c): Constraint that 3x3 subgrids cannot contain the same number
block_constraints = [
one_hot(sum([q[i + m // 3, j + m % 3, k] for m in range(9)]))
for i in range(0, 9, 3)
for j in range(0, 9, 3)
for k in range(9)
]
```

We have now obtained all the constraints, and we can add up all these constraints together to obtain
the
final constraint object. Here, since only `block_constraints`

is still in a list, they are
added together using the `amplify.sum`

imported in the cell above, converted to the
`amplify.ConstraintList`

type, and then added together with the other constraints.

```
constraints = (
row_constraints + col_constraints + num_constraints + sum(block_constraints)
)
```

This completes the formulation of the Sudoku problem. If the Ising machine can find the combination of variables that satisfies all the constraints, this allows us to determine the numbers to be placed for all the cells, giving the solution for Sudoku.

## Execution of the Ising machine¶

```
from amplify import solve, FixstarsClient
from datetime import timedelta
client = FixstarsClient()
client.parameters.timeout = timedelta(milliseconds=1000) # timeout is 1000 ms
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" # If you use Amplify in a local environment, enter the Amplify API token.
# Solve the problem
result = solve(constraints, client)
if len(result) == 0:
raise RuntimeError("Some of the constraints are not satisfied.")
values = result.best.values
```

The result of the execution is stored in `values`

. By using the `evaluate`

member
function of `q`

, the solution is arranged in the original variable shape.
Then, by searching for `k`

such that `q[i][j][k] = 1`

for all `i`

and
`j`

, we can get the number `k + 1`

as the solution for the indicated cell.

```
q_values = q.evaluate(values)
answer = np.array([np.where(np.array(q_values[i]) != 0)[1] + 1 for i in range(9)])
```

Finally, the `print_sudoku`

function is used to output the solution.

```
print_sudoku(answer)
```

## Generalization of Sudoku¶

So far we have dealt with a $9 \times 9$ grid with $3 \times 3$ subgrids, but the Ising machine can easily handle Sudoku puzzles with extended problem sizes, such as $16 \times 16$ and $25 \times 25$. We let the number of cells in a Sudoku puzzle be $N\times N,\,(N\in\mathbb{Z})$ and the corresponding subgrids be in the form of $n\times n$ (where $N=n^2,\,(n\in\mathbb{Z})$). For example, a basic $9\times9$ Sudoku would have $N=9$ and $n=3$.

We generalize the previous code using $N$ and $n$ and solve the $16\times16$ mass Sudoku as an example.

```
from amplify import VariableGenerator, solve, one_hot, FixstarsClient
import numpy as np
n = 4 # block size
N = n * n # total number of blocks
# Initial configuration for n = 4 (N = 16)
# Reference: https://www.free-sudoku-puzzle.com/puzzle_fours/solve/3/238
initial = np.array(
[
[0, 0, 0, 7, 0, 0, 0, 1, 5, 0, 3, 16, 4, 0, 15, 0],
[0, 11, 0, 0, 0, 0, 5, 0, 0, 2, 12, 6, 0, 0, 7, 14],
[4, 0, 0, 0, 7, 8, 9, 0, 11, 0, 1, 15, 0, 0, 10, 0],
[10, 0, 0, 0, 0, 0, 0, 15, 13, 0, 9, 7, 8, 0, 0, 1],
[13, 0, 0, 16, 15, 0, 4, 9, 0, 0, 14, 0, 11, 0, 1, 0],
[8, 0, 5, 0, 0, 0, 10, 0, 0, 0, 0, 0, 15, 0, 0, 0],
[0, 0, 0, 11, 0, 0, 0, 8, 16, 7, 0, 9, 0, 0, 0, 0],
[0, 0, 0, 14, 0, 0, 3, 0, 4, 0, 0, 5, 13, 0, 0, 0],
[0, 0, 0, 0, 3, 0, 0, 14, 0, 0, 4, 0, 9, 12, 8, 15],
[0, 0, 0, 0, 0, 1, 7, 10, 0, 15, 8, 11, 0, 0, 0, 0],
[0, 0, 0, 0, 11, 12, 0, 0, 0, 0, 16, 0, 3, 5, 0, 0],
[0, 0, 0, 0, 0, 16, 0, 0, 0, 5, 0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 6, 16, 0, 15, 1, 5, 0, 14, 2, 0, 0],
[0, 0, 0, 3, 0, 0, 0, 0, 9, 0, 0, 14, 0, 1, 0, 4],
[2, 0, 12, 0, 0, 0, 0, 0, 0, 16, 13, 0, 6, 0, 3, 5],
[1, 0, 0, 0, 0, 15, 0, 0, 2, 11, 6, 12, 7, 9, 0, 10],
]
)
# Display a sudoku grid with numbers
def print_sudoku(sudoku):
width = len(str(N))
for i in range(len(sudoku)):
line = ""
if i % n == 0 and i != 0:
print("-" * ((width + 1) * n * n + 2 * (n - 1)))
for j in range(len(sudoku[i])):
if j % n == 0 and j != 0:
line += "| "
line += str(sudoku[i][j]).rjust(width) + " "
print(line)
q = VariableGenerator().array("Binary", N, N, N)
for i, j in zip(*np.where(initial != 0)):
k = initial[i, j] - 1
q[i, :, k] = 0 # Constraint (a)
q[:, j, k] = 0 # Constraint (b)
q[i, j, :] = 0 # Constraint (d)
for m in range(N):
q[(n * (i // n) + m // n), (n * (j // n) + m % n), k] = 0 # Constraint (c)
q[i, j, k] = 1 # Substitute 1
# (a): Constraints that each row cannot contain the same number
row_constraints = one_hot(q, axis=1)
# (b): Constraints that each column cannot contain the same number
col_constraints = one_hot(q, axis=0)
# (d): Constraints that only one number can be in a cell
num_constraints = one_hot(q, axis=2)
# (c): Constraints that a nxn subgrid cannot contain the same number
block_constraints = [
one_hot(sum([q[i + m // n, j + m % n, k] for m in range(N)]))
for i in range(0, N, n)
for j in range(0, N, n)
for k in range(N)
]
constraints = (
row_constraints + col_constraints + num_constraints + sum(block_constraints)
)
client = FixstarsClient()
client.parameters.timeout = timedelta(milliseconds=10000) # timeout is 10 seconds
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" # If you use Amplify in a local environment, enter the Amplify API token.
# Solve the problem
result = solve(constraints, client)
if len(result) == 0:
raise RuntimeError("Some constraints are unsatisfied.")
values = result.best.values
q_values = q.evaluate(values)
answer = np.array([np.where(np.array(q_values[i]) != 0)[1] + 1 for i in range(N)])
```

```
print_sudoku(initial)
```

```
print_sudoku(answer)
```