Urychlování výpočtů

Antonín Slavík
Matematicko-fyzikální fakulta UK

Příklad: Mandelbrotova množina

n-tá iterace funkce:  (f^n)(z)=(f ◦ ⋯ ◦ f)(z)
Mandelbrotova funkce:  (Subscript[f, c])(z)=(z^2)+c   (c, z jsou komplexní čísla)
Mandelbrotova množina:  M={c∈ℂ; posloupnost (Subscript[{(Subscript[f, c]^n)(0)}, n=1]^∞) je omezená}

Jak vypadají iterace (Subscript[f, c]^n)(0)?

f[z_]:=(z^2)+c

NestList[f,0,5]//TableForm

Jak vypadá Mandelbrotova množina?

Zvolíme výřez komplexní roviny: {x+i y∈ℂ; (Subscript[x, min])≤x≤(Subscript[x, max])& (Subscript[y, min])≤y≤(Subscript[y, max])}

Určíme rozlišení obrázku ve vodorovném směru ((Subscript[x, res])) a ve svislém směru ((Subscript[y, res]))

Ve zvoleném výřezu umístíme (Subscript[x, res])× (Subscript[y, res]) bodů, které odpovídají různým hodnotám parametru c

Pro každou hodnotu c rozhodneme, zda (Subscript[{(Subscript[f, c]^n)(0)}, n=1]^∞) je omezená

P\[RHacek]ibli\[ZHacek]ná metoda: Testujeme, zda b\[EHacek]hem prvních 100 iterací opustí (Subscript[f, c]^n)(0) kruh se st\[RHacek]edem 0 a polom\[EHacek]rem 2. 
(V\[EHacek]ta: Pokud pro n\[EHacek]jaké k∈ℕ platí \[LeftBracketingBar](Subscript[f, c]^k)(0)\[RightBracketingBar]>2, pak (lim_(n→∞))(Subscript[f, c]^n)(0)=∞.)

iCount[c_] := Module[{cnt = 0, z = 0}, While[Abs[z] ≤2 && cnt < 100, cnt++; z = (z^2)+c]; cnt]

mandelbrot[{{xmin_, xmax_}, {ymin_, ymax_}}, {xres_, yres_}] := 
    Table[iCount[x + I*y], {y, ymin, ymax, (ymax - ymin)/(yres-1)}, {x, xmin, xmax, (xmax - xmin)/(xres-1)}]

showMatrix[matrix_] := ArrayPlot[matrix, Mesh -> False, Frame -> False, 
   ColorFunction -> (If[#1 == 1, Black, Hue[(5*#1)/6]] & )]

showMatrix[mandelbrot[{{-2.1,0.5},{-1.3,1.3}},{100,100}]]

Timing[A=mandelbrot[{{-2.1,0.5},{-1.2,1.2}},{400,400}];]

showMatrix[A]

Poznámka: V nových verzích Mathematiky je součástí jádra funkce MandelbrotSetPlot pro zobrazování Mandelbrotovy množiny.

MandelbrotSetPlot[{-2.1-1.2I,0.5+1.2I},ColorFunction->Hue]

Sledování průběhu výpočtu

Během výpočtu chceme informovat uživatele o tom, kolik procent je hotovo.

hotovo (v procentech) = ((y-ymin)/(ymax-ymin))·100

K zobrazování informace použijeme funkci Monitor.

?Monitor

mandelbrot[{{xmin_, xmax_}, {ymin_, ymax_}}, {xres_, yres_}] := 
   Monitor[Table[Table[iCount[x + I*y], 
       {x, xmin, xmax, (xmax - xmin)/(xres-1)}], {y, ymin, ymax, (ymax - ymin)/(yres-1)}], 
    StringForm["`` % done.",Round[(100*(y - ymin))/(ymax - ymin)]]];

A=mandelbrot[{{-2.1,0.5},{-1.2,1.2}},{400,400}];

ProgressIndicator je grafický prvek určený k zobrazování průběhu výpočtu. Můžeme jej použít společně s Monitor.

?ProgressIndicator

ProgressIndicator[8,{0,10}]

mandelbrot[{{xmin_, xmax_}, {ymin_, ymax_}}, {xres_, yres_}] := 
Monitor[Table[Table[iCount[x + I*y], {x, xmin, xmax, (xmax - xmin)/(xres-1)}], 
      {y, ymin, ymax, (ymax - ymin)/(yres-1)}], ProgressIndicator[y,{ymin,ymax}]];

A=mandelbrot[{{-2.1,0.5},{-1.2,1.2}},{400,400}];

Paralelní výpočty

Možnosti paralelních výpočtů v Mathematice:

Jeden počítač, vícejádrový procesor (Local Kernels)

Více počítačů v síti (Lightweight Grid, Cluster Integration, Remote Kernels)

Viz Parallel Kernel Configuration v menu Evaluation a dokumentaci v Parallel Computing Tools User Guide.

$ProcessorCount vrací počet jader, které jsou k dispozici v procesoru počítače.

Kernels[] vrací seznam aktuálně běžících jader.

LaunchKernels[] spustí všechna dostupná jádra.

$ProcessorCount

Kernels[]

LaunchKernels[]

Parallelize[výraz] vyhodnotí zadaný výraz, přičemž se pokusí výpočet paralelizovat.

Příklad: Hledání Mersennových prvočísel, tj. prvočísel tvaru (2^k)-1.
Pro měření délky paralelního výpočtu se místo Timing používá AbsoluteTiming (Timing nezapočítává čas, kdy běží paralelní jádra).

AbsoluteTiming[Select[Range[2000], PrimeQ[2^#-1]&]]

AbsoluteTiming[Parallelize[Select[Range[2000], PrimeQ[2^#-1]&]]]

Příklad: Paralelní výpočet Mandelbrotovy množiny

mandelbrot[{{xmin_, xmax_}, {ymin_, ymax_}}, {xres_, yres_}] := 
    Parallelize[Table[iCount[x + I*y], {y, ymin, ymax, (ymax - ymin)/(yres-1)}, {x, xmin, xmax, (xmax - xmin)/(xres-1)}]]

Uživatelem definované funkce, které jsou k paralelnímu výpočtu zapotřebí, je nutné distribuovat jednotlivým jádrům.

DistributeDefinitions[iCount]

AbsoluteTiming[mandelbrot[{{-2.1,0.5},{-1.2,1.2}},{400,400}];]

(Využití jednotlivých jader je vidět v Parallel Kernel Status v menu Evaluation.)

Kompilace funkcí do Wolfram Virtual Machine

Potřebujeme-li často vyhodnocovat nějakou funkci, je dobré ji „zkompilovat“. Kompilace je přeložení kódu v Mathematice do sledu instrukcí zapsaných ve vnitřním jazyce Mathematiky (Wolfram Virtual Machine).

Nemá smysl kompilovat zabudované funkce (ty jsou již optimalizovány), pouze uživatelem definované funkce sestavené z následujících funkcí:

Compile`CompilerFunctions[]//Sort

Výhody a nevýhody:

Kompilované funkce jsou často mnohonásobně rychlejší

Kompilované funkce mohou pracovat pouze s celými čísly, s reálnými a komplexními čísly v machine precision nebo se seznamy uvedených typů čísel

Je třeba předem specifikovat datové typy použitých proměnných (kompilace se přizpůsobuje datovým typům)

?Compile

Funkce pro výpočet hodnoty polynomu:

f[x_]=(x^5)+3 (x^4)+2 x+1;

Jak dlouho trvá výpočet 100 000 funkčních hodnot?

Timing[Table[f[x],{x,0,1,(1/(10^5))}];]

Stejná funkce v kompilované podobě. Pozor, funkci je třeba definovat pomocí = (nikoliv pomocí :=).

g=Compile[{x},(x^5)+3 (x^4)+2 x+1];

Timing[Table[g[x],{x,0,1,(1/(10^5))}];]

Při kompilaci funkce potřebuje Mathematica znát datové typy všech proměnných. Pokud není řečeno něco jiného, předpokládá, že jsou reálné.

Do funkce f lze dosazovat i komplexní čísla:

f[1+I]

Pokus o použití kompilované funkce s komplexním argumentem selže, je vyhodnocen nezkompilovaný funkční předpis.

g[1+I]

Deklarace x jako komplexní proměnné (povolené datové typy jsou _Integer, _Real, _Complex nebo seznamy těchto hodnot):

g=Compile[{{x,_Complex}},(x^5)+3 (x^4)+2 x+1];

g[1+I]

I pro reálné argumenty je nyní výsledek v komplexním tvaru:

g[0]

Zkompilovaná funkce neumí počítat s přesnými čísly:

g[Sqrt[2]]

f[Sqrt[2]]

Kompilace funkce, jejímž argumentem je seznam:

body=RandomReal[1,{(10^6),2}];(* (10^6) náhodných bodů ve čtverci [0,1]x[0,1] *)

Timing[Map[{-#[[2]],#[[1]]}&,body];](* rotace bodů o 90 stupňů *)

rotace=Compile[{{pt,_Real,1}},{-pt[[2]],pt[[1]]}]

Timing[Map[rotace,body];]

iCount – původní verze:

iCount[c_] := Module[{cnt = 0, z = 0}, While[Abs[z]≤2 && cnt < 100, cnt++; 
z = z^2 + c]; cnt]

Quit[]

iCount – kompilovaná verze:

iCount=Compile[{{c,_Complex}},Module[{cnt = 0,z=0.0I},While[Abs[z]≤2 && cnt < 100, cnt++; z = (z^2)+c]; cnt]];
mandelbrot[{{xmin_, xmax_}, {ymin_, ymax_}}, {xres_, yres_}] := 
    Table[iCount[x + I*y], {y, ymin, ymax, (ymax - ymin)/(yres-1)}, {x, xmin, xmax, (xmax - xmin)/(xres-1)}]
showMatrix[matrix_] := ArrayPlot[matrix, Mesh -> False, Frame -> False, 
   ColorFunction -> (If[#1 == 1, Black, Hue[(5*#1)/6]] & )]

Poznámka: Přiřazením z=0.0I  oznamujeme, že z je komplexní číslo.

Timing[A=mandelbrot[{{-2.1,0.5},{-1.2,1.2}},{400,400}];]

showMatrix[A]

Mandelbrotova množina – paralelní výpočet pomocí kompilované funkce:

mandelbrot[{{xmin_, xmax_}, {ymin_, ymax_}}, {xres_, yres_}] := 
    Parallelize[Table[iCount[x + I*y], {y, ymin, ymax, (ymax - ymin)/(yres-1)}, {x, xmin, xmax, (xmax - xmin)/(xres-1)}]]

LaunchKernels[]

DistributeDefinitions[iCount]

AbsoluteTiming[A=mandelbrot[{{-2.1,0.5},{-1.2,1.2}},{400,400}];]

Clear[iCount,mandelbrot]

Kompilace funkcí do jazyka C

Funkce je možné kompilovat nejen do Wolfram Virtual Machine, ale i do jazyka C (a následně do strojového kódu) pomocí volby CompilationTarget->"C". Takové funkce obvykle počítají několikanásobně rychleji, než při kompilaci do Wolfram Virtual Machine.

Na počítači musí být nainstalován vhodný překladač C (ve Windows např. MS Visual Studio Compiler, v Linuxu gcc).

Needs["CCompilerDriver`"]

CCompilers[]

DefaultCCompiler[]

iCount=Compile[{{c,_Complex}},Module[{cnt = 0,z=0.0I},While[Abs[z]≤2 && cnt < 100, cnt++; z = (z^2)+c]; cnt],CompilationTarget->{"C","WorkingDirectory"->NotebookDirectory[]}];
mandelbrot[{{xmin_, xmax_}, {ymin_, ymax_}}, {xres_, yres_}] := 
    Table[iCount[x + I*y], {y, ymin, ymax, (ymax - ymin)/(yres-1)}, {x, xmin, xmax, (xmax - xmin)/(xres-1)}]

Timing[A=mandelbrot[{{-2.1,0.5},{-1.2,1.2}},{400,400}];]

showMatrix[A]

Poznámka. Mathematica umí komunikovat s jinými aplikacemi pomocí Wolfram Symbolic Transfer Protocol (WSTP). Ty části výpočtu, které chceme urychlit, můžeme naprogramovat v jiném jazyce (C/C++, Java, .NET). Podrobněji viz WSTP and External Program Communication.

Další poznámky k efektivitě programů

Ukládání hodnot funkcí („dynamické programování“)

Uvažujme definici funkce v následujícím tvaru:

f[x_] := f[x] = definice funkce

Při prvním vyhodnocení f[x0] s konkrétní hodnotou x0 se zároveň vytvoří nové globální pravidlo tvaru f[x0]=hodnota (uloží se tedy vlastně hodnota f[x0]). Bude-li později potřeba opět vyhodnotit f[x0], použije se toto nové pravidlo, neboť má přednost před obecnějším pravidlem f[x_]:=...

Tato technika se někdy nazývá dynamické programování.

fib[0]=1;
fib[1]=1;
fib[n_]:=fib[n-1]+fib[n-2]

Timing[fib[30]]

Výpočet je velmi pomalý, protože hodnoty se počítají opakovaně:

Trace[fib[8],fib[_]]

Dynamické programování:

fib2[0]=1;
fib2[1]=1;
fib2[n_]:=fib2[n]=fib2[n-1]+fib2[n-2]

Timing[fib2[30]]

Pravidla uložená v paměti:

?fib2

Clear[fib,fib2]

Práce se seznamy

Funkce pracující s celými seznamy jsou v Mathematice velmi rychlé. Doporučuje se co nejvíce využívat funkcí jako Map, MapThread, Select, Cases, funkcí s atributem Listable apod. Naopak je vhodné vyvarovat se provádění operací po jednotlivých prvcích pomocí cyklů, jako je např. For.

Chceme-li např. vykreslit seznam bodů v rovině, můžeme použít buď {Point[{a,b}],Point[{c,d}],...}, nebo Point[{{a,b},{c,d},...}]; druhá varianta se vykresluje výrazně rychleji.

pts=RandomReal[1,{10000,2}];(* 10000 náhodných bodů ve čtverci [0,1]x[0,1] *)

pts1=Point[pts];(* Point[{{a,b},{c,d},...}] *)

pts2=Map[Point,pts]; (* {Point[{a,b}],Point[{c,d}],...} *)

Timing měří pouze čas, který potřebovalo jádro k provedení výpočtu. Chceme měřit i čas potřebný k vykreslení, proto použijeme SessionTime[] (vrací čas v sekundách, který uplynul od spuštění jádra).

s=SessionTime[];
Graphics[{PointSize[0.006],pts1},Frame->True]
SessionTime[]-s

s=SessionTime[];
Graphics[{PointSize[0.006],pts2},Frame->True]
SessionTime[]-s

Máme-li různobarevné body, můžeme použít volbu VertexColors.

colors=Map[Hue[#/2]&,Range[0,1,1/(Length[pts]-1)]];
Graphics[{PointSize[0.006],Point[pts,VertexColors->colors]}]

Clear[s,pts,pts1,pts2]

Podobně se chovají také Line a Polygon:

Line[seznam čar] je rychlejší než {Line[...],Line[...],...}.
Polygon[seznam polygonů] je rychlejší než {Polygon[...],Polygon[...],...}.

Opakované provádění operací, které mění délku seznamu (Append, AppendTo, Prepend, PrependTo, Insert, Delete) je pomalé.

Generování náhodných čísel tak dlouho, dokud součet nepřekročí 100 000:

SeedRandom[0];Timing[a = {};sum = 0;While[sum < 100000,r = RandomInteger[9];sum += r;AppendTo[a, r]];
Length[a]]

Rychlejší verze pomocí Sow a Reap:

SeedRandom[0];Timing[sum = 0;{r,{a}} = Reap[While[sum < (10^5), r = RandomInteger[9]; sum += r; Sow[r]]];Length[a]]

Clear[a]

Úsporná reprezentace seznamů

SparseArray

SparseArray slouží k úsporné reprezentaci řídkých seznamů nebo matic, ve kterých má většina prvků stejnou (typicky nulovou) hodnotu.

SparseArray[{2->a,5->b}]

Normal převádí SparseArray na klasický seznam.

Normal[%]

Zadání rozměru:

SparseArray[{2->a,5->b},10]

Normal[%]

Matice:

SparseArray[{{1,2}->a,{3,2}->b,{3,3}->c}]

Normal[%]

SparseArray[{3->a,5->b},7,1](* nespecifikované prvky nebudou nuly, ale jedničky *)

Normal[%]

Pravidla mohou obsahovat také vzory:

Normal[SparseArray[i_->i^2,10]]

Normal[SparseArray[{{_,3}->a,{i_,i_}->1},{4,4}]]

Konverze obyčejné matice na řídkou:

dense=DiagonalMatrix[Range[1000]];

sparse=SparseArray[dense]

Úspora paměti:

ByteCount[dense]-ByteCount[sparse]

Clear[dense,sparse]

SparseArray se chová stejně jako seznam:

s=SparseArray[{2->a,4->b,5->c},10]

s[[2]]

Length[s]

s.s

Clear[s]

m=SparseArray[{i_,j_}/;Abs[i-j]<=1:>RandomInteger[{1,10}],{100,100}]

MatrixPlot[m]

Det[m]

m10=MatrixPower[m,10]

MatrixPlot[m10]

Clear[m,m10]

PackedArray

PackedArray je úsporná reprezentace seznamů, ve kterých mají všechny prvky stejný datový typ (Integer, Real nebo Complex). Mathematica ji často používá automaticky. Developer`PackedArrayQ[seznam] vrací True, pokud je daný seznam uložen jako PackedArray.

realVec = Table[RandomReal[], {10000000}];
intVec = Table[RandomInteger[], {10000000}];
list={1,2,5,6};

{Developer`PackedArrayQ[realVec], Developer`PackedArrayQ[intVec],Developer`PackedArrayQ[list]}

Developer`ToPackedArray[seznam] a Developer`FromPackedArray[seznam] převádějí mezi obyčejnou a úspornou reprezentací.

ByteCount[realVec] (* spotřeba paměti v úsporné reprezentaci: 8 bajtů na jedno desetinné číslo *)

ByteCount[Developer`FromPackedArray[realVec]] (* spotřeba paměti ve standardní reprezentaci: 24 bajtů na jedno desetinné číslo *)

Výpočty s PackedArray jsou rychlejší:

Norm[realVec] // AbsoluteTiming

realVec2=Developer`FromPackedArray[realVec];

Norm[realVec2] // AbsoluteTiming

Přidání celého čísla na konec seznamu znemožní úspornou reprezentaci:

ByteCount[Append[realVec, 1]]

Pokud požádáme o PackedArray, převede se celé číslo na desetinné:

appendedRealVec = Developer`ToPackedArray[Append[realVec, 1], Real];
ByteCount[appendedRealVec]

Last[appendedRealVec]

Cvičení

Uvažujme rovinnou oblast M={(x,y)∈(ℝ^2):f(x,y)<=0}, kde f(x,y)=(((4 (x^2)+4 (y^2)-1))^3)-32(x^2)(y^3). Tato množina je obsažena ve čtverci C=[-1,1]×[-1,1] (lze se o tom přesvědčit pomocí RegionPlot). Vypočítejte přibližnou hodnotu obsahu množiny M následujícím způsobem:

i) Vygenerujte 1 000 000 náhodných bodů ve čtverci C.
ii) Určete, kolik z těchto bodů leží v množině M.
iii) Odhadněte obsah množiny M jako (počet bodů ležících v M)/(počet všech bodů)*(obsah čtverce C).

Změřte, jak dlouho trvá druhý krok a k jakému urychlení dojde, pokud funkci f zkompilujete.

Součástí Mathematiky je slovník DictionaryLookup[] obsahující více než 90 000 anglických slov a dále některé ukázkové texty v angličtině i v jiných jazycích; jejich seznam lze získat pomocí ExampleData["Text"]. 
Načtěte text obsahující Shakespearovy sonety (ShakespearesSonnets), rozdělte jej na seznam jednotlivých slov a najděte v něm všechna slova, která se nevyskytují ve slovníku DictionaryLookup[]. Pokuste se najít co nejrychlejší způsob řešení a změřte, jak dlouho výpočet trvá. Umíte-li úlohu vyřešit více způsoby, porovnejte jejich rychlosti.
Návod. Ukázkový text lze načíst pomocí ExampleData[{"Text","ShakespearesSonnets"}] a poté jej rozdělit na slova pomocí TextWords. Při porovnávání slov v textu a ve slovníku neberte v úvahu velikost písmen (převeďte všechna slova v textu i ve slovníku na malá písmena pomocí ToLowerCase).