#!/usr/bin/python3

import numpy
import math
import matplotlib.pyplot as plt
from operator import mul
from functools import reduce

polynomial_power = 5

function = math.sin
function_domain = [0.0, 2.0*math.pi]


def compute_vandermonde_L_inv(x):
    size = (len(x), len(x))
    L_inv = numpy.matrix([[0.0 for i in range(size[0])]
        for j in range(size[1])])

    for i in range(size[0]):
        for j in range(size[1]):
            if i == j == 0:
                L_inv.itemset(i, j, 1.0)
            elif i < j:
                L_inv.itemset(i, j, 0.0)
            else:
                ks = list(range(i + 1))
                del(ks[j])
                product_elements = [1.0 / (x[j] - x[k]) for k in ks]
                product = reduce(mul, product_elements)
                L_inv.itemset(i, j, product)
    return L_inv

def compute_vandermonde_U_inv(x):
    size = (len(x), len(x))
    U_inv = numpy.matrix([[0.0 for i in range(size[0])]
        for j in range(size[1])])

    for i in range(size[0]):
        for j in range(size[0]):
            print(i, j)
            if i == j:
                U_inv.itemset(i, j, 1.0)
            elif j == 0:
                U_inv.itemset(i, j, 0.0)
            else:
                up_left = 0 if i == 0 else U_inv.item(i-1, j-1)
                U_inv.itemset(i, j, up_left -\
                    U_inv.item(i, j-1) * x[j-1])
    return U_inv


# compute N 2d-space points of given function
function_samples_count = polynomial_power + 1
function_samples_x = numpy.linspace(function_domain[0],
                    function_domain[1], num=function_samples_count)
function_samples_y = [function(x) for x in function_samples_x]


# solve system of linear equations

# convert list of function values to matrix
Y = numpy.matrix([[y] for y in function_samples_y])

# decompose Vandermonde matrix to upper right and lower left matrices
# and get their inverse
U_inv = compute_vandermonde_U_inv(function_samples_x)
L_inv = compute_vandermonde_L_inv(function_samples_x)

#find coefficients
C = U_inv * L_inv * Y

# convert 1-column matrix to list
C = [C.item(i, j) for i in range(C.shape[0]) for j in range(C.shape[1])]

# show calculated polynomial coefficients
print('coefficients:', C)

# plot
plot_x = numpy.linspace(function_domain[0], function_domain[1], num=1000)
plot_y = [sum([float(c)*(x**p) for c, p in zip(C, range(0, polynomial_power + 1))])
            for x in plot_x]

plt.plot(function_samples_x, function_samples_y,'--')
plt.plot(plot_x, plot_y)
plt.show()