Gradient

import matplotlib.pyplot as plt
import numpy as np
import math as math


class Vec2:
    x = 0
    y = 0

    def __init__(self, x, y):
        self.x = x
        self.y = y

    def __add__(self, other):
        return Vec2(self.x + other.x, self.y + other.y)

    def __sub__(self, other):
        return Vec2(self.x - other.x, self.y - other.y)

    def __mul__(self, other):
        return Vec2(self.x * other, self.y * other)

    def __truediv__(self, other):
        return Vec2(self.x / other, self.y / other)

    def __iadd__(self, other):
        self.x += other.x
        self.y += other.y
        return self

    def __isub__(self, other):
        self.x -= other.x
        self.y -= other.y
        return self

    def __imul__(self, other):
        self.x *= other
        self.y *= other
        return self

    def __repr__(self):
        return "{} {}".format(self.x, self.y)

    def SquaredLength(self):
        return self.x * self.x + self.y * self.y

    def Length(self):
        return math.sqrt(self.SquaredLength())

    def __getitem__(self, key):
        if key == 0:
            return self.x;
        return self.y;


# Отобразить линии уровня функции двух переменных и траекторию работы алгоритма
def plot_3d_with_trajectory(
        f, trajectory, x_left, x_right, y_left, y_right, resol=10, title="F(X, Y)"
):
    # f - отрисовываемая функция
    # trajectory - последовательность точек траектории
    # x_left - левая граница оси X
    # x_right - правая граница оси X
    # y_left - левая граница оси Y
    # y_right - правая граница оси Y
    # resol - разрешение, число отрисовываемых точек на единицу длины
    # title - заголовок графика

    fig = plt.figure()
    ax = fig.add_subplot(111, projection="3d")

    ax.set_xlabel("X")
    ax.set_ylabel("Y")
    ax.set_title(title)

    x = np.linspace(x_left, x_right, int((x_right - x_left) * resol))
    y = np.linspace(y_left, y_right, int((y_right - y_left) * resol))
    x, y = np.meshgrid(x, y)
    z = f(x, y)

    ax.contour3D(x, y, z)

    points_x = [point.x for point in trajectory]
    points_y = [point.y for point in trajectory]
    points_z = [f(point.x, point.y) for point in trajectory]

    ax.scatter(points_x, points_y, points_z, c="red")

    plt.show()


def derivative(func, point, h=0.001):
    res = func(point.x, point.y)
    return Vec2(func(point.x + h, point.y) - res, func(point.x, point.y + h) - res) / h;


def gradient(func, startPoint, method, eps=0.00001):
    newPoint = startPoint
    result = []
    while True:
        lastPoint = newPoint
        newPoint = lastPoint - derivative(func, lastPoint) * method(lastPoint, func)
        result.append(lastPoint)
        if np.absolute(func(newPoint.x, newPoint.y) - func(lastPoint.x, lastPoint.y)) < eps:
            return result
    return result


class constantMethod():
    value = 0.01

    def __init__(self, value=0.01):
        self.value = value

    def __call__(self, point, func):
        return self.value


class splitMethod():
    value = 0.01
    eps = 1
    delta = 0.5

    def __init__(self, value=0.01, eps=1, delta=0.5):
        self.value = value
        self.eps = eps
        self.delta = delta

    def __call__(self, point, func):
        div = derivative(func, point)
        fval = func(point.x, point.y)
        while func(point.x - self.value * div.x, point.y - self.value * div.y) > (
                fval - self.eps * self.value * div.SquaredLength()):
            self.value *= self.delta
        return self.value


class goldMethod():
    a = 0.0001
    b = 1
    eps = 0.00001

    def __init__(self, a=0.0001, b=1, eps=0.00001):
        self.a = a
        self.b = b
        self.eps = eps
        if a > b:
            self.a, self.b = self.b, self.a

    def __call__(self, point, func):
        fi = 1.6180339887
        grad = derivative(func, point)
        x1 = self.b - (self.b - self.a) / fi
        x2 = self.a + (self.b - self.a) / fi
        resa = self.a
        resb = self.b
        y1 = func(point.x - grad.x * x1, point.y - grad.y * x1)
        y2 = func(point.x - grad.x * x2, point.y - grad.y * x2)
        while np.absolute(resa - resb) > self.eps:
            if y1 <= y2:
                resb = x2
                x2 = x1
                x1 = resb - (resb - resa) / fi
                y2 = y1
                y1 = func(point.x - grad.x * x2, point.y - grad.y * x2)
            else:
                resa = x1
                x1 = x2
                x2 = resa + (resb - resa) / fi
                y1 = y2
                y1 = func(point.x - grad.x * x1, point.y - grad.y * x1)
        return (resb + resa) * 0.5


class fibonacciMethod():
    a = 0.0001
    b = 1
    n = 50

    def __init__(self, a=0.0001, b=1, n=50):
        self.a = a
        self.b = b
        self.n = n
        if a > b:
            self.a, self.b = self.b, self.a

    def fib(self, n):
        if n == 0 or n == 1:
            return 1
        a = 1
        b = 1
        res = 0
        for i in range(1, n):
            res = a + b
            a = b
            b = res
        return res

    def __call__(self, point, func):
        grad = derivative(func, point)
        resa = self.a
        resb = self.b
        for k in range(self.n):
            fibo = self.fib(self.n - k + 1)
            lambd = resa + self.fib(self.n - k - 1) / fibo * (resb - resa)
            mu = resa + self.fib(self.n - k) / fibo * (resb - resa)
            if func(point.x - grad.x * lambd, point.y - grad.y * lambd) > func(point.x - grad.x * mu,
                                                                               point.y - grad.y * mu):
                resa = lambd
            else:
                resb = mu
        return (resb + resa) * 0.5


def dixotomMin(func, point, grad, eps=0.001, delta=0.001, a=-100, b=100):
    t = Vec2(0, 0)
    x1 = (a + b - delta) * 0.5
    x2 = (a + b + delta) * 0.5
    total = 0
    while b - a > eps * 2:
        total += 1
        if total > 1000:
            break
        t = point + grad * x1
        y1 = func(t.x, t.y)
        t = point + grad * x2
        y2 = func(t.x, t.y)
        if y1 > y2:
            a = x1
        else:
            b = x2
        x1 = (a + b - delta) * 0.5
        x2 = (a + b + delta) * 0.5
    return (a + b) * 0.5


def gradientFletcherReeves(func, startPoint, eps=0.00001):
    newPoint = startPoint
    result = []
    beta = 0
    k = 0
    m = 100
    bFin = False
    while k <= m:
        lastPoint = newPoint
        grad = derivative(func, lastPoint)
        result.append(lastPoint)
        gradLen2 = grad.SquaredLength()
        if gradLen2 < eps * eps:
            return result
        if (k == 0):
            grad *= -1
        else:
            beta = gradLen2 / prevGradLen2
            grad *= beta - 1
        prevGradLen2 = gradLen2
        minn = dixotomMin(func, lastPoint, grad)
        newPoint = lastPoint + grad * minn
        norm = (newPoint - lastPoint).Length()
        if np.absolute(func(newPoint.x, newPoint.y) - func(lastPoint.x, lastPoint.y)) < eps and norm < eps:
            if bFin == False:
                bFin = True
            else:
                return result
        k += 1
    return result


# def function(x, y):
# return 10 * x * x + y * y
def function(x, y):
    return x * x + y * y + 2 * x * y


trajectory = gradient(function, Vec2(1, 5), constantMethod(0.001))
print(trajectory[len(trajectory) - 1])
print(len(trajectory))
plot_3d_with_trajectory(function, trajectory, -20, 20, -20, 20, 10, "Constant (X, Y)")

trajectory = gradient(function, Vec2(1, 5), splitMethod(0.1, 0.1, 0.1))
print(trajectory[len(trajectory) - 1])
print(len(trajectory))
plot_3d_with_trajectory(function, trajectory, -20, 20, -20, 20, 10, "Split (X, Y)")

trajectory = gradient(function, Vec2(1, 5), goldMethod())
print(trajectory[len(trajectory) - 1])
print(len(trajectory))
plot_3d_with_trajectory(function, trajectory, -20, 20, -20, 20, 10, "Gold (X, Y)")

trajectory = gradient(function,Vec2(1, 5), fibonacciMethod())
print(trajectory[len(trajectory) - 1])
print(len(trajectory))
plot_3d_with_trajectory(function, trajectory, -20, 20, -20, 20, 10, "Fibonacci (X, Y)")

trajectory = gradientFletcherReeves(function, Vec2(1, 5))
print(trajectory[len(trajectory) - 1])
print(len(trajectory))
plot_3d_with_trajectory(function, trajectory, -20, 20, -20, 20, 10, "FletcherReeves (X, Y)")