Pyomo/maxpasssum_visual_multi.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
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
"""
Developed by Alireza Soroudi
Constraint programming in ORTools
alireza.soroudi@gmail.com
puzzle source: project 18,67
https://projecteuler.net/problem=67

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

data = [
    [75],
    [95, 64],
    [17, 47, 82],
    [18, 35, 87, 10],
    [20, 4, 82, 47, 65],
    [19, 1, 23, 75, 3, 34],
    [88, 2, 77, 73, 7, 63, 67],
    [99, 65, 4, 28, 6, 16, 70, 92],
    [41, 41, 26, 56, 83, 40, 80, 70, 33],
    [41, 48, 72, 33, 47, 32, 37, 16, 94, 29],
    [53, 71, 44, 65, 25, 43, 91, 52, 97, 51, 14],
    [70, 11, 33, 28, 77, 73, 17, 78, 39, 68, 17, 57],
    [91, 71, 52, 38, 17, 14, 91, 43, 58, 50, 27, 29, 48],
    [63, 66, 4, 68, 89, 53, 67, 30, 73, 16, 69, 87, 40, 31],
    [4, 62, 98, 27, 23, 9, 70, 98, 73, 93, 38, 53, 60, 4, 23]
]

rows = [r for r in range(len(data))]
cols = [len(data[r]) for r in rows]
col_max = max(cols)
n = 0
dic_data = {}
for r in rows:
    for c, v in enumerate(data[r]):
        n += 1
        dic_data[n] = (r, c, v)
dic_by_loc = {(r, c): n for n, (r, c, v) in dic_data.items()}

nodes = [i for i in range(1, n + 1)]
neighbour = {}
for n in nodes:
    (r, c, v) = dic_data[n]
    neighbour[n] = []
    if r + 1 in rows:
        cols = [j for j in range(len(data[r + 1]))]
        if c + 1 in cols:
            neighbour[n].append(dic_by_loc[r + 1, c + 1])
        if c in cols:
            neighbour[n].append(dic_by_loc[r + 1, c])

print(neighbour)
rowmax = max(rows)


def opt_path(direction: int, minlim: int, maxlim: int):
    model = cp_model.CpModel()
    x = {(i, j): model.new_bool_var(f"flow_{i}_{j}") for i in nodes
         for j in nodes if i in neighbour and j in neighbour[i]}

    fn = {i: model.new_bool_var(f"fn_{i}") for i in nodes if dic_data[i][0] == rowmax}
    select = {i: model.new_bool_var(f"select_{i}") for i in nodes}

    arcs = [(i, j, v) for (i, j), v in x.items()] + [(j, 1, v_link) for j, v_link in fn.items()] + [(i, i, v.Not()) for
                                                                                                    i, v
                                                                                                    in select.items()]
    model.add_circuit(arcs)

    for (i, j), v in x.items():
        model.add(v <= select[i])
        model.add(v <= select[j])

    expr_fn = [v_link for j, v_link in fn.items()]
    model.add_exactly_one(expr_fn)

    expr_fn_select = [select[j] for j, v_link in fn.items()]
    model.add_exactly_one(expr_fn_select)
    for j, v_link in fn.items():
        model.add(v_link <= select[j])

    expr = [v * dic_data[i][2] for i, v in select.items()]
    if direction > 0:
        model.maximize(sum(expr))
    elif direction < 0:
        model.minimize(sum(expr))
    else:
        model.add(cp_model.LinearExpr.Sum(expr) <= int(0.8 * maxlim))
        model.add(cp_model.LinearExpr.Sum(expr) >= int(1.5 * minlim))

    solver = cp_model.CpSolver()
    results = solver.Solve(model)
    print(solver.status_name(results))
    print('OF = ', solver.best_objective_bound)
    out = {(i, j) for (i, j), v in x.items() if solver.value(v) > 0}
    return out, int(solver.best_objective_bound)


xp, ofmax = opt_path(1, 0, 20000)
xn, ofmin = opt_path(-1, 0, ofmax)
xr, of = opt_path(0, ofmin, ofmax)

plt.figure(figsize=(5, 5))

for n, (r, c, v) in dic_data.items():
    bias = 0.5 * col_max + 0.5 - 0.5 * r
    plt.scatter(bias + c, rowmax - r, s=220, c='grey', alpha=0.4, zorder=0)
    # plt.text(bias + c - 0.2, rowmax - r - 0.7, s=str(n))
    plt.text(bias + c - 0.2, rowmax - r - 0.1, s=str(v), fontweight='bold')

plt.figure(figsize=(5, 5))

for n, (r, c, v) in dic_data.items():
    bias = 0.5 * col_max + 0.5 - 0.5 * r
    plt.scatter(bias + c, rowmax - r, s=220, c='grey', alpha=0.4, zorder=0)
    plt.text(bias + c - 0.2, rowmax - r - 0.1, s=str(v), fontweight='bold')
plt.axis('off')
plt.tight_layout()

for r in rows:
    for (i, j) in xp:
        r0, c0, v0 = dic_data[i]
        r1, c1, v1 = dic_data[j]
        bias0 = 0.5 * col_max + 0.5 - 0.5 * r0
        bias1 = 0.5 * col_max + 0.5 - 0.5 * r1

        x0, y0 = bias0 + c0, rowmax - r0
        x1, y1 = bias1 + c1, rowmax - r1
        if r0 == r:
            plt.plot([x0, x1], [y0, y1], lw=3, color='r', zorder=1)
        if r == r0 == rowmax - 1:
            plt.plot([x0, x1], [y0, y1], lw=3, color='r', zorder=1, label='max')
    for (i, j) in xn:
        r0, c0, v0 = dic_data[i]
        r1, c1, v1 = dic_data[j]
        bias0 = 0.5 * col_max + 0.5 - 0.5 * r0
        bias1 = 0.5 * col_max + 0.5 - 0.5 * r1

        x0, y0 = bias0 + c0, rowmax - r0
        x1, y1 = bias1 + c1, rowmax - r1
        if r0 == r:
            plt.plot([x0, x1], [y0, y1], lw=3, color='g', zorder=1)
        if r == r0 == rowmax - 1:
            plt.plot([x0, x1], [y0, y1], lw=3, color='g', zorder=1, label='min')
    for (i, j) in xr:
        r0, c0, v0 = dic_data[i]
        r1, c1, v1 = dic_data[j]
        bias0 = 0.5 * col_max + 0.5 - 0.5 * r0
        bias1 = 0.5 * col_max + 0.5 - 0.5 * r1

        x0, y0 = bias0 + c0, rowmax - r0
        x1, y1 = bias1 + c1, rowmax - r1
        if r0 == r:
            plt.plot([x0, x1], [y0, y1], lw=3, color='gold', zorder=1)
        if r == r0 == rowmax - 1:
            print(i)
            plt.plot([x0, x1], [y0, y1], lw=3, color='gold', zorder=1, label='random')

    plt.legend()
    plt.savefig(f"row {r}.png")
plt.show()