-
Notifications
You must be signed in to change notification settings - Fork 0
/
MC.py
97 lines (89 loc) · 2.69 KB
/
MC.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
import random
def validate(iters, typ=0):
ctrl = 0
try:
if typ == 0:
ctrl = int(iters)
else:
ctrl = float(iters)
return False
except ValueError:
print "No ingresaste bien el numero."
return True
def get_data(get_message,retry_message,typ=0):
data = raw_input(get_message)
while validate(data, typ):
data = raw_input(retry_message)
if typ == 0:
return int(data)
else:
return float(data)
def get_range():
a = get_data("Por favor ingresa el limite inferior de integracion: ",\
"Por favor reingresa el numero: ", 1)
b = get_data("Por favor ingresa el limite superior de integracion: ",\
"Por favor reingresa el numero: ", 1)
return [a,b]
def get_iters():
n = get_data("Por favor ingresa el numero de iteraciones: ",\
"Por favor reingresa el numero: ", 0)
return n
def find_max(f, rang):
lim = 10.**(-7)
pp = rang[0]
np = rang[1]
while abs(np - pp) > lim:
mid = (np - pp) / 2 + pp
if f(np) > f(pp) and f(mid) > f(pp):
if f(np) > f(mid):
pp = mid
elif f(np) < f(mid):
pp = np
np = mid
else:
return np
elif f(pp) > f(np) and f(mid) > f(np):
if f(pp) > f(mid):
np = mid
elif f(pp) < f(mid):
np = pp
pp = mid
else:
return pp
else:
if f(pp) == f(np) and f(mid) > f(pp):
left = find_max(f, [pp, mid])
right = find_max(f, [mid,np])
if left > right:
return left
elif left < right:
return right
elif left == right:
return left
else:
return mid
else:
return pp
return np
def monte_carlo(f, iters, rang, f_max, box_x):
inner_points = 0
for i in range(iters):
x = random.random() * box_x + rang[0]
y = random.random() * f_max
v = f(x)
if ( v > 0 and y < v ) or ( v < 0 and y > v ):
inner_points += 1
box_area = abs(box_x) * abs(f_max)
point_relation = inner_points * 1./iters
area = box_area * point_relation
return area
def mmc(f):
area = 0
iters = get_iters()
rang = get_range()
f_max = f(find_max(f, rang))
box_x = rang[1] - rang[0]
for i in range(100):
area += monte_carlo(f, iters, rang, f_max, box_x)
area /= 100
return area