#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
Naive Python implementation of Romberg's method for computing 1D integrals, with a recursive function.

The recursive implementation is simpler to understand, but really less efficient that my first implementation (using a dictionary to store the previous values of R(n, m)).

Note: thanks to the student Jithendra Sai V (group) for sharing his initial program.


@reference: https://en.wikipedia.org/wiki/Romberg%27s_method#Method
@date: Sun Apr 05 15:00:27 2015.
@author: Lilian Besson for CS101 course at Mahindra Ecole Centrale 2015.
@licence: MIT Licence (http://lbesson.mit-license.org).
"""

# from __future__ import division
# No issue with integer division ! / is now the mathematically correct /.

import math

def romberg(f, xmin, xmax, n=8, m=8):
    r""" Compute the R(n, m) value recursively, to approximate \int_a^n f(x) dx.

    - The integrand f must be of class C^(2n+2) for the method to be correct.
    - Time complexity is O(n × m).
    - Memory complexity is also O(n × m).
    """
    assert n >= m
    if n == 0 and m == 0:
        return ((xmax-xmin)/2.0)*(f(xmin)+f(xmax))
    elif m == 0:
        h = (xmax-xmin)/float(2**n)
        N = (2**(n-1))+1
        term = math.fsum(f( xmin + ((2*k)-1)*h ) for k in xrange(1, N))
        return (term*h) + (0.5)*romberg(f, xmin, xmax, n-1, 0)
    else:
        return (1.0/((4**m)-1))*( (4**m)*romberg(f, xmin, xmax, n, m-1) - romberg(f, xmin, xmax, n-1, m-1) )


# %% Examples
if __name__ == '__main__':
    print "Performing some examples with this Romberg method (romberg recursive function):"
    n = m = 5
    a = 0
    b = 1

    def f1(x):
        """ Quadrant of a circle."""
        return (1-x**2)**0.5

    print "For the function f1: x --> (1-x**2)**0.5, its integral from", a, "to", b, "is:"
    int1 = romberg(f1, a, b, n, m)
    print "approximatively", int1, "with a Romberg method with n =", n, "and m =", m

    def f2(x):
        """ Upward parabola."""
        return x**2

    print "For the function f2: x --> x**2, its integral from", a, "to", b, "is:"
    int2 = romberg(f2, a, b, n, m)
    print int2, "(expected value is 1/3)"
    print "approximatively with a Romberg method with n =", n, "and m =", m

    f3 = math.sin
    b = math.pi
    print "For the function f3: x --> sin(x), its integral from", a, "to", b, "is:"
    int3 = romberg(f3, a, b, n, m)
    print int3, "(expected value is 2)"
    print "approximatively with a Romberg method with n =", n, "and m =", m

    f4 = math.exp
    n, m = 5, 5
    a, b = -4, 19
    print "For the function f4: x --> exp(x), its integral from", a, "to", b, "is:"
    int4 = romberg(f4, a, b, n, m)
    print int4, "(expected value is exp(19) - exp(-4) ~=", f4(b) - f4(a), ")"
    print "approximatively with a Romberg method with n =", n, "and m =", m

    n, m = 10, 10
    a, b = -4, 19
    print "For the function f4: x --> exp(x), its integral from", a, "to", b, "is:"
    int4 = romberg(f4, a, b, n, m)
    # 10 loops, best of 3: 41.2 ms per loop
    print int4, "(expected value is exp(19) - exp(-4) ~=", f4(b) - f4(a), ")"
    print "approximatively with a Romberg method with n =", n, "and m =", m

    a, b = -1000, 20
    print "For the function f4: x --> exp(x), its integral from", a, "to", b, "is:"
    int4 = romberg(f4, a, b, n, m)
    # 10 loops, best of 3: 71.4 ms per loop
    print int4, "(expected value is exp(20) - exp(-1000) ~=", f4(b) - f4(a), ")"
    print "approximatively with a Romberg method with n =", n, "and m =", m

# End