#! /usr/bin/env python
# -*- coding: utf-8 -*-
"""
Simple Monte-Carlo method for a 1D integral.

Note: thanks to Anu Apurvaa Sindol and Devavrat Yargop (A1 group) for their initial program.


@date: Sat Mar 18 12:21:29 2015.
@author: Lilian Besson for CS101 course at Mahindra Ecole Centrale 2015.
@licence: MIT Licence (http://lbesson.mit-license.org).
"""

# WARNING: To pick real number in an interval [a, b], we must use uniform(a, b), not randrange(a, b) !
from random import uniform


def montecarlo(f, xmin, xmax, ymin, ymax, n):
    """ Compute an approximation of the integral of f(x) for x from xmin to xmax.

    - n is the number of random points to pick (should be big, like 1 lack at least).
    - WARNING: ymin should be smaller than the minimum of f(x), for x in [xmin, xmax], and ymax should be bigger that the maximum of f(x).
    """
    # Here we are cautious about the arguments
    if n <= 0 or not isinstance(n, int):
        raise ValueError("montecarlo() the argument n has to be an integer.")
    if not (xmax >= xmin and ymax >= ymin):
        raise ValueError("montecarlo() invalid arguments xmax < xmin or ymax < ymin.")

    # This will count the number of points below the curve
    nb_below = 0

    area_rectangle = (xmax - xmin) * (ymax - ymin)

    for _ in xrange(0, n):
        x = uniform(xmin, xmax)
        if not(ymin <= f(x) <= ymax):
            raise ValueError("montecarlo() ymin and ymax are not correct: for x = {}, f(x) = {} is NOT in the interval [ymin, ymax] = [{}, {}].".format(x, f(x), ymin, ymax))
        y = uniform(ymin, ymax)
        # Here we need to be cautious, if f(x) is positive or not !
        if 0 <= y <= f(x):
            nb_below += 1
        elif f(x) <= y <= 0:
            nb_below -= 1

    observed_probability = float(nb_below)/float(n)

    return area_rectangle * observed_probability
# End of the function montecarlo


def yminmax(f, xmin, xmax, n, threshold = 0.02):
    """ Experimental guess of the values for ymin, ymax, by randomly picking 2*n points.

    - Note: there are more efficient and trustworthy methods.
    - threshold is here to increase a little bit the size of the window, to be cautious. Default is 5%.
    """
    ymin, ymax = f(xmin), f(xmax)  # initial guess!
    for _ in xrange(0, 2*n):
        x = uniform(xmin, xmax)
        if f(x) < ymin:
            ymin = f(x)
        if ymax < f(x):
            ymax = f(x)
    # Now we just increase a little bit the size of the window [ymin, ymax]
    if ymin < 0:
        ymin *= (1+threshold)  # ymin becomes bigger!
    else:
        ymin *= (1-threshold)  # ymin becomes smaller!
    if ymax < 0:
        ymax *= (1-threshold)  # ymax becomes smaller!
    else:
        ymax *= (1+threshold)  # ymax becomes bigger!

    return ymin, ymax


# %% Here, we perform some examples
if __name__ == '__main__':
    n = 100000

    # Example 1
    xmin, xmax = 1, 6

    def f(x):
        return x
    # f = lambda x: x  # quicker!

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

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

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

    print "Formally, we compute the integral as 12.5 (25/2)."


    # Example 2
    xmin, xmax = 0., 1.

    f = lambda x: x**3  # quicker!

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

    print "\nFor this function f: 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(f, xmin, xmax, ymin, ymax, n)
    print "This leads to an approximated integral value of {}.".format(intf)

    ymin, ymax = yminmax(f, xmin, xmax, n)  # guess!
    print "But experimentally, I found a possible ymin, ymax to be {}, {}.".format(ymin, ymax)
    intf = montecarlo(f, xmin, xmax, ymin, ymax, n)
    print "This leads to a second approximated integral value of {}.".format(intf)
    print "Formally, we compute the integral as 0.25 (1/4)."


    # Example 2
    xmin, xmax = 0.8, 3.

    import math
    f = lambda x: 1./(1+(math.sinh(2*x))*(math.log(x)**2))  # quicker!

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

    print "\nFor this function f: 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(f, xmin, xmax, ymin, ymax, n)
    print "This leads to an approximated integral value of {}.".format(intf)

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