# Cmput 455 sample code
# Numerical Integration using Monte Carlo method
# FB - 201006137
# Code by author "FB36", file recipe-577263-1.py from  http://code.activestate.com/recipes/577263-numerical-integration-using-monte-carlo-method/
# Adapted to Python3 and slightly modified for Cmput 455 by Martin Mueller

import math
import random

def f(x):
    return x*x - 1

# find ymin-ymax range
def findYRange(function, xmin, xmax):
    numSteps = 1000000 # bigger the better but slower!
    ymin = function(xmin)
    ymax = ymin
    for i in range(numSteps):
        x = xmin + (xmax - xmin) * float(i) / numSteps
        y = function(x)
        if y < ymin: ymin = y
        if y > ymax: ymax = y
    return ymin, ymax

# Monte Carlo
numPoints = 1000 # bigger the better but slower!

def monteCarloIntegrate(function, xmin, xmax, ymin, ymax):
    sampleSum = 0
    for _ in range(numPoints):
        x = xmin + (xmax - xmin) * random.random()
        y = ymin + (ymax - ymin) * random.random()
        if math.fabs(y) <= math.fabs(function(x)):
            if function(x) > 0 and y > 0 and y <= function(x):
                sampleSum += 1 # area over x-axis is positive
            if function(x) < 0 and y < 0 and y >= function(x):
                sampleSum -= 1 # area under x-axis is negative
    return sampleSum

def numericalIntegrate(function, xmin, xmax):
    ymin, ymax = findYRange(function, xmin, xmax)
    count = monteCarloIntegrate(function, xmin, xmax, ymin, ymax)
    rectArea = (xmax - xmin) * (ymax - ymin)
    fnArea = rectArea * float(count) / numPoints
    print("Numerical integral =", fnArea)

#numericalIntegrate(math.sin, 0.0, math.pi) # True value = 2
numericalIntegrate(f, 0.0, 10.0) # True value = 323.33333...