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 np.sin(4 * x)


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 m(i, )


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))

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.empty(1000)
spline = np.empty(1000)
for i in range(1000):
    points[i] = i
    for m in range(N-1):
        if ((i <= X[m+1]) && (i >= X[m]))
        {
            spline[i] = _s(m+1, i)
            break;
        }
        else
            if (i > X[n - 1])
            {
                spline[i] = _s(N-1, i)
                break;
            }
            else
                if (i < X[0])
                {
                    spline[i] = _s(1, i)
                    break;
                }

plt.plot(points, spline)
plt.show()

# for i in range(100):
#    print(i, ',', np.sin(4 * i), sep='')