1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71

#!/usr/bin/env python

from itertools import product
from zibopt import scip

# 0 indicates a cell value is not given
problem = [
    [0, 0, 0,   6, 9, 2,   0, 4, 0],
    [7, 0, 0,   0, 0, 0,   8, 9, 0],
    [0, 0, 0,   0, 0, 0,   0, 0, 6],

    [0, 0, 9,   0, 1, 7,   0, 0, 3],
    [0, 0, 7,   0, 8, 0,   5, 0, 0],
    [8, 0, 0,   4, 6, 0,   1, 0, 0],

    [5, 0, 0,   0, 0, 0,   0, 0, 0],
    [0, 8, 6,   0, 0, 0,   0, 0, 1],
    [0, 3, 0,   7, 2, 8,   0, 0, 0]
]

if __name__ == '__main__':
    # Indexes for creating variables and constraints
    rows = range(len(problem))
    cols = rows
    vals = range(1, 10)
    groups = (0, 3, 6)

    # 9x9x9 structure for storing all binary variables
    # If x[i][j][k] == 1, then problem[i][j] == k.
    x = dict((i, dict((j, {}) for j in cols)) for i in rows)

    # Initialize one binary variable per cell value
    solver = scip.solver()
    for i, j, k in product(rows, cols, vals):
        x[i][j][k] = solver.variable(
            vartype = scip.BINARY,
            lower   = 1 if problem[i][j] == k else 0
        )
     
    # Each cell takes on exactly one value
    for i, j in product(rows, cols):
        solver += sum(x[i][j][k] for k in vals) == 1

    # Each value occurs in each row once
    for i, k in product(rows, vals):
        solver += sum(x[i][j][k] for j in cols) == 1

    # Each value occurs in each column once
    for j, k in product(cols, vals):
        solver += sum(x[i][j][k] for i in rows) == 1

    # Each 3x3 group has all unique values    
    for m, n, k in product(groups, groups, vals):
        solver += sum(x[i][j][k] for i, j in product(range(m,m+3), range(n,n+3))) == 1

    # All variables have 0 coefficients, so we can either max or min            
    solution = solver.maximize()
    if solution:
        # Convert solution values to a nice solution matrix
        for i, j in product(rows, cols):
            for k in vals:
                if solution[x[i][j][k]]:
                    problem[i][j] = k
                    break
                
        for row in problem:
            print(row)

    else:
        print('invalid')