Programowanie funkcyjne/Strumienie

Wstęp

Przedstawiając konstrukcje imperatywne, powiedzieliśmy, że w programowaniu imperatywnym upływ czasu w świecie rzeczywistym jest modelowany przez upływ czasu w programie komputerowym. W "czystym" programowaniu funkcyjnym nie ma zmiennych, ani zmian stanów. Trudno też mówić o upływie czasu w programie. Wszak program składa się z implementacji pojęć matematycznych. Jedyna zależność związana z czasem polega na tym, że pojęcia wykorzystywane muszą być zdefiniowane przed pojęciami korzystającymi z nich.

Strumienie to po prostu ciągi wartości, być może nieskończone. Programowanie strumieniowe pozwala jednak na utworzenie związku między czasem i obliczeniami.

W fizyce mamy pojęcie czasoprzestrzeni. Sama czasoprzestrzeń nie zmienia się w czasie -- jest to twór matematyczny ujmujący upływ czasu i historię całej przestrzeni. Historia pojedynczej cząstki jest reprezentowana przez jej trajektorię w czasoprzestrzeni, nazywaną linią świata cząstki.

Jak zobaczymy, strumienie mogą być użyte do przedstawiania historii obliczeń iteracyjnych. Strumień może reprezentować całą, być może nieskończoną historię pewnych obliczeń. Pozwala to "jednym rzutem oka" spojrzeć na całą historię obliczenia. Jak zobaczymy, pozwoli to również na implementowanie pewnych ciekawych operacji na "liniach świata obliczeń".

Wiele z podanych w tym wykładzie przykładów pochodzi z [AS], przy czym zostały przeniesione z języka Scheme do Ocamla.

Definiowanie form specjalnych

Jak już pisaliśmy wcześniej, w Ocamlu argumenty procedur są obliczane gorliwie, tzn. najpierw są obliczane ich wartości, a dopiero potem przystępujemy do obliczania wartości procedury. Wyjątkiem są formy specjalne, np. instrukcja warunkowa if-then-else --- w zależności od wartości warunku, tylko jeden z pozostałych członów jest obliczany.

Czy można takie formy specjalne zaimplementować samemu? Tak. Potrzeba do tego dwóch rzeczy. Musimy umieć odroczyć moment obliczania wartości i wprowadzać własne makrodefinicje.

O uleniwianiu i odraczaniu obliczeń mówiliśmy w poprzednim wykładzie. W Ocamlu możemy nie tylko definiować własne makrodefinicje, ale wręcz zmieniać składnię języka (sic!). Za analizę składniową odpowiedzialny jest preprocesor P4. Nie będziemy tu mówić o nim, gdyż wykracza to poza zakres tego kursu. Jednak raz użyjemy definicji wprowadzającej prostą makrodefinicję, która realizuje formę specjalną.

Implementacja strumieni

Formalnie strumień to ciąg elementów --- być może nieskończony. Pomysł polega jednak na tym, żeby był on "leniwy", tzn. obliczane były tylko potrzebne wartości. Przyjmujemy przy tym, że wartość stają się potrzebne zgodnie z ich kolejnością występowania w strumieniu. Tak więc w każdej chwili strumień składa się ze skończonego ciągu obliczonych (potrzebnych) wartości oraz odroczonego (być może nieskończonego) ciągu pozostałych elementów. Implementując strumienie korzystamy z uleniwiania (zaimplementowanego przez moduł Lazy). Strumień, jeśli nie jest pusty, to para: wartość pierwszego elementu i odroczony strumień pozostałych wartości.

type 'a stream = Nil | Cons of 'a * 'a stream Lazy.t;;
type 'a stream = Nil | Cons of 'a * 'a stream Lazy.t

let empty s = s = Nil;;
val empty : 'a stream -> bool = <fun>

Podstawowe "cegiełki" do budowy strumieni są dokładnie takie same jak do budowy list: pusty strumień oraz dołączenie elementu na początku strumienia. Dwie podstawowe operacje służące do dekonstrukcji strumieni, to procedury zwracające głowę i ogon danego strumienia.

let head s = 
  match s with
    Nil -> failwith "Empty" |
    Cons (h, _) -> h;;

let tail s = 
  match s with
    Nil -> failwith "Empty" |
    Cons (_, t) -> force t;; 

Tak zdefiniowany typ danych jest trochę niewygodny w użyciu, gdyż zawsze gdy używamy Cons, powinniśmy go używać razem z lazy. Aby ułatwić sobie życie, zdefiniujemy Cons jako makro, które przed drugim swoim parametrem dodaje słowo kluczowe lazy.

$\text{<mrow data-mjx-texclass="ORD">Cons</mrow>}(x,y)\equiv \text{<mrow data-mjx-texclass="ORD">Cons</mrow>}(x,\text{<mrow data-mjx-texclass="ORD">lazy</mrow>}y)$

Oto niezbędne "zaklęcie", które podajemy bez wyjaśnień. Zainteresowanych odsyłamy do dokumentacji preprocesora P4.

EXTEND
  Pcaml.expr: LEVEL "simple"
  [ [ UIDENT "Cons"; param=SELF ->
        let h,l =
          match param with <:expr< ($list:[h;l]$) >> -> (h,l)
                           | _ -> raise Not_found
        in
          <:expr< Cons ($h$) (lazy $l$)>> ] ]
  ;
END;;

Mając do dyspozycji takie konstruktory i selektory możemy zdefiniować kilka pomocniczych operacji na strumieniach:

let rec const_stream c = 
  Cons (c, const_stream c);; 

Definiuje nieskończony strumień, którego wszystkie elementy są równe c.

let ones = const_stream 1;;

To jest nieskończony strumień jedynek.

let rec filter_stream p s = 
  if empty s then 
    Nil
  else 
    let h = head s 
    and t = tail s
    in 
      if p h then 
        Cons (h, filter_stream p t)
      else
        filter_stream p t;;

Działa podobnie jak procedura filter dla list, tworzy strumień złożony wyłącznie z tych elementów danego strumienia, które spełniają podany predykat. Jeśli strumień będący parametrem jest nieskończony, to nieskończenie wiele jego elementów powinno spełniać predykat, tak żeby strumień wynikowy był również nieskończony.

let rec stream_ref s n = 
  if n = 0 then 
    head s
  else
    stream_ref (tail s) (n - 1);;

Zwraca element strumienia z podanej pozycji. Pozycje w strumieniu są indeksowane od 0.

let rec stream_map f s = 
  if empty s then 
    Nil
  else
    Cons (f (head s), stream_map f (tail s));;

Działa podobnie jak procedura map dla list, stosuje dane przekształcenie f do wszystkich elementów danego strumienia i tworzy strumień wyników.

Uwaga

W tym wykladzie pokazujemy, jak zaimplementować strumienie, choć Ocaml zawiera moduł (Streams) implementujący strumienie. Jednak są to dwie różne struktury danych. Strumienie zaimplementowane w module Streams odpowiadają bardziej pojęciu strumienia w C++.

Przykłady strumieni nieskończonych

Zdefiniujmy strumień liczb naturalnych:

let rec integers_from x = 
  Cons(x, integers_from (x+1));;

let nats = integers_from 0;;

Podobnie możemy zdefiniować strumień liczb Fibonacciego:

let fibs = 
  let rec fibo a b = 
    Cons (a, fibo b (a+b))
  in fibo 0 1;;

Jak widać, jeżeli jesteśmy w stanie skonstruować iteracyjną procedurę wyliczającą kolejne elementy strumienia, to tym samym jesteśmy w stanie zdefiniować strumień tych wartości. Procedura taka stanowi ukrytą w strumieniu "maszynerię", która na żądanie podaje kolejne elementy strumienia. Takie definicje strumieni nazywamy nie uwikłanymi (w odróżnieniu od uwikłanych, które poznamy dalej). Charakterystyczne dla definicji nie uwikłanych jest to, że bytem definiowanym rekurencyjnie jest procedura, a nie strumień. Strumień jest zdefiniowany dopiero w oparciu o rekurencyjnie zdefiniowaną procedurę.

Strumienie to tylko jeden ze sposobów konstruowania złożonych struktur danych. Można więc go łączyć z dowolnymi innymi sposobami konstruowania danych i tworzyć strumienie list, listy strumieni, czy nawet strumienie strumieni. Oto przykład, strumień strumieni powstałych z danego (nieskończonego) strumienia przez usunięcie coraz większej liczby początkowych elementów (poczynając od zera elementów):

let rec tails s = 
  Cons (s, tails (tail s));;

Strumień liczb pierwszych -- sito Eratostenesa

Spróbujmy skonstruować nieskończony strumień liczb pierwszych, generowanych metodą sita Eratostenesa. Sposób konstrukcji takiego strumienia możemy przedstawić w postaci schematu przypominającego schematy blokowe układów przetwarzających sygnały. Strumienie są na nim zaznaczone liniami ciągłymi, a pojedyncze wartości przerywanymi.

Poniższy diagram przedstawia budowę procedury sito implementującej pojedyncze sito. Zauważmy, że jest on rekurencyjny, tzn. procedura ta pojawia się również na tym schemacie. Odpowiada to temu, że liczby przesiane przez jedno sito trafiają na kolejene sito, itd.

Schemat ten przekłada się na następujący program:

let sitko p s = 
  filter_stream (function x -> (x mod p) <> 0) s;;

let rec sito s  = 
  Cons (head s, sito (sitko (head s) (tail s)));;
 
let primes = 
  sito (integers_from 2);;

Dzięki leniwości strumieni, w zależności od tego ile liczb odczytamy ze strumienia, powstanie odpowiednia liczba sit i odpowiednia ilość liczb zostanie przez nie przesiana.

Definicje uwikłane

Definiując nieskończone strumienie nie musimy tego zawsze robić poprzez podanie odpowiedniej procedury rekurencyjnej (integers-from, fibgen, sieve).Możemy użyć elementów strumienia do zdefiniowania jego samego.

let rec ones =
  Cons (1, ones);;
      

Do budowy bardziej skomplikowanych strumieni potrzebujemy dodatkowych operacji budujących strumienie. Na przykład, do utworzenia strumienia liczb naturalnych czy liczb Fibonacciego potrzebujemy dodawania strumieni.

let rec add_int_streams s1 s2 = 
  if empty s1 or empty s2 then 
    Nil
  else
    Cons (head s1 + head s2, add_int_streams (tail s1) (tail s2));;
  
let rec nats =
  Cons (0, add_int_streams ones nats);;
      

Sposób obliczania kolejnych elementów strumieni zdefiniowanych w sposób uwikłany najlepiej przedstawić sobie w postaci tabelki:

<applet code="PSViewer" archive="images/d/dd/Psviewer.jar" align="center" width="600" height="250"> <param name="DIR" value="images/tab1/"> </applet>

Zamiast dodawania strumienia jedynek, możemy użyć operacji zwiększania o jeden:

let succ x = x + 1;;
let rec nats = 
  Cons(0, stream_map succ nats);;

<applet code="PSViewer" archive="images/d/dd/Psviewer.jar" align="center" width="600" height="250"> <param name="DIR" value="images/tab2/"> </applet>

let rec fibs = 
  Cons (0, Cons (1, add_int_streams fibs (tail fibs)));;
  
let tails s = 
  let rec w = 
    Cons (s, stream_map tail w)
  in w;;

<applet code="PSViewer" archive="images/d/dd/Psviewer.jar" align="center" width="600" height="250"> <param name="DIR" value="images/tab3/"> </applet>

Możemy też w uwikłany sposób zdefiniować strumień liczb pierwszych. Użyjemy do tego predykatu sprawdzającego, czy liczba jest pierwsza:

let rec primes = 
  Cons (2, filter_stream prime (integers_from 3))
      

Natomiast predykat prime zdefiniujemy używając ... strumienia liczb pierwszych:

and prime n = 
  let rec iter ps = 
    if square (head ps) > n then 
      true
    else if divisible n (head ps) then false
    else iter (tail ps)
  in 
    iter primes;;
     

Całość działa poprawnie, ponieważ strumień jest leniwą strukturą danych. Pierwszy element strumienia liczb pierwszych, 2, jest dany explicite. Natomiast sprawdzając kolejne liczby, czy są pierwsze, zawsze mamy już obliczone potrzebne liczby pierwsze.

Konkretnie, największa obliczona liczba pierwsza $p$ wystarcza do przetestowania pierwszości liczb do $p^2$ włącznie. Natomiast kolejna liczba pierwsza $p^{\prime }$ spełnia $p^{\prime }<2p\le p^2$.

<applet code="PSViewer" archive="images/d/dd/Psviewer.jar" align="center" width="600" height="250"> <param name="DIR" value="images/tab4/"> </applet>

Iteracje jako strumienie

W przypadku definicji nie uwikłanych używaliśmy procedur iteracyjnych do zdefiniowania strumieni. Możemy ten mechanizm odwrócić i użyć strumieni do opisania procesów iteracyjnych - kolejne elementy strumienia mogą reprezentować kolejne kroki iteracji. Oto przykład, strumień kolejnych przybliżeń pierwiastka, metodą Newtona:

let sqrt_improve x g = (g +. x /. g) /. 2.0;;
  
let sqrt_stream x = 
  let rec guesses = 
    Cons (1.0, 
      stream_map (sqrt_improve x) guesses)
  in guesses;;
     

Spróbujmy przybliżyć liczbę $\pi$. Użyjemy do tego celu szeregu zbieżnego do $\frac{\pi }{4}$:

$\frac{\pi }{4}=1-\frac{1}{3}+\frac{1}{5}-\frac{1}{7}+\cdots$

let scale_stream c s = 
  stream_map (function x -> x *. c) s;;
  
let pi_summands = 
  let succ x = 
    if x > 0.0 then 
      -.x -. 2.0
    else
      -.x +. 2.0
  in
    let rec s = 
      Cons (1.0, stream_map succ s)
    in
      stream_map (fun x -> 1.0 /. x) s;;
  
let partial_sums s = 
  let rec ps = 
      Cons (head s, add_float_streams (tail s) ps)
  in ps;;
     
let pi_stream = scale_stream 4.0 (partial_sums (pi_summands 1.0));;
      

Taki strumień jest co prawda zbieżny, ale bardzo wolno. Po zsumowaniu pierwszych 1000 elementów ustalone są dopiero trzy pierwsze cyfry: $3.14$. Jest to dosyć oczywiste, jeżeli zauważymy, że suma pierwszych $n$ elementów jest obarczona błędem rzędu $\frac{1}{n}$

Można "przyspieszyć" zbieżność tego szeregu stosując "akcelerator Eulera". Jeśli $S_n$ jest sumą pierwszych $n$ wyrazów szeregu, to szereg zakcelerowany ma postać.

$S_{n+1}-\frac{(S_{n+1}-S_n{)}^2}{S_{n-1}-2S_n+S_{n+1}}$

Działa on dobrze dla ciągów o malejących modułach błędów przybliżenia.

Dlaczego on działa? Przedstawmy kolejne sumy częściowe szeregu w postaci granica plus błąd: $S_n=x+e_n$, $|e_n|>|e_n+1|>0$. Wówczas elementy przyspieszonego szeregu mają postać:

$\begin{array}{c}x+e_{n+1}-\frac{(x+e_{n+1}-x-e_n{)}^2}{x+e_{n-1}-2x-2e_n+x+e_{n+1}}=x+e_{n+1}-\frac{(e_{n+1}-e_n{)}^2}{e_{n-1}-2e_n+e_{n+1}}=\\ =x+\frac{e_{n+1}{e}_{n-1}-2e_{n+1}{e}_n+e_{n+1}^2-e_{n+1}^2+2e_n{e}_{n+1}-e_n^2}{e_{n-1}-2e_n+e_{n+1}}=x+\frac{e_{n+1}{e}_{n-1}-e_n^2}{e_{n-1}-2e_n+e_{n+1}}\end{array}$

Jeżeli np. $|e_n|=e\cdot c^n$ i znaki $e_n$ są takie same lub naprzemienne, to $e_{n+1}{e}_{n-1}-e_n^2=0$, czyli ciąg natychmiast osiąga granicę.

.... sprawdzić, kiedy można strumień wiele razy przyspieszać ...

Oto implementacja akceleratora Eulera:

let euler_transform st =
  let transform s = 
    let s0 = stream_ref s 0
    and s1 = stream_ref s 1
    and s2 = stream_ref s 2
    in (s2 -. square (s2 -. s1) /. (s0 -. 2.0 *. s1 +. s2))
  in
    stream_map transform (tails st);;
      

Taki strumień jest już zbieżny w rozsądnym czasie. Przeliczenie 1000 elementów przyspieszonego szeregu daje przybliżenie $3.141592653$. Jeszcze lepsze wyniki daje wielokrotne przyspieszanie. Skonstruujmy strumień kolejnych przyspieszeń strumienia sum częściowych i wybierzemy z niego strumień pierwszych elementów:

let rec pi_table = 
  Cons (pi_stream, stream_map euler_transform pi_table);;
  
let pi_stream_acc = stream_map head pi_table;;
  
let pi = stream_ref pi_stream_acc 8;;
      

Taki strumień już w 8 elemencie osiąga granice precyzji arytmetycznej.

Przykład: Trójkąt Pascala można sobie przedstawić jako strumień strumieni liczb całkowitych. Zdefiniuj w sposób uwikłany taki strumień.

<applet code="PSViewer" archive="images/d/dd/Psviewer.jar" align="center" width="600" height="300"> <param name="DIR" value="images/tab5/"> </applet>

Strumienie par - liczby Ramanujana

Przedstawimy teraz jako przykład operowania strumieniami, konstrukcję strumienia liczb Ramanujana. Liczby Ramanujana to takie liczby naturalne, które można przedstawić jako sumy sześcianów dwóch liczb naturalnych na dwa różne sposoby. Pierwsze z tych liczb, to: 1729, 4104, 13832, 20683. Niech waga pary liczb naturalnych to będzie suma sześcianów jej składowych. Pomysł konstrukcji strumienia liczb Ramanujana polega na skonstruowaniu strumienia (nieuporządkowanych) par liczb naturalnych uporządkowanego wg. niemalejących wag, wychwyceniu par o powtarzających się wagach i wypisaniu tych wag.

Następująca procedura łączy dwa strumienie uporządkowane wg. niemalejących wag elementów w jeden uporządkowany strumień. Waga jest określona w postaci funkcji w

let rec merge_weighted w s1 s2 = 
  if empty s1 then s2 else
  if empty s2 then s1 else 
    let h1 = head s1 
    and h2 = head s2
    in 
      if w h1 < w h2 then 
        Cons (h1, merge_weighted w (tail s1) s2)
      else if w h1 > w h2 then 
        Cons (h2, merge_weighted w s1 (tail s2))
      else if h1 = h2 then 
        Cons (h1, merge_weighted w (tail s1) (tail s2))
      else 
        Cons (h1, Cons (h2, merge_weighted w (tail s1) (tail s2)));;
      

Następujące definicje określają wagi par oraz funkcję tworzącą strumień par uporządkowanych wg. wag.

let cube x = x * x * x;;
  
let weight (x, y) = cube x + cube y;;
   
let rec pairs s = 
  Cons ((head s, head s), 
        merge_weighted 
          weight
          (stream_map (function x -> (head s, x)) (tail s))
          (pairs (tail s)));;
      

Następująca procedura pozostawia jedynie reprezentantów sekwencji elementów o powtarzających się wagach.

let rec non_uniq w s = 
  let rec skip we s =
    if empty s then Nil 
    else if we = w (head s) then skip we (tail s)
    else s
  in
    'if empty s then Nil 
    else if empty (tail s) then Nil 
    else
      let h1 = head s
      and h2 = head (tail s)
      in 
        if w h1 = w h2 then 
          Cons (w h1, non_uniq w (skip (w h1) s))
        'else
          non_uniq w (tail s);;
     

Strumień liczb Ramanujana możemy zdefiniować następująco:

let ramanujan = 
  in non_uniq weight (pairs nats);;

Co jeszcze ...

Strumienie jako reprezentacja danych, np. szeregów potęgowych. Operacje na funkcjach reprezentowanych w taki sposób.

Laboratoria i ćwiczenia

Laboratoria i ćwiczenia do wykładu