Cvičení č. 4: Fraktál

Připravíme program pro otestování, zda bod komplexní roviny je (přesněji zda není) prvkem Mandelbrotovy množiny. Prověříme pak mnoho bodů, na cvičení volených náhodně, za domácí úkol vybíraných systematicky, a zhotovíme z nich obrázek podobný obrázku z Wikipedie). Pro pěkné rozlišení obrázku potřebujeme prověřit řádově miliony bodů a pro každý bod to budou až tisíce operací; v Pythonu proto použijeme balíček Numba, abychom se přiblížili fortranské rychlosti běhu.

Mandelbrotova množina

Populární fraktál zvaný Mandelbrotova množina je množinou komplexních čísel \(c\), pro které posloupnost komplexních čísel \(z_n\),

\[z_{n+1}=z_n^2+c,\quad z_0=0,\]

pro \(n\to\infty\) zůstane omezená. Neomezenost posloupnosti pro pevné \(c\) lze po úvaze ztotožnit s existencí únikového indexu \(N\), pro který (poprvé) nastane \(|z_N|>2\). Bod komplexní roviny do Mandelbrotovy množiny buď patří nebo nepatří, obrázky jako tento navíc probarvují vyloučené body barvou přiřazenou únikovému indexu \(N\).

V následujícím pythonském skriptu se zvolený obdélník v komplexní rovině pokrývá náhodně vybíranými body, pro které se provede až nmax kroků uvedené posloupnosti a testuje se splnění únikové podmínky. Na výstup se posílají souřadnice testovaného bodu a menší z čísel \(N\) a nmax+1. Body volíme náhodně pomocí funkce random.random(), která vrací (při každém volání jiné) pseudonáhodné reálné číslo z intervalu \(\langle 0,1)\).

# Python: příprava bodů pro kreslení Mandelbrotovy množiny
import random
xmin=-2; xmax=0.6; ymin=-1.2; ymax=1.2  # obdélník v komplexní rovině
nmax=100                     # max. počet kroků testovací posloupnosti
random.seed()
while True:
  x=random.random()          # náhodné číslo z intervalu <0,1)
  y=random.random()
  x=xmin+(xmax-xmin)*x       # přepočet na <xmin,xmax)
  y=ymin+(ymax-ymin)*y       # přepočet na <ymin,ymax)
  c=complex(x,y)             # náhodný bod v obdélníku
  z=0+0j                     # výchozí z_0
  nstop=0
  for n in range(nmax):
    z=z*z+c
    if abs(z)>2:
      nstop=n+1              # únikový index nalezen
      break
  if nstop==0: nstop=nmax+1  # únikový index nenalezen, volíme nmax+1
  print('%8.5f%9.5f%5d'%(x,y,nstop))
  # print('{:8.5f} {:8.5f} {:4d}'.format(x,y,nstop))
  # print(f'{x:8.5f} {y:8.5f} {nstop:4d}')

Ve vnitřním cyklu for n se mnohokrát provádí zhruba desítka elementárních operací (real sčítání a násobení) a 1 real odmocnina. Obojí může být svou časovou náročností srovnatelné a můžeme zkusit eliminovat zbytnou odmocninu ve výrazu abs(z)>2 nahrazením ekvivalentní podmínkou z.real**2+z.imag**2>4. Pro rozumné rozlišení výstupních obrázků však potřebujeme podstatně větší zrychlení skriptu, pro což se opět obrátíme k balíčku Numba. Jak jsme viděli minule, Numba umí překládat funkce, náš skript tedy upravíme v procedurálním stylu: test příslušnosti bodu k Mandelbrotově množině shrneme do funkce testMandelbrot(x,y), vracející únikový index, resp. nmax+1 při jeho nenalezení. Vypisovat data x, y a nstop budeme vně funkce, Numba nerada překládá příkaz print. (Přeloží jednoduchý print, nepřeloží formátovaný print nebo print do souboru.) Ostatně, jak si ještě na přednášce vyložíme, výpisy do funkce nepatří.

Ustoupíme také od nekonečného cyklu (přerušit lze klávesovou kombinaci Ctrl+C) k cyklu běžícímu MaxTime sekund. K tomu použijeme modul time s funkcí time.time(), který vrací reprezentaci aktuálního času v sekundách. Obvykle se uloží startovací čas, tic=time.time(), a koncový čas, toc=time.time(), a doba běhu měřeného úseku je pak toc-tic.

Pro rozšíření našeho repertoáru ještě rozlišíme směrování výpisů: souřadnice bodů a únikový index budeme vypisovat na standardní výstup stdout, zatímco souhrnnou statistiku (počet bodů cnt1 v Mandelbrotově množině, celkový počet prověřených bodů cnt2 a jejich procentuální podíl) pošleme ven alternativním kanálem zvaným standardní chybový výstup stderr (zpřístupněný modulem sys). Standardní výstup pak z příkazového řádku (Windows nebo Linuxu) přesměrujeme pomocí operátoru > do souboru, stderr ponecháme na obrazovce.

# Python: příprava bodů pro kreslení Mandelbrotovy množiny
# Balíčky Numba pro zrychlení, time pro měření času a sys pro stderr

from numba import njit
xmin=-2; xmax=0.6; ymin=-1.2; ymax=1.2  # obdélník v komplexní rovině
nmax=100                    # max. počet kroků testovací posloupnosti
MaxTime=10                  # doba výpočtu [sec]
output=1                    # 0: jen statistika, 1: data na obrazovku

@njit
def testMandelbrot(x,y):    # výpočet únikového indexu pro bod (x,y)
  c=complex(x,y)
  z=0+0j
  nstop=0
  for n in range(nmax):
    z=z*z+c
    if abs(z)>2:
    # if z.real**2+z.imag**2>4:
      nstop=n+1
      break                 # únikový index nalezen, jinak nmax+1
  return nstop if nstop>0 else nmax+1

import random               # pro náhodná čísla
import time                 # pro měření času
import sys                  # pro výpis na stderr
cnt1=0; cnt2=0              # čítače bodů
tic=time.time()             # startovací čas
while time.time()-tic<MaxTime:        # dokud je doba výpočtu menší než MaxTime
  x=xmin+(xmax-xmin)*random.random()  # náhodné číslo v <xmin,xmax)
  y=ymin+(ymax-ymin)*random.random()  # náhodné číslo v <ymin,ymax)
  nstop=testMandelbrot(x,y)           # únikový index
  if output>0: print('%8.5f%9.5f%5d'%(x,y,nstop))
  # print('{:8.5f} {:8.5f} {:4d}'.format(x,y,nstop))
  # print(f'{x:8.5f} {y:8.5f} {nstop:4d}')
  if nstop==nmax+1: cnt1+=1 # čítač bodů v Mandelbrotově množině
  cnt2+=1                   # čítač všech prověřených bodů
sys.stderr.write('Mandelbrot: '+str(cnt1)+' / '+str(cnt2)+' = '+f'{100*cnt1/cnt2:5.1f}'+' %\n')
                            # výpis na stderr pomocí write(string), včetně \n pro odřádkování

Soubor s výstupními daty bychom mohli získat obdobně jako na minulém cvičení. Teď si vyzkoušíme přesměrování standardního výstupu skriptu z obrazovky do souboru příkazem:

python mandelbrot.py > mandelbrot.dat

Pozor – pokud přesměrování výstupu provedeme ve VS Code v terminálu s PowerShellem (ve výzvě příkazového řádku je PS), výstupní soubor bude kódován v UTF-16 a nebude pro gnuplot čitelný. Rozumné je změnit Terminal: Select Default Profile (po F1 v Command Palette) na Command Prompt, kde problém s kódováním nevzniká. Terminál PowerShellu se pak uzavře pomocí Kill (Delete) neboli ikonkou koše vpravo u ikony powershell a otevře se nový terminál, nyní už Command Prompt.

Fortran

Fortranská verze je obdobná. Funkci testMandelbrot(x,y) vnoříme do programu za slůvko contains. Náhodná čísla z <0,1) generujeme voláním podprogramu call random_number(x). Pro výpis na stdout používáme příkaz print *,data, zatímco pro výpis na stderr použijeme (obdobně jako minule pro výpis do souboru) příkaz write (0,*) data, kde hvězdičku případně nahradíme formátovou specifikací. Měření doby běhu provádíme s užitím matlabovského názvosloví tic toc, simulovaného níže přiloženým fortranským modulem mTicToc.

! Fortran: příprava bodů pro kreslení Mandelbrotovy množiny
! Modul mTicToc pro měření času v matlabovském stylu tic-toc

include 'mTicToc.f90'           ! vložení obsahu souboru s modulem mTicToc

program Mandelbrot
use mTicToc
implicit none
real :: xmin=-2,xmax=0.6,ymin=-1.2,ymax=1.2 ! obdélník v komplexní rovině
integer :: nmax=100             ! max. počet kroků testovací posloupnosti
real :: MaxTime=10              ! doba výpočtu [sec]
integer:: output=1              ! 0: jen statistika, 1: data na obrazovku
real :: x,y
integer :: nstop,cnt1,cnt2

cnt1=0; cnt2=0                  ! čítače bodů
call tic()                      ! uložení startovacího času
do while (toc()<MaxTime)        ! dokud je doba výpočtu menší než MaxTime
  call random_number(x)         ! náhodné číslo v <0,1)
  call random_number(y)
  x=xmin+(xmax-xmin)*x          ! přepočet do <xmin,xmax)
  y=ymin+(ymax-ymin)*y
  nstop=testMandelbrot(x,y)     ! únikový index
  if (output>0) print '(f8.5,f9.5,i5)',x,y,nstop
  if (nstop==nmax+1) cnt1=cnt1+1 ! čítač bodů v Mandelbrotově množině
  cnt2=cnt2+1                   ! čítač všech prověřených bodů
enddo
write (0,'(a,i0,a,i0,a,f5.1,a)') 'Mandelbrot: ',cnt1,' / ',cnt2,' = ',100.*cnt1/cnt2,' %'
                                ! výpis statistiky na stderr
contains

! Výpočet únikového indexu pro bod (x,y)
integer function testMandelbrot(x,y) result (nstop)
real x,y
complex c,z
integer n
c=cmplx(x,y)
z=(0.,0.)
nstop=0
do n=1,nmax
  z=z*z+c
  if (abs(z)>2) then
  ! if (z%Re**2+z%Im**2>4) then
    nstop=n
    exit                        ! únikový index nalezen
  endif
enddo
if (nstop==0) nstop=nmax+1      ! únikový index nenalezen, vrátí se nmax+1
end function

end program

Pro měření času použijeme a budeme používat modul mTicToc inspirovaný matlabovským měřením času. Voláním podprogramu call tic() zafixujeme startovací čas, voláním funkce toc() získáme počet sekund uplynulých od startovacího času.

! Wallclock timer
! Usage: call tic(); call task(); print *,toc()

MODULE mTicToc

implicit none
private
public tic,toc
integer,parameter :: wp=8
real(wp) :: time

CONTAINS

subroutine tic()
integer count,count_rate
call system_clock(count,count_rate)
time=real(count,wp)/real(count_rate,wp)
end subroutine

real(wp) function toc()
integer count,count_rate
call system_clock(count,count_rate)
toc=real(count,wp)/real(count_rate,wp)-time
end function

END MODULE

Program přeložíme s optimalizací a při nenulové hodnotě proměnné output spustíme příkazem:

gfortran -O -o mandelbrot mandelbrot.f90 && mandelbrot > mandelbrot.dat

Benchmark

Porovnáme rychlost skriptů v interpretovaném Pythonu, v Numbou překládaném Pythonu a v optimalizovaném Fortranu (gfortran -O). Zahrneme varianty s pouhým výpočtem (bez pomalého výpisu všech dat), s ukládáním do souboru a s vypisováním na obrazovku. Ve variantách bez výpisu zkusíme nahradit test abs(z)>2 matematicky ekvivalentním z.real**2+z.imag**2>4, ovšem snad rychlejším, neboť neobsahujícím odmocninu. Programy necháme běžet 10 sekund a shrneme počty prověřených bodů z obdélníka <-2,0.6>x<-1.2,1.2>.

Tabulka počtu prověřených bodů pro 10sekundový běh (typ procesoru: Intel i5-11600K, Windows, zdrojové soubory)

varianta

milionů bodů

poznámka

Python no sqrt

1.1

bez výpisů do souboru stejně pomalé jako s nimi

Python abs

2.8

bez výpisů do souboru zrychlení 2x

Python file

1.2

pomalý výpis do souboru

Python screen

0.1

ještě pomalejší výpis na obrazovku

Numba no sqrt

18.5

Numba: zrychlení 17x

Numba abs

16.7

Numba file

1.2

při pomalých výpisech není zrychlení výpočtů znát

Numba screen

0.1

gfortran no sqrt

95.0

výpočty 5x rychlejší než Numba

gfortran abs

39.0

gfortran file

5.8

gfortran screen

0.1

výpis na obrazovku pokaždé stejně pomalý

Vyplývají odsud následující poznatky: výpisy na obrazovku bývají o řád pomalejší než výpisy do souboru (na SSD disku) a vycházejí zde stejně pomalu pro Python, Numbu i Fortran; při vypnutí výpisů běží interpretovaný Python stejně rychle nebo 2krát rychleji (podle použité podmínky), zatímco Python s Numbou běží rychleji o řád a gfortran je ještě 5krát rychlejší; eliminací odmocniny zrychlila Numba mírně a Fortran více než 2krát. Celkově je mezi pomalým výpisem na obrazovku a rychlými fortranskými počty bez výpisů třířádový rozdíl (100 tisíc vs. 100 milionů prověřených bodů za 10 sekund).

DÚ: Fraktál v gnuplotu

Pro vykreslení obrázku použijeme gnuplot. Následující skript pro gnuplot mandelbrot.gp definuje spektrum (set palette) a rozsah (set cbrange) barevné palety, nastaví PNG výstup o zvoleném rozlišení (set term pngc size) do souboru (set out file) a volá příkaz plot file palette pro datový soubor, v jehož prvních dvou sloupcích se očekávají souřadnice kresleného bodu a ve třetím sloupci index barvy. Příkazem unset key jsme vypnuli popisek se jménem souboru, pomocí příkazů unset xtics; unset ytics bychom vypnuli popisy os.

# gnuplot: vykreslení Mandelbrotovy množiny
set palette defined \
  (0 'dark-blue',6 'dark-blue',30 'white',70 'gold',100 'brown',100.001 'black',101 'black')
set cbrange [0:101]
set terminal pngcairo size 1024,768
# set term pngcairo size 2048,1536
set output 'mandelbrot.png'
unset key
# unset xtics; unset ytics
plot [-2.:.6][-1.2:1.2] 'mandelbrot.dat' using 1:2:3 with dots palette
# plot 'mandelbrot.dat' using 1:2:3 with dots palette
set output

Skript se spustí buď v gnuplotu příkazem load 'mandelbrot.gp' nebo mimo gnuplot příkazem gnuplot mandelbrot.gp. Tento příkaz můžeme vložit do konfigurace extenze Code Runner editoru VS Code a skript spouštět z něj, včetně prohlížení výstupního PNG obrázku. Při systematickém průchodu obdélníkem by výsledný obrázek mohl vypadat podobně jako tento:

_images/mandelbrot.png

Domácí úkol č. 4: uložte si výstupní data do souboru a pomocí gnuplotu a přiloženého skriptu vytvořte PNG obrázek. Obrázky získané výše uvedenými skripty budou znehodnoceny bílými tečkami, neboť náhodně vybíranými body se v rozumném čase celý obdélník nepokryje. Upravte tedy (součást domácího úkolu) program tak, aby se zvolený obdélník pokrýval systematicky a s předepsaným rozlišením, tedy projděte po řadě dva vnořené indexované cykly pro reálné a imaginární souřadnice testovaných bodů. Výpočetní rozlišení bude vhodné volit jemnější než zobrazovací rozlišení, řekněme 2krát v každém směru, jinak budou gnuplotovské obrázky mřížkované. A také zoomujte obdélníkem blíže k hranici Mandelbrotovy množiny, ať vidíme nějaký fraktální detail. Při hlubokém zoomování pozor na vhodnou volbu nmax a také na formátovací specifikaci v příkazu print; obojí může být potřeba zvětšit.

Nápověda k vnořeným cyklům:

# Python a průchod skrz uzly obdélníka <0,xmax> x <0,ymax> s výpisem vzdálenosti od uzlu (0,0)
ixmax=3; iymax=4
for ix in range(ixmax+1):
  for iy in range(iymax+1):
    print(ix,iy,f'{abs(complex(ix,iy)):5.3f}',end='/')
print('stop')
! Fortran a průchod skrz uzly obdélníka
ixmax=3; iymax=4
do ix=0,ixmax
  do iy=0,iymax
    print '(2(i0,1x),f5.3,a$)',ix,iy,abs(cmplx(ix,iy)),'/'
  enddo
enddo
print '(a)','stop'
end

Výkonnější algoritmy pro diagnostiku Mandelbrotovy množiny popisuje stránka Plotting algorithms for the Mandelbrot set. Můžeme tam nalézt i postup, jak dodat obrázkům spojité barevné přechody mezi diskrétními hodnotami únikového indexu. Posloupnost \(z_n\) vedoucí k Mandelbrotově množině lze zobecnit a kreslit např. Multibrot sets. K více dimenzím směřuje zobecnění zvané Mandelbulb.

Také jsme dříve na cvičení někdy stihli i jiný typ úchvatných fraktálů: logistická mapa a Markusovy-Ljapunovovy fraktály (skripty v Pascalu).

Py a f90 zdrojové soubory s dnešními ukázkami: zip.