generators_ready

from ROOT import TH1F,TF1, gRandom, TStopwatch
from math import sin, cos, fabs
################ const ##########################################3
numb = 100000
x1 = -1
x2 = 1
coef = 1

############ normirovla functsii ###########################

def norm(a,b):
    return 1.0 / 10.0 * (2 * b ** 5 + 15 * b + 5 * sin(b) * cos(b) - 2 * a ** 5 - 15 * a - 5 * sin(a) * cos(a))

############### normirovka slagaemih #############

def norm3(a,b):
    return  1.0 * (b ** 5 - a ** 5) / 5.0

def norm2(a,b):
    return 1.0 * (b - a)

def norm1(a,b):
    return (1 / 2.0) * (b + sin(b) * cos(b)) - (1 / 2.0) * (a + sin(a) * cos(a))

################## beta ###################

def beta1(a, b):
    return 1.0 * norm1(a,b) / (norm1(a,b) + norm2(a,b) + norm3(a, b))

def beta2(a, b):
    return 1.0 * norm2(a,b) / (norm1(a,b) + norm2(a,b) + norm3(a, b))

def beta3(a, b):
    return 1.0 * norm3(a,b) / (norm1(a,b) + norm2(a,b) + norm3(a, b))

####################### functions of composition ###################

def func2(xm, p):
    x = xm[0]
    a = p[1]
    b = p[2]
    k = p[0]
    if norm2(a,b) <= 0:
        raise 'Wrong arguments in function: func1 norm1'
    return k * norm2(a,b) * (1)

def func3(xm, p):
    x = xm[0]
    a = p[1]
    b = p[2]
    k = p[0]
    if norm3(a,b) <= 0:
        raise 'Wrong arguments in function: func3 norm1'
    return k * norm3(a,b) * (x ** 4)

def func1(xm, p):
    x = xm[0]
    a = p[1]
    b = p[2]
    k = p[0]
    if norm1(a,b) <= 0:
        raise 'Wrong arguments in function: func1 norm1'
    return k * norm1(a,b) * ((cos(x)) ** 2)

######################## main function ###################################
def func(xm, p):
    #print a, b
    x = xm[0]
    a = p[1]
    b = p[2]
    k = p[0]
    if norm(a, b) <= 0:
    raise 'Wrong arguments in function: norm'
    return k * norm(a,b) * (x ** 4 + 1 + (cos(x)) ** 2)
################################### methods ##################################
def neumann(a, b, c, d):
    if d == 0:
    if b <= a:
            raise 'Wrong number in neumann'
        fmax = func([ b ], [c, a, b])
        while True:
        mu = gRandom.Uniform(0, fmax)
        r = gRandom.Uniform(a,b)
        if mu <= func( [ r ], [c, a, b]):
        return r
        raise 'Could never reach this point'
    else:
    if b <= a:
            raise 'Wrong number in neumann'
        fmax = 1 * c * norm1(a,b)
        while True:
        mu = gRandom.Uniform(0, fmax)
        r = gRandom.Uniform(a,b)
        if mu <= func1( [ r ], [c, a, b]):
        return r
        raise 'Could never reach this point'
###############do the const#######################
b1 = beta1(x1,x2)
b2 = beta2(x1,x2)
b3 = beta3(x1,x2)
#######################################################
def regression(a,b,c):
    r = gRandom.Uniform(0,1)
    if r < b1:
    return neumann(a,b,c,1)
    if (r < (b1 + b2) and r >= b1):
    return fndirect(a,b,c,2)
    if (r >= (b1 + b2) and r <= 1):
    return fndirect(a,b,c, 3)

def fndirect(a, b, c, d):
    if d == 3:
    ua = 1.0 * a ** 5 / 5
    ub = 1.0 * b ** 5 / 5
    u = gRandom.Uniform(ua,ub)
    if u < 0:
        t = fabs(u)
        r = -(5.0 * t) ** (1.0 / 5)
    else:
        r = (5.0 * u) ** (1.0 / 5)
        return r
    if d == 2:
    u = gRandom.Uniform(a,b)
    r = u
    return r
####################################### main #############################
f = TF1("f", func, x1, x2, 3)
f.SetParameter(0,x1)
f.SetParameter(1,x2)
f.SetParameter(2,coef)
f.Draw()
raw_input()
h1 = TH1F("h1", "neumann method",100,x1,x2 )
h2 = TH1F("h2", "regression", 100, x1, x2)
#zapolnenie neumann
sw = TStopwatch()
sw.Start()
for i in range(0, numb):
    r = neumann(x1, x2, coef, 0)
    h1.Fill(r)
sw.Stop()
sw.Print()
#print norm1(x1,x2), norm2(x1,x2), norm3(x1,x2), norm1(x1,x2) + norm2(x1,x2) + norm3(x1,x2)
h1.Draw()
print "Mean is ", h1.GetMean()
print "RMS is " , h1.GetRMS()
raw_input()
f.FixParameter(1,x1)
f.FixParameter(2,x2)
h1.Fit(f)
raw_input()
sw1 = TStopwatch()
sw1.Start()
for i in range(0, numb):
    r =regression(x1,x2,coef)
#    r = fndirect(x1,x2,coef,3)
    h2.Fill(r)
sw1.Stop()
sw1.Print()
h2.Draw()
print "Mean2 is ", h2.GetMean()
print "RMS2 is " , h2.GetRMS()
raw_input()
h2.Fit(f)
raw_input()