import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

data = pd.read_csv('input.csv')
x = data['x'].tolist()
f = data['f(x)'].tolist()
start = float(data['a'][0])
end = float(data['b'][0])
df0 = 0
dfn = -0.939692621


def func(x):
    return x**3 + x**2 + x + 1


def h(i):
    return x[i] - x[i - 1]


def coeffs_vector(A, B):
    alpha = np.empty(shape=(N))
    beta = np.empty(shape=(N))
    cv = np.empty(shape=(N))

    y = A[0][0]
    alpha[0] = -A[0][1] / y
    beta[0] = B[0] / y

    for i in range(1, N - 1):
        y = A[i][i] + A[i][i - 1] * alpha[i - 1]
        alpha[i] = -A[i][i + 1] / y
        beta[i] = (B[i] - A[i][i - 1] * beta[i - 1]) / y

    cv[N - 1] = (B[N - 1] - A[N - 1][N - 2] * beta[N - 2]) / (A[N - 1][N - 1] + A[N - 1][N - 2] * alpha[N - 2])
    for i in range(N - 2, 0, -1):
        cv[i] = alpha[i] * cv[i + 1] + beta[i]



    return cv


def s(i, arg):
    cv = coeffs_vector(A, B)
    return (cv[i - 1] * func(x[i] - arg)) / (6 * (x[i] - x[i - 1])) + (cv[i] * func(arg - x[i - 1])) / (6 * (x[i] - x[i - 1])) + (f[i - 1] - (cv[i - 1] * (x[i] - x[i - 1]) * (x[i] - x[i - 1])) / 6) * ((x[i] - arg) / (x[i] - x[i - 1])) + (f[i] - (cv[i] * (x[i] - x[i - 1]) * (x[i] - x[i - 1]) / 6)) * ((arg - x[i - 1]) / (x[i] - x[i - 1]))


def _s(i, arg):
    cv = coeffs_vector(A, B)
    a = M[i-1] * (x[i] - arg)**3 / 6 * h(i)
    b = M[i] * (arg - x[i-1])**3 / 6 * h(i)
    c = (f[i-1] - M[i]*h(i)**2 / 6) * ((x[i] - arg) / h(i))
    d = (f[i] - M[i]*h(i)**2 / 6) * ((arg - x[i-1]) / h(i))
    return a + b + c + d

# print('Input data\n~~~~~~~~~~~~~')
# print(data)


N = len(x)
A = np.zeros(shape=(N,N))
B = np.empty(shape=(N))
M = np.empty(shape=(N))


def main():
    # ================ Заполнение матриц ===============


    A[0][0] = (x[1] - x[0]) / 3
    A[0][1] = (x[1] - x[0]) / 6
    A[N - 1][N - 2] = (x[N - 1] - x[N - 2]) / 6
    A[N - 1][N - 1] = (x[N - 1] - x[N - 2]) / 3
    for i in range(1, N - 1):
        for j in range(0, N):
            if j < i - 1 or j > i + 1:
                A[i][j] = 0
            elif j == i - 1:
                A[i][j] = (x[i] - x[i - 1]) / 6
            elif j == i:
                A[i][j] = (x[i] - x[i - 1] + x[i+1] - x[i]) / 3
            elif j == i + 1:
                A[i][j] = (x[i + 1] - x[i]) / 6


    for i in range(N - 1):
        B[i] = (f[i + 1] - f[i]) / (x[i + 1] - x[i]) - (f[i] - f[i - 1]) / (x[i] - x[i - 1])

    M[0] = ((A[1][1] - A[0][1]) / h(0) - df0 - h(0) * B[0] / 6) * (3 / h(0))
    for i in range(N - 2):
        M[i + 1] = B[i]
    M[N - 1] = dfn
    # print(A)
    # print(B)

    points = np.linspace(-1, 2, 0.01)
    spline = np.empty(300)
    k = 0
    for i in range(0, 300):
        points[k] = i
        for m in range(N-1):
            if (i <= x[m+1]) and (i >= x[m]):
                spline[i] = _s(m+1, i)
                break
            else:
                if i > x[N - 1]:
                    spline[k] = _s(N-1, i)
                    break
                else:
                    if i < x[0]:
                        spline[k] = _s(1, i)
                        break
        k += 1

    #plt.plot(np.linspace(0, 100, 10), [func(i) for i in range(0, 100, 10)], color='red')
    plt.plot(points, spline)
    plt.show()


main()
#for i in range(100):
#    print(i, ',', func(i), sep='')