Second algorithm for approximating Pareto set.

second_report.md

Պարետոյի բազմությունը մոտարկելու երկրորդ ալգորիթմ։

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
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 visualize(polygon, f1, f2, pareto, plot_bad, name):
    n = polygon.shape[0]

    fig, ax = plt.subplots(ncols=2, figsize=(12, 6))

    ax[0].set_xlabel("$x$")
    ax[0].set_ylabel("$y$")

    ax[1].set_xlabel("$f_1$")
    ax[1].set_ylabel("$f_2$")

    colors = matplotlib.cm.rainbow(np.linspace(0, 1, len(pareto)))

    dict_x = {}
    dict_f = {}

    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')

    border = []
    for i, A in enumerate(pareto):
        x, y = A
        ax[0].scatter([x], [y], c=colors[i])
        k = "{:.5f} {:.5f}".format(x, y)
        if k not in dict_x:
            dict_x[k] = []
        dict_x[k].append(i+1)

        f1_val = x * f1[0] + y * f1[1]
        f2_val = x * f2[0] + y * f2[1]
        ax[1].scatter([f1_val], [f2_val], c=colors[i])
        k = "{:.5f} {:.5f}".format(f1_val, f2_val)
        if k not in dict_f:
            dict_f[k] = []
        dict_f[k].append(i+1)

        border.append((pareto[i], f1_val))

    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[0].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[1].plot([p1_x, p2_x], [p1_y, p2_y], 'r')

    for k, v in dict_x.iteritems():
        ax[0].text(float(k.split(' ')[0]), float(k.split(' ')[1]), "  {}".format(v))

    for k, v in dict_f.iteritems():
        ax[1].text(float(k.split(' ')[0]), float(k.split(' ')[1]), "  {}".format(v))

    fig.savefig(name + ".png")


def main():
    (polygon, f1, f2) = read_problem()
    A, b = get_constraints(polygon)

    pareto = []
    visualize(polygon, f1, f2, pareto, plot_bad=True, name='iter0')

    alphas = np.linspace(0, 1, 10)
    for iteration, alpha in enumerate(alphas):
        print "Solving for coefficients [{}, {}]".format(alpha, 1-alpha)

        # Make linear programming problem and solve
        sol = linprog(c=-alpha*np.array(f1) - (1-alpha)*np.array(f2),
                       A_ub=A,
                       b_ub=b,
                       bounds=(None, None))

        sol_x = getattr(sol, "x")
        print "\tsolution: {}".format(sol_x)
        pareto.append(sol_x)

        visualize(polygon, f1, f2, pareto, plot_bad=True, name="iter{}".format(iteration + 1))

    visualize(polygon, f1, f2, 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  

Աշխատանքի ընթացքում ծրագիրը կտպի՝

Solving for coefficients [0.0, 1.0]
        solution: [-1. -4.]
Solving for coefficients [0.111111111111, 0.888888888889]
        solution: [-1. -4.]
Solving for coefficients [0.222222222222, 0.777777777778]
        solution: [-1. -4.]
Solving for coefficients [0.333333333333, 0.666666666667]
        solution: [-1. -4.]
Solving for coefficients [0.444444444444, 0.555555555556]
        solution: [-2. -2.]
Solving for coefficients [0.555555555556, 0.444444444444]
        solution: [-2. -2.]
Solving for coefficients [0.666666666667, 0.333333333333]
        solution: [-2. -2.]
Solving for coefficients [0.777777777778, 0.222222222222]
        solution: [-2.  3.]
Solving for coefficients [0.888888888889, 0.111111111111]
        solution: [-2.  3.]
Solving for coefficients [1.0, 0.0]
        solution: [ 0.  5.]

Վերջնական պատասխանը կընդունի հետևյալ տեսքը՝