#! /usr/bin/env python
# -*- coding: utf-8 -*-
"""
Demo of different kinds of Gaussian quadrature to compute a definite integral in 1D.

Reference is: https://en.wikipedia.org/wiki/Gaussian_quadrature.

Note: thanks to Samhita M. from A4 group for sharing her program with me (used as a starting point to write all this).

@date: Mon Apr 13 11:40:55 2015.
@author: Lilian Besson for CS101 course at Mahindra Ecole Centrale 2015.
@licence: MIT Licence (http://lbesson.mit-license.org).
@author:
"""

from math import cos, sin, sqrt, pi
import math


# %% Different strategies

import numpy as np

# Here we are cheating a little bit...
from scipy.special.orthogonal import p_roots

def legendre(f, n):
    """ Easy case for integral of f(x) from -1 to 1, with n points (weight is 1).

    - We use Gauss-Legendre quadrature,
    - Reference is https://en.wikipedia.org/wiki/Gaussian_quadrature#Gauss.E2.80.93Legendre_quadrature.
    """
    # Thanks to http://stackoverflow.com/a/27116609
    xi, wi = p_roots(n+1)  # pylint: disable=W0632
    return np.sum(wi * f(xi))


def chebyshev_1(f, n):
    """ Easy case for integral of f(x)/sqrt(1-x**2) from -1 to 1, with n points.

    - We use Chebyshev-Gauss quadrature (first case)
    - Reference is https://en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature.
    """
    wi = pi / float(n)
    list_xi = [cos(pi * (2.0*i - 1.0) / (2.0*n)) for i in xrange(1, n+1)]
    return sum(wi * f(xi) for xi in list_xi)


def chebyshev_2(f, n):
    """ Easy case for integral of f(x)*sqrt(1-x**2) from -1 to 1, with n points.

    - We use Chebyshev-Gauss quadrature (second case)
    - Reference is https://en.wikipedia.org/wiki/Chebyshev%E2%80%93Gauss_quadrature.
    """
    list_wi = [(pi / (n+1)) * (sin(i * pi / (n+1)))**2 for i in xrange(1, n+1)]
    list_xi = [cos(i * pi / (n+1)) for i in xrange(1, n+1)]
    return sum(wi * f(xi) for xi, wi in zip(list_xi, list_wi))


# %% Unified in one function
def gaussian_quadrature(g, a, b, n, family='chebyshev_1'):
    """ Compute an approximation of the integral of g(x) for x from a to b, with n points.

    - Family specifies what family of points and weigts should be used (legendre, chebyshev_1, chebyshev_2 as far as now).
    - Reference: https://en.wikipedia.org/wiki/Gaussian_quadrature#Change_of_interval.
    """
    assert a <= b
    size = (b - a)/2.0
    aver = (a + b)/2.0
    # We need to integrate the linear change of variable
    # and the weight into the function

    if family == 'legendre':
        f = lambda x: g(size * x + aver)
        return size * legendre(f, n)
    elif family == 'chebyshev_1':
        f = lambda x: g(size * x + aver) * (sqrt(1 - x**2))
        return size * chebyshev_1(f, n)
    elif family == 'chebyshev_2':
        f = lambda x: g(size * x + aver) / (sqrt(1 - x**2))
        return size * chebyshev_2(f, n)
    else:
        raise ValueError("gaussian_quadrature(): unknown family of points and weights {}".format(family))


# %% Examples
if __name__ == '__main__':

    a, b = -4, 13
    print r"Computing \int_{%g}^{%g} \exp(x) dx = ?" % (a, b)

    def f1(x):
        """ Exponential function, aware of numpy arrays."""
        return np.exp(x)

    N = 10
    print "Number of points:", N
    print "Approximative value (Gauss-Legendre quadrature):\n", gaussian_quadrature(f1, a, b, N, family='legendre')


    # Second example
    def f2(x):
        """ Exponential function, for float numbers."""
        return math.exp(x)

    N = 10000
    print "Number of points:", N
    print "Approximative value (first case of Chebyshev-Gauss quadrature):\n", gaussian_quadrature(f2, a, b, N, family='chebyshev_1')
    print "Approximative value (second case of Chebyshev-Gauss quadrature):\n", gaussian_quadrature(f2, a, b, N, family='chebyshev_2')
    print (r"The exact value was: \int_{%g}^{%g} \exp(x) dx = \exp(%g) - \exp(%g) =" % (a, b, b, a)), math.exp(b) - math.exp(a)