############################################################################
## CALCULUS - BIOINFORMATICS, 2018-19 T3
## Authors: Toni Guillamon, Alex Haro
############################################################################

import numpy as np

  #RombergScheme(f,a,b,n):
  #    r        is the Romberg integration array; most accurate
  #             answer should be at bottom of right-most column,
  #    f        is the function to integrate,
  #    [a, b]   is the interval of integration, and
  #    n        is the number of levels of recursion to use
  #
  # Adapted from the Octave code by Jonathan R. Senning <jonathan.senning@gordon.edu>
  # Gordon College
  #
  # Based on an algorithm in "Numerical Mathematics and Computing"
  # 4th Edition, by Cheney and Kincaid, Brooks-Cole, 1999.
  # ----------------------------------------------------------------------
  
def RombergScheme(f,a,b,n):    
    r=np.zeros((n+1,n+1))
    h=b-a
    r[0,0]=0.5*h*(f(a)+f(b))
    for i in range(0,n):
        # Compute the halved stepsize and use this to sum the function at
        # all the new points (in between the points already computed)
        h=0.5*h
        sum=0.0
        for k in range(1,2**(i+1),2):
            sum = sum+f(a+k*h)
        # Compute the composite trapezoid rule for the next level of
        # subdivision.  Use Richardson extrapolation to refine these values
        # into a more accurate form.
        r[i+1,0]=0.5*r[i,0]+sum*h
        for j in range(0,i+1):
            print(i,j)
            r[i+1,j+1]=r[i+1,j]+(r[i+1,j]-r[i,j])/(4**(j+1)-1)
    return r

#usage
RombergScheme(lambda x: np.sin(x**2),0.,2.,4)
  
############################################################################
## Octave version
############################################################################

function r = romberg( f, a, b, n )
  % ----------------------------------------------------------------------
  % Jonathan R. Senning <jonathan.senning@gordon.edu>
  % Gordon College
  % February 17, 1999
  %
  % Usage: r = romberg(@f, a, b, n)
  %
  %    where r        is the Romberg integration array; most accurate
  %                   answer should be at bottom of right-most column,
  %          f        is the handle (precede with @) of the function to
  %                   integrate,
  %          [a, b]   is the interval of integration, and
  %          n        is the number of levels of recursion to use
  %
  % Based on an algorithm in "Numerical Mathematics and Computing"
  % 4th Edition, by Cheney and Kincaid, Brooks-Cole, 1999.
  %
  % The algorithm in the text indexes R(n,m) from 0,0 but Octave starts
  % indices at 1 so that is reflected in the code below.
  % ----------------------------------------------------------------------

  r = zeros( n + 1, n + 1 );
  h = b - a;
  r(1,1) = 0.5 * h * ( f( a ) + f( b ) );

  for i = 1 : n

    % Compute the halved stepsize and use this to sum the function at
    % all the new points (in between the points already computed)
    h = 0.5 * h;
    sum = 0.0;
    for k = 1 : 2 : 2^i - 1
      sum = sum + f( a + k * h );
    endfor

    % Compute the composite trapezoid rule for the next level of
    % subdivision.  Use Richardson extrapolation to refine these values
    % into a more accurate form.
    r(i+1,1) = 0.5 * r(i,1) + sum * h;
    for j = 1 : i
      r(i+1,j+1) = r(i+1,j) + ( r(i+1,j) - r(i,j) ) / ( 4^j - 1 );
    endfor

  endfor

endfunction