Պարետոյի բազմությունը մոտարկելու առաջին ալգորիթմ։
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import shapely.geometry as geom
from scipy.optimize import linprog
def read_problem():
fin = open("input.txt", "r")
# import sys
# fin = sys.stdin
print "Enter number of vertices: ",
n = int(fin.readline())
print "Enter the coordinates of vertices in clockwise or counter-clockwise order"
polygon = []
for i in range(n):
x, y = map(float, fin.readline().strip().split(' '))
polygon.append([x, y])
polygon = np.array(polygon)
# There are two objectives to maximize: f_1(x) and f2(x)
print "Enter coefficients of f_1(x): ",
mas = map(float, fin.readline().strip().split(' '))
assert len(mas) == 2
f1 = np.array(mas)
print "Enter coefficients of f_2(x): ",
mas = map(float, fin.readline().strip().split(' '))
assert len(mas) == 2
f2 = np.array(mas)
print "Problem is successfully inputed !!!"
return polygon, f1, f2
def get_constraints(polygon):
""" Takes polygon and returns the hyperplanes: Ax <= b
"""
mx = np.mean(polygon[:, 0])
my = np.mean(polygon[:, 1])
n = len(polygon)
A = []
bs = []
for i in range(n):
j = (i + 1) % n
ax, ay = polygon[i]
bx, by = polygon[j]
a = by - ay
b = ax - bx
c = -a * ax - b * ay
# ax + by + c = 0
val = np.sign(a * mx + b * my + c)
c = -c
assert val != 0
if val == 1:
a, b, c = -a, -b, -c
A.append([a, b])
bs.append(c)
return np.array(A), np.array(bs)
def is_inside(x, y, polygon):
polygon = geom.polygon.Polygon(list(polygon))
return polygon.contains(geom.Point(x, y))
def sample_point(polygon):
lx = np.min(polygon[:, 0])
rx = np.max(polygon[:, 0])
ly = np.min(polygon[:, 1])
ry = np.max(polygon[:, 1])
while True:
x = np.random.uniform(lx, rx)
y = np.random.uniform(ly, ry)
if is_inside(x, y, polygon):
return (x, y)
def visualize(polygon, f1, f2, points, pareto, plot_bad, name):
n = polygon.shape[0]
fig, ax = plt.subplots(ncols=3, figsize=(18, 6))
ax[0].set_xlabel("$x$")
ax[0].set_ylabel("$y$")
ax[1].set_xlabel("$x$")
ax[1].set_ylabel("$y$")
ax[2].set_xlabel("$f_1$")
ax[2].set_ylabel("$f_2$")
colors = matplotlib.cm.rainbow(np.linspace(0, 1, len(points)))
for i in range(n):
j = (i + 1) % n
xs = [polygon[i][0], polygon[j][0]]
ys = [polygon[i][1], polygon[j][1]]
ax[0].plot(xs, ys, 'b')
ax[1].plot(xs, ys, '--b')
for i, A in enumerate(points):
x, y = A
ax[0].scatter([x], [y], c=colors[i])
ax[0].text(x, y, " {}".format(i + 1))
f1_val = x * f1[0] + y * f1[1]
f2_val = x * f2[0] + y * f2[1]
ax[2].scatter([f1_val], [f2_val], c=colors[i])
ax[2].text(f1_val, f2_val, " {}".format(i + 1))
px, py = pareto[2*i]
p_f1 = px * f1[0] + py * f1[1]
p_f2 = px * f2[0] + py * f2[1]
ax[2].plot([f1_val, p_f1], [f2_val, p_f2], '--', color=colors[i])
px, py = pareto[2*i + 1]
p_f1 = px * f1[0] + py * f1[1]
p_f2 = px * f2[0] + py * f2[1]
ax[2].plot([f1_val, p_f1], [f2_val, p_f2], '--', color=colors[i])
border = []
for i, A in enumerate(pareto):
x, y = A
f1_val = x * f1[0] + y * f1[1]
f2_val = x * f2[0] + y * f2[1]
bad = False
eps = 1e-7
for j, B in enumerate(pareto):
if i == j:
continue
f1_val_b = B[0] * f1[0] + B[1] * f1[1]
f2_val_b = B[0] * f2[0] + B[1] * f2[1]
if f1_val_b > f1_val - eps and f2_val_b > f2_val - eps:
bad = True
if not bad:
border.append((pareto[i], f1_val))
if plot_bad or not bad:
ax[0].scatter([x], [y], c=colors[i//2])
ax[0].text(x, y, " ${}_{}$".format(i // 2 + 1, i % 2 + 1))
ax[2].scatter([f1_val], [f2_val], c=colors[i//2])
ax[2].text(f1_val, f2_val, " ${}_{}$".format(i // 2 + 1, i % 2 + 1))
border = sorted(border, key=lambda item: item[1])
for i in range(len(border)-1):
x1, y1 = border[i][0]
x2, y2 = border[i+1][0]
ax[1].plot([x1, x2], [y1, y2], 'r')
p1_x = f1[0] * x1 + f1[1] * y1
p1_y= f2[0] * x1 + f2[1] * y1
p2_x = f1[0] * x2 + f1[1] * y2
p2_y = f2[0] * x2 + f2[1] * y2
ax[2].plot([p1_x, p2_x], [p1_y, p2_y], 'r')
fig.savefig(name + ".png")
def main():
(polygon, f1, f2) = read_problem()
A, b = get_constraints(polygon)
points = []
pareto = []
visualize(polygon, f1, f2, points, pareto, plot_bad=True, name='iter0')
n_iterations = 10
for iteration in range(n_iterations):
print "Starting iteration {} ... ".format(iteration + 1)
x, y = sample_point(polygon)
points.append([x, y])
print "Selected point: {} {}".format(x, y)
# Make two linear programming problems and solve
f1_val = x * f1[0] + y * f1[1]
f2_val = x * f2[0] + y * f2[1]
# hold f1, maximize f2
sol1 = linprog(c=-np.array(f2),
A_ub=A,
b_ub=b,
A_eq=np.array(f1).reshape((1,2)),
b_eq=np.array([f1_val]),
bounds=(None, None))
pareto.append(getattr(sol1, "x"))
# hold f2, maximize f1
sol2 = linprog(c=-np.array(f1),
A_ub=A,
b_ub=b,
A_eq=np.array(f2).reshape((1, 2)),
b_eq=np.array([f2_val]),
bounds=(None, None))
pareto.append(getattr(sol2, "x"))
visualize(polygon, f1, f2, points, pareto, plot_bad=True, name="iter{}".format(iteration + 1))
visualize(polygon, f1, f2, points, pareto, plot_bad=False, name="final")
if __name__ == "__main__":
main()Աշխատեցնենք հետևյալ օրինակի վրա՝
Enter number of vertices: 9
Enter the coordinates of vertices in clockwise or counter-clockwise order
-1 -4
-2 -2
-2 3
-1 4
0 5
4 3
6 0
4 -2
2 -3
Enter coefficients of f_1(x): -1 1
Enter coefficients of f_2(x): -2 -2
Աշխատանքի ընթացքում ծրագիրը կտպի՝
Starting iteration 1 ...
Selected point: -0.6729050969 3.11166326347
Starting iteration 2 ...
Selected point: 0.982273905471 -1.46971433859
Starting iteration 3 ...
Selected point: 0.546581372964 1.62623339798
Starting iteration 4 ...
Selected point: -0.0519679063657 0.807525924481
Starting iteration 5 ...
Selected point: 2.42075415604 0.0976427951623
Starting iteration 6 ...
Selected point: 2.08518380468 0.660652821881
Starting iteration 7 ...
Selected point: 0.962018890124 0.508583870988
Starting iteration 8 ...
Selected point: -1.12124948861 -1.51789923326
Starting iteration 9 ...
Selected point: 2.11942995637 -0.120783055278
Starting iteration 10 ...
Selected point: 4.9275614714 -0.688260368964
Վերջնական պատասխանը կընդունի հետևյալ տեսքը՝
