import math
from matplotlib import pyplot as plot


class CauchySolver:
    def __init__(self, phi, derivative, x0, xn, k):
        self.phi = phi
        self.derivative = derivative
        self.x0 = x0
        self.xn = xn
        self.k = k

    def get_setka(self, n):
        h = (self.xn - self.x0) / n
        x = [self.x0 + i * h for i in range(n + 1)]
        n = n + 1
        return x, n, h

    def f(self, x, y):
        return self.derivative(x) + self.k * (y - self.phi(x))

    def y_euler(self, n):
        x, n, h = self.get_setka(n)
        y = []
        for i in range(n):
            if i == 0:
                y.append(self.phi(x[i]))
            else:
                y.append(y[i - 1] + h * self.f(x[i - 1], y[i - 1]))
        return y


    def y_adams(self, n):
        y = self.y_runge_cutte_b(n)[:4]
        x, n, h = self.get_setka(n)
        for i in range(4, n):
            f1 = self.f(x[i - 1], y[i - 1])
            f2 = self.f(x[i - 2], y[i - 2])
            f3 = self.f(x[i - 3], y[i - 3])
            f4 = self.f(x[i - 4], y[i - 4])
            y.append(y[i - 1] + h / 24 * (55 * f1 - 59 * f2 + 37 * f3 - 9 * f4))
        return y

    def y_pred_corr(self, n, l):
        x, n, h = self.get_setka(n)
        y = []
        for i in range(n):
            if i == 0:
                y.append(self.phi(x[i]))
            else:
                yp = y[i - 1] + h * self.f(x[i - 1], y[i - 1])
                if l is None:
                    y.append(yp)
                else:
                    while True:
                        yc = y[i - 1] + (self.f(x[i - 1], y[i - 1]) + self.f(x[i], yp)) * h / 2
                        if abs(yc - yp) / abs(yp) <= l:
                            y.append(yc)
                            break
                        yp = yc

        return y

    def y_abm(self, n):
        y = self.y_runge_cutte_a(n)[:4]
        x, n, h = self.get_setka(n)
        for i in range(4, n):
            f1 = self.f(x[i - 1], y[i - 1])
            f2 = self.f(x[i - 2], y[i - 2])
            f3 = self.f(x[i - 3], y[i - 3])
            f4 = self.f(x[i - 4], y[i - 4])
            yc = y[i - 1] + h / 24 * (55 * f1 - 59 * f2 + 37 * f3 - 9 * f4)
            y.append(y[i - 1] + h  / 24 * (9 * self.f(x[i], yc) + 19 * f1 - 5 * f2 + f3))
        return y

    def y_runge_cutte_a(self, n):
        x, n, h = self.get_setka(n)
        y = []
        for i in range(n):
            if i == 0:
                y.append(self.phi(x[i]))
            else:
                k1 = self.f(x[i - 1], y[i - 1])
                k2 = self.f(x[i - 1] + h / 2, y[i - 1] + k1 * h / 2)
                k3 = self.f(x[i - 1] + h, y[i - 1] - k1 * h + 2 * k2 * h)
                y.append(y[i - 1] + h * (k1 + 4 * k2 + k3) / 6)
        return y

    def y_runge_cutte_b(self, n):
        x, n, h = self.get_setka(n)
        y = []
        for i in range(n):
            if i == 0:
                y.append(self.phi(x[i]))
            else:
                k1 = self.f(x[i - 1], y[i - 1])
                k2 = self.f(x[i - 1] + h / 2, y[i - 1] + k1 * h / 2)
                k3 = self.f(x[i - 1] + h / 2, y[i - 1] +  k2 * h / 2)
                k4 = self.f(x[i - 1] + h, y[i - 1] + k3 * h)
                y.append(y[i - 1] + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6)
        return y

    def y_gira(self, n, eps=1e-3):
        y = self.y_runge_cutte_b(n)[:4]
        x, n, h = self.get_setka(n)
        for i in range(4, n):
            yy = 4 * h * self.f(x[i - 1], y[i - 1]) + (y[i - 4] - 10 * y[i - 1]) / 3 - 2 * y[i - 3] + 6 * y[i - 2]
            y.append((48 * y[i - 1] - 36 * y[i - 2] + 16 * y[i - 3] - 3 * y[i - 4] + 12 * h * self.f(x[i], yy)) / 25)
        return y

    def local_regarding_eps(self, y):
        x, n, h = self.get_setka(len(y) - 1)
        eps = []
        for i in range(n):
            phi_i = self.phi(x[i])
            eps.append(abs(phi_i - y[i]) / max(abs(phi_i), abs(y[i])))
        return eps

    def absolute_eps(self, y):
        x, n, h = self.get_setka(len(y) - 1)
        eps = []
        for i in range(n):
            eps.append(abs(self.phi(x[i]) - y[i]))
        return eps

    def valuation_runge(self, yh1, yh2, p):
        diff = []
        for i in range(len(yh1)):
            diff.append((2 ** p / (2 ** p - 1)) * abs(- yh2[2 * i] + yh1[i]))
        return diff

    def plot(self, y):
        x, n, h = self.get_setka(len(y) - 1)
        phi_x = [self.x0 + i * (self.xn - self.x0) / 1000 for i in range(1000)]
        phi_y = [self.phi(xi) for xi in phi_x]
        plot.plot(phi_x, phi_y, color='red', label='phi')
        plot.plot(x, y, color='blue', label='y')
        plot.show()

def plot_emax(f, error, p, **kwargs):
    a = 100
    b = 1000
    step = 1

    label = kwargs['label']
    del kwargs['label']

    c = 20
    x = []
    e_max = []
    sf = []
    for n in range(a, b, step):
        x.append(n)
        e_max.append(max(error(f(n=n, **kwargs))))
        # sf.append(c * n ** (-p))

    plot.xscale('log')
    plot.yscale('log')
    plot.plot(x, e_max, label = label)
    # plot.plot(x, sf, color='red')
    # plot.show()


def phi(x):
    return math.cos(x ** 3)

def phi_derivative(x):
    return -math.sin(x ** 3) * 3 * x ** 2

def main():
    n = 10
    p = 2

    solver = CauchySolver(phi=phi, derivative=phi_derivative, x0=-math.pi, xn=math.pi, k=0.8)

    # y = solver.y_runge_cutte_b(n=n)
    # y = solver.y_adams(n=n)
    # y2= solver.y_adams(n=2*n)
    # y = solver.y_pred_corr(l = 0.01, n=n)
    # y2 = solver.y_pred_corr(l = 0.01, n=n)
    # y = solver.y_abm(n=n)
    # y2 = solver.y_abm(n=2* n)
    # y = solver.y_euler(n=n)
    # y2 = solver.y_euler(n = 2*n)
    # y = solver.y_runge_cutte_a(n=n)
    # y = solver.y_gira(n=n)
    # y2 = solver.y_gira(n=2*n)
    # y2 = solver.y_runge_cutte_a(n=2 * n)
    # pz = math.log( abs(max(solver.local_regarding_eps(y)) / max(solver.local_regarding_eps(y2))), 2)
    # pz2 = math.log( abs(max(solver.absolute_eps(y)) / max(solver.absolute_eps(y2))), 2)


    # n = 10
    # e_max = None

    # while True:
        # y = solver.y_euler(n=n)
        # e = max(solver.absolute_eps(y))
        # if e_max is not None:
            # if abs(e - e_max) < 10e-3:
                # print(n, (solver.xn - solver.x0) / n)
                # break

        # n = n * 2
        # e_max = e
# графики е_max
    # aa = [(solver.y_euler, 1, 'эйлер'),(solver.y_adams, 4, 'адамс'),(solver.y_pred_corr, 2, 'итерац'),(solver.y_abm, 4, 'АБМ'),(solver.y_runge_cutte_a, 3, 'РК(3)'),(solver.y_runge_cutte_b, 4, 'РК(4)'), (solver.y_gira, 4, 'Гира)')]
    # for f, p, l in aa:
        # plot_emax(f=f, error=solver.absolute_eps, p=p, label=l)
    # plot.legend()
    # plot.show()
# итерац.метод
    plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=None, label='10000')
    plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-10, label='-10')
    plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-8, label='-8')
    plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-6, label='-6')
    plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-4, label='-4')
    plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-2, label='-2')
    plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-1, label='-1')
    plot_emax(f=solver.y_pred_corr, error=solver.absolute_eps, p=p, l=1e-0, label='-0')
    plot.legend()
    plot.show()

    # for i in range(20):
        # y = solver.y_abm(n=n)
        # y2 = solver.y_abm(n=2*n)
        # y = solver.y_gira(n=n)
        # y2 = solver.y_gira(n=2*n)
        # y = solver.y_runge_cutte_b(n=n)
        # y2 = solver.y_runge_cutte_b(n=2*n)
        # y = solver.y_runge_cutte_a(n=n)
        # y2 = solver.y_runge_cutte_a(n=2*n)
        # y = solver.y_euler(n=n)
        # y2 = solver.y_euler(n= 2 * n)
        # y = solver.y_adams(n=n)
        # y2 = solver.y_adams(n=2*n)
        # y = solver.y_pred_corr(l = 10 ** (-1), n=n)
        # y2 = solver.y_pred_corr(l = 10 ** (-1), n=2 * n)
        # print(solver.local_regarding_eps(y))
        # pz = math.log( abs(max(solver.local_regarding_eps(y)) / max(solver.local_regarding_eps(y2))), 2)
        # pz2 = math.log( abs(max(solver.absolute_eps(y)) / max(solver.absolute_eps(y2))), 2)
        # print('p* (loc) =  %.2E, p* (abs) = %.2E' % (pz, pz2))
        # print('max локальная относительная %.2E, max абсолютная %.6E' % (max(solver.local_regarding_eps(y)), max(solver.absolute_eps(y))))
        # print('%d оценка по правилу Рунге %.2E, %.2E' % (n, max(solver.valuation_runge(y, y2, p)), max(solver.absolute_eps(y))))
        # plot_emax(f=solver.y_gira, error=solver.absolute_eps, p=p)
        # n = n * 2
    # print('оценка по правилу Рунге %.6E' % (max(solver.valuation_runge(y, y2, p))))
    # print('p* (loc) =  %.6f, p* (abs) = %.6f' % (pz, pz2))
    # print('max локальная относительная %.2E, max абсолютная %.2E' % (max(solver.local_regarding_eps(y)), max(solver.absolute_eps(y))))

    # solver.plot(y)
    # solver.plot(solver.absolute_eps)
    # plot_emax(f=solver.y_gira, error=solver.absolute_eps, p=p)

    # l10 = [3.14856056E+02, 9.61512887E+02, 3.60834374E+02, 6.84264503E+01, 1.60062067E+01, 3.93924985E+00, 9.81012935E-01, 2.45017209E-01, 6.12395731E-02, 1.53089731E-02, 3.82718587E-03, 9.56793005E-04, 2.39198143E-04, 5.97994648E-05, 1.49499250E-05, 3.73755566E-06, 9.34497514E-07, 2.33762336E-07, 5.83910977E-08, 1.47224043E-08]
    # l6 = [3.14855847E+02, 9.61511494E+02, 3.60834328E+02, 6.84264467E+01, 1.60062119E+01, 3.93924970E+00, 9.81014023E-01, 2.45016802E-01, 6.12410550E-02, 1.53098508E-02, 3.82874286E-03, 9.57928022E-04, 2.38187413E-04, 6.02441899E-05, 1.44482460E-05, 3.61183274E-06, 9.02953146E-07, 2.25740056E-07, 5.64369144E-08, 1.41134897E-08]
    # l4 = [3.14815323E+02, 9.61416445E+02, 3.60816867E+02, 6.84247258E+01, 1.60066843E+01, 3.93913479E+00, 9.81308365E-01, 2.45172731E-01, 6.12491721E-02, 1.54368763E-02, 3.90645617E-03, 9.24641114E-04, 2.31157243E-04, 5.77892692E-05, 1.44472601E-05, 3.61180901E-06, 9.02952284E-07, 2.25740057E-07, 5.64368925E-08, 1.41134897E-08]
    # l2 = [3.12414690E+02, 9.54277666E+02, 3.58308448E+02, 6.84714427E+01, 1.60458342E+01, 3.93507879E+00, 9.65446205E-01, 2.26064204E-01, 5.91838519E-02, 1.47948051E-02, 3.69854349E-03, 9.24627030E-04, 2.31156215E-04, 5.77890156E-05, 1.44472457E-05, 3.61180909E-06, 9.02952284E-07, 2.25740057E-07, 5.64368925E-08, 1.41134897E-08]

    # x = [x for x in range(len(l10))]

    # plot.plot(x, l10)
    # plot.plot(x, l6)
    # plot.plot(x, l4)
    # plot.plot(x, l2)
    # plot.show()

main()