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 (Q_table == None):
    Q_table = np.zeros((n + 1, n + 1));

  for i in np.arange(0, n + 1):
    Q_table[i][0] = vec_f[i];
  
  for j in np.arange(1, n + 1):
    for i in np.arange(j, n + 1):
      # compute Q_{i,j}
      Q_table[i][j]  = 0.0;
      Q_table[i][j] += (x - vec_x[i - j])*Q_table[i][j - 1];
      Q_table[i][j] -= (x - vec_x[i])*Q_table[i - 1][j - 1];
      Q_table[i][j] /= (vec_x[i] - vec_x[i - j]);
                   
  return Q_table[n][n], 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');
  xx = np.linspace(-math.pi,math.pi,int(1e2));  
  yy = np.cos(xx);
  plt.plot(xx, yy, '-', linewidth=1.0, markersize=12, 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 "================================================================================"