---

from collections import defaultdict
import math
import numpy as np

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


# n - число точек! Не степень!
COMPARE_N = 101  # число точек для сравнения по норме (значения записываются в файл)
PLOT_N = 101  # число точек для построения графика
A, B = -3, 3  # границы
# N_LIST = [3, 10, 20]  # число узлов для сравнений
N_LIST = [i for i in range(3, 31)]  # число узлов для сравнений
DO_PLOT = False  # если True, откроет matplotlib окна с графиками


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


def f2(x):
    return abs(abs(x)-1)


def chebyshev_nodes_on_interval(a, b, n):
    k = np.arange(1, n + 1)
    nodes = (a + b) / 2 + ((b - a) / 2) * np.cos((2*k - 1) * np.pi / (2*n))
    return list(nodes)


def newton_interpolation(points):
    """
    Интерполирует функцию методом Ньютона по массиву координат
    Возвращает функцию для приблизительного рассчета значения в точке
    """
    x, y = list(zip(*points))  # unzipping points
    x, y = list(x), list(y)
    n = len(x)
    table = y.copy()  # строим таблицу разделенных разностей

    # считаем по 1 столбцу
    # храним элементы "прижатыми" к концу массива
    # в начале массива сохраняем коэффициенты для будущего многочлена
    for j in range(1, n):
        for i in range(n-1, j-1, -1):
            table[i] = (table[i] - table[i-1]) / (x[i] - x[i-j])

    # noinspection PyShadowingNames
    def func(t):
        res = table[-1]
        # считаем значение многочлена
        for i in range(n-2, -1, -1):
            res = res * (t - x[i]) + table[i]
        return res

    return func


def count_func(f, a, b, n, numbers_by_interval_func=np.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 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)
    if DO_PLOT:
        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'
        })
        cur_plt.legend(prop={'size': 6})
    return eps


def write_appr_func_to_file(f, func_name, a, b, n):
    func_points = count_func(f, a, b, COMPARE_N, numbers_by_interval_func=np.linspace)
    with open(f'output/{func_name}_{n}_appr.txt', 'w') as f:
        for x, y in func_points:
            f.write(f'{x} {y}\n')


def main(func, func_name):
    if DO_PLOT:
        figure, axis = plt.subplots(len(N_LIST), 2)
        plt.subplots_adjust(left=0.05, bottom=0.05, right=0.95, top=0.95, hspace=0.5)
    else:
        axis = defaultdict(lambda: plt)

    # 101 точное значений функции на равностоящих узлах (запись в файл)
    func_prec_points = count_func(func, A, B, COMPARE_N, numbers_by_interval_func=np.linspace)
    with open(f'output/{func_name}_{COMPARE_N}_prec.txt', 'w') as f:
        for x, y in func_prec_points:
            f.write(f'{x} {y}\n')

    # равностоящие узлы
    # генерируем значения в точках, в которых будем строить интерполяцию
    # n - число точек! Не степень!
    for i, n in enumerate(N_LIST):
        # P1/2
        base_points = count_func(func, A, B, n, numbers_by_interval_func=np.linspace)
        func_approx = newton_interpolation(base_points)
        eps_p = compare_and_plot(func, func_approx, A, B, plot_n=PLOT_N, compare_n=COMPARE_N,
                                 cur_plt=axis[i, 0])
        write_appr_func_to_file(func_approx, f'P{func_name[1:]}', A, B, n)
        if DO_PLOT:
            # noinspection PyUnresolvedReferences
            axis[i, 0].set_title(f'Равноотстоящие, n={n}')
        # C1/2
        base_points = count_func(func, A, B, n, numbers_by_interval_func=chebyshev_nodes_on_interval)
        func_approx = newton_interpolation(base_points)
        eps_c = compare_and_plot(func, func_approx, A, B, plot_n=PLOT_N, compare_n=COMPARE_N,
                                 cur_plt=axis[i, 1])
        write_appr_func_to_file(func_approx, f'C{func_name[1:]}', A, B, n)
        if DO_PLOT:
            # noinspection PyUnresolvedReferences
            axis[i, 1].set_title(f'Чебышевские, n={n}')
        print(f'n={n} |dP|={eps_p:.8f} |dC|={eps_c:.8f}')

    if DO_PLOT:
        plt.show()


if __name__ == '__main__':
    main(f1, 'f1')
    # main(f2, 'f2')