Programování pro fyziky 2023/24 – přednášky č. 7 a 8

Dnešní téma nás vyzbrojí k procedurálnímu programování. Chápejme jej jako vítanou možnost, jak pomocí funkcí a procedur rozšířit možnosti imperativního strukturovaného programování. Funkcionální programování jde ještě dále, možnost tlačí spíš až k nutnosti, snaží se nahradit i samotný imperativní strukturovaný styl.

Funkce a procedury

Členění zdrojového kódu do menších úseků (jednotek, program units) má řadu výhod: oddělenost jednotek přináší lepší srozumitelnost a testovatelnost; zřetelné hranice jednotek (lokalita, scoping) umožňují přesné definování vstupních a výstupních dat a zapouzdření lokálních dat; dobře navržené jednotky lze používat opakovaně, mohou se sdružovat do knihoven (libraries). Funkce (functions) se podobně jako v matematice používají pro získání jedné návratové hodnoty, procedury (procedures, subprograms, subroutines) jsou určeny pro ostatní účely; některé jazyky procedury nemají a simulují je pomocí funkcí. Voláním jména funkce ve výrazu, resp. příkazem volání procedury se provede odskok na začátek funkce či procedury, po jejím vykonání se chod programu vrací zpět na místo volání. S volající programovou jednotkou komunikují funkce/procedury jednak prostřednictvím svého rozhraní, tj. svých argumentů (parametrů), jednak prostřednictvím dat známých oběma jednotkám (globální data). Data deklarovaná ve funkci/proceduře jsou zapouzdřená, vně nedostupná (lokální data).

Python deklaruje funkce příkazem def. Tělo funkce je vyznačeno, jak je u pythonských konstrukcí zvykem, odsazením. Návratovou hodnotu funkce definuje příkaz return, po jehož dosažení funkce končí. Nenarazí-li při průchodu funkce na příkaz return, funkce skončí po provedení posledního příkazu těla funkce s prázdnou návratovou hodnotou None. (Lze psát i explicitně return None.) Funkci lze ukončit i po doporučeníhodných testech chybových stavů příkazem assert, nebo kdykoliv – spolu s celým skriptem – voláním funkce exit(). Je dobrým zvykem poznamenat k záhlaví funkce komentář o jejím účelu, v Pythonu obvykle formou dokumentačního řetězce (docstring) v trojitých apostrofech nebo uvozovkách – ten je pak (pro funkce f) dostupný jako atribut f.__doc__ nebo voláním funkce help(f). V Pythonu se žije bez procedur, ale snadno, procedury se simulují funkcemi bez návratové hodnoty čili s return None.

# Python: funkce s jedním argumentem, docstring, assert, exit(), return
def f(x):
  '''funkce f(x): test assert, exit() a return'''
  assert x>=0, 'Funkce f nechce x<0.'
  if x==0: exit('exit z f(x): x nulové')
  return x

# Funkce bez návratové hodnoty (ekvivalent procedury)
def p(x,y):
  '''funkce p(x) bez return'''
  print('p(x,y):',x+y)

print(f.__doc__)         # výpis docstring
# help(f)                # alternativa
print('f(-1) = ',f(-1))  # akce assert
print('f(0) = ',f(0))    # akce exit()
print('f(1) =',f(1))     # volání funkce a výpis návratové hodnoty
p(1,1)                   # volání bez užití návratové hodnoty

Funkce mohou být řazeny sekvenčně. Mohou mít rozličný počet argumentů (včetně žádného argumentu nebo jejich volitelného počtu). Návratová hodnota může být jediná – v Pythonu to ovšem může být n-tice (tuple), dokonce bez závorek, tak to pak ani jako jedna návratová hodnota nemusí vypadat.

# Python: funkce s různým počtem argumentů a návratových hodnot, volitelné argumenty
def f0():
  return 0

def f1(x):
  return x

def f2(x,y):
  return x+y

# funkce s volitelnými argumenty (s předepsanou default hodnotou)
def f3(x=0,y=0,z=0):
  return x+y+z

def f4(x,y,z):
  return x,y,z                              # totéž co return (x,y,z)

print('f0:',f0())
x=1
print('f1:',f1)                             # volání bez závorek: f1: <function f1 at ...>
# print('f1:',f1())                         # TypeError: missing argument
print('f1:',f1(0),f1(x),f1(x+1),f1(f1(3)))  # f1: 0 1 2 3
print('f2:',f2(1,1),f2(f2(1,1),f2(1,1)))    # f2: 2 4
print('f3:',f3(),f3(1),f3(1,1),f3(1,1,1))   # f3: 0 1 2 3
print('f4:',f4(1,1,1),type(f4(1,1,1)))      # f4: (1, 1, 1) <class 'tuple'>

V objektově orientovaném programování (OOP) se mohou funkce a procedury (zde: metody) vázat na strukturovaný typ (třídu) a potažmo na proměnnou tohoto typu (objekt). Objekt je typicky argumentem metod, což inspirovalo alternativní syntaxi volání: místo method(object,other_arguments) lze volat object.method(other_arguments). V Pythonu je často v nabídce obojí volání, ale někdy jen procedurální, jindy jen objektové, anebo je omezena jejich použitelnost. Např. dotaz na velikost NumPy polí size má obě podoby, ale spolu s balíčkem Numba pro zrychlení funkcí jejich jednorázovým překladem lze použít jen objektový zápis A.size:

# Python: OOP syntaxe
import numpy as np
from numba import njit
def test1():
  A=np.zeros((2,2))
  print(np.size(A),A.size,np.shape(A),A.shape)  # výpis: 4 4 (2, 2) (2, 2)
  B=np.ones(A.shape)
  print(sum(sum(B)),B.sum())                    # výpis: 4.0 4.0
@njit
def test2():
  A=np.zeros((2,2))
  print(A.size,np.shape(A),A.shape)  # njit nekompatibilní s np.size(A)
  B=np.ones(A.shape)
  print(B.sum())                     # njit nekompatibilní se sum(B)
test1()
test2()

Předávání dat

Úkolem je zajistit jednak jednosměrné předávání vstupních argumentů (z místa volání do funkce/procedury) a výstupních argumentů (z funkce/procedury ven), jednak obousměrné předávání argumentů, které tak jsou vstupně-výstupní. Snahou též bývá zajistit efektivní předávání rozměrnějších dat, tj. polí a struktur. V klasických programovacích jazycích se to děje tzv. předáváním hodnotou nebo odkazem (call by value/by reference). Předání hodnotou znamená vytvoření nového paměťového místa, do kterého je zkopírována hodnota argumentu z volající jednotky (tím může být proměnná nebo hodnota vypočteného výrazu); při návratu se hodnota zpět nekopíruje. Předání odkazem předpokládá, že předávaný argument má své paměťové místo (musí být adresovatelný). Do volané jednotky se předá adresa argumentu; případné změny argumentu prováděné uvnitř procedury se pak realizují vně, v paměťovém místě předaného argumentu. Pokud se vstupní argumenty v proceduře nemění, je vhodnější je výslovně označit jako konstantní. Zvláště vhodné je předávat odkazem paměťově náročné argumenty, např. velká pole.

Lokální proměnné se alokují v paměti až po vstupu do funkce a zanikají při odchodu z ní. Pythonská funkce vidí do vnějších dat existujících při vstupu do funkce; pokud se do nich pokusí zapsat, buď se (jako default) vytvoří nová lokální proměnná a překryje vnější data, nebo (při uvedení proměnné v příkazu global) se vnější proměnná updatuje – vnější proměnná je pak z pohledu funkce globální proměnnou. V Pythonu můžeme předávání argumentů shrnout tak, že všechny argumenty jsou vstupní a výstup je realizován návratovou hodnotou, která může být skalárem i kolekcí dat různého typu. Pokud se ocitne argument funkce na levé straně přiřazovacího příkazu, vytvoří se nová lokální proměnná a překryje argument.

# Python: argumenty a lokální a globální proměnné
def f(a,b):       # argumenty a,b
  global e
  print('f: a,b,c,e =',a,b,c,e) # f: a,b,c,e = 1 1 1 1 (lokální d zatím vypsat nelze)
  b=2             # update b: nová lokální proměnná překrývající argument
  d=4             # update d: nová lokální proměnná překrývající vnější proměnnou
  e=5             # update globální proměnné
  return a,b,c,d,e

a=1; b=1; c=1; d=1; e=1
print('a,b,c,d,e =',a,b,c,d,e)  # a,b,c,d,e = 1 1 1 1 1
print('f =',f(a,b))             # f = (1, 2, 1, 4, 5)
print('a,b,c,d,e =',a,b,c,d,e)  # a,b,c,d,e = 1 1 1 1 5

V Pythonu není předurčena cesta pro definování globálních konstant nebo statických lokálních proměnných funkcí. S jistou dávkou programátorovy kázně je lze simulovat.

# Python: globální konstanta, statická lokální proměnná funkce
nmax=10    # globální konstanta
counter=0  # statická lokální proměnná funkce

def allocate(string='A',n=nmax):
  global counter
  A=string*n
  counter+=1
  return A

print(allocate('A'))      # AAAAAAAAAA
print(allocate('ABC',3))  # ABCABCABC
print(allocate('XYZ',1))  # XYZ
print('citac:',counter)   # citac: 3

Tvrzení, že lokální update argumentů se nepropíše vně funkce, platí jen pro neměnné (immutable) argumenty, což jsou především skaláry standardních typů a n-tice (tuples). Lze však updatovat obsah argumentů modifikovatelných (mutable) datových typů, např. seznamů, slovníků nebo NumPy polí.

# Python: modifikovatelnost argumentu funkce
def f(arg):
  arg+=arg
  return arg

import numpy as np
arg=1;    print(f(arg),arg)           # 2 1
arg=(1,); print(f(arg),arg)           # (1, 1) (1,)
arg=[1];  print(f(arg),arg)           # [1, 1] [1, 1]
arg=np.array([1]); print(f(arg),arg)  # [2] [2]

Ještě příklad s updatem pole v argumentu pomocí přiřazení do sekce a v cyklech:

# Python: update argumentů typu pole ve funkcích
import numpy as np

def f(a,b,c,d):
  a=[1,2,3]          # update celého argumentu: nový seznam
  b[:]=[1,2,3]       # update sekce: přenos hodnot do prvků pole
  for n in range(len(c)): c[n]=n+1  # update prvků pole indexovaným cyklem: přenos hodnot do prvků pole
  for x in d: x=x+1  # update prvků pole cyklem "foreach" neúčinný
  return a,b,c,d

a=np.zeros(3)
b=np.ones(3)
c=np.full(3,2.)
d=np.full(3,3.)
print('a,b,c,d =',a,b,c,d)  # a,b,c,d = [0. 0. 0.] [1. 1. 1.] [2. 2. 2.] [3. 3. 3.]
print('f:',f(a,b,c,d))      # f: ([1, 2, 3], array([1., 2., 3.]), array([1., 2., 3.]), array([3., 3., 3.]))
print('a,b,c,d =',a,b,c,d)  # a,b,c,d = [0. 0. 0.] [1. 2. 3.] [1. 2. 3.] [3. 3. 3.]

Předávat lze i odkazy na funkce, funkce může být argumentem funkce. Argumenty jsou volitelné, mají-li v hlavičce funkce přiřazenu default hodnotu. Operátor * rozbaluje posloupnosti (n-tice, seznamy, NumPy pole) na jednotlivé prvky, např. uvedením u skutečného argumentu (ve volání funkce). Totéž u formálního argumentu (v hlavičce funkce) vyznačuje možnost proměnného počtu argumentů.

# Python: funkce v argumentu funkce, volitelné argumenty, proměnný počet argumentů
def f1(x): return x
def f2(x): return x*2

def g(f,x,d=0):                   # první argument funkce, třetí argument volitelný
  result=f(x+d)
  return result

def h(init,*args):                # proměnný počet argumentů
  result=init
  for arg in args: result+=arg
  return result
                                  # výpisy: 1, 4, 6 10, 6 10
print(g(f1,1))                    # volání bez třetího volitelného argumentu
print(g(f=f2,d=1,x=1))            # argumenty přiřazeny pomocí formálních jmen
print(h(1,2,3),h(1,2,3,4))        # volání s proměnným počtem argumentů
# print(h(range(4)),h(range(5)))  # nelze: TypeError
print(h(*range(4)),h(*range(5)))  # rozbalení posloupnosti ve skutečném argumentu

Přiřazovat skutečné argumenty k formálním argumentům lze pozičně nebo pomocí jmen formálních argumentů, jako v ukázce ve volání funkce g. Přiřazení argumentu jménem je možné i pro standardní funkce jazyka. Jaká výhoda, když je editor (např. VS Code) poučený a popis rozhraní standardních i nově definovaných funkcí při editaci nabízí.

Numba: překládání funkcí

Už na druhé přednášce jsme si ukázali na algoritmu sčítání aritmetické posloupnosti srovnání, že přeložený optimalizovaný fortranský program běží řádově rychleji než ekvivalentní pythonský skript (50krát pro indexovaný cyklus, 100krát pro cyklus s podmínkou). Příčinou je kontinuální interpretace pythonských cyklů – příkazy v jejich těle jsou stále znovu a znovu překládány, což trvá mnohem déle než následné provedení přeloženého kódu. Pomoc nabízí balíček Numba, jehož hlavní schopností je překlad pythonských funkcí. Ten je podmíněn řadou předpokladů (pro Python i NumPy), ovšem pokud nevybočujeme, užití Numby je snadné: před každou funkci určenou k překladu předepíšeme direktivu (dekorátor) @njit pro just-in-time překlad funkce v nopython režimu.

Přepišme naše staré ukázky pro součet aritmetické posloupnosti do tvaru funkcí, přidejme import Numby a dekorátory a zvyšme horní mez aritmetické posloupnosti na miliardu, stejně jako dříve ve Fortranu:

# Python & Numba: součet aritmetické posloupnosti, resp. harmonická čísla
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+i
  # s=s+1/i
  return s

@njit
def sum_for(n):
  s=0.
  for i in range(1,n+1):  # indexovany cyklus
    s=s+i
  # 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)

Srovnáme-li s dobou běhu ekvivalentního fortranského kódu, jsme tamtéž:

! Fortran: soucet aritmeticke posloupnosti, resp. harmonicka cisla
function sum_for(n) result (s)
  real(8) s
  s=0
  do i=1,n
    s=s+i
  ! 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

Když už máme před sebou tyto ukázky, zmíníme i další možnost, jak počítat v Pythonu rychle: můžeme volat z pythonského skriptu fortranské funkce. Slouží k tomu součást NumPy zvaná F2PY. V našem případě bychom ve Windows nebo Linuxu voláním příkazu

f2py -c -m f90 s.f90      # nebo f2py3 ...

přeložili (předpokladem je přítomný gfortran) naši fortranskou ukázku uloženou ve zdrojovém souboru s.f90 do pythonského modulu f90 (ve Windows vznikne soubor s příponou dll v podadresáři f90\.libs, který je třeba přikopírovat k našim skriptům) a ten pak užili tímto skriptem:

# Python & f2py: 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)

Doba běhu bude samozřejmě totožná s dobou běhu originálního fortranského kódu.

Rekurzivní funkce

Rekurzivní neboli na sebe se odvolávající zápisy jsou v matematice běžné s účelem převést úlohu na obdobnou, ale jednodušší. Např. pro n-faktoriál píšeme:

\[\begin{split}\begin{array}{rcl} 0! & = & 1\\ N! & = & N(N-1)! \quad\mbox{ pro }\quad N=1,\dots \end{array}\end{split}\]

V leckterých případech (i zde) existuje sice i alternativní nerekurzivní zápis, \(N!=\prod_{n=1}^N n\), ne však vždy. Rekurzivní algoritmy se často používají např. při práci s dynamickými datovými strukturami nebo při řešení úlohy půlením. (Př. Setřídit pole? Setřiď každou polovinu zvlášť a vhodně je sluč.) Programovací jazyky rekurzivní volání funkcí umožňují, přímo (jednotka volá sama sebe) i nepřímo (první jednotka volá druhou, druhá volá první). Nezbytnou součástí rekurzivní funkce je podmíněný příkaz s testem ukončení (zde: nerekurzivní větev pro 0-faktoriál).

Ukázka: nerekurzivní a rekurzivní implementace výpočtu součtu řady (raději než součinu ve faktoriálu, který brzy přeteče). V nerekurzivní funkci je počet aritmetických operací úměrný nmax, jak vidíme přímo z mezí cyklu. V rekurzivní funkci to bude také tak, a navíc se přičítá režie pro opakované volání (zanořování se do) funkce. To stojí nejen čas, ale i paměť, a zanoříme-li se hlouběji, snadno narazíme na limity systému. Python s Numbou a GNU Fortran (ve Windows) jsou ochotny se zanořit řádově 10000×, pak běh tiše spadne, Python bez Numby končí u zanoření řádově 1000×.

# Python: součet řady v cyklu vs. rekurzivně
import sys
from numba import njit
sys.setrecursionlimit(1000000)
imax = 10000                       # počet opakování výpočtu
nmax = 10000                       # argument pro výpočet (s Numbou lze více než bez ní)

# sčítání v cyklu
@njit
def Sum(nmax):
  result=0
  for n in range(1,nmax+1): result+=n
  return result

# rekurzivní funkce
@njit
def SumRec(n):
  if n<=0:
    result=0                       # nerekurzivní větev
  else:
    result=n+SumRec(n-1)           # rekurzivní volání
  return result

print('V cyklu...   ',end='')
for i in range(1,imax+1): result=Sum(nmax)
print(nmax,result)
print('Rekurzivne...',end='')
for i in range(1,imax+1): result=SumRec(nmax)
print(nmax,result)

Oblíbenou, ale tragicky neefektivní ukázkou bývá rekurzivní výpočet Fibonacciho čísel:

\[\begin{split}\begin{array}{rcl} F_0 & = & 0\\ F_1 & = & 1\\ F_n & = & F_{n-1}+F_{n-2} \quad\mbox{ pro }\quad n>1 \end{array}\end{split}\]

Pythonské řešení pomocí cyklu i rekurze může vypadat takto:

# Python: Fibonacciho posloupnost v cyklu vs. rekurzivně
from numba import njit
nmax=33
counter=0

# výpočet v cyklu
# @njit
def FibFor(nmax):
  if   nmax<=0: result=0
  elif nmax==1: result=1
  else:
    fm2=0; fm1=1
    for n in range(nmax-1):
      result=fm1+fm2
      fm2,fm1=fm1,result
  return result

# rekurzivní funkce
# @njit
def FibRec(n):
  global counter                        # nelze s Numbou
  counter+=1
  if   n<=0: result=0                   # nerekurzivní větve
  elif n==1: result=1
  else: result=FibRec(n-1)+FibRec(n-2)  # rekurzivní volání
  return result

print('V cyklu...   ',end=' ')
print(nmax,FibFor(nmax))
print('Rekurzivne...',end=' ')
print(nmax,FibRec(nmax),counter)

a pro nmax=40 si na rekurzivní řešení počkáme kolem minuty (proměnná counter dokládá, že rekurzivní funkce je volána 331_160_281krát). Sčítání pomocí cyklu je samozřejmě bleskové. S Numbou (a zakomentovanými řádky týkajícími se globální proměnné counter) bychom sice dobu běhu rekurze stáhli k sekundě, ale už pro nmax=50 bychom byli zpátky nad minutou – to bychom těch padesát součtů provedli rychleji ručně. Háček je zjevně v množství duplicitních výpočtů při dvojím rekurzivním volání v těle funkce, např. FibRec(48) je voláno pro hodnotu argumentu n=50 a nezávisle na tom i pro n=49, a tak pořád dál.

Lambda funkce

Krátké, tzv. lambda neboli anonymní funkce nabízejí úspornější syntaxi pro případ, kdy je tělo funkce místo sekvence příkazů tvořeno pouhým výrazem. S argumenty se zachází totožně jako u běžných (globálních) funkcí, stejně tak návratovou hodnotou může být skalár i kolekce. Bývají jako tzv. first-class functions argumentem funkcí vyššího řádu, např. populární trojice Map-Filter-Reduce:

# Python: lambda funkce v přiřazení a jako argumenty funkcí map-filter-reduce
f=lambda x: x**2
g=lambda x,y: x**2+y**2
h=lambda x: (x,x**2,x**3)
print('lambda:',f(2),g(3,4),h(5))                  # lambda: 4 25 (5, 25, 125)

from functools import reduce
a=[1,3,5,2,4]
print('list  :',a)                                 # list  : [1, 3, 5, 2, 4]
print('map   :',list(map(lambda x: x**2,a)))       # map   : [1, 9, 25, 4, 16]
print('filter:',list(filter(lambda x: x%2!=0,a)))  # filter: [1, 3, 5]
print('reduce:',reduce(lambda x,y: x+y,a))         # reduce: 15

Modularizace

Funkce sdružené v souboru vytvářejí modul; modul module je třeba před použitím importovat příkazem import module a jméno modulové funkce, např. f(), při volání kvalifikovat jménem modulu, module.f(), nebo zkráceně po import module as m jen m.f(). Jednotlivé funkce z modulu lze importovat pomocí from m import f, všechny funkce také pomocí from m import *; taková jména pak není třeba kvalifikovat. Z adresáře s moduly lze vytvořit balíček. Součástí distribuce Pythonu je The Python Standard Library.

# Python: modul s funkcemi
f=lambda x: x**2
g=lambda x,y: x**2+y**2
h=lambda x: (x,x**2,x**3)
def ff(x):
  return x**2
def gg(x,y):
  return x**2+y**2
def hh(x):
  return x,x**2,x**3
# Python: volání funkcí z modulu m
import m
print('module m:',m.f(2), m.g(3,4), m.h(5))   # module m: 4 25 (5, 25, 125)
print('module m:',m.ff(2),m.gg(3,4),m.hh(5))  # module m: 4 25 (5, 25, 125)
from m import *
print('module m:',ff(2),  gg(3,4),  hh(5))    # module m: 4 25 (5, 25, 125)