import math
import numpy as np
import matplotlib.pyplot as plt

def neville(x, vec_x, vec_f, Q_table = None, i0 = -1, j0 = -1):
    n = np.size(vec_x) - 1;  # x0, x1, ..., xn.

    if (i0 == -1):
        i0 = n;

    if (j0 == -1):
        j0 = n;

    if (Q_table == None):
        Q_table = np.zeros((n + 1, n + 1));

    if (n <= 0):
        Q_table[i0][j0] = vec_f[0];
    else:
        for i in range(1, (n) + 1):
            for j in range(1, (i) + 1):
                # compute Q_{i,j-1}
                vec_x_tilde = vec_x[i - (j - 1):(i) + 1];
                vec_f_tilde = vec_f[i - (j - 1):(i) + 1];
                Q_i_jm1, Q_table = neville(x, vec_x_tilde, vec_f_tilde, Q_table, i, j - 1);
                
                # compute Q_{i-1,j-1}
                vec_x_tilde = vec_x[(i - 1) - (j - 1):(i - 1) + 1];
                vec_f_tilde = vec_f[(i - 1) - (j - 1):(i - 1) + 1];
                Q_im1_jm1, Q_table = neville(x, vec_x_tilde, vec_f_tilde, Q_table, i - 1, j - 1);
            
                # compute Q_{i,j}
                Q_table[i][j]  = 0.0;
                Q_table[i][j] += (x - vec_x[i - j])*Q_i_jm1;
                Q_table[i][j] -= (x - vec_x[i])*Q_im1_jm1;
                Q_table[i][j] /= (vec_x[i] - vec_x[i - j]);
                     
    return Q_table[i0][j0], Q_table;
        
if __name__ == "__main__":
    # run the module as script

    print("================================================================================");
    print("Neville's Method");
    print("Date: November, 2014.");
    print("Author: Paul J. Atzberger.");
    print("--------------------------------------------------------------------------------");
    vec_x = np.linspace(-math.pi,math.pi,10);  
    vec_f = np.cos(vec_x);
    x     = math.pi/4.0;

    print(" ");
    print("Interpolating polynomial will be determined using the data:");
    np.set_printoptions(precision=4)
    print("vec_x: ");
    print(vec_x);
    print("vec_f: ");
    print(vec_f);
    print("x: ");
    print("%.4e"%x);

    # Compute the interpolation
    print(" ");
    print("Computing the interpolating polynomial using Neville's Method.");
    P_x, Q_table = neville(x, vec_x, vec_f);
    #P_x = 1.3;

    print(" ");
    print("Interpolating polynomial P(x) has value:");
    print("P(%.4e) = %.4e"%(x,P_x));
    print(" ");
    print("Q_table has value:");
    print(Q_table);

    # Plot the results
    print(" ");
    print("Plottting the results.")
    plt.figure(1, facecolor='white');
    plt.clf();
    plt.plot(vec_x, vec_f, '.', linewidth=1.0, markersize=12, color='blue');
    plt.plot(vec_x, vec_f, '-', linewidth=1.0, color='blue');
    plt.plot(x, P_x, '.', markersize=15, color='red');
    plt.xlabel('x');
    plt.ylabel('y');
    plt.title("Neville's Method");
    plt.show();

    print("Done.")
    print "================================================================================"