5.Bardziej zaawansowane tematy.

Programowanie rekursywne i dynamiczne

Programy rekursywne to programy które "odwołują się do siebie samych".  W Mathematice bardzo łatwo jest programować rekursywnie choć w przeciweństwie do języków programistycznych takich jak LISP (do których język Mathematiki ma pewne podobieństwo) Mathematica nie jest zoptymalizowana do tego celu i nie uważnie napisane programy rekursywne mogą być bardzo powolne.  Jednym ze sposobów przyspieszenia ich jest  metoda zwana "programowaniem dynamicznym",  lub, bardziej precyzyjnie, metoda "funkcji które zapamiętują swoje wartości".   Jest to bardzo pożyteczna metoda którą najłatwiej zrozumieć rozpatrując przykłady jej użycia.

Szybkie obliczanie liczb Fibonacciego.

Liczby Fibonacciego są rozwiązaniami równania rekursywnego

({{a(n)==a(n-2)+a(n-1)}, {a(0)==0}, {a(1)==1}})

Funkcja RSolve rozwiązuje wiele równań tego typu:

a[n]/.RSolve[{a[n]==a[n-1]+a[n-2],a[0]==0,a[1]==1},a[n],n][[1]]

Fibonacci[n] to po prostu zakodowana n-ta liczba Fibonacciego. Aby zobaczyć definicję użwyamy FunctionExpand:

FunctionExpand[Fibonacci[n]]

Teraz spróbujemy sami zdefiniować liczby Fibonacciego. Zrobimy to na dwa sposoby. Pierwszy, to "zwyczajna rekursja".

Fib1[0]=0;Fib1[1]=1;

Fib1[n_]:=Fib1[n-1]+Fib1[n-2]

Timing[Fib1[30]]

Niestety, to podejście jest bardzo powolne i już  obliczenie 30tej liczby Fibonacciego zabiera zauważalną ilość czasu. Następnie użyjemy "programowania dynamicznego". Na pierwszy rzut oka definicja wygląda dziwnie.  Należy zwrócić uwagę na równoczesne użycie := i =.

Clear[Fib2]

Fib2[0]=0;Fib2[1]=1;

Fib2[n_]:=Fib2[n]=Fib2[n-1]+Fib2[n-2]

Timing[Fib2[30]]

Tym razem działa to błyskawicznie.

Co się dzieje podczas wykonywanie tych dwóch definicji-programów można prześledzić za pomocą funkcji Trace. Rożnicę widać wyraźnie:

Trace[Fib1[5]]

Widzimy że Fib1 wielkorotnie powtarza te same obliczenia.

W przypadku Fib2 pierwsze oblicznie i następne są zupełnie inne. Najpierw usuńmy jesze raz Fib2 i definiujmy funkcję od początku:

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

Trace[Fib2[5]]

Teraz wywołajmy Fib2 jeszcze raz, tym razem dla n=6:

Trace[Fib2[6]]

Widzimy wyraźnie że Mathematica zapamiętała wartości Fib2[n] dla n≤5 i więcej nie musiała ich obliczać.

Możemy także sprawdzić bezpośrednio co Matheamtica wie o Fib1 i Fib2

?Fib1

?Fib2

Oczywiście to duży zysk w szybkości w obliczeniach dzięki użyciu programowania dynamicznego ma pewien koszt w użytej pamięci. Pamięć możemy odzyskać przez:

Clear[Fib2]

Imperatywne (proceduralne) i funkcyjne  programowanie

Dotychczas rozważaliśmy dwa paradygmaty programowania używane w  Mathematice: programowanie oparte na regułach oraz funkcyjne.  Teraz przejdziemy do programowania proceduralnego zwanego także imperatywnym. Jest to oczywiście najbardziej powszechny paradygmat programowania, odający dosyć wiernie sposób działania większości współczesnych komputerów.  Typowe dla tego paradygmatu jest przypisywanie wartości zmiennym i używanie pętli w celu zmiany tych wartości.  Zaczniejmy od przykładu.

Clear[x]

x=1;x=x+1;x=x+1;x=x+1;x=x+1;x=x+1;x=x+1;

Zauważmy że nie otrzymaliśmy żadengo wyniku. Tymniemniej x ma oczekiwaną wartość.

x

Przyjrzyjmy się formie wewnętrznej (FullForm) obliczanego wyrażenia:

FullForm[Hold[x=1;x=x+1;x=x+1;x=x+1;x=x+1;x=x+1;x=x+1;]]

Zauważmy że Mathematica automatycznie daje nam wartość ostatniego wyrażenia w ComposedExpression. Jeśli więc nie chcemy tego wyrażenia otrzymać, używamy ostatniego argumentu Null. 
W wielu proceduralnych językach wartości zmiennej x nie zwracana  o ile nie użyje się  Return[x] lub Print[x]. Ponieważ Mathematica zawsze zwraca wartość obliczanego wyrażenia, Return[x] służy do zupełnie innego celu. Jeśli na końcu pętli w programie napisanym w Mathematice widzimy Return[x] (zamiast prostego x) możemy być prawie pewni że autor nauczył się programować w języku Fortran (lub C) i nadal pisze
programy w tym języku. Nie jest to najbardziej wydajne i eleganckie podejście do programowania w Mathematice ale działa! Teraz powtórzymy to samo używając pętli Do. Zauważmy że sama pętla Do nic nie zwraca - dla tego na końcu kodu umieściliśmy x:

Clear[x]

x=1;Do[x=x+1,{6}];x

To samo można zapisać  nieco krócej:

x=1;Do[x++,{6}];x

Mathematica ma także inne pętle, np. While i For. Ogólnie lepiej starać się ich unikać ponieważ zazwyczaj można uzyskać w Mathematice znacznie lepsze wyniki pod względem szybkości używając funkcyjnych konstrukcji Nest, Fold, FixedPoint i Accumulate.  Jest to, oczywiście, zwiazane z naturą Mathematiki, nie z programowaniem w ogólności. Nawet w Mathematice różnica w szybkości programów napisanych proceduralnie i funkcyjnie znika w sytuacjach w których udaje się użyć tzw. kompilacji.

Zmienne Lokalne

Jeśli w programie przypisujemy zmiennym wartości i zmieniamy je,  musimy zwracać szczegułną uwagę aby przypdakowo nie zmienić wartości zmiennych które chcialiśmy zachować.  Najprostszym sposobem zabezpieczenia się przed taką ewentualnością jest używanie lokalnych zmiennych.  Mathematica ma trzy podstawowe konstrukcje lokalizujące zmienne: Block, Module and With. Każda z nich działa inaczej i ma swoje silne i słabe strony.

?Block

Najpierw prosty przykład lokalizacji przez Block. Przypisujemy zmiennej x wartość 3. Wewnątrz Block zmieniamy jej wartość na 1. Widzimy że "na zewnątrz" wartość x nie uległa zmianie.

x=3;Block[{x},x=1]

x

Dokładnie tak samo zachowa się Module

?Module

x=3;Module[{x,y,z=1},x=5;y=x+z]

x

Block and Module działają zupełnie inaczej. Kiedy lokalizujemy zmienną używając Block jej wartość jest najpierw zapisana, potem tymczasowo odebrana wraz z wszystkimi atrybutami zmiennej, i na koniec, po wyjściu z  Block oryginalna wartość zostaje przywrócona zmiennej. W przypadku Module nazwy zmiennych są zmienione wewnątrz Module,  tak że nie konfliktują z zewnętrznymi. W przeciwieńswie do Module i Block,   wszystkie lokalne zmienne w With muszą mieć przypisane wartości początkowe:

With[{x=1,y=2},x+y]

Następnie, zauważmy jedną cechę która odróżnia Module i With od Block. W tych pierwszych inicializacja (przypisanie wartości początkowych) dla wszystkich zmiennych odbywa się niezależnie:

ClearAll[f]

g[x_List]:=Module[{u=Length[x],v=u+1},v]

g[{1,2,3}]

ClearAll[g]

g[x_List]:=With[{u=Length[x],v=u+1},v]

g[{1,2,3}]

Block zachowuje się inaczej:

g[x_List]:=Block[{u=Length[x],v=u+1},v]

g[{1,2,3}]

A oto inny ważny przykład ilustrujący różnicę między Block i Module:

foo:=x

x=1;foo

Block[{x=2},foo]

Chociaż  foo było zdefiniowane na zewnątrz struktury Block, jego wartość zmieniała się wewnątrz Block. Module  zachowuje się inaczej:

foo

Module[{x=2},foo]

A więc,  ewaluacja wewnatrz Module zależna jest  tylko od oryginalnej definicji zmiennej czy funkcji, co nie jest prawdą w przypadku Blocku.

Pętle i funkcyja iteracji

Programy napisane według imperatywnego paradygmatu zmieniają wartości przypisane jakiejś zmiennej, po czym, aby otrzymać wartość zmiennej należy ją jawnie ewaluować. Rozważmy prosty przykład używający pęteli Do:

Timing[x=1;Do[x++,{20000}];x]

Teraz zrobimy to samo w paradygmie funkcyjnej.  Programując w tym stylu używamy funkcji które zwracają swoje wartości zamiast zmieniać wartość jakiejś zmiennej.  Zamiast pęteli używamy funkcji których argumentami są funkcje. Przykładem takiej funkcji jest Nest:

?Nest

Nest[f,a,4]

Jest także blisko spokrewniona funkcja NestList

?NestList

NestList[f,a,4]

Zamiast używać pętli Do jak powyżej,  uzyskujemy ten sam wynik programując funkcyjnie:

Nest[#+1&,1,20000]//Timing

Zauważmy dużą różnicę w szybkości.

FoldList[f,a,{b,c,d,e}]

Pierwszy argument FoldList musi być funkcją dwóch argumentów. Poniższe przykłady ilustrują działanie FoldList:

FoldList[Plus,0,{a,b,c,d}]

FoldList[Times,1,{a,b,c,d}]

Zwróćmy uwagę na jeszcze jedną funkcję tego typu:

?Accumulate

Accumulate[{a,b,c,d,e}]

Oczywiście ten sam wynik można uzyskać przy pomocy bardziej ogólnej funkcji FoldList:

FoldList[Plus,0,{a,b,c,d,e}]

Jednakże, bardziej wyspecializowana funkcja Accumulate, jest szybsza w działaniu na numerycznych argumentach:

ls=RandomInteger[{1,100},{1000}];

a=(Accumulate[ls];//Timing)

b=(Rest[FoldList[Plus,0,ls]];//Timing)

First[b]/First[a]

Last[a]==Last[b]

Większa szybkość wyspecializowanych funkcji w porównaniu z bardziej ogólnymi jest ważnym elementem w programowaniu w Mathematice.

MouseOver

Block i zmienne globalne.

Block jest najczęściej używany w celu tymczasowej zmiany wartości zmiennych globalnych. Na przykład, globalna zmienna $RecursionLimit ma domyślną wartość:

$RecursionLimit

Oznacza to że pętla w której długość rekursji przekroczy 256 kroków zostanie autmatycznie zatrzymana. Oczywiście ma to ułatwić wyłapywanie błędów programistycznych które bez tego ograniczenia prowadziłyby do nieskończonych pętli. Czasem jednak to ogranicznie jest nie wygodne. Dla przykładu, wróćmy jeszcze raz do  dobrze znanej nam definicji:

Clear[Fib]

Fib[1]=1;Fib[2]=1;Fib[n_]:=Fib[n]=Fib[n-1]+Fib[n-2];

Spróbyjmy

Fib[3000]

Ograniczenie zmiennej  $RecursionLimit do 256 powoduje że nasz kod nie chce działać. Mogłibyśmy zmienić $RecursionLimit do wiekszej liczby albo całkowicie go usunąć robiąc go równym ∞, ale robiąc to narażamy się na nieprzyjemne konsekwencje w przyszłości. Dużo lepiej jest zmienić tymczasowo wartość $RecursionLimit  za pomocą Block. Zanim jednak to zrobimy musimy usunąć wszystkie definicje Fib bo Mathematica zapamiętała "zatrzymane" wartości.

Clear[Fib]

Fib[1]=1;Fib[2]=1;Fib[n_]:=Fib[n]=Fib[n-1]+Fib[n-2];

Block[{$RecursionLimit=∞},Fib[3000]]

Zauważmy że  globalna wartość $RecursionLimit pozostaje nie zmieniona:

$RecursionLimit

Przykład: Symulacja Ruchu Browna

Jako przykład użycia funkcji FoldList, NestList  i Accumulate podajemy konstrukcję symulacji ruchu Browna (procesu Wienera). W Mathematice 9 konstrukcja staje się dużo łatwiejsza dzięki nowym wbudowanym funkcjom RandomFunction  i WienerProcess. Kod który podajemy poniżej działa w wersjach wyższych niż 6.

Jedna ścieżka

Ścieżkę ruchu Browna (przebytą w czasie 1)  przybliżamy przez losowe błądzenie z n krokami, gdzie każdy krok jest liczbą rzeczywistą o rozkładzie normalnym z średnią 0 i standardowym odchyleniem (Sqrt[(1/n)]). Jest szereg sposobów zaprogramowanie tego w Mathematice. Podajemy trzy

Accumulate  (najszybsza metoda)

BrownianMotion[n_]:=Accumulate[Prepend[RandomReal[NormalDistribution[0,Sqrt[1/n]],{n}],0]]

FoldList (powolniejsza)

BrownianMotion[n_]:=FoldList[Plus,0,RandomReal[NormalDistribution[0,Sqrt[1/n]],{n}]]

NestList (jeszcze powolniejsza)

BrownianMotion[n_]:=NestList[#+RandomReal[NormalDistribution[0,Sqrt[1/n]]]&,0,n]

ListLinePlot[BrownianMotion[2000],DataRange->{0,1}]

Wiele ścieżek

Clear[BrownianMotion]

BrownianMotion[time_,steps_,paths_]:=Transpose[Accumulate[Join[{ConstantArray[0,paths]},Transpose[RandomReal[NormalDistribution[0,Sqrt[time/steps]],{paths,steps}]]]]]

ListLinePlot[BrownianMotion[1,100,10],DataRange->{0,1},PlotRange->All]

Manipulate[BlockRandom[SeedRandom[r];ListLinePlot[BrownianMotion[time,steps,paths],DataRange->{0,time},PlotRange->{-2,2}]],{{steps,100,"number of steps"},10,300,1,Appearance->"Labeled"},{{paths,10,"number of paths"},1,50,1,Appearance->"Labeled"},{{time,1,"time"},0.5,10,Appearance->"Labeled"},{{r,0,""},Button["randomize",r=RandomInteger[2^64-1]]&},SaveDefitition->True,Initialization:>(BrownianMotion[time_,steps_,paths_]:=Transpose[Accumulate[Join[{ConstantArray[0,paths]},RandomReal[NormalDistribution[0,Sqrt[time/steps]],{steps,paths}]]]])]

MouseOver