Pyomo/snake_egg_modular.py at main · OptimizationExpert/Pyomo · GitHub

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
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
"""
Developed by Alireza Soroudi
alireza.soroudi@gmail.com
The Snake: Draws a continuous path moving up, down, left, or right between cells.
The head and tail start at marked circular cells.
The snake cannot touch itself straight-on (orthogonally), but can touch diagonally,
 and cannot revisit squares.
The Eggs: The empty cells left over outside the snake form distinct, connected groups (eggs).
Sizes: These groups match a precise set of sizes sequentially
(such as one of each size from 1 to N-1 in standard large grids).
"""

import random

import matplotlib.pyplot as plt
from ortools.sat.python import cp_model

st, fn, N = 98, 58, 10
st, fn, N = 98, 58, 10
st, fn, N = 16, 19, 6
st, fn, N = 25, 13, 6


def create_grid(N):
    grid = {}
    c = 0
    for i in range(N):
        for j in range(N):
            c += 1
            grid[c] = (i, j)
    return grid


grid = create_grid(N)
nodes = [c for c in grid.keys()]
connect = {(c1, c2) for c1 in nodes for c2 in nodes if
           abs(grid[c1][0] - grid[c2][0]) + abs(grid[c1][1] - grid[c2][1]) == 1}


def opt_model(st, fn, N):
    model = cp_model.CpModel()
    solver = cp_model.CpSolver()
    numbers = [n for n in range(N)]
    non_snake_numbers = [n for n in range(1, N)]
    cell_assign = {(i, n): model.new_bool_var(f"assign_to_number_{i}_{n}") for i in nodes for n in numbers}
    U = {(i, j): model.new_bool_var(f"connection_{i}_{j}") for i in nodes for j in nodes if (i, j) in connect}
    snake = {i: cell_assign[i, 0] for i in nodes}

    flow_n = {(i, j, n): model.new_int_var(0, n - 1, f"flow_{i}_{j}_{n}") for i in nodes for j in nodes for n in
              non_snake_numbers if
              (i, j) in connect}
    source_n = {(i, n): model.new_bool_var(f"source_{i}_{n}") for i in nodes for n in non_snake_numbers}

    arcs = [(i, j, v) for (i, j), v in U.items()] + [(i, i, v.Not()) for i, v in snake.items()]
    arcs += [(st, fn, True)]
    model.add_circuit(arcs)

    for (i, n), v in cell_assign.items():
        if n in non_snake_numbers:
            model.add(source_n[i, n] <= v)
    for i in nodes:
        expr = [cell_assign[i, n] for n in numbers]
        model.add_exactly_one(expr)

    for n in non_snake_numbers:
        expr_assign = [cell_assign[i, n] for i in nodes]
        model.add(sum(expr_assign) == n)
        expr_src = [source_n[i, n] for i in nodes]
        model.add_exactly_one(expr_src)
    for i in nodes:
        for n in non_snake_numbers:
            out_flow = sum([flow_n[i, j, n] for j in nodes if (i, j, n) in flow_n])
            in_flow = sum([flow_n[j, i, n] for j in nodes if (j, i, n) in flow_n])
            model.add(n * source_n[i, n] - cell_assign[i, n] == out_flow - in_flow).only_enforce_if(cell_assign[i, n])

    for (i, j, n), v in flow_n.items():
        model.add(v <= n * cell_assign[i, n])
        model.add(v <= n * cell_assign[j, n])

    model.add(cell_assign[fn, 0] == True)
    model.add(cell_assign[st, 0] == True)
    for n in non_snake_numbers:
        model.add(cell_assign[fn, n] == False)
        model.add(cell_assign[st, n] == False)

    # Maximize x
    expressions = [v for (i, j), v in U.items()]
    model.minimize(sum(expressions))

    status = solver.Solve(model)
    print(solver.status_name(status))
    return U, cell_assign, source_n, numbers, solver, status


def visualize(U, cell_assign, source_n, numbers, solver):
    KOLOR = ['grey', 'red', 'blue', 'green', 'gold', 'brown',
             'dodgerblue', 'cyan', 'magenta', 'orange', 'purple', 'black']
    fig, ax = plt.subplots(figsize=(5, 5))
    for n in numbers:
        for i in nodes:
            if solver.value(cell_assign[i, n]) > 0:
                x0, y0 = grid[i]
                square = plt.Rectangle((x0 + .5, y0 + 0.5), 1, 1, fill=True, linewidth=2, facecolor=KOLOR[n],
                                       edgecolor='purple', alpha=0.3)
                ax.add_patch(square)

    for (i, j), v in U.items():
        if solver.value(v) > 0:
            (x0, y0) = grid[i]
            (x1, y1) = grid[j]
            plt.plot([x0 + 1, x1 + 1], [y0 + 1, y1 + 1], c='k', lw=5)

    for i in [fn, st]:
        x0, y0 = grid[i]
        plt.scatter(x0 + 1, y0 + 1, s=200, c='r', alpha=1, zorder=2)
    ax.set_xlim(0.5, N + 0.5)
    ax.set_ylim(0.5, N + 0.5)
    ax.set_aspect('equal')
    ax.grid(True)
    plt.xticks(range(1, 1 + N))
    plt.yticks(range(1, 1 + N))
    plt.title("Square")
    plt.tight_layout()
    plt.savefig(f"square{fn}-{st}.png")
    plt.show()


seen = []
for r in range(100):
    (st, fn) = random.sample(range(1, 1 + N ** 2), 2)
    if (st, fn) not in seen and (fn, st) not in seen:
        print(st, fn)
        seen.append((st, fn))
        U, cell_assign, source_n, numbers, solver, status = opt_model(st, fn, N)
        if solver.status_name(status) == 'OPTIMAL':
            visualize(U, cell_assign, source_n, numbers, solver)