2011/05/07

Term algebra

Jeremy 在昨天的 AoP meeting 重講 maximum segment sum,最後提出一個 datatype-generic version。(他寫在他最新的 blog post Horner's Rule。)中間我們拐到 T-algebras:令 T 是 monad,那麼我們稱 f : T a → aT-algebra 的意思是它滿足 f . μ = f . T ff . η = id。(μ : TT → Tη : Id → T 是伴隨 monad T 來的兩個 natural transformations,在 Haskell 裡面是 joinreturn 這兩個 polymorphic functions。)這時 Nick 問道:他一直看不出這裡說的 algebra 和我們一般講的 algebra(monoids, groups, rings, fields, ...)有什麼關係。Well asked! Jeremy 於是給了下面這一段 monad 循循善誘版介紹。

我們知道一個 monoid 是某個集合上面定義一個 associative binary operator 並且有一個 unit。比方說,List A 是一個 monoid,用的 associative binary operator 是 concatenation "++",unit 則是 empty list []。它們的 type 是

[]   : List A
(++) : List A × List A → List A
稍微複雜一點的例子像是架在某個 field K 上的 vector space V,這時我們有 vector addition、unit、和 scalar product:
0   : V
(+) : V × V → V
(·) : K × V → V
這些 operators 必須滿足某些額外的條件。但無論那些 operators 的 signature 有多複雜,我們發現它們的 result type 都是那個底層的集合。Category theorists 於是用 sum types 把這些 operators 收集起來,list monoid 的話是
1 + List A × List A  →  List A
vector space 則是
1 + V × V + K × V  →  V
最後我們把箭號左邊的東西抽象成某個 "signature" functor F,上述兩個 type 就都能寫成 F X → X 的型式,其中 X 是底層的集合(List AV)。(對於 list monoid 我們讓 F X := 1 + X × X,vector space 則是 F X := 1 + X × X + K × X。)於是 F X → X 這種型式的 arrows 我們就叫它做 F-algebra,只要適當定義 F 就能表現各式各樣的 operators,至於這些 operators 應該滿足的條件就還要另外陳述。以 monoids 為例,令 f : 1 + X × X → X,裡面藏的那個 binary operator (++) = f . inr 的 associativity 是
(x ++ y) ++ z = x ++ (y ++ z)
這可以繼續改寫成 point-free 型式:
(x ++ y) ++ z = x ++ (y ++ z)
≡ ((++) . ((++) × id)) ((x, y), z) = ((++) . (id × (++))) (x, (y, z))
≡   { 令 assoc ((x, y), z) := (x, (y, z)) }
  ((++) . ((++) × id)) ((x, y), z) = ((++) . (id × (++)) . assoc) ((x, y), z)
≡   { extensionality }
  (++) . ((++) × id) = (++) . (id × (++)) . assoc
願意的話可以繼續畫成 commutative diagrams 之類的。如果描述的性質足夠泛化,甚至可以寫成「不須依賴 F 實際定義」的型式。

所以給一個 F-algebra f : F A → A,小學生可以用它列出算式然後求算 normal form。但國中生開始處理未知數的時候,事情就變得複雜一些,因為「列算式」和「求值」這兩個階段更明確地分開了。假設我們考慮整數加法和乘法,小學生看到的算式都是立刻可以求算出值的,在 Haskell 裡面就相當於用 +* 直接列式,比方說 3 + 2 * 5,那直接是個 Int,小學生做的事情只是把它化簡成 normal form 而已。但國中生看到的式子含有未知數,那些式子沒辦法求算得一個整數值,所以更精確的說法是他們先定義一個 datatype

data Expr X = Var X | Add (Expr X) (Expr X) | Mult (Expr X) (Expr X)
其中 X 是未知數 x, y, z, ... 這些符號的集合(請忽略 Haskell 的 naming conventions XD),然後用 Var, Add, Mult 這些 constructors 列式,如 Var 'x' `Add` (Var 'y' `Mult` Var 'z'),另外再有一個求值函式把 X 上的賦值擴充到 Expr X 上:
eval : (X → Int) → Expr X → Int
eval σ (Var x)    = σ x
eval σ (Add a b)  = eval σ a + eval σ b
eval σ (Mult a b) = eval σ a * eval σ b
如果我們要把加法和乘法這兩個 operators 表示成 F-algebra,那麼我們會定義 F Y := Y × Y + Y × Y。顯然 Expr XF 是有關聯的 — 前者乃導自後者,Expr X 其實是 μY. X + F Y!一個算式 Expr X 可能是個未知數,或是一個 F (Expr X),也就是某一個 operator 下繼續裝更多算式。現在有兩件事情應該滿容易能接受:首先,給一個未知數 x : X,我們可以把它變成一個算式 Var x : Expr X。再來,Expr XX 可以是任意的集合,也就是說我們所謂的未知數其實可以是各種奇怪的東西,就像高年級小學生可以寫 "□ + ▵ * ◯" 一樣。而其中一種我們可以當作未知數的東西就是 Expr X,像 "[x + y] + [y * z] * [z + x]",方括號的意思是說括起來的部份其實我們看作是一個「未知數」,這個算式的 type 是 Expr (Expr X)。但絕大部份人看到這個算式都會覺得我們只是寫一個單純的算式,裡面有未知數 x, y, z,也就是說,他們看到的算式型別是 Expr X。因此給一個 Expr (Expr X),我們可以把中括號抹去,讓它變成一個 Expr X。這兩件事情正是 monad 附帶的那兩個 natural transformations return : X → T Xjoin : T (T X) → T X,它們互動的時候會滿足某些看起來很自然的條件,比方說把一個算式看作是未知數(加上中括號)再把中括號抹掉會得到原來的算式(join . return = id)之類的。所以我們剛得知 Expr 是個 monad。再想一想,「抹括號」這件事情其實和 substitution 息息相關。如果我們有一個 Expr Xx + y * z,然後對於每個未知數我們都指定一個要取代它的 Expr Y,也就是說我們有一個 function X → Expr Y,那麼代換後我們會得到 [u + v] + [v * w] * [w + u] 這種東西,其中 u, v, wY 裡的未知數,算式型別是 Expr (Expr Y),抹去括號後我們就得到一個 Expr Y。整個過程的型別是 Expr X → (X → Expr Y) → Expr Y,也就是 bind operator。

List monad 用算式觀點來看應該十分自然(因為正好只是 free monoid),比方說 concat : List (List A) → List A 確實是把第二層的中括號抹掉。IO monad 我相信一樣可以用算式來看,只是它的 operators 複雜一點。想法大概是 IO a 只是一段程式碼(syntax tree — 所以也是某種算式),我們寫 Haskell 程式去產生那些程式碼。至於程式碼的執行是高一層的事情,就好像外界有個 eval 在跑一樣。Expreval 只是把 Add 翻譯成 +,但 IO 的 eval 翻譯出來的東西比較複雜,得把某個 operand 的結果餵給別的 operand 讓那個 operand 算出更多程式碼,所以情況會像是跑一段程式碼把結果丟給 Haskell 程式產出更多程式碼再繼續跑一樣。Haskell 程式負責的永遠只是產出程式碼,那當然都是 pure data/computations 嘍。

--
先還個債,以後想到再寫清楚一點⋯ XD

Labels: ,

2010/07/17

Concurrent Haskell

最近做的東西是 encoding session types in Agda,我這邊目前最終的希望能把 Agda terms 編譯成 Concurrent Haskell 程式。這樣的話,總得要先會寫 Concurrent Haskell 吧!以下是我設想應該要翻譯出來的程式(的加糖版本):

 1: {-# OPTIONS -XScopedTypeVariables #-}
 2: 
 3: import Control.Concurrent
 4:   (forkIO, MVar, newEmptyMVar, putMVar, takeMVar)
 5: import Data.Function
 6: import Unsafe.Coerce
 7: 
 8: -- untyped synchronous channels
 9: 
10: data Chan = Chan (MVar ()) (MVar ())
11: 
12: newChan :: IO Chan
13: newChan =
14:   do d <- newEmptyMVar
15:      a <- newEmptyMVar
16:      putMVar a ()
17:      return (Chan d a)
18: 
19: readChan :: Chan -> IO a
20: readChan (Chan d a) =
21:   do v <- takeMVar (unsafeCoerce d)
22:      putMVar a ()
23:      return v
24: 
25: writeChan :: Chan -> a -> IO ()
26: writeChan (Chan d a) v =
27:   do takeMVar a
28:      putMVar d (unsafeCoerce v)
29: 
30: -- service (channel) provider
31: 
32: data Service = Service (MVar (MVar Chan)) (MVar (MVar Chan))
33: 
34: accept :: Service -> IO Chan
35: accept (Service acc req) =
36:   do m <- newEmptyMVar
37:      putMVar acc m
38:      c <- takeMVar m
39:      return c
40: 
41: request :: Service -> IO Chan
42: request (Service acc req) =
43:   do m <- newEmptyMVar
44:      putMVar req m
45:      c <- takeMVar m
46:      return c
47: 
48: newService :: IO Service
49: newService =
50:   do acc <- newEmptyMVar
51:      req <- newEmptyMVar
52:      forkIO $ fix (\s ->
53:        do m1 <- takeMVar acc
54:           m2 <- takeMVar req
55:           c <- newChan
56:           putMVar m1 c
57:           putMVar m2 c
58:           s)
59:      return (Service acc req)
60: 
61: -- test program
62: 
63: server :: Service -> IO ()
64: server s =
65:   do c1 <- accept s
66:      n :: Integer <- readChan c1
67:      c2 :: Chan <- newChan
68:      writeChan c1 c2
69:      xs :: [Integer] <-
70:        sequence (take (fromInteger n) (repeat (readChan c2)))  -- simplified
71:      writeChan c1 (sum xs)
72:      return ()
73: 
74: client :: Service -> IO Integer
75: client s =
76:   do c1 <- request s
77:      writeChan c1 (2 :: Integer)
78:      c2 :: Chan <- readChan c1
79:      writeChan c2 (512 :: Integer)
80:      writeChan c2 (256 :: Integer)
81:      v :: Integer <- readChan c1
82:      return v
83: 
84: main :: IO ()
85: main =
86:   do s <- newService
87:      forkIO (server s)
88:      v <- client s
89:      print v

不過這是剛好可以動的版本,如果有兩組以上的 server/client 就會當掉。持續實驗中⋯

--
因為在 Agda 已經過了 typecheck,翻譯成 Haskell 理論上就不用堅持 channels 必須是 typed,hence the uses of unsafeCoerce。


去辦台胞證,捷運上想到這樣做其實還弄不出 synchronous channels⋯ 寫的人可能一寫完離開就又進去自己把東西拿走了,但正確行為應該要卡住,等讀的人拿走才離開。


改成這樣似乎可以:

-- untyped synchronous channels

data Chan = Chan (MVar ()) (MVar ())

newChan :: IO Chan
newChan =
  do d <- newEmptyMVar
     a <- newEmptyMVar
     return (Chan d a)

readChan :: Chan -> IO a
readChan (Chan d a) =
  do v <- takeMVar (unsafeCoerce d)
     putMVar a ()
     return v

writeChan :: Chan -> a -> IO ()
writeChan (Chan d a) v = 
  do putMVar d (unsafeCoerce v)
     takeMVar a

另外「main thread 一停、整個程式就停」的行為很討厭,所以我寫了一個 "Join" 等那些分岔出去的 threads 結束。

newtype Join = Join (MVar [MVar ()])

newJoin :: IO Join
newJoin = newMVar [] >>= return . Join

join :: Join -> IO ()
join (Join j) = takeMVar j >>= mapM_ takeMVar >> putMVar j []

forkIO' :: Join -> IO () -> IO ThreadId
forkIO' (Join j) p =
  do m <- newEmptyMVar
     modifyMVar_ j (return . (m:))
     forkIO (p >> putMVar m ())

現在的 main program 就可以寫

main :: IO ()
main =
  do j <- newJoin
     s <- newService
     forkIO' j (server s)
     forkIO' j (client s)
     forkIO' j (client s)
     forkIO' j (client s)
     forkIO' j (server s)
     forkIO' j (server s)
     join j

--
這樣的翻譯應該是忠實的?挺難證的樣子,但不證這類東西感覺上就沒有 contribution 啦!

Labels:

2009/05/01

GHCi

用 Haskell 寫程式確實相當順暢!我把 128 個 iterations 的狀態全部輸出到一個檔案(約 5M),稍後要檢視裡面的資料。只要善用 deriving (Show, Read)unsafePerformIO + strict IO,GHCi 立刻就變成十分專業的檢視器:把資料拉進 GHCi 當作一個正常變數,然後就可以用任意的 Haskell expressions 操作、顯示這些資料。方便極了!

--
幸好這次用 Haskell XD。

Labels:

2009/04/23

hmatrix

因為要用 Haskell 寫作業(某種 fuzzy clustering),所以先弄一套算矩陣的 library 備用,挑上的是 hmatrix-0.5.1.1。想不到其實不好裝:在 Mac 上只要裝 GSL 而不用裝 ATLAS,後者讓 PB 白編了一個晚上(而且還編不出來!);然後 ghci 即使加了 -framework vecLib 選項,跑 runTests 20 時仍必跳出 segmentation fault,看來是連不到 LAPACK 的函式,幸好編譯好的程式就沒這問題。

--
Library 本身看起來是挺好用的啦 XD。

Labels:

2008/12/28

農夫過河

解題動機來自 scm 老師在嵐達網上貼的 forum post,摘錄其中的問題描述:

大家都知道這個農夫過河問題:農夫、羊、狼、與包心菜準備過河,河中有一條小船。只有農夫能撐船,而農夫每次只能帶一樣東西過河。當農夫不在時,狼會吃羊,羊會吃包心菜。怎麼讓他們安然過河呢?

雖然不太了解農夫為什麼要帶著狼,但單純的 searching 就可以輕鬆解決這道題目,而 searching 在 Haskell 裡面可以很乾淨地分成兩段,第一段造出 search tree,第二段對這棵 search tree 做 breadth-first search。

先引入 module Data.List。

> import Data.List  -- using (unfoldr)
我們的 search tree 採用 "rose tree":這種樹有一個 root,其下是一群 subtrees,收納在一個 list 裡面。
> data Tree a = Node a [Tree a]  deriving Show
在 rose tree 上一樣可以定義 unfold,拿一個種子和「展開一步」的 function funfoldt 會用 f 把種子變成 root element 和新一列種子,然後遞迴地把那一列種子長成 subtrees。
> unfoldt :: (b -> (a, [b])) -> b -> Tree a
> unfoldt f b = let (a, bs) = f b in Node a (map (unfoldt f) bs)
接下來我們要寫一個函式
> bfs :: Tree a -> [a]
對一棵 rose tree 做 breadth-first traversal。傳統的 breadth-first search 用一個 queue 記錄尚未拜訪的 nodes,每次從這個 queue 拿出一個 node,然後把這個 node(尚未拜訪過)的鄰居放進 queue 裡面,直到 queue 全部清空。bfs 實作的是同一個想法,直接把尚未巡訪過的 (sub-) trees 記錄在一個充當 queue 的 list 裡面,每次拿出一棵樹就輸出它的 root,然後把其下的 list of subtrees 串在 queue 後面。因此 bfs 可以寫成一個 unfoldr,用 queue 當作種子:
> bfs :: Tree a -> [a]
> bfs = unfoldr f . (:[])
>   where
>     f [] = Nothing
>     f (Node a ts : xs) = Just (a, xs ++ ts)

接下來才真正開始解農夫過河問題。我們用一個 4-tuple 代表過河的狀態:

> type State = (Bool, Bool, Bool, Bool)
四個 components 分別是農夫、狼、羊、包心菜的狀態,False 代表還沒過河,True 代表已經過河。我們希望造一棵 Tree State,滿足
  • root 是 (False, False, False, False)
  • 當某個 state s 可經過一步抵達 state t 時,就讓 t 成為 s 的 child,而且
  • 出現在樹中的 state 都沒有慘劇發生,即不會有哪個生物被吞進別人的肚子裡。
實作上為了方便,我們把「從 root 走到某個 state」中間採取的所有行動記錄在一個字串裡面,並直接把這個字串和對應的 state 綁在一起。
> type Trace = String
於是 search tree 就寫成
> stateSpace :: Tree (State, Trace)
> stateSpace = unfoldt f ((False, False, False, False), "") 
>   where
其中 f 拿到一個 state(和對應的 trace)時應該告訴 unfoldt 那個 state 是新的元素,並從那個 state 算出所有可能到達的 feasible states 當作新的種子。所謂 "feasible" 就是「沒有慘劇發生」,而「慘劇發生」就是「狼羊同處、農夫卻不在一旁」或「羊菜同處、農夫卻不在一旁」,因此
>     feasible :: State -> Bool
>     feasible (a, w, g, c) = not (w == g && a /= w || g == c && a /= g)
為了方便表達「當農夫和狼同處時,農夫可帶狼過河」,我們寫一個 helper function
>     ccons :: Bool -> a -> [a] -> [a]
>     ccons b x = if b then (x:) else id
ccons 是 "conditional cons" 的縮寫:ccons b x xsbTrue 時是 x : xs,否則是 xs。如果舊的 state 是 (a, w, g, c),一步可抵達的 states 就是
ccons (a == w) (not a, not w,     g,     c) $
ccons (a == g) (not a,     w, not g,     c) $
ccons (a == c) (not a,     w,     g, not c) [ (not a, w, g, c) ]
例如 ccons (a == w) (not a, not w, g, c) 講的是「如果農夫和狼在一塊,就把『農夫和狼一起過河』後的 state 加入 list」。有了這個 list,再用 filter feasible 就可篩出一步可抵達的 feasible states。實際上我們還需要處理 traces,所以 unfoldt 用的 (coalgebra) f 寫成
>     f s @ ((a, w, g, c), tr) = (s, filter (feasible . fst)
>       $ ccons (a == w) ((not a, not w,     g,     c), 'W' : tr)  -- Wolf
>       $ ccons (a == g) ((not a,     w, not g,     c), 'G' : tr)  -- Goat
>       $ ccons (a == c) ((not a,     w,     g, not c), 'C' : tr)  -- Cabbage
>       $ (:[])          ((not a,     w,     g,     c), 'A' : tr)) -- Alone
至此就展開了一棵 search tree。最後我們只要對這棵樹做 BFS,找到那些 (True, True, True, True) states 並收集對應的 traces 就結束了。因為當初記錄的 traces 是反序的,所以要套個 reverse 把順序調整好。
> solutions :: [Trace]
> solutions = [ reverse tr | ((True, True, True, True), tr) <- bfs stateSpace ]
在 ghci 下打 head solutions 就可以看到最短的一組解:
*Main> head solutions
"GAWGCAG"

對於沒看過 Haskell 的人而言,其中一個令人好奇的現象應該是「stateSpace 是一棵 infinite search tree」。在「正常」的語言裡面,bfs stateSpace 會先試圖把 stateSpace 算完,但 Haskell 採用 lazy evaluation,除非真有需要否則不會把值算出來,於是程式執行時實際發生的事情會很類似用「正常」語言寫的 BFS,每次看一個 state 找出一步可抵達的 feasible states 放進 queue 裡面,然後重複。(算是自動做了某種程度的 deforestation 吧?)如果只要 solutions 的第一個元素,找到對應的 state 時程式就會停止,因為不需要算出其餘的解。

--
個人認為像這樣「拆成兩部份寫」是非常優雅的寫法 XD。

Labels:

2008/08/03

Regular Expressions and Basic Computability Theory in Haskell

這篇 blogpost 討論的 regular expressions 是用 formal language theory 裡面常見的定義。

定義 令 Σ 是給定的 alphabet,則 regular expressions 是以下列規則建造:
  1. x ∈ Σ 是個 regular expression;
  2. ε 是個 regular expression;
  3. 是個 regular expression;
  4. RS 是 regular expressions,則 R ∪ S 也是個 regular expression;
  5. RS 是 regular expressions,則 RS 也是個 regular expression;
  6. R 是個 regular expression,則 R* 也是個 regular expression。

為了方便,以下簡稱 regular expressions 為 regexes。在 Haskell 裡面我們定義一個新的 datatype RegExp a,其中 a 是 regex 的 alphabet type。

data RegExp a = Literal a
              | Epsilon
              | Empty
              | Union (RegExp a) (RegExp a)
              | Append (RegExp a) (RegExp a)
              | Kleene (RegExp a)
              deriving (Show, Read)
基本上就只是把數學定義照抄一遍而已。

Regular Expression Matching

讓我們照一般對 regex 的直覺寫一個 match function,拿一個 RegExp a[a],傳回一個 Bool 表示那個 [a] 是否匹配那個 RegExp a。這可以輕鬆地透過 Haskell 的 pattern matching 機制完成:對於每種可能的 regex,指定對應的判斷方式。

Literal x 這樣的 regex 只能匹配「恰有一個元素 x 的 list」。其他的 lists 都不對。

match (Literal x) [y] | y == x = True
match (Literal x) _ = False
一個 list xs 成功匹配 Epsilon 的充要條件是「xs 是 empty list」。
match Epsilon xs = null xs
不可能有任何 list 和 Empty 匹配。
match Empty _ = False
想匹配 r ∪ s,只要和其中一個成功匹配就行了。
match (r `Union` s) xs = match r xs || match s xs
先跳到 Kleene star。我們知道 R* = ε ∪ RR*,所以匹配左側 regex 就等於匹配右側 regex。因為 Haskell 的 laziness,這樣等同於測試是否匹配 ε 然後 R 然後 RR 然後 RRR 然後一直下去,發現匹配結果成功或失敗就停止。
match (Kleene r) xs = match (Epsilon `Union` (r `Append` Kleene r)) xs
最後是 r `Append` s 的情況。我們知道一個 list xs 要想匹配這個 regex,必須把 xs 分割成兩段 yszs(使得 xs = ys ++ zs),而且 yszs 分別和 rs 成功匹配,但我們不知道該怎麼分割。幸好有個簡單的做法:試過每種可能的分割就行了。Haskell Prelude 有個函式 splitAt,給定一個長度 n,把一個 list 分割成長度 n 的 prefix 和剩下的部份,形成一個 pair。例如
splitAt 4 [1..10] = ([1,2,3,4],[5,6,7,8,9,10])
對於一個 list xs,所有可能的分割就是
[splitAt 0 xs, splitAt 1 xs, ..., splitAt (length xs) xs]
flip 改寫一下
[flip splitAt xs 0, flip splitAt xs 1, ..., flip splitAt xs (length xs)]
現在變成把 flip splitAt xs 套用到 [0..length xs] 的每個元素,而這就是 map 做的事情:
map (flip splitAt xs) [0..length xs]
到這裡我們就得到一個 list of pairs,是 xs 所有可能的分割方式。接下來對於每個 pair,我們要試著匹配 sublists 和對應的 regex。這再用一次 map 就能達成:
map (\(xs, ys) -> match r xs && match s ys)
或者用 Control.Arrow.*** 可以寫成
map (uncurry (&&) . (match r *** match s))
至此我們會得到一個 [Bool],每一個 Bool 都是一個可能分割的匹配結果。只要有一個匹配成功就行了,所以最後用 or 取這些 Bools 的 disjunction。至此就完成 r `Append` s 情況的定義。
match (r `Append` s) xs =
  or $ map (\(xs, ys) -> match r xs && match s ys) $
  map (flip splitAt xs) [0..length xs]

A Semantic Function

讓 GHCi 告訴我們 match 的型別:

match :: (Eq a) => RegExp a -> [a] -> Bool
它說 alphabet a 的元素必須可以比較相等與否,然後 match 接收一個 RegExp a 和一個 [a],傳回一個 Bool。對這種 function type 表示法稍微熟一點的人就知道這樣的描述是簡化過的。-> 是 right-associative,所以 match 的型別其實是
match :: (Eq a) => RegExp a -> ([a] -> Bool)
也就是說,match 是接收一個 RegExp a 的函式,結果是 [a] -> Bool 型別的函式。而 [a] -> Bool 其實是 set of [a]'s 的一種表述方式,即一個集合的 characteristic (indicator) function。令 p :: [a] -> Bool,若我們採納「p xsTruexs 在集合 S 裡面」的詮釋,那麼 p 就唯一決定了 S(cf. axiom of specification,即對於任意集合 S,{ x | x ∈ S and P(x) } 是個 well-defined set)。因此 match 其實是 RegExp a 的一個 semantic function,把任何的 regex 映射到它的語意,即一個 set of [a]'s。這也正是數學上我們對於 regex 語意的定義:一個 regex 的語意是 Σ* 的子集合。

Another Semantic Function

既然如此,我們何不另外寫一個更直接的 semantic function,用 list 表示集合,直接把 RegExp a 映射到 [[a]] 呢?

enumerate :: RegExp a -> [[a]]
一個 regex 對應的集合可能是無窮集,但 Haskell 的 infinite lists 很好操作,所以不構成(大)問題。接著我們就再次對 regex 每個可能的結構定義對應的語意。前三種情形很簡單:
enumerate (Literal x) = [[x]]
enumerate Epsilon = [[]]
enumerate Empty = []
Kleene star 也用和 match 相同的伎倆。
enumerate (Kleene r) = enumerate (Epsilon `Union` (r `Append` Kleene r))
到了 Union,如果允許集合元素重複出現在 list 裡面的話,我們可能會定義
enumerate (r `Union` s) = enumerate r ++ enumerate s
把表示 set 的兩個 lists 接在一起,就是兩個集合的聯集,看起來很自然。不過這在計算上有個嚴重瑕疵,考慮
enumerate (Kleene (Literal '0') `Union` Literal '1')
這個結果其實等同於
enumerate (Kleene (Literal '0'))
因為 "1"(計算上)永遠不會出現在 list 裡面!用 computability theory 的術語,我們希望 enumerate r 對應的集合是 recursively enumerable,亦即每個元素都會在有窮時間內被列出來。但照剛才的定義,regex 0* ∪ 1 對應的 list 前段是 0* 對應的 infinite list,enumeration 於是就會永遠困在這個 infinite list 裡面走不到盡頭,也就不可能在有窮時間內列出 "1" 了。

一種解決方式是把 finite list 擺在前面,但一般而言我們不知道一個 list 是不是 finite。然而這裡我們需要的只是特例:判斷 enumerate r 是不是 finite,而這可以從 r 的結構得知。

finite (Literal _) = True
finite Epsilon = True
finite Empty = True
finite (Union r s) = finite r && finite s
finite (Append r s) = finite r && finite s
finite (Kleene _) = False
finite rfinite s 有一個成立時,就把有窮的那一個 list 放在前面;如果兩個 lists 都是 infinite,我們就交替走訪這兩個 lists。如此便可達成 enumerability,讓每一個元素在有窮時間內出現。
enumerate (r `Union` s)
  | finite r  = enumerate r ++ enumerate s
  | finite s  = enumerate s ++ enumerate r
  | otherwise = alternate (enumerate r) (enumerate s)
  where alternate (x : xs) (y : ys) = x : y : alternate xs ys
或者其實可以直接交替走訪,如果有一個 list 耗盡了就接另一個 list 的剩餘部份。亦即
enumerate (r `Union` s) = alternate (enumerate r) (enumerate s)
  where alternate (x : xs) (y : ys) = x : y : alternate xs ys
        alternate xs [] = xs
        alternate [] ys = ys
最後一個情況是 r `Append` s。初步嘗試是直接算 enumerate renumerate s 的 "Cartesian product"
enumerate (r `Append` s) = [xs ++ ys | xs <- enumerate r, ys <- enumerate s]
然而這一樣會使得原本是 recursively enumerable 的集合無法成功列舉。List comprehension 會先從 enumerate r 拿第一個元素出來、稱之為 xs,然後對於 enumerate s 的每個元素 ys 計算 xs ++ ys。如果 enumerate s 是無窮的,我們就永遠無法前進到 enumerate r 的第二個元素(如果有的話)。

enumerate renumerate s 視為兩個座標軸,我們的目標是讓平面上的任一點在有窮時間內出現。這正是有理數的列舉方式,即「有理數集可數」證明的核心。然而這樣的 traversal 在 Haskell 該怎麼寫呢?

觀察下面這張示意圖

走訪順序是紅、橙、黃、綠,一路往外。這些路線對於座標軸的投影就是我們觀察一個 list 的順序。因此觀察 [1..] 這個 list 的順序是

[1,1,2,1,2,3,1,2,3,4,...]
觀察 ['A'..] 這個 list 的順序則是
['A','B','A','C','B','A','D','C','B','A',...]
zip 把這兩個 lists 拉在一起就是我們要的 enumeration。如果兩個 lists 都是 infinite lists,這兩個觀察順序就很單純。對於一個 list xs,第一種觀察順序可以寫成
concat $ map (flip take xs) [1..]
第二種則可以寫成
concat $ map (reverse . flip take xs) [1..]
或是用 unfoldr 寫得(似乎)更有效率一點
concat $ unfoldr f (xs, [])
  where f (y : ys, zs) = Just (y : zs, (ys, y : zs))
對於 finite lists,這兩個觀察順序就比較麻煩,特別是當兩個 lists 都有窮的時候。但只要有一個 finite list,我們就可以把它放在內層迴圈而輕鬆達成 enumerability。所以整個 r `Append` s 的情況就寫成
enumerate (r `Append` s)
  | finite s  = [xs ++ ys | xs <- enumerate r, ys <- enumerate s]
  | finite r  = [xs ++ ys | ys <- enumerate s, xs <- enumerate r]
  | otherwise = zigzag (enumerate r) (enumerate s)
  where zigzag xs ys = zipWith (++) (concat $ map (flip take xs) [1..])
                                    (concat $ unfoldr f (ys, []))
        f (y : ys, zs) = Just (y : zs, (ys, y : zs))
至此便完成 enumerate 的定義。

Matching via Enumeration

有了 enumerate,我們可以定義另一個 match' function:

match' r xs = xs `elem` enumerate r
因為我們剛剛確保 enumerate r 確實列出 r 對應的集合,所以當 xs 匹配 r 時,match' r xs 一定會停。然而當 xs 不匹配 r 而且 enumerate r 是 infinite 時,match' r xs 就不會停了。因此 match' 只證明了 regular expression matching problem 是 semidecidable(只有答案是「對」的時候才保證會停),要證明 decidability 必須靠原本的 match

--
好長的一篇通識文(花了四小時多)XD。


習題 上面說 match 證明了 regex matching problem 是 decidable,其實是錯的。請給一個反例使 match 不會終止。如何加以修改使 match 遇到每種輸入都終止呢?enumerate 有類似的問題嗎?

--
沒有人抓包,只好自己動手 XD。

Labels: ,

2008/08/02

起點

以下描述一個需要 heavy optimisation 的 Haskell 模擬程式。這是 Weijin 最近加入中研院物理所李副所長旗下接到的第一個習題,用電腦模擬某種 Ising model 並作圖。對我而言這是練習 practical Haskell 的大好機會 XD。

一開始定義「晶格」為 n * n 的矩陣,每一格的值可能是 +1 或 -1。這裡直接用 array of Bools 表示。

type Lattice = IOUArray (Int, Int) Bool
可以把 lattice 印成 0 和 1 的矩陣。
showLattice lat n = getElems lat >>= putStr . concat .
  zipWith (flip (++)) (cycle $ replicate (n-1) "" ++ ["\n"]) .
  map (show . fromEnum)
給定一個 lattice,可以計算某一格的能量:把這一格的值和它上下左右四個鄰居分別相乘後加起來取負值。Lattice 的左右、上下邊界是黏在一起的,例如左邊界上的格子的左邊鄰居就是同一列右邊界的格子。物理上叫它作 "periodic boundary condition",拓樸上這個 lattice 就變成一個甜甜圈,日常生活的例子就是不撞牆的貪食蛇 XD。下面這個函式取得某個座標和它的鄰居的值。
neighbors :: Lattice -> Int -> (Int, Int) ->
             IO (Bool, Bool, Bool, Bool, Bool)
neighbors lat n (i, j) =
  do c <- readArray lat (i, j)
     let dec x = let xdec = x - 1 in if xdec <  0 then n - 1 else xdec
         inc x = let xinc = x + 1 in if xinc >= n then 0     else xinc
     u <- readArray lat (dec i, j)
     d <- readArray lat (inc i, j)
     l <- readArray lat (i, dec j)
     r <- readArray lat (i, inc j)
     return (c, u, d, l, r)
然後是計算能量的函式。
energy :: Lattice -> Int -> (Int, Int) -> IO Int
energy lat n (i, j) =
  do (c, u, d, l, r) <- neighbors lat n (i, j)
     let f b1 b2 = if b1 == b2 then 1 else -1
     return $! -(f c u + f c d + f c l + f c r)
可以印出整個 lattice 的平均能量。
showAverage lat n =
  do r <- getBounds lat
     xs <- sequence (map (energy lat n) (range r))
     print (fromIntegral (sum xs) / fromIntegral (n^2))
這個模擬程式主要要做的事情是這樣:隨機選一格,算出它目前的能量和翻轉它(1 和 -1 互換)之後的能量。如果能量差小於(等於)0 就立刻翻轉;如果能量差大於 0 就以某個機率翻轉它。這個機率由某個溫度參數 t 和能量差 e 決定。
flipProb :: Double -> Double -> Double
flipProb t e = exp (-e / t)
下面的函式實作上述行為。
touch :: Lattice -> Int -> Double -> IO ()
touch lat n t =
  do (i, j) <- newStdGen >>= return . choose n
     (nv, e) <- energyDiff lat n (i, j)
     let flip = writeArray lat (i, j) nv
     if e <= 0
        then flip
        else do f <- newStdGen >>= return . fst . random
                if f < flipProb t (fromIntegral e)
                   then flip else return ()
(其實我不知道可不可以這樣仰賴 newStdGen。)其中 energyDiffenergy 的簡單改版。
energyDiff :: Lattice -> Int -> (Int, Int) -> IO (Bool, Int)
energyDiff lat n (i, j) =
  do (c, u, d, l, r) <- neighbors lat n (i, j)
     let f b1 b2 = if b1 == b2 then 1 else -1
         g c = -(f c u + f c d + f c l + f c r)
     return $! let nc = not c in (nc, g nc - g c)
choose 則是一個 pure function。
choose :: RandomGen g => Int -> g -> (Int, Int)
choose n g = let (x, _) = randomR (0, n^2 - 1) g
             in  (x `div` n, x `mod` n)
一個計算單位是 n^2 次 touch,也就是平均每格碰一次。聽說要從 n = 4 算到 n = 64,每個 n 各做兩億個單位,所以計算量驚人… 總之 main 先讀入 n、溫度、和要模擬的單位數,亂數初始化 lattice,每做完一個單位就輸出目前的平均能量,最後輸出 lattice 的狀態。
main = do n <- readLn
          t <- readLn
          cnt <- readLn
          lat <-
            newStdGen >>= newListArray ((0, 0), (n-1, n-1)) . randoms
          sequence_ $ flip map [1..cnt]
            (\i -> do replicateM_ (n^2) (touch lat n t)
                      putStr (show i ++ "\t")
                      showAverage lat n)
          showLattice lat n
目前大概慢 Weijin 的 C++ 版本至少幾百倍甚至千倍,所以這也是探索 Haskell 程式最佳化的大好機會 XD。

--
另外想寫一篇 regex 與基本 computability theory 的通識文章 XD。

Labels:

2008/08/01

好慢

我覺得「把 Haskell 程式寫對」相對容易,可是「把 Haskell 程式寫快」很難… 以下是一個 2-dimensional mutable array of booleans 的型別:

import Data.Array.IO

type Lattice = IOUArray (Int, Int) Bool
因為稍後要頻繁地做 indexing 和 state change(是某種模擬程式),所以用 plain list 應該不是好選擇。我現在想把 n * n 個元素初始化成 random booleans,一種寫法是
{-# OPTIONS -XPatternSignatures #-}
import System.Random

main = do n <- readLn
          lat :: Lattice <-
            newStdGen >>= newListArray ((0, 0), (n-1, n-1)) . randoms
          return ()
可是隨便給個 n = 1000 就要跑好幾秒,如果把 return () 換成 getElems lat >>= print 還會 stack overflow!

--
啊,我不會寫啦!

Labels:

2008/07/31

The Evolution of a Haskell Programmer

這篇〈The Evolution of a Haskell Programmer〉是模仿另一篇很有名的〈The Evolution of a Programmer〉寫的,好幾項都深得我心啊!XD 以下是我對各條的感想:

Freshman Haskell programmer
中規中矩。
Sophomore Haskell programmer, at MIT (studied Scheme as a freshman)
Operators 全部改成 prefix form 有笑點 XD。
Junior Haskell programmer (beginning Peano player)
嗯,我支持 (n+k) pattern!XD
Another junior Haskell programmer (... joined the "Ban n+k patterns"-movement)
這程式碼才叫噁心好不好 XD。
Senior Haskell programmer/Another senior Haskell programmer/Yet another senior Haskell programmer
標準的幾個 fold 技巧。
Memoizing Haskell programmer
唔,我跟 scan 不太熟 XD。
Pointless (ahem) "Points-free" Haskell programmer (studied at Oxford)
耶,pointfree 萬歲!(原來 pointfree 是 Oxford 的標誌喔 XD。)
Iterative Haskell programmer (former Pascal programmer)
For loop 移植得算相當成功 XD。
Iterative one-liner Haskell programmer (former APL and C programmer)
懶得看 XD。
Accumulating Haskell programmer
這個好一點 XD。
Continuation-passing Haskell programmer
XD
Boy Scout Haskell programmer
Y combinator,太炫了 XD。
Combinatory Haskell programmer
我跟 combinatory logic 完全不熟 XD。
List-encoding Haskell programmer
這啥啊?XD
Interpretive Haskell programmer
這太勁爆了,直接讓我想到 Oleg XD。
Static Haskell programmer
我以前在 C++ 可能會玩這種 ad hoc type-level programming XD。
Beginning graduate Haskell programmer
這個我最近才剛好無聊又重寫一次 XD。(其實是下下下面那一個 XD。)
Origamist Haskell programmer (always starts out with the "basic Bird fold")
不寫成 fold 是罪過啊!XD 不過 unfoldr 我喜歡用 Maybe 定的 XD。
Cartesianally-inclined Haskell programmer
喔耶,莫斐生大家庭登場!
Ph.D. Haskell programmer (ate so many bananas that his eyes bugged out, now he needs new lenses!)
最近看到 datatype-generic programming 剛好拿出來複習 XD。Natural numbers 直接用 Maybe functor 定就好啦!XD
Post-doc Haskell programmer
進入未知領域了… Comonad 和 zygomorphism 是啥啊?
Tenured professor (teaching Haskell to freshmen)
scm 老師的中年寫照?XD
我看得懂的附註(像 "ate so many bananas that his eyes bugged out, now he needs new lenses!")就留下來,反之捨棄 XD。

--
今天沒 meeting 真傷心…

Labels:

2008/07/27

Line Numbers

練習用 State monad 寫「為一個檔案加上行號」的程式。想法是這樣:先用 lines 把 input string 切割成 list of lines,然後把每個 string 變成 State monad 之後接成一個超大的 State monad,最後加以核算。可是結果看起來不太令人滿意(以下是用 "runghc lineno.hs < lineno.hs" 指令產生):

 1: import Data.List
 2: import Control.Arrow
 3: import Control.Monad.State
 4: 
 5: main = getContents >>=
 6:          putStr . uncurry ($) . (id *** g) . flip runState 0 .
 7:          foldr (combine . monadise) (return (const "")) . lines
 8:   where g :: Integer -> Integer -> String
 9:         g n = (++ ": ") .
10:               until ((>= 1 + (floor . logBase 10 . fromInteger) n) . length)
11:                     (' ':) . show
12: 
13: monadise str = do State (const () &&& (1+))
14:                   n <- get
15:                   return (\g -> g n ++ str ++ "\n")
16: 
17: combine s1 s2 = do f <- s1
18:                    g <- s2
19:                    return (\h -> f h ++ g h)
無論是程式複雜度和執行速度都比兩年前寫的 Ruby script 還要慢。當然編譯過後會比 Ruby 快啦,可是拿 compiled program 和 interpreted program 比速度完全是作弊。

--
改用 StateT monad transformer 寫看看好了…


用 StateT 沒有比較好寫的樣子,所以只稍微把上面的版本改寫一下,好像比較順眼:

 1: import Control.Monad.State
 2: 
 3: main =
 4:   do xs <- getContents
 5:      let (ys, n) = runState (monadise n (lines xs)) 0
 6:      foldr ((>>) . putStrLn) (return ()) ys
 7: 
 8: monadise n = foldr (combine . makeLine (width n)) (return [])
 9: 
10: width = (1+) . floor . logBase 10 . fromInteger
11: 
12: makeLine w s =
13:   do modify (1+)
14:      n <- get
15:      let lineNoStr = until ((>= w) . length) (' ' :) (show n)
16:      return (lineNoStr ++ ": " ++ s)
17: 
18: combine s1 s2 = do x  <- s1
19:                    xs <- s2
20:                    return (x : xs)

--
我果然還不太會寫 Haskell…


果然是 imperative programs 寫太久了,看到這問題就想用 state。其實用個 zipWith 把每一行和 [1..] 的對應元素拉在一起,再做點轉換就好了。

 1: main =
 2:   do xs <- getContents
 3:      let xss = lines xs
 4:      let toPrefix = (++ ": ") . padding (width (length xss)) ' ' . show
 5:      let ys = zipWith ((++) . toPrefix) [1..] xss
 6:      foldr ((>>) . putStrLn) (return ()) ys
 7: 
 8: width = (1+) . floor . logBase 10 . fromIntegral
 9: 
10: padding w c = until ((>= w) . length) (c :)

--
有更好一點 XD。

Labels:

2008/07/23

Monad Transformer

繼續 Haskell 效能測試,想試試看 memoisation 能不能進一步加快執行速度,但這麼一來就需要結合兩個 State monads ─ 於是就輪到「魔那變形金剛」上場啦!

 1: import Control.Monad.Identity
 2: import Control.Monad.State
 3: import qualified Data.IntMap as IM
 4: 
 5: cycleLength :: Monad m => Int -> StateT (IM.IntMap Int) m Int
 6: cycleLength 1 = return 1
 7: cycleLength n =
 8:   do calculated <- gets (IM.member n)
 9:      if calculated
10:         then do x <- gets (IM.lookup n)
11:                 return (runIdentity x)
12:         else do x <- if even n then cycleLength (n `div` 2)
13:                                else cycleLength (3 * n + 1)
14:                 modify (IM.insert n (x + 1))
15:                 return (x + 1)
16: 
17: maxCycleLength :: Monad m => [Int] -> StateT (IM.IntMap Int) m Int
18: maxCycleLength [] = return 0
19: maxCycleLength (x : xs) = do a <- cycleLength x
20:                              b <- maxCycleLength xs
21:                              return (max a b)
22: 
23: driver :: StateT (IM.IntMap Int) (State [Int]) [Int]
24: driver =
25:   do eof <- lift $ gets null
26:      if eof
27:         then return []
28:         else do i <- lift $ gets head
29:                 lift $ modify tail
30:                 j <- lift $ gets head
31:                 lift $ modify tail
32:                 x <- maxCycleLength [min i j .. max i j]
33:                 xs <- driver
34:                 return (x : xs)
35:                           
36: main = getContents >>= foldr ((>>) . print) (return ()) .
37:                        evalState (evalStateT driver IM.empty) .
38:                        map read . words
寫得已經很像 imperative code 了。可惜這個版本遇上測資 "1 100000" 的時候會 stack overflow XD。

--
是沒辦法做最佳化的關係嗎?

Labels:

分開

承接上次的 Haskell 效能測試,我突然想到應該把 IO code 和 pure code 再分清楚一點,讓 State monad 單純傳回 list of results,最後用個 foldr 把這個 list 轉成一串 IO () 再全部接起來執行。因為 Haskell 的 laziness,這樣寫的執行期表現和上次沒什麼差別。

 1: import Control.Monad.State
 2: 
 3: cycleLength 1 k = k
 4: cycleLength n k | even n    = cycleLength (n `div` 2) $! (k + 1)
 5:                 | otherwise = cycleLength (3 * n + 1) $! (k + 1)
 6: 
 7: driver :: State [Int] [Int]
 8: driver = do eof <- gets null
 9:             if eof then return []
10:                    else do this <- commit
11:                            remains <- driver
12:                            return (this : remains)
13: 
14: commit :: State [Int] Int
15: commit = do i <- nextInt
16:             j <- nextInt
17:             return
18:               (maximum . map (flip cycleLength 1) $ [min i j .. max i j])
19:          where nextInt = State $ \(x : xs) -> (x, xs)
20: 
21: main = getContents >>= foldr ((>>) . print) (return ()) .
22:                        evalState driver . map read . words
驚人的是單單這樣改一下就讓執行時間下降到 0.6--0.7 秒,和 Java 不相上下!現在我的偏心得到 justification 了 XD。

--
而且現在我也覺得 IO monad 和 State monad 相處得比較融洽了 XD。

Labels:

2008/07/21

競速 Reprise

承接兩年前的 C/C++/Java/Ruby 大亂鬥,這次我寫了 Haskell 版本:


 1: import Control.Monad.State
 2: 
 3: cycleLength 1 k = k
 4: cycleLength n k | even n    = cycleLength (n `div` 2) $! (k + 1)
 5:                 | otherwise = cycleLength (3 * n + 1) $! (k + 1)
 6: 
 7: driver :: State [Int] (IO ())
 8: driver = do eof <- gets null
 9:             if eof then return (return ())
10:                    else do this <- commit
11:                            remains <- driver
12:                            return (this >> remains)
13: 
14: commit :: State [Int] (IO ())
15: commit = do i <- nextInt
16:             j <- nextInt
17:             return $ print
18:               (maximum . map (flip cycleLength 1) $ [min i j .. max i j])
19:          where nextInt = State $ \(x : xs) -> (x, xs)
20: 
21: main = getContents >>= evalState driver . map read . words
用 ghc 編譯並開 -O 最佳化,執行時間大約是 1.1 秒,是 C++ 版本的 5 倍多。其實還勉強可以接受啦 XD。(偏心?XD)至於程式碼我覺得有一點點彆扭,總覺得 IO monad 和 State monad 沒有搭配得很好。

--
現在在讀 Real World Haskell 的線上 beta 版

Labels:

2008/07/20

IO Monad

我終於克服對於 (IO) monad 的恐懼,開始能夠寫簡單的 Haskell programs with IO 了 XD。以下是一段 GHCi 的 session,亂數產生長 100 的 list,測試 insertion sort:

Prelude> :module System.Random
Prelude System.Random> let split f g x = (f x, g x)
Prelude System.Random> let cross f g = split (f . fst) (g . snd)
Prelude System.Random> let randomList n g = if n == 0 then ([], g) else let (x, 
g') = next g in cross (x:) id $ randomList (n - 1) g'
Loading package old-locale-1.0.0.0 ... linking ... done.
Loading package old-time-1.0.0.0 ... linking ... done.
Loading package random-1.0.0.0 ... linking ... done.
Prelude System.Random> let insert x xs = case xs of [] -> [x]; y : ys -> if x <= 
y then x : y : ys else y : insert x ys
Prelude System.Random> do {xs <- newStdGen >>= return . map (`mod` 100) . fst . 
randomList 100; print $ foldr insert [] xs}
[0,0,2,6,6,7,7,8,8,9,11,12,13,13,13,13,14,15,17,18,19,19,21,21,23,25,25,25,27,28
,33,33,34,34,35,35,35,35,37,38,39,43,44,44,45,45,46,46,47,47,48,49,50,51,53,55,5
6,56,58,58,60,61,62,62,63,63,64,65,65,66,67,67,69,70,71,73,75,75,77,78,80,81,81,
82,83,83,86,87,89,90,92,92,92,92,94,95,96,96,97,99]

--
馬上要面對的就是 Haskell 的一大堆 libraries XD。


寫成一串好像比較好看 XD:

Prelude System.Random> newStdGen >>= print . foldr insert [] . map (`mod` 100) .
 fst . randomList 100
[0,1,2,5,6,9,9,12,12,12,12,13,13,13,15,17,21,22,23,24,25,28,29,29,30,31,32,33,34
,34,36,36,36,39,42,43,44,44,46,46,48,48,48,51,52,53,53,53,54,54,55,57,58,58,60,6
1,61,62,63,65,67,67,68,68,68,69,70,70,70,70,71,72,74,74,74,74,75,76,76,77,78,80,
80,80,80,82,87,87,87,88,93,93,94,94,94,94,94,95,97,99]

Labels:

2008/04/09

一百萬

突然想測試一下 Haskell 跑得有多快(或是有多慢),所以進 GHCi:

Prelude> [1..] !! 1000000
*** Exception: stack overflow
啥?XD 為了確定不是 interpreter 太嫩,我用 ghc 編譯了我的第一份 native code:
main = putStr (show ([1..] !! 1000000))
執行結果是
Stack space overflow: current size 8388608 bytes.
Use `+RTS -Ksize' to increase it.
所以還真的滿出來了!(!!) 的正常定義應該是
(!!) : [a] -> Int -> a
(x : xs) !! 0       = x
(x : xs) !! (1 + n) = xs !! n
看起來是 tail recursion 呀,為什麼會 stack overflow 勒?是 [1..] 的關係?

--
ghc -C 產出來的 C code 看不懂 XD。

Labels:

2008/04/03

Unfold Fusion

一直沒辦法感受 coinduction 的直覺意義。在豁然貫通的那一個時刻來臨以前,最簡單的方式大概就是用 unfold fusion 嘍,因為 category theory 裡面的 duality 很明確很好操作,只要把 inductive datatypes 那一套的箭頭全部反過來就行了。

第一個 coinductive proof 當然是

map f . iterate f = iterate f . f        (*)
我們需要一些基本定義(in Haskell)。
-- constructor
-- Maybe (A, B) is isomorphic to 1 + A * B
alpha :: Maybe (a, [a]) -> [a]
alpha Nothing       = []
alpha (Just (a, t)) = a : t

-- deconstructor
ahpla :: [a] -> Maybe (a, [a])
ahpla []       = Nothing
ahpla (x : xs) = Just (x, xs)

-- products
split f g a = (f a, g a)
cross f g   = split (f . fst) (g . snd)

-- base functor
listrF f g Nothing       = Nothing
listrF f g (Just (a, t)) = Just (cross f g (a, t))

接著是一些 coinductive lists 的相關定義。

unfoldr :: (b -> Maybe (a, b)) -> b -> [a]
unfoldr f = alpha . listrF id (unfoldr f) . f

-- type functor
mapana :: (a -> b) -> [a] -> [b]
mapana f = unfoldr (listrF f id . ahpla)

iterate :: (a -> a) -> a -> [a]
iterate f = unfoldr (Just . split id f)

為了對稱,我們把 inductive lists 的相關定義重寫一遍。

foldr :: (Maybe (a, b) -> b) -> [a] -> b
foldr f = f . listrF id (foldr f) . ahpla

mapcata :: (a -> b) -> [a] -> [b]
mapcata f = foldr (alpha . listrF f id)
以上其實就是莊老師在 FLOLAC '07 講的一部份東西。真恐怖,FLOLAC 花了快一年時間都還沒消化完 XD。

Inductive lists 有 fusion law

h . foldr f = foldr g    <=    h . f = g . listrF id h
以及雖為特例但很好用的 type functor fusion law
foldr f . mapcata g = foldr (f . listrF g id)
現在把 function composition 的方向顛倒,就直接得到 coinductive lists 的 fusion law 和 type functor fusion law!
unfoldr f . h = unfoldr g    <=    f . h = listrF id h . g

mapana g . unfoldr f = unfoldr (listrF g id . f)

我沒預期要鋪這麼多路 XD。Anyway,(*) 式左邊可直接從 type functor fusion law 得到:

  map f . iterate f
=   { definition of iterate }
  map f . unfoldr (Just . split id f)
=   { type functor fusion }
  unfoldr (listrF f id . Just . split id f)
=   { definition of listF }
  unfoldr (Just . cross f id . split id f)
=   { product functor }
  unfoldr (Just . split f f)
右邊很自然地湊湊看 fusion condition:
  (Just . split id f) . f
=   { split fusion }
  Just . split f (f . f)
=   { product functor }
  Just . cross id f . split f f
=   { definition of listrF }
  listrF id f . Just . split f f
=   { let g = Just . split f f }
  listrF id f . g
所以依 fusion law 就得到
iterate f . f = unfoldr (Just . split f f)
於是 (*) 式兩邊都等於同一個 unfold,證畢。Coinduction 大概就先這樣弄吧 XD。

--
寫完才發現 this post should have been written in English XD。

Labels: ,

2008/03/30

Parser Combinators in Haskell

(以下大部份說明文字亂譯自 scm 老師的原文 XD。)

考慮這種 binary tree:

data Tree = A | B | Bin (Tree, Tree)
     deriving Show
我們可以用比較簡短的語法印出這種樹,把 A 印為 aB 印為 bBin (t, u) 印為 (t,u)。例如 (Bin (A,Bin (B,A)) 就印為 (a,(b,a))。這樣的語法可用以下的 context-free grammar 描述:
S ::= ( S , S ) | a | b
在一般的程式語言之內寫一個 recursive-descent parser 並非易事。但若有了 higher-order functions,我們就可以如下為上述語言建造出一個 parser:
pTree :: Parser Char Tree
pTree =  (lit '(' *> pTree <*> lit ',' *> pTree <* lit ')') `using` Bin
      <|> lit 'a'                                           `using` const A
      <|> lit 'b'                                           `using` const B
pTree 的定義和 grammar S 面貌相似:每個 terminal(小括號、逗號、ab)都套上 lit,concatenation (sequencing) 以 *><*、和 <*> 明確表示,各種可能情況之間則以 <|> 分隔。'using' 子句說明如何建造出對應的樹。

這是怎麼做的呢?

Parser a b 的輸入是一連串的 tokens of type a,它將這串 tokens 的某個 prefix 解析為某個型別為 b 的東西。這個型別該怎麼定義呢?很顯然,Parser a b 接收一個 list of a 並傳回一個 b-structure:

type Parser a b = [a] -> b
但 parser 不一定會耗盡所有的 input tokens,這樣的設計在我們把 parsers 串起來的時候會有用。因此我們除了讓 parser 傳回一個 b-structure 以外,也傳回尚未消耗的 list of tokens:
type Parser a b = [a] -> (b, [a])
最後,parser 可能會因為輸入不合法而解析失敗,或者它可能傳回多個可能的解析結果。因此我們讓 parser 傳回「可能的解析結果」所形成的序列(list):
type Parser a b = [a] -> [(b, [a])]
如果解析失敗,parser 就傳回 []

我們定義另一個函式 parse「執行」一個 parser。這個函式選出第一個「耗盡整個輸入字串」的解析結果:

parse :: Parser a b -> [a] -> b
parse p = fst . head . filter (null . snd) . p
fail 這個 parser 的解析永遠是失敗的:
fail : Parser a b
fail xs = []
與之相反的是 succeed,它總是立即解析成功,不消耗任何輸入:
succeed :: Parser a ()
succeed xs = [((), xs)]
換句話說,這個 parser 接受空字串。

lit x 這個 parser 檢查輸入的第一個 token 是否等於 x,如果真的是這樣就傳回 x(因此它的型別是 Parser a a)以及尚未解析的剩餘輸入。這是唯一可行的解析方式。若不然,它就傳回一個 empty list 代表解析失敗,因為沒有以 x 起頭的解析方式。

lit :: Eq a => a -> Parser a a
lit x (y : xs) | x == y = [(y, xs)]
lit x _ = []

p1 <|> p2 這個 parser(p1 OR p2)試著用 p1p2 分別進行解析,然後合併兩者結果。那就單純是把兩個 parsers 傳回的 lists 串接在一起。

infixr 6 <|>
(<|>) :: Parser a b -> Parser a b -> Parser a b
(p1 <|> p2) xs = p1 xs ++ p2 xs

p1 <*> p2 這個 parser(p1 followed by p2)先用 p1 解析輸入字串,再用 p2 解析剩餘的字串。假設 p1 造出一個 b-structure、p2 造出一個 c-structure,那麼 p1 <*> p2 會把兩者的結果放進一個 pair (b, c) 裡面當作自己的結果。這個 parser combinator 用 list comprehension 會比較好寫:輸入 xs,對於 p1 xs 的每個解析結果 (b, ys),用 p2 解析 ys 產生一群結果 (c, zs),然後把 (b, c) 和剩餘的輸入 zs 組合成一個 pair 放進 resulting list。

(<*>) :: Parser a b -> Parser a c -> Parser a (b, c)
(p1 <*> p2) xs = [ ((b, c), zs) | (b, ys) <- p1 xs, (c, zs) <- p2 ys ]

給定一個 parser p :: Parser a b 和一個函式 f :: b -> c(p `using` f) 是一個新的 parser,把 f 套用到 p 的結果:

using :: Parser a b -> (b -> c) -> Parser a c
p `using` f = map (cross f id) . p
其中 cross 是標準的 product functor:
split f g a = (f a, g a)
cross f g   = split (f . fst) (g . snd)

p1 *> p2 很類似 p1 <*> p2,但它會丟棄 p1 的結果。這個 combinator 可用 <*>using 定義:

infixr 8 *>
(*>) :: Parser a b -> Parser a c -> Parser a c
p1 *> p2 = (p1 <*> p2) `using` snd
<**> 的對偶:
infixr 8 <*
(<*) :: Parser a b -> Parser a c -> Parser a b
p1 <* p2 = (p1 <*> p2) `using` fst

至此我們已經定義了一個內嵌於 Haskell 的 domain-specific embedded language,專門用來解析語言。讓我們試用這些 combinators。以下這個簡單的 grammar 描述「反二進位數字」:

B ::= 0 | 1 | 0 B | 1 B
例如 "11010" 的值是 11(因為數字的順序是反的)。欸,選這種奇怪的例子只是為了練習啦 XD。以下就是它的 parser:
pB :: Parser Char Int
pB =     lit '0'         `using` const 0
     <|> lit '1'         `using` const 1
     <|> (lit '0' *> pB) `using` (*2)
     <|> (lit '1' *> pB) `using` ((+1) . (*2))
一些執行例子:pB "1101" = [(1,"101"),(3,"01"),(3,"1"),(11,"")]parse pB "1101" = 11

如果我們改用下面這個 grammar 可以嗎?

B ::= 0 | 1 | B 0 | B 1
很可惜,這個 grammar 是個 left-recursive grammar,而我們的 recursive-descent parser 只能處理 LL grammar。LL 的第二個 L 代表解析方式是 leftmost derivation,所以如果有 left-recursion 就會遞迴個沒完。

回到我們一開始定義的 parser pTree。如我們一再提及,用這樣的方法寫出的 parsers 都是 recursive-descent (hence top-down) parsers,也就是從 start symbol 逐漸展開一棵樹(unfold?!)。詳情可參考 dragon book 2/e, section 4.4,其中 p.218 有精美插圖 XD。單純的 recursive-descent parsing 可能需要 backtrack,因為展開某個 nonterminal 的時候可能選錯 production rule。這裡的精神是 "nondeterministically" 選一個 production rule,實作上用 subset construction 把所有可能的結果收集起來,到最後才挑出成功解析的結果。

另一種印樹的形式是把 Bin (t, u) 印為 (t)u。例如 Bin (Bin (A,B), Bin (A,B)) 印為 ((a)b)(a)b。對應的 grammar 是

S ::= ( S ) S | a | b
Parser 則可寫為
pTree2 :: Parser Char Tree
pTree2 =
      lit 'a'                                   `using` const A
  <|> lit 'b'                                   `using` const B
  <|> (lit '(' *> pTree2 <*> lit ')' *> pTree2) `using` Bin

最後讓我們試一個稍微實際一點的例子,處理「反二進位數字」的加法與乘法算式。算式的 grammar 是

Expr ::= Expr + Term | Term
Term ::= Term * Factor | Factor
Factor ::= B | ( Expr )
是個超標準的版本。(我高中看 C++ grammar 的時候就看過啦 XD。)不過這個 grammar 是 left-recursive,所以我們必須做點變換,消除 left-recursion(ref. dragon book, section 4.3.3)。這過程其實就是把 non-empty snoc-lists 換成對應的 cons-lists。轉換後的 grammar 變成
Expr ::= Term Expr'
Expr' ::= + Term Expr' | ε
Term ::= Factor Term'
Term' ::= * Factor Term' | ε
Factor ::= B | ( Expr )
其實差不多就是 dragon book 的 example 4.27。

定義 Token 型別為

data Token = Number Int | Plus | Multiply | LeftP | RightP
     deriving (Eq, Show)
pToken :: Parser Char [Token] 會把輸入字串裂解為一連串的 tokens。為了讓效率好一點,我們定義 munch :: Parser a b -> Parser a b 把一個 parser 轉換為 "greedy" parser ─ resulting parser 只挑出剩餘字串長度最短的那個結果。
munch :: Parser a b -> Parser a b
munch p = cond null id ((:[]) . foldr1 f) . p
  where f p q = if sndlen p < sndlen q then p else q
        sndlen = length . snd
其中 cond 是 McCarthy conditional form:
cond p f g a = if p a then f a else g a
於是我們就可以定義 pToken
pToken :: Parser Char [Token]
pToken =
      succeed               `using` const []
  <|> (lit ' ' *> pToken)
  <|> (lit '+' *> pToken)   `using` (Plus :)
  <|> (lit '*' *> pToken)   `using` (Multiply :)
  <|> (lit '(' *> pToken)   `using` (LeftP :)
  <|> (lit ')' *> pToken)   `using` (RightP :)
  <|> (munch pB <*> pToken) `using` (uncurry (:) . cross Number id)
其中第二條是忽略空格的規則。令
s = "001101 + 101 * (011 + 1 + 10)"
那麼 parse pToken s = [Number 44,Plus,Number 5,Multiply,LeftP,Number 6,Plus,Number 1,Plus,Number 1,RightP]

接下來我們定義 abstract syntax tree 的型別。

data ExprTree = Operand Int
              | PlusBin ExprTree ExprTree
              | MultiplyBin ExprTree ExprTree
     deriving Show
對應的 fold operator 是
foldet f g h (Operand n) = h n
foldet f g h (PlusBin t u) = f (foldet f g h t) (foldet f g h u)
foldet f g h (MultiplyBin t u) = g (foldet f g h t) (foldet f g h u)
想對一個算式求值,可用以下函式:
evalExpr :: [Char] -> Int
evalExpr = foldet (+) (*) id . parse pExpr . parse pToken
其中 pExpr :: Parser Char ExprTreeExpr grammar 的 parser。例如 parse pExpr (parse Token s) = PlusBin (Operand 44) (MultiplyBin (Operand 5) (PlusBin (PlusBin (Operand 6) (Operand 1)) (Operand 1)))。定義 pExpr 會碰上一點麻煩,因為經過 left-recursion elimination 的 grammar 和我們想建造的樹形不一樣了!更精確地講,pExpr'(注意是 Expr "prime")傳回的是個 partial syntax tree,裡面有個洞是要放 + 左邊的樹。But hey, we are working in a functional language! 我們就讓 pExpr' 的結果是 ExprTree -> ExprTree,這個 function 拿一個 ExprTree 塞到那棵 partial syntax tree 的洞裡面,傳回一棵完整的樹。這樣我們就可以定義出全部的 parsers:
pExpr :: Parser Token ExprTree
pExpr = (pTerm <*> pExpr') `using` uncurry (flip ($))

pExpr' :: Parser Token (ExprTree -> ExprTree)
pExpr' =
      succeed
          `using` const id
  <|> (lit Plus *> pTerm <*> pExpr')
          `using` (\(t, f) -> \u -> f (PlusBin u t))

pTerm :: Parser Token ExprTree
pTerm = (pFactor <*> pTerm') `using` (uncurry (flip ($)))
  
pTerm' :: Parser Token (ExprTree -> ExprTree)
pTerm' =
      succeed
          `using` const id
  <|> (lit Multiply *> pFactor <*> pTerm')
          `using` (\(t, f) -> \u -> f (MultiplyBin u t))

pFactor :: Parser Token ExprTree
pFactor = lit LeftP *> pExpr <* lit RightP
      <|> pNumber `using` Operand

pNumber :: Parser Token Int
pNumber (Number x : ys) = [(x, ys)]
pNumber _ = []

有了傳回 syntax tree 的 parser,要轉 postfix form 也很容易:

postExpr :: [Char] -> [Char]
postExpr = foldet f g h . parse pExpr . parse pToken
           where f s t = s ++ t ++ "+"
                 g s t = s ++ t ++ "*"
                 h x   = "[" ++ show x ++ "]"
例如 postExpr s = "[44][5][6][1]+[1]+*+"

--
好長 Orz。

Labels: ,

2008/03/29

Parsing Combinators in Haskell

今天排定的進度算結束了,唸編譯器是其中一項,所以就順勢來做 scm 老師這次的 Haskell exercises。目標是寫出一組可合成 context-free parser 的 higher-order functions,精神和 Boost spirit library 應該是一樣的。後者我是上了大學才聽聞,沒怎麼研究過。

看來兩者產生的都是 recursive-descent (hence LL) parsers。然後 Haskell 版處理 nondeterminism 的方式是用標準的 subset construction,其中 sets 以 lists 實作。

--
先貼一篇佔位,之後會補完 XD。


(3:06AM) 寫了一個很簡單的算式求值器,支援 "reversed binary numbers" 的加法與乘法,可用小括號表示子算式。架構是標準的 tokenisation → parsing → tree flattening。為了讓加和乘是 left-associative,我寫出來的 grammar 是 left-recursive,所以還用上才剛讀到的 left-recursion elimination。

--
顯然又是寫篇 literate Haskell 的時候了,待續 XD。


我直接在 blog 上寫一篇好了,parsing combinator 的符號和 lhs2TeX 不太有默契 XD。

--
請見《Parser Combinators in Haskell》XD。

Labels: ,

2008/03/20

Modelling Digital Logic in Haskell

scm 老師這次出的 Haskell 練習題是用 Haskell 寫 digital logic components。這正是 OUCL 大學部教 digital logic 的方式,真是太浪漫了 XD。一條接線是以一個 list 表示,其中每個元素是那條線在各個離散時間點的訊號值。

type Wire = [Bool]
也就是說,Wire 這個 list 的 index 是時間(單位是 infinitesimal 吧,我們暫時不考慮 physical delays XD)。另一種正交觀點是在同一個時間點看多條線(一般就是一個 bus)的 states,我稱之為 time slice(從 OS 借來的詞)。
type TimeSlice = [Bool]
此時 list 的 index 是線的編號。後面會發現區分這兩個型別是有意義的。

我們先從簡單的開始。以下的函式 binary c 把一個整數轉為 binary representation:

binary :: Int -> Integer -> TimeSlice
binary c = take c . map ((/= 0) . (`mod` 2)) . iterate (`div` 2)
這產生出線數為 c 的 time slice。binary c 有反函式 numeric c
numeric :: Int -> TimeSlice -> Integer
numeric c = sum . tri (*2) . map bit2num . take c
tri f(就是出現在 Horner's rule 裡面的那個函式)和 bit2num 的定義是:
tri f = foldr (\x xs -> x : map f xs) []
bit2num False = 0 
bit2num True  = 1
我們也可以產生承載常數的 bus:
constant :: Int -> Integer -> [Wire]
constant c = map repeat . binary c
其中 repeat a = a : repeat a 是 Haskell 內建的函式(吧)。

接著我們試試加法器的核心功能。以下函式 add c 把兩個 time slices 視為二進位數字加在一起,產生一個新的 time slice:

add :: Int -> TimeSlice -> TimeSlice -> TimeSlice
add c = curry (binary (c + 1) . uncurry (+) . parallel (numeric c))
add c 把它收到的兩個 time slices 轉為數值加起來,然後轉換回長度 c + 1 的 time slice,多一位是最高位的進位。這其實是加法電路的 specification,以後我希望能從這裡導出 ripple adder, carry look-ahead adder 等等的。(Hope it's possible. XD)

加法器的型別是

adder :: Int -> [Wire] -> [Wire] -> [Wire]
[Wire] 是一條一條的線,但 add c 處理的是一個一個的 time slice,我們必須寫個函式在這兩個觀點間轉換。稍微想一想就會發覺,這其實就是 matrix transposition,所以我們可以寫
transpose = foldr (zipWith (:)) (repeat [])
然後我們就可以寫出加法器:
adder c = curry (transpose . uncurry (zipWith (add c)) . parallel transpose)
其中 parallel 是 "square" functor。
split f g a = (f a, g a)
cross f g = split (f . fst) (g . snd)
parallel f = cross f f
寫完了很高興,我們測試一下,在 GHCi 下輸入
map head (adder 5 (constant 5 7) (constant 5 8))
卻發現 adder is unproductive!(意思是什麼東西都吐不出來。)追蹤一下就會發現,add 產生出來的是個 infinite list of time slices,而 transpose 是個 foldr,找不到 base case 可以停下來。因此我們被迫寫另一個 transpose,功能一樣但用 unfoldr 實作:
wirewise :: [TimeSlice] -> [Wire]
wirewise = unfoldr f
           where f ([] : xss) = Nothing
                 f xss        = (Just . split (map head) (map tail)) xss
其中 unfoldr f 是標準的定義:
unfoldr f a = case f a of
                Nothing -> []
                Just (x, b) -> x : unfoldr f b
本來的 transpose 也改名一下:
timewise :: [Wire] -> [TimeSlice]
timewise = foldr (zipWith (:)) (repeat [])
adder 於是就變成
adder :: Int -> [Wire] -> [Wire] -> [Wire]
adder c = curry (wirewise . uncurry (zipWith (add c)) . parallel timewise)
這樣我們就有一個把兩條二進位數字流加起來的加法器,而且是 productive。例如剛剛在 GHCi 下輸入的算式現在就能正確輸出結果了:
*Main> map head (adder 5 (constant 5 7) (constant 5 8))
[True,True,True,True,False,False]

我刻意把各個部份寫得很 compositional,希望可以多做到一點 derivations ─ 雖然我覺得規模有點恐怖了,而且還有 infinite lists 和 unfold XD。

--
還有得玩,不過要過一陣子 XD。

Labels:

2008/03/04

Haskell Exercises --- Done

寫完了!真正的程式只有二十幾行,產出的 pdf 卻有六頁 XD。

--
後面就比較懶得正式證明 XD。

Labels: , ,