#! /usr/bin/env python
# -*- coding: utf-8 -*-
"""
Examples of several algorithm for numerical integration problems.
See the `integrals.py <integrals.html>`_ file for more details.


- *Date:* Saturday 18 June 2016, 18:59:23.
- *Author:* `Lilian Besson <https://bitbucket.org/lbesson/>`_ for the `CS101 course <http://perso.crans.org/besson/cs101/>`_ at `Mahindra Ecole Centrale <http://www.mahindraecolecentrale.edu.in/>`_, 2015.
- *Licence:* `MIT Licence <http://lbesson.mit-license.org>`_, © Lilian Besson.
"""

from __future__ import division, print_function, absolute_import  # Python 2/3 compatibility
import math
from integrals import *


# Do some quick tests (it should be in the 'tests.py' file):
if __name__ == '__main__':
    print("\nSome examples of numerical integrations:\n")

    # # %% Example 1
    f = math.cos
    print("Plotting and computing the integral of cos(x) from 0 to pi.")
    print("Formally, the value is about 0 (= sin(pi) - sin(0)).")
    print("With 15 left rectangles:", riemann_left(f, 0, math.pi, n=15))
    print("With 15 right rectangles:", riemann_right(f, 0, math.pi, n=15))
    print("With 15 center rectangles:", riemann_center(f, 0, math.pi, n=15))
    print("With 3 trapezoides:", trapez(f, 0, math.pi, n=3))

    # print("With 15 left rectangles:", plot_riemann_left(f, 0, math.pi,
    #     namef=r"$x \mapsto \cos(x)$",n=15, figname="riemannleft.png"))
    # print("With 15 right rectangles:", plot_riemann_right(f, 0, math.pi,
    #     namef=r"$x \mapsto \cos(x)$",n=15, figname="riemannright.png"))
    # print("With 15 center rectangles:", plot_riemann_center(f, 0, math.pi,
    #     namef=r"$x \mapsto \cos(x)$",n=15, figname="riemanncenter.png"))
    # print("With 3 trapezoides:", plot_trapez(f, 0, math.pi,
    #     namef=r"$x \mapsto \cos(x)$",n=3, figname="trapezoides.png"))

    # # %% Example 2
    # f = lambda x: x**2
    # print("Plotting and computing the integral of x**2 from 0 to 3.")
    # print("Formally, the value is about 9 (= 3**3 / 3).")
    # print("With 8 left rectangles:", plot_riemann_left(f, 0, 3,
    #     namef="$x \mapsto x^2$", n=8, figname="riemannleft2.png"))
    # print("With 8 right rectangles:", plot_riemann_right(f, 0, 3,
    #     namef="$x \mapsto x^2$", n=8, figname="riemannright2.png"))
    # print("With 8 center rectangles:", plot_riemann_center(f, 0, 3,
    #     namef="$x \mapsto x^2$", n=8, figname="riemanncenter2.png"))
    # print("With 3 trapezoides:", plot_trapez(f, 0, 3,
    #     namef="$x \mapsto x^2$", n=3, figname="trapezoides2.png"))

    # %% Example 3 with Monte-Carlo
    n = 10000

    xmin, xmax = 1, 6

    def f1(x):
        """ f1(x) = x"""
        return x

    # For this particular example, its easy:
    ymin, ymax = xmin, xmax

    print("\nFor this function f1: x → x, on I = [{}, {}], we will pick {} random points.".format(xmin, xmax, n))
    print("Manually, I chose ymin = {}, ymax = {}.".format(ymin, ymax))

    random.seed(1)
    intf = montecarlo(f1, xmin, xmax, n, ymin, ymax)
    print("This leads to an approximated integral value of {}.".format(intf))
    print("Error % is", 100 * abs(17.5 - intf))

    ymin, ymax = yminmax(f1, xmin, xmax, n)  # guess!
    print("But experimentally, I found a possible ymin, ymax to be {}, {}.".format(ymin, ymax))
    intf = montecarlo(f1, xmin, xmax, n, ymin, ymax)
    print("This leads to a second approximated integral value of {}.".format(intf))
    print("Error % is", 100 * abs(17.5 - intf))

    # plot_montecarlo(f, xmin, xmax, 1500, ymin, ymax, namef=r"$x \mapsto x$", figname="montecarlo.png")

    print("Formally, we compute the integral as 17.5 (36/2 - 1/2).")

    # %% Example 4 with Monte-Carlo
    xmin, xmax = 0.0, 1.0

    def f2(x):
        """ f1(x) = x"""
        return x ** 3

    # For this particular example, its easy:
    ymin, ymax = xmin, xmax

    print("\nFor this function f2: x → x^3, on I = [{}, {}], we will pick {} random points.".format(xmin, xmax, n))
    print("Manually, I chose ymin = {}, ymax = {}.".format(ymin, ymax))
    intf = montecarlo(f2, xmin, xmax, n, ymin, ymax)
    print("This leads to an approximated integral value of {}.".format(intf))
    print("Error % is", 100 * abs(0.25 - intf))

    ymin, ymax = yminmax(f2, xmin, xmax, n)  # guess!
    print("But experimentally, I found a possible ymin, ymax to be {}, {}.".format(ymin, ymax))
    intf = montecarlo(f2, xmin, xmax, n, ymin, ymax)
    print("This leads to a second approximated integral value of {}.".format(intf))
    print("Error % is", 100 * abs(0.25 - intf))

    # plot_montecarlo(f2, xmin, xmax, ymin, ymax, 1500, namef=r"$x \mapsto x^3$", figname="montecarlo2.png")
    print("Formally, we compute the integral as 0.25 (1/4).")

    # %% Example 5 with Monte-Carlo
    xmin, xmax = 1e-8, 3.

    def f3(x):
        """ f3(x) = x"""
        return 1.0 / (1 + (math.sinh(2 * x)) * (math.log(x)**2))

    # For this particular example, its easy:
    ymin, ymax = 0, 1

    print("\nFor this function f3: x → 1/(1+sinh(2x)log(x)^2) on I = [{}, {}], we will pick {} random points.".format(xmin, xmax, n))
    print("Manually, I chose ymin = {}, ymax = {}.".format(ymin, ymax))
    intf = montecarlo(f3, xmin, xmax, n, ymin, ymax)
    print("This leads to an approximated integral value of {}.".format(intf))

    ymin, ymax = yminmax(f3, xmin, xmax, n)  # guess!
    print("But experimentally, I found a possible ymin, ymax to be {}, {}.".format(ymin, ymax))
    intf = montecarlo(f3, xmin, xmax, n, ymin, ymax)
    print("This leads to a second approximated integral value of {}.".format(intf))

    # plot_montecarlo(f3, xmin, xmax, ymin, ymax, 1500, namef=r"$x \mapsto 1/(1+\sinh(2x)\log(x)^2)$", figname="montecarlo3.png")

    # plot_montecarlo(f3, xmin, xmax, ymin, ymax, 10000, namef=r"$x \mapsto 1/(1+\sinh(2x)\log(x)^2)$", figname="montecarlo4.png")


# End of tests.py