Sumace 1/3: cykly, Numba a f2py, redukční funkce

Budeme počítat harmonická čísla:

\[H_n=\sum_{i=1}^n 1/i\,.\]

Zkusíme dvě varianty cyklů v Pythonu a porovnáme s obdobným řešením ve Fortranu. V Pythonu nás zklame rychlost běhu, způsobené kontinuálním překládáním (interpretováním) skriptů. Vložíme proto cykly do funkcí a použijeme pythonský balíček Numba pro zrychlení běhu funkcí jednorázovým překladem. Zkusíme také volat fortranskou funkci z Pythonu. Prověříme i další způsob, jak tlačit Python k rychlosti – počítat na NumPy polích spíše než s pythonskými seznamy. Pro kontrolu správnosti výpočtů můžeme použít vzorec podle Wikipedie:

\[H_n=\ln n+\gamma+\frac1{2n}-\varepsilon_n \,,\]

kde \(\varepsilon_n\) klesá rychleji než \(1/8n^2\). Eulerova konstanta \(\gamma \approx\) 0.577 215 664 901 533.

Cykly while a for

Dvě varianty pythonských cyklů:

# Harmonická čísla & Python: cyklus s podmínkou (cyklus while)
# https://en.wikipedia.org/wiki/Harmonic_series_(mathematics)

# n=int(input('Zadejte n: '))  # vstup dat z klávesnice
n=100_000_000
i=0; s=0.
while i<n:                     # cyklus s podmínkou
  i+=1; s+=1/i                 # blok cyklu: update proměnných i a s
print('n,s =',n,s)
# Harmonická čísla & Python: indexovaný cyklus (cyklus for)
# n=int(input('Zadejte n: '))
n=100_000_000
s=0.
for i in range(1,n+1):  # indexovaný cyklus od 1 do n (bez koncove hodnoty)
  s+=1/i                # blok cyklu: update proměnné s explicitně, update indexu i implicitně
print('n,s =',n,s)

Fortranská obdoba. S GNU Compiler Collection (download z equation.com) přeložíme a spustíme ve Windows i Linuxu pomocí příkazů:

gfortran -O -o file file.f90 && file
! Harmonická čísla & Fortran: cyklus s podmínkou
! https://en.wikipedia.org/wiki/Harmonic_series_(mathematics)

real(8) s
integer(8) n,i
! print '(a$)','Zadejte n: '; read *,n  ! vstup dat z klávesnice
n=100000000_8
i=0; s=0
do while (i<n)       ! cyklus s podminkou
  i=i+1; s=s+1._8/i  ! blok cyklu: update proměnných i a s
end do
print *,'n,s =',n,s,log(real(n,8))+.577215664901533_8+.5_8/n
end program
! Harmonická čísla & Fortran: indexovaný cyklus
real(8) s
integer(8) n,i
! print '(a$)','Zadejte n: '; read *,n
n=100000000_8
call cpu_time(tic)
s=0
do i=1,n  ! indexovany cyklus od 1 do n vcetne, krok implicitne 1
  s=s+1._8/i   ! blok cyklu: update promenne s explicitne, update indexu i implicitne
end do
call cpu_time(toc)
print *,'n,s =',n,s,log(real(n,8))+.577215664901533_8+.5_8/n,toc-tic
end program

Numba a zrychlování funkcí

Balíček Numba přináší možnost překladu pythonských skriptů do strojového kódu a přiblížení se rychlostem běhu programů v jiných překládaných jazycích, z nichž za chvíli vytěžíme Fortran. Překlad Numbou proběhne automaticky po spuštění skriptu (just-in-time compilation v nopython módu). Úspěch překladu je podmíněn splněním řady omezení. S našimi skripty to půjde snadno, přidáme řádek from numba import njit, cykly vložíme do funkcí a ty označíme pro překlad pomocí direktivy (metapříkazu, dekorátoru) @njit:

# Harmonická čísla & Python s funkcemi a Numba
from numba import njit
from time import time
n=1_000_000_000

@njit
def sum_while(n):
  i=0; s=0.
  while i<n:              # cyklus s podminkou
    i=i+1
    s=s+1/i
  return s

@njit
def sum_for(n):
  s=0.
  for i in range(1,n+1):  # indexovany cyklus
    s=s+1/i
  return s

tic=time(); print(n,sum_while(n)); toc=time(); print(toc-tic)
tic=time(); print(n,sum_for(n));   toc=time(); print(toc-tic)

f2py a fortranské funkce v Pythonu

Fortranský program s funkcí:

! Harmonická čísla & Fortran s funkcí
function sum_for(n) result (s)
  real(8) s
  s=0
  do i=1,n
    s=s+1._8/i
  end do
end function

real(8) sum_for
n=1000000000
call cpu_time(tic); print *,n,sum_for(n); call cpu_time(toc); print *,toc-tic
end program

Příkazem f2py, který je součástí balíčku NumPy, fortranskou funkci přeložíme (-c) a připravíme pythonský modul (-m):

f2py -c -m f90 file.f90

Import modulu a volání fortranské funkce z pythonského skriptu:

# Harmonická čísla & Python a f2py pro import fortranských funkcí
import f90
from time import time
n=1_000_000_000
tic=time(); print(n,f90.sum_for(n)); toc=time(); print(toc-tic)

Redukční funkce pro pole a seznamy

V Pythonu sčítáme v cyklu, sčítáme pole a sčítáme seznam, bez Numby i s Numbou:

# Harmonická čísla & Python a Numba: sčítání v cyklu, součet prvků pole, součet prvků seznamu
import numpy as np
from numba import njit
from time import time
n=1_000_000_000

# sčítání v cyklu for
@njit
def sum_for_loop(n):
  s=0.
  for i in range(1,n+1):
    s=s+1/i
  return s

# sčítání prvků pole redukční funkcí np.sum
@njit
def sum_for_array(n):
  a=np.empty(n)
  for i in range(1,n+1):
    a[i-1]=1/i
  s=np.sum(a)
  return s

# sčítání prvků seznamu redukční funkcí sum
@njit
def sum_for_list(n):
  a=[]
  for i in range(1,n+1):
    a.append(1/i)
  s=sum(a)
  return s

tic=time(); print(n,sum_for_loop(n));  toc=time(); print(toc-tic)
tic=time(); print(n,sum_for_array(n)); toc=time(); print(toc-tic)
tic=time(); print(n,sum_for_list(n));  toc=time(); print(toc-tic)

Ve Fortranu sečteme cyklem a v poli:

! Harmonická čísla & Fortran: sčítání v cyklu a součet prvků pole
function sum_for_loop(n) result (s)
  real(8) s
  s=0
  do i=1,n
    s=s+1._8/i
  end do
end function

function sum_for_array(n) result (s)
  real(8) s
  real(8) a(n)
  do i=1,n
    a(i)=1._8/i
  end do
  s=sum(a)
end function

real(8) sum_for_loop,sum_for_array
n=1000000000
call cpu_time(tic); print *,n,sum_for_loop(n), log(real(n))+.577215; call cpu_time(toc); print *,toc-tic
call cpu_time(tic); print *,n,sum_for_array(n),log(real(n))+.577215; call cpu_time(toc); print *,toc-tic
end program

Domácí úkol bez časového limitu: Nahraďte v některém z našich (nebo vašich) skriptů harmonickou řadu některou řadou pro výpočet \(\pi\). Inspirovat se můžete naším výběrem vzorců nebo vzorci z Wikipedie.