Sledování a urychlování
výpočtů

Antonín Slavík
MFF 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 Mathematice 10 a novější je součástí jádra funkce MandelbrotSetPlot pro zobrazování Mandelbrotovy množiny.

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

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

Jednoduché řešení: Při dokončení výpočtu každého řádku vypíšeme informaci pomocí Print nebo PrintTemporary.

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

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

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

Výstup nevypadá hezky – lepší je použít 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

Viz Parallel Kernel Configuration a Parallel Kernel Status 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 a pokusí se 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.

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.)

Počítání v síti – více možností komunikace: Lightweight Grid, Cluster Integration, Remote Kernels

Asi nejsnazší je Lightweight Grid:

Na server je potřeba nainstalovat Lightweight Grid Manager (licence na MFF)

Na klientském počítači vyhledáme server pomocí Parallel Kernel Configuration, Lightweight Grid, Discover More Kernels

LaunchKernels[] spouští všechna dostupná jádra i na vzdálených počítačích

Kompilované funkce

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

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 čísly v machine precision nebo s komplexními čísly v machine precision

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+3x^4+2x+1;

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

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

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

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

Timing[Table[g[x],{x,0,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+3x^4+2x+1];

g[1+I]

Počítání reálných funkčních hodnot se mírně zpomalilo:

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

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]]

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}];]

Komunikace s externím programem

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).

Příklad: Implementace funkce pro výpočet Mandelbrotovy množiny v jazyce C.
Místo komplexních čísel počítáme s dvojicemi reálných: Vztah z=(z^2)+c má reálnou složku x=(x^2)-(y^2)+(Subscript[c, x]) a imaginární složku y=2x y+(Subscript[c, y]).

#include "wstp.h"
 
void mandelbrot(double xmin, double xmax, double ymin, double ymax, int xres, int yres)
{
   double x, y, dist, cx, cy, temp;
   int column, row, counts;
     
   WSPutFunction(stdlink, "List", yres);
   for(row = 0; row < yres; row++)
  	{
   	  cy = ymin + (ymax - ymin)*row/(yres - 1.0);
   	  WSPutFunction(stdlink, "List", xres);
   	  for(column = 0; column < xres ; column++)
    		{
     		  x = 0;
     		  y = 0;
     		  cx = xmin + (xmax - xmin)*column/(xres - 1.0);	  
     		  dist = 0;
     		  for(counts = 0; counts < 100 && dist < 4; counts++)
      			{
       			  temp = x*x - y*y + cx;
       			  y = 2*x*y + cy;
       			  x = temp;
       			  dist = x*x + y*y;
       		  }
      		  WSPutInteger(stdlink, counts);
     	  }
     }
 }

int main(int argc, char *argv[])
{ 
 	return WSMain(argc, argv);
}

Zkompilovaný zdrojový kód se nachází v mandelbrot/mandelbrot.exe.

SetDirectory[NotebookDirectory[]];

?Install

Clear[mandelbrot]

Install["mandelbrot/mandelbrot.exe"]

?mandelbrot

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

showMatrix[A]

Technické poznámky

Jak zkompilovat program v C tak, aby s ním Mathematica uměla komunikovat?
(Podrobněji viz WSTP and External Program Communication)

Vyrobíme pomocný soubor mandelbrot.tm (template file), kde je popsáno, které funkce implementované v C mají být dostupné v Mathematice (podrobněji viz Setting Up External Functions to Be Called from the Wolfram Language):

:Begin:
:Function:   mandelbrot
:Pattern: mandelbrot[{{a_?NumberQ,b_?NumberQ},{c_?NumberQ,d_?NumberQ}},{e_?NumberQ,f_?NumberQ}] 
:Arguments:  {N[a],N[b],N[c],N[d],N[e],N[f]}
:ArgumentTypes: {Real, Real, Real, Real, Integer, Integer}
:ReturnType:   Manual
:End:

Zpracujeme mandelbrot.tm programem mprep (viz dokumentaci k mprep), výsledkem je kód v C (např. mandelbrottm.c)

Použijeme vhodný překladač C na soubory mandelbrot.c a mandelbrottm.c k sestavení mandelbrot.exe (navíc je potřeba hlavičkový soubor wstp.h a knihovny dodávané s Mathematicou)

Celý proces lze automatizovat pomocí Makefile souboru.

Vhodné překladače: 
gcc pod Unixem nebo ve Windows (prostředí Cygwin), MS Visual Studio

Automatická kompilace funkcí do C

Funkce je možné automaticky kompilovat do C pomocí Compile a volby CompilationTarget->"C".

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

Needs["CCompilerDriver`"]

CCompilers[]

DefaultCCompiler[]

Clear[iCount,mandelbrot]

iCount=Compile[{{c,_Complex}},Module[{cnt = 0,z=0.0I},While[Abs[z]≤2 && cnt < 100, cnt++; z = z^2 + c]; cnt],CompilationTarget->"C"];
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]

Průlet Mandelbrotovou množinou

Přibližování k bodu ((Subscript[x, 0]),(Subscript[y, 0])):

{(Subscript[x, 0]),(Subscript[y, 0])}={-1.408,0.136};Manipulate[showMatrix[mandelbrot[{{(Subscript[x, 0])-0.8^n,(Subscript[x, 0])+0.8^n},{(Subscript[y, 0])-0.8^n,(Subscript[y, 0])+0.8^n}},{300,300}]],{n,-3,32,.5}]

Efektivita programů – obecné poznámky

Dynamické programování

Vytváření globálních přepisovacích pravidel za běhu programu se nazývá dynamické programování. Typicky jde o definice funkcí v následujícím tvaru:

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

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

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]

Point, Line a Polygon

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

Příklad – atraktor Duffingovy rovnice. Najdeme numerické řešení rovnice pro t∈[0,50000] a chceme zobrazit hodnotu a derivaci pro t=k·2π, k∈{0,1,2,…}.

sol=NDSolveValue[{-x[t]+((x[t])^3)+0.15 x'[t]+x''[t]==0.3 Cos[t],x[0]==-1,x'[0]==1},x,{t,0,50000},MaxSteps->Infinity];
dsol=D[sol[t],t][[0]];
duffingPts=Map[({sol[#],dsol[#]})&,Range[0,50000,2Pi]]

Length[duffingPts]

Porovnání rychlostí zobrazení: Point[{{a,b},{c,d},...}] vs. {Point[{a,b}],Point[{c,d}],...}

pts1=Point[duffingPts]

pts2=Map[Point,duffingPts]

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[duffingPts]-1)]];
Graphics[{PointSize[0.006],Point[duffingPts,VertexColors->colors]},Background->GrayLevel[0.4]]

Clear[sol,dsol,s,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[...],...}.

Přidávání prvků do seznamu

Přidávání prvku na začátek nebo na konec dlouhého seznamu pomocí Append, AppendTo, Prepend, PrependTo je pomalé.

SeedRandom[0];Timing[a = {};sum = 0;While[sum < 10^5,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]]