2015-06-12 22 views
5

Napisałem poniższy kod w Pythonie, aby oszacować wartość Pi. Nazywa się to metodą Monte Carlo. Oczywiście przez zwiększenie liczby próbek kod staje się wolniejszy i zakładam, że najwolniejsza część kodu znajduje się w części próbkowania. Jak mogę przyspieszyć?Jak zwiększyć wydajność do szacowania `Pi`in Python

from __future__ import division 
import numpy as np 

a = 1 
n = 1000000 

s1 = np.random.uniform(0,a,n) 
s2 = np.random.uniform(0,a,n) 

ii=0 
jj=0 

for item in range(n): 
    if ((s1[item])**2 + (s2[item])**2) < 1: 
     ii = ii + 1 

print float(ii*4/(n)) 

Czy sugerujesz inne (prawdopodobnie szybsze) kody?

+2

Jakiś powód, dla którego musiałbyś to obliczyć? Dokonano tego już tysiąc razy (np. Przeczytaj z pliku lub zakoduj go). – Caramiriel

+3

Powolność jest nieodłącznym elementem algorytmu - za każdym razem, gdy potrzebna jest jeszcze jedna cyfra precyzji, czas działania wzrasta o 10x. Sprawdź [Pi - nowoczesne poszukiwanie więcej cyfr] (https://en.wikipedia.org/wiki/Pi#Modern_quest_for_more_digits) w celu szybszego podejścia. – Kevin

+0

Czy jest jakikolwiek powód do zbudowania dwóch tablic z góry? Dlaczego nie dostać dwóch liczb 'random.random()' wewnątrz * pętli? Jednak, jak zaznacza @Kevin, ten algorytm jest w zasadzie "O (n)", więc dla dużych 'n' jakiekolwiek zmiany w precyzyjnej implementacji będą jedynie minimalnymi różnicami w stosunku do ogólnego środowiska wykonawczego. – jonrsharpe

Odpowiedz

8

Wąskim gardłem jest tutaj twoja pętla for. Pętle Python for są stosunkowo powolne, więc jeśli potrzebujesz powtórzyć ponad milion pozycji, możesz zyskać dużą szybkość, unikając ich w ogóle. W tym przypadku jest to dość łatwe. Zamiast tego:

for item in range(n): 
    if ((s1[item])**2 + (s2[item])**2) < 1: 
     ii = ii + 1 

to zrobić:

ii = ((s1 ** 2 + s2 ** 2) < 1).sum() 

To działa, ponieważ numpy posiada wbudowane wsparcie dla optymalizacji operacji tablicowych. Pętla występuje w języku c zamiast pythona, więc jest znacznie szybsza. Zrobiłem szybki test, dzięki czemu można zobaczyć różnicę:

>>> def estimate_pi_loop(x, y): 
...  total = 0 
...  for i in xrange(len(x)): 
...   if x[i] ** 2 + y[i] ** 2 < 1: 
...    total += 1 
...  return total * 4.0/len(x) 
... 
>>> def estimate_pi_numpy(x, y): 
...  return ((x ** 2 + y ** 2) < 1).sum() 
... 
>>> %timeit estimate_pi_loop(x, y) 
1 loops, best of 3: 3.33 s per loop 
>>> %timeit estimate_pi_numpy(x, y) 
100 loops, best of 3: 10.4 ms per loop 

Oto kilka przykładów rodzajów działalności, które są możliwe, tak więc masz poczucie jak to działa.

Kwadratura tablicy:

>>> a = numpy.arange(5) 
>>> a ** 2 
array([ 0, 1, 4, 9, 16]) 

Dodawanie tablic:

>>> a + a 
array([0, 2, 4, 6, 8]) 

Porównanie tablic:

>>> a > 2 
array([False, False, False, True, True], dtype=bool) 

Sumowanie wartości logiczne:

>>> (a > 2).sum() 
2 

Jak zapewne zdajesz sobie sprawę, istnieją szybsze sposoby oszacowania Pi, ale przyznaję, że zawsze podziwiałem prostotę i skuteczność tej metody.

2

Masz przydzielone niezliczone tablice, więc powinieneś używać ich na swoją korzyść.

for item in range(n): 
    if ((s1[item])**2 + (s2[item])**2) < 1: 
     ii = ii + 1 

staje

s1sqr = s1*s1 
s2sqr = s2*s2 
s_sum = s1sqr + s2sqr 
s_sum_bool = s_sum < 1 
ii = s_sum_bool.sum() 
print float(ii*4/(n)) 

Gdzie jesteś kwadratury tablic, ich zsumowanie, sprawdzenie, czy suma ta jest mniejsza niż 1, a następnie zsumowanie wartości logicznych (false = 0, prawda = 1), aby uzyskać całkowita liczba spełniająca kryteria.

1

I upvoted odpowiedź senderle „s, ale w przypadku, gdy nie chcesz, aby zmienić swój kod zbyt wiele:

numba jest biblioteką przeznaczone do tego celu.

Wystarczy zdefiniować algorytm jako funkcja i dodać @jit dekorator:

from __future__ import division 
import numpy as np 
from numba import jit 

a = 1 
n = 1000000 

s1 = np.random.uniform(0,a,n) 
s2 = np.random.uniform(0,a,n) 

@jit 
def estimate_pi(s1, s2): 
    ii = 0 
    for item in range(n): 
     if ((s1[item])**2 + (s2[item])**2) < 1: 
      ii = ii + 1 
    return float(ii*4/(n)) 

print estimate_pi(s1, s2) 

na moim laptopie, zyskuje około 20x dla n = 100000000 SpeedUp oraz SpeedUp 3x dla n = 1000000.