---

from math import sin, cos
import numpy as np
from scipy.interpolate import CubicSpline

import matplotlib
matplotlib.use('TkAgg')
import matplotlib.pyplot as plt


# n - число точек! Не степень!
BUILD_N = 15 + 1  # число точек для построения сплайна
COMPARE_N = 100 + 1  # число точек для сравнения по норме (значения записываются в файл)
A, B = -3, 3  # границы


def f1(x):
    return sin(cos(x))


def d_f1(x):
    return -sin(x) * cos(cos(x))


def d2_f1(x):
    return -sin(cos(x))*(sin(x)**2) - cos(cos(x))*cos(x)


def linspace(a, b, n):
    """
    Возвращает n равностоящих узлов на отрезке [a, b]
    """
    return np.linspace(a, b, n)


def count_func(f, a, b, n, numbers_by_interval_func=linspace):
    """
    Считает значение функции на отрезке [a, b] в каждой из n точек
    Использует функцию-генератор узлов (numbers_by_interval_func)
        по умолчанию np.linspace - равностоящие узлы
    Возвращает массив пар вида [(x1, y1), (x2, y2), ...]
    """
    x_values = list(numbers_by_interval_func(a, b, n))
    y_values = list(map(f, x_values))
    return list(zip(x_values, y_values))


def plot_func(f, a, b, n, cur_plt=plt, plot_kwargs=None):
    """
    Строить функцию f на отрезке a, b по n равностоящим точкам в matplotlib
    """
    if plot_kwargs is None:
        plot_kwargs = {}
    x_values = np.linspace(a, b, n)
    y_values = list(map(f, x_values))
    cur_plt.plot(x_values, y_values, **plot_kwargs)


def build_spline(f, a, b, n=BUILD_N):
    """
    Строит функцию - кубический сплайн по отрезку и числу точек
    Возвращает функцию для приблизительного рассчета значения в точке
    """
    # считаем значения в n точках
    func_prec_points = count_func(f, a, b, n, numbers_by_interval_func=linspace)
    # в левой границе совпадает первая производная
    d_f1_a = d_f1(a)
    # в правой границе совпадает вторая производная
    d2_f1_b = d2_f1(b)
    # unzipping points
    x, y = list(zip(*func_prec_points))
    x, y = list(x), list(y)
    # считаем сплайн
    return CubicSpline(x, y, bc_type=((1, d_f1_a), (2, d2_f1_b)))


def build_matr(points, a, b, n):
    # unzipping points
    x, y = list(zip(*points))
    x, y = list(x), list(y)

    # в левой границе совпадает первая производная
    d_f1_a = d_f1(a)
    # в правой границе совпадает вторая производная
    d2_f1_b = d2_f1(b)

    h = (b-a) / n
    # совпадение по первой производной
    ldiag = []
    cdiag = [h/3]
    rdiag = [h/6]
    vals = []
    vals.append((y[1] - y[0])/h - d_f1_a)

    for i in range(1, n):
        ldiag.append(h/6)
        cdiag.append(2*h/3)
        rdiag.append(h/6)
        vals.append(
            (y[i+1] - 2 * y[i] + y[i-1]) / h
        )
    # совпадение по второй производной
    ldiag.append(0)
    cdiag.append(1)
    vals.append(d2_f1_b)
    return [ldiag, cdiag, rdiag, vals]


def get_moments(matr, n):
    x_i, r_i = 0, 0
    coeff = []
    # прямой ход
    for i in range(n+1):
        a = 0 if i == 0 else matr[0][i - 1]
        b = -1 * matr[1][i]
        c = 0 if i == n else matr[2][i]

        delta = b - a * x_i
        x_i = c / delta
        r_i = (a * r_i - matr[3][i]) / delta
        coeff.append([x_i, r_i])
    # обратный ход
    x = 0
    res = [0]*(n+1)
    for i in range(n, -1, -1):
        x = coeff[i][0] * x + coeff[i][1]
        res[i] = x
    return res


def my_build_spline(f, a, b, n=BUILD_N-1):
    """
    Строит функцию - кубический сплайн по отрезку и степени(!!!)
    Возвращает функцию для приблизительного рассчета значения в точке
    """
    func_prec_points = count_func(f, a, b, n+1, numbers_by_interval_func=linspace)

    matr = build_matr(func_prec_points, a, b, n)
    moments = get_moments(matr, n)

    # unzipping points
    x_prec, y_prec = list(zip(*func_prec_points))
    x_prec, y_prec = list(x_prec), list(y_prec)

    def func(x):
        pos = 0  # промужеток, которому принадлежит x
        h = (b-a) / n
        for i in range(1, n+1):
            if x_prec[i - 1] <= x <= x_prec[i]:
                pos = i
                break
        res = 0
        res += moments[pos-1] * pow((x_prec[pos]-x), 3) / (6*h)
        res += moments[pos] * pow(x-x_prec[pos-1], 3) / (6*h)
        res += (y_prec[pos-1] - (h*h/6)*moments[pos-1]) * (x_prec[pos]-x)/h
        res += (y_prec[pos] - (h*h/6)*moments[pos]) * (x-x_prec[pos-1])/h
        return res
    return func


def compare_funcs_and_count_eps(base_func, approx_func, a, b, n):
    """
    Считает норму двух функций по числу равностоящих точек
    """
    base_points = count_func(base_func, a, b, n)
    base_values = [e[1] for e in base_points]  # только y
    approx_points = count_func(approx_func, a, b, n)
    approx_values = [e[1] for e in approx_points]  # только y
    eps = max([abs(base_value - approx_value)
               for base_value, approx_value
               in zip(base_values, approx_values)])
    return eps


def compare_and_plot(base_func, approx_func, a, b, plot_n, compare_n, cur_plt=plt):
    """
    Считает норму для двух функций на отрезке с заданным числом равностоящих узлов
    После чего строит обе функции на равностоящих узлах
    :param base_func: точная функция
    :param approx_func: приближенная функция
    :param a: левая граница
    :param b: правая граница
    :param plot_n: число узлов для построения графика
    :param compare_n: число узлов для расчёта нормы
    :param cur_plt: matplotlib.pyplot на котором строить графики
    :return: максимум-норму на compare_n равностоящих узлах
    """
    eps = compare_funcs_and_count_eps(base_func, approx_func, a, b, compare_n)
    plot_func(base_func, a, b, plot_n, cur_plt, plot_kwargs={
        'label': 'base'
    })
    plot_func(approx_func, a, b, plot_n, cur_plt, plot_kwargs={
        'label': 'approximate'
    })
    # plot_func(lambda e: abs(approx_func(e) - base_func(e)), a, b, plot_n, cur_plt, plot_kwargs={
    #     'label': 'eps'
    # })
    cur_plt.legend(prop={'size': 6})
    return eps


def main():
    # строим сплайн
    # func_spline = build_spline(f1, A, B, BUILD_N)
    func_spline = my_build_spline(f1, A, B, BUILD_N-1)
    # сравниваем в равностоящих узлах
    eps_p = compare_and_plot(f1, func_spline, A, B, plot_n=COMPARE_N, compare_n=COMPARE_N,
                             cur_plt=plt)
    print(f'Погрешность: {eps_p}')
    plt.show()


if __name__ == '__main__':
    main()