igor's_montecarlo

from ROOT import TF2, gRandom, gROOT, TStopwatch
from math import exp, log, sqrt
import sys

def func(xx):
    x = xx[0]
    y = xx[1]
    return exp(x ** 2 + y) / (1 + x * sqrt(y) / 2)

def func_g(xx, pp):
    x = xx[0]
    y = xx[1]
    norm = pp[0]
    return exp(y) / 2 * norm

def func_fg(xx):
    x = xx[0]
    y = xx[1]
    return func ([x, y]) / func_g([x, y], [norm])

def fndirect(x1, x2, y1, y2):
    if x1 < 0 or x2 > 0.5 or y2 > 3 or y1 < 1:
        raise 'Wrong arguments in fndirect'

    u = gRandom.Uniform(umin, umax);

    ry = log(u)
    rx = gRandom.Uniform(x1, x2)
    return rx, ry



def direct_integrate(fun, left, right, number):
    xmin = left[0] ; xmax = right[0]
    ymin = left[1] ; ymax = right[1]
    sum_f = 0
    sum_fSquared = 0
    for n in range(0,number):
        x = gRandom.Uniform(xmin, xmax)
    y = gRandom.Uniform(ymin, ymax)
        f = fun([x, y])
        sum_f += f
        sum_fSquared += f ** 2
    mean = sum_f * (xmax - xmin) * (ymax - ymin) / number
    disp = sum_fSquared * ((xmax - xmin) * (ymax - ymin)) ** 2 / number
    return mean, disp

def signSampl_integrate(fun_f, fun_g, left, right, number):
    xmin = left[0] ; xmax = right[0]
    ymin = left[1] ; ymax = right[1]
    s = 0
    sum_fSquared = 0
    for n in range(0,number):
        x, y = fndirect(xmin, xmax, ymin, ymax)
        f = fun_f([x, y]) / fun_g([x, y], [norm])
        s += f
        sum_fSquared += f ** 2
    mean = s * (x2 - x1) / number
    disp = sum_fSquared * (x2 - x1) ** 2 / number
    return mean, disp

x1 = 0
x2 = 0.5
y1 = 1
y2 = 3

norm = 2 / (exp(3) - exp(1))
f = TF2("f", func, x1, x2, y1, y2, 0)
f.Draw("surf")
raw_input()

fg = TF2("fg", func_fg, x1, x2, y1, y2, 0)
fg.Draw("surf")
raw_input()

step = 1000
totalMean = 0
totalDisp = 0
print 'mean method\n'
for i in range(1, 11):
    result = direct_integrate(func, [x1, y1], [x2, y2], step)
    totalMean += result[0]
    totalDisp += result[1]
    print str(i * step) +  '\tmean: ' +  str(result[0]) + '\trms: ' + str(sqrt((result[1] - result[0] ** 2) / step)) + \
       '\ntotal mean: ' +   str(totalMean / i) + '\ttotal rms: ' + str(sqrt((totalDisp / i - (totalMean / i) ** 2) / (step * i)))
totalMean = 0
totalDisp = 0
umin = exp(1)
umax = exp(3)
print '\nsignificant sampling\n'
for i in range(1, 11):
    result = signSampl_integrate(func, func_g, [x1, y1], [x2, y2], step)
    totalMean += result[0]
    totalDisp += result[1]
    print str(i * step) +  '\tmean: ' +  str(result[0]) + '\trms: ' + str(sqrt((result[1] - result[0] ** 2) / step)) + \
       '\ntotal mean: ' +   str(totalMean / i) + '\ttotal rms: ' + str(sqrt((totalDisp / i - (totalMean / i) ** 2) / (step * i)))

raw_input()