# Интегрирование заданной функции F(x) на отрезке [A; B]
# Методы:
# - левых прямоугольников
# - средних прямоугольников
# - осреднения (Монте-Карло)
# - геометрический (Монте-Карло)
#

#===================
# Интегрируемая функция F(x)
#===================
#
def F(x):
    return x**2-x

# Подключение библиотеки random для генерации псевдослучайных чисел.
# (Псевдо-)случайные числа используются в методах Монте-Карло.
# В данной программе используюется только функция uniform(),
# которая генерирует случайное число с равномерным распределением
# на заданном отрезке.
#
from random import uniform

#===================
#   Integrate_MC_averaging(F, A,B, N)
# Интегрирование методом осреднения (Монте-Карло)
# ПАРАМЕТРЫ:
#   F   - интегрируемая функция
#   A,B - начало и конец отрезка интегрирования
#   N   - количество выбираемых случайных точек
# РЕЗУЛЬТАТ:
#   - значение интеграла
#===================
#
def Integrate_MC_averaging(F, A,B, N):
    sum_Fi = sum( F(uniform(A,B)) for i in range(N) )
    average = sum_Fi/N          # математическое ожидание для F(случайный-X)
    return (B-A) *average       # вычисление значения интеграла

#===================
#   Sign(x)
# Возвращает знак вещественного числа: 
#   +1 для x>0
#    0 для x==0
#   -1 для x<0
# Нужна исключительно для упрощения кода геометрического интегрирования.
#===================
#
def Sign(x):                           # функция "знак числа"
    return (x>0) - (x<0)

#===================
#   Integrate_MC_geometric(F, A,B, N)
# Интегрирование геометрическим методом (Монте-Карло)
# ПАРАМЕТРЫ:
#   F   - интегрируемая функция
#   A,B - начало и конец отрезка интегрирования
#   Ymin,Ymax - границы по оси OY, требуется Ymin <= min(F(x) <= max(F(x) <= Ymax
#   N   - количество выбираемых случайных точек
# РЕЗУЛЬТАТ:
#   - значение интеграла
#===================
#
def Integrate_MC_geometric(Func, A,B, Ymin,Ymax, N):
    K = 0                              # счётчик попаданий
    for i in range(N):
        x = uniform(A,B)               # выбрали случайный X
        y = uniform(Ymin,Ymax)         # выбрали случайный Y
        Fx = Func(x)                   # вычислили значение функции
        if y*(Fx-y)>=0: K+= Sign(y)    # проверка точки на попадание между графиком и осью OX
    average = K/N                      # частота попаданий
    return (B-A)*(Ymax-Ymin)*average   # вычисление значения интеграла

#===================
#   Integrate_LeftRect(F, A,B, N)
# Интегрирование методом "левых" прямоугольников.
# ПАРАМЕТРЫ:
#   F   - интегрируемая функция
#   A,B - начало и конец отрезка интегрирования
#   N   - количество выбираемых случайных точек
# РЕЗУЛЬТАТ:
#   - значение интеграла
#===================
#
def Integrate_LeftRect(F, A,B, N):
    h = (B-A)/N                     # шаг интегрирования
    sum_Fi = sum( F( h*i ) for i in range(N) ) # формула суммирования по методу левых пр-ков
    return sum_Fi*h                 # вычисление значения интеграла
    
#===================
#   Integrate_MidRect(F, A,B, N)
# Интегрирование методом "средних" прямоугольников.
# ПАРАМЕТРЫ:
#   F   - интегрируемая функция
#   A,B - начало и конец отрезка интегрирования
#   N   - количество выбираемых случайных точек
# РЕЗУЛЬТАТ:
#   - значение интеграла
#===================
#
def Integrate_MidRect(F, A,B, N):
    h = (B-A)/N                     # шаг интегрирования
    sum_Fi = sum( F( h*i+h/2 ) for i in range(N) ) # формула суммирования по методу средних пр-ков
    return sum_Fi*h                 # вычисление значения интеграла



A = 0               # начало отрезка интегрирования
B = 3               # конец отрезка интегрирования

print("Интегрирование функции {} на отрезке [{}; {}] различными способами".format("F(x)=x**2-x", A,B))

# Выполняется вычисление интеграла функции F(x)=x**2-x на отрезке [0; 3]
# несколькими методами, с разным количеством точек/узлов:  N = 10, 1000, 100000.
# При интегрировании по методу Монте-Карло выполняется 
#
for N in [100, 10000, 1000000]:
    print("\n N=", N)       # N = количество узлов / случайных точек
    print("метод  левых  пр-ков:  ", Integrate_LeftRect(F, A,B, N))
    print("метод средних пр-ков:  ", Integrate_MidRect (F, A,B, N))

    m = 10                  # количество 
    MC_avg = [Integrate_MC_averaging(F, A,B, N//m) for i in range(m)]
    print("Монте-Карло осреднения:", sum(MC_avg)/len(MC_avg) )

    MC_geom = [Integrate_MC_geometric(F, A,B, -1,9, N//m) for i in range(m)]
    print("Монте-Карло геометрич.:", sum(MC_geom)/len(MC_geom) )


#============================================
# Дополнительно: способы вычисления доверительного интервала
#
# Computing the confidence interval
#   1)
##  from functools import reduce
##  variance1 = reduce(lambda s,x: s+(x-average)**2 , Fi, 0)/(len(Fi)-1)
##  2)
##  tmp = [(x-average)**2 for x in Fi]
##  variance2 = sum(tmp)/(len(Fi)-1)
##  3)
##  variance3 = sum(map(lambda x: (x-average)**2, Fi))/(len(Fi)-1)
#
# Вычисление верхней и нижней оценки по уровню доверия 95%
#
##  variance = (sum(map(lambda x: (x-average)**2, Fi))/(N-1))
##  deviance = variance**0.5
##
##  lower_eval = (B-A)*(average-(deviance/N**0.5)*1.96)
##  upper_eval = (B-A)*(average+(deviance/N**0.5)*1.96)
##
##  print("N =", N)
##  print("random integral =", integral_eval)
##  print("bounds: {} to {}".format(lower_eval, upper_eval))