mopt

import sys
import math
import numpy

class function:
    a = [1.11, 0.6404, 0.0844, 1.064, -0.5425, 0.826, -0.7579, -0.3042, -0.4513]
    x0 = [1, 1, 1]
    epsilon = 0.000001

    def calculate(self, x1, x2, x3):
        return self.a[0] * x1 * x1 + 2 * self.a[1] * x1 * x2 + 2 * self.a[2] * x1 * x3 + self.a[3] * x2 * x2 + 2 * self.a[4] * x2 * x3 + self.a[5] * x3 * x3 + 2 * self.a[6] * x1 + 2 * self.a[7] * x2 + 2 * self.a[8] * x3

    def dFdx1(self, x):
        return 2 * self.a[0] * x[0] + 2 * self.a[1] * x[1] + 2 * self.a[2] * x[2] + 2 * self.a[6]
    def dFdx1dx1(self):
        return 2 * self.a[0]
    def dFdx1dx2(self):
        return 2 * self.a[1]
    def dFdx1dx3(self):
        return 2 * self.a[2]
    def dFdx2(self, x):
        return 2 * self.a[1] * x[0] + 2 * self.a[3] * x[1] + 2 * self.a[4] * x[2] + 2 * self.a[7]
    def dFdx2dx2(self):
        return 2 * self.a[3]
    def dFdx2dx3(self):
        return 2 * self.a[4]
    def dFdx3(self, x):
        return 2 * self.a[2] * x[0] + 2 * self.a[4] * x[1] + 2 * self.a[5] * x[2] + 2 * self.a[8]
    def dFdx3dx3(self):
        return 2 * self.a[5]

    def Jacobi(self, x):
        return [-self.dFdx1(x), -self.dFdx2(x), -self.dFdx3(x)]
    def CountAlpha(self, x, h):
        f2 = [[self.dFdx1dx1(), self.dFdx1dx2(), self.dFdx1dx3()],
              [self.dFdx1dx2(), self.dFdx2dx2(), self.dFdx2dx3()],
              [self.dFdx1dx3(), self.dFdx2dx3(), self.dFdx3dx3()]]
        h_copy = h[:]
        res = [0, 0, 0]
        for i in range(0, 3):
            for j in range(0, 3):
                res[i] += h_copy[j] * f2[i][j]
        return -(self.dFdx1(x) * h[0] + self.dFdx2(x) * h[1] + self. dFdx3(x) * h[2]) / (h[0] * res[0] + h[1] * res[1] + h[2] * res[2])

    def MethodNewtonSubFunction(self, h):
        res = [0, 0, 0]
        iJacobi = [[self.dFdx1dx1(), self.dFdx1dx2(), self.dFdx1dx3()],
                   [self.dFdx1dx2(), self.dFdx2dx2(), self.dFdx2dx3()],
                   [self.dFdx1dx3(), self.dFdx2dx3(), self.dFdx3dx3()]]
        inverse = numpy.linalg.inv(iJacobi)
        for i in range(0, 3):
            for j in range(0, 3):
                res[i] += -h[j] * inverse[i][j]
        return res

    def CountH(self, x0, x1, h):
        res = [-self.dFdx1(x1), -self.dFdx2(x1), -self.dFdx3(x1)]
        c = (self.dFdx1(x1) * self.dFdx1(x1) + self.dFdx2(x1) * self.dFdx2(x1) + self.dFdx3(x1) * self.dFdx3(x1)) / (self.dFdx1(x0) * self.dFdx1(x0) + self.dFdx2(x0) * self.dFdx2(x0) + self.dFdx3(x0) * self.dFdx3(x0))
        for i in range(0, 3):
            res[i] += h[i] * c
        return res


    def Spusk(self):
        k = 0
        if (len(self.x0) > 3):
            raise IOError("Исключение")
        res_x = self.x0[:]
        while(True):
            ak = self.CountAlpha(res_x, self.Jacobi(res_x))
            nx1 = res_x[0] - ak * self.dFdx1(res_x)
            nx2 = res_x[1] - ak * self.dFdx2(res_x)
            nx3 = res_x[2] - ak * self.dFdx3(res_x)
            if math.sqrt(math.pow(res_x[0] - nx1, 2) + math.pow(res_x[1] - nx2, 2) + math.pow(res_x[2] - nx3, 2)) < self.epsilon:
                break
            res_x[0] = nx1
            res_x[1] = nx2
            res_x[2] = nx3
            k += 1
        print('Наискорейший спуск (итерации): ', k)
        return res_x

    def MethodNewton(self):
        l = 0
        res_x = self.x0[:]
        i = 0
        while(i != 1):
            multiplying = self.MethodNewtonSubFunction(self.Jacobi(res_x))
            for j in range(0, 3):
                res_x[j] -= multiplying[j]
            i += 1
            l += 1
        print('Метод Ньютона (итерации): ', l)
        return res_x

    def MethodGradienov(self):
        j = 0
        h0 = self.Jacobi(self.x0)
        res_x = self.x0[:]
        a = self.CountAlpha(self.x0, h0)
        for i in range(0, 3):
            x1 = [res_x[0] + a * h0[0], res_x[1] + a * h0[1], res_x[2] + a * h0[2]]
            h0 = self.CountH(res_x, x1, h0)
            res_x = x1[:]
            a = self.CountAlpha(self.x0, h0)
            j += 1
        print('Метод градиентов (итерации)', j)
        return res_x

if __name__ == "__main__":
    f = function()
    pointsSp = f.Spusk()
    pointsNewton = f.MethodNewton()
    pointsGradient = f.MethodGradienov()

    print('Наискорейший спуск: ')
    print(pointsSp[0], pointsSp[1], pointsSp[2])
    print('F = ', f.calculate(pointsSp[0], pointsSp[1], pointsSp[2]))

    print('Метод Ньютона: ')
    print(pointsNewton[0], pointsNewton[1], pointsNewton[2])
    print('F = ', f.calculate(pointsNewton[0], pointsNewton[1], pointsNewton[2]))

    print('Сопряженные градиенты: ')
    print(pointsGradient[0], pointsGradient[1], pointsGradient[2])
    print('F = ', f.calculate(pointsGradient[0], pointsGradient[1], pointsGradient[2]))