２次元配列を scan する ― Haskellプログラミングの一風景

この記事は、Haskell Advent Calendar 2024 13日目の記事です。

はじめに

プログラミングをしている様子を傍から見るのは面白いものです。少くとも、この記事の筆者はそう感じています。いろいろな気付きがあって楽しいものです。
ここでは、シンプルなプログラミングのお題をめぐって、あれこれプログラミングしてみます。Haskellプログラミングのエキスパートにとっては他愛もないものでしょう。それでも、「えぇー、それ、自分とは違う」ところがあれば、あるいは、「あらら、自分ならもっとうまくやる」なんてことがあれば、それはそれで、ツッコんで楽しめると思います（たぶん）。

プログラミングのお題

$10\times 10$ の行列が与えられ、$(1,1)$成分から右か下に向って成分を足しながら進んでいくとき、$(i,j)$成分へ至る道のうち最も和が小さいものを求める（以下略）
（haskell-jp Slack に投稿されたコメントより）

素朴な実装

minPath0

このお題を素朴にminPath0として実装してみましょう。
まずは仕様を型で表現しましょう。

minPath0 :: Array (Int,Int) Int -- 与えられた行列
         -> (Int,Int)           -- 指定された到達点
         -> Int                 -- 到達点までの最小和

この型シグネチャーを見つめながら、実装を考えてみましょう。

与えられた行列はsa、指定された到達点は(i,j)としましょう。

minPath0 sa (i,j) = undefined

点(i,j)に到達する直前には、点(i-1,j)または点(i,j-1)にいたはず。到達点までの最小和は、直前点までの最小和の大きくない方の値と到達点の値との和になります。

minPath0 sa (i,j) = min 
                        (minPath0 sa (i-1,j)) -- 点(i-1,j)までの最小和
                        (minPath0 sa (i,j-1)) -- 点(i,j-1)までの最小和
                   + 
                        sa ! (i,j)            -- 到達点の値

ただし、出発点(1,1)の場合は直前点は上にも左にもなく、上辺にある点(1,j)の場合は直前点は左だけ、左辺にある点(i,1)の場合は直前点は上だけです。これを定義に追加すると以下のとおりになります。

minPath0 sa (1,1) = sa ! (1,1)
minPath0 sa (1,j) = minPath0 sa (1,j-1) + sa ! (1,j)
minPath0 sa (i,1) = minPath0 sa (i-1,1) + sa ! (i,1)
minPath0 sa (i,j) = minPath0 sa (i-1,j) `min` minPath0 sa (i,j-1) + sa ! (i,j)

この実装は、問題をそのまま素朴に表現したものですね。少くとも、これで求まる値は正しい値であると容易に想像できると思います。

実際の計算

では、実際に具体的に計算してみましょう。minPath0に与える行列をsampleとして定義しましょう。
件のコメントでは以下のようなものを例として挙げてありました。

sample :: Array (Int,Int) Int
sample = listArray ((1,1),(10,10)) [i+j | i <- [1..10], j <- [1..10]]

中身を表示してみましょう。

stack exec -- ghc -e 'elems sample' src/MinPath.hs

[2,3,4,5,6,7,8,9,10,11,3,4,5,6,7,8,9,10,11,12,4,5,6,7,8,9,10,11,12,13,5,6,7,8,9,10,11,12,13,14,6,7,8,9,10,11,12,13,14,15,7,8,9,10,11,12,13,14,15,16,8,9,10,11,12,13,14,15,16,17,9,10,11,12,13,14,15,16,17,18,10,11,12,13,14,15,16,17,18,19,11,12,13,14,15,16,17,18,19,20]

見やすいように、ここでは手動で^^;10要素ずつ改行すると、

[2,3,4,5,6,7,8,9,10,11
,3,4,5,6,7,8,9,10,11,12
,4,5,6,7,8,9,10,11,12,13
,5,6,7,8,9,10,11,12,13,14
,6,7,8,9,10,11,12,13,14,15
,7,8,9,10,11,12,13,14,15,16
,8,9,10,11,12,13,14,15,16,17
,9,10,11,12,13,14,15,16,17,18
,10,11,12,13,14,15,16,17,18,19
,11,12,13,14,15,16,17,18,19,20]

これを見るとわかるように（中身を見るまでもなく、配列の作り方をみればわかりますが）、これは、正しい経路であれば、どの経路をとおっても、その軌跡は同じリストになり、和は同じになってしまいます。１つの例としてはなりたちますが、特例中の特例なので、もうすこし、散けたものを作りましょう。

散けた要素をもつ配列の生成

とりあえず、[1 .. n]をシャッフルしたリストをn個用意する、shuffledを定義しよう。

Shuffle.hs

module Shuffled where

import Data.List
import System.Random -- random パッケージ
import List.Shuffle  -- list-shuffle パッケージ

shuffled :: Int -> IO [[Int]]
shuffled n = take n . unfoldr phi <$> getStdGen where
    bs  = [1 .. n]
    phi = Just . shuffle bs

これで、$10\times 10$の数字が散けた行列を表示して、それをそのまま、コードに貼り付けてしまいましょう。

stack exec -- ghc -e 'shuffled 10' src/Shuffled.hs

[[8,9,6,10,5,4,2,7,3,1],[1,8,9,5,2,4,3,6,7,10],[3,7,2,5,4,8,6,10,1,9],[5,1,7,6,9,4,3,8,2,10],[6,9,2,3,1,8,10,4,7,5],[10,5,8,3,1,9,2,7,4,6],[3,7,9,6,2,4,1,8,10,5],[9,3,2,1,8,7,6,5,4,10],[6,9,2,1,8,10,5,7,4,3],[2,6,3,10,7,1,4,9,8,5]]

これをもとにsampleを作成しましょう。

sample :: Array (Int,Int) Int 
sample = listArray ((1,1),(10,10)) (concat bss) where
    bss = [[8,9,6,10,5,4,2,7,3,1]
          ,[1,8,9,5,2,4,3,6,7,10]
          ,[3,7,2,5,4,8,6,10,1,9]
          ,[5,1,7,6,9,4,3,8,2,10]
          ,[6,9,2,3,1,8,10,4,7,5]
          ,[10,5,8,3,1,9,2,7,4,6]
          ,[3,7,9,6,2,4,1,8,10,5]
          ,[9,3,2,1,8,7,6,5,4,10]
          ,[6,9,2,1,8,10,5,7,4,3]
          ,[2,6,3,10,7,1,4,9,8,5]]

では、これを使って、右下の点(10,10)での最小和を求めてみましょう。ついでに簡単に計算時間も測ることにしましょう。

% time stack exec -- ghc -e 'minPath0 sample (10,10)' src/MinPath.hs

結果は、以下のとおり。

66
--
Program : stack exec -- ghc -e 'minPath0 sample (10,10)' src/MinPath.hs
CPU     : 98%
user    : 0.383s
system  : 0.114s
total   : 0.506s
--

もうすこし、大きいサイズ$100\times 100$の配列を用意して実験しましょう。

genSample :: Int -> IO ()
genSample n = writeFile "sample.txt"
            . (++"\n")
            . show
            . listArray ((1,1),(n,n)) . concat =<< shuffled n

getSapmle :: IO (Array (Int,Int) Int)
getSample = read <$> readFile "sample.txt"

genSample 100 で $100\times 100$のサンプル行列を作成してファイルにテキスト保存します。

% stack exec -- ghc -e 'genSample 100' src/MinPath.hs

これを getSample で読み込んで使用します。

% for i in {10..15}
for> time stack exec -- ghc -e "flip minPath0 (${i},${i}) <$> getSample" src/MinPath.hs

832
--
Program : stack exec -- ghc -e "flip minPath0 (${i},${i}) <$> getSample" src/MinPath.hs
CPU     : 97%
user    : 0.412s
system  : 0.098s
total   : 0.521s
--
922
--
Program : stack exec -- ghc -e "flip minPath0 (${i},${i}) <$> getSample" src/MinPath.hs
CPU     : 98%
user    : 0.577s
system  : 0.092s
total   : 0.682s
--
989
--
Program : stack exec -- ghc -e "flip minPath0 (${i},${i}) <$> getSample" src/MinPath.hs
CPU     : 99%
user    : 1.242s
system  : 0.097s
total   : 1.351s
--
1152
--
Program : stack exec -- ghc -e "flip minPath0 (${i},${i}) <$> getSample" src/MinPath.hs
CPU     : 99%
user    : 3.788s
system  : 0.109s
total   : 3.904s
--
1149
--
Program : stack exec -- ghc -e "flip minPath0 (${i},${i}) <$> getSample" src/MinPath.hs
CPU     : 100%
user    : 13.356s
system  : 0.110s
total   : 13.460s
--
1282
--
Program : stack exec -- ghc -e "flip minPath0 (${i},${i}) <$> getSample" src/MinPath.hs
CPU     : 100%
user    : 51.578s
system  : 0.152s
total   : 51.667s
--

さすがに遅いですね。計算時間が、指数オーダーになっているようです。

原因は、minPath0に、ツリー構造再帰がある（下線部）からですね。

minPath0 sa (1,1) = sa ! (1,1)
minPath0 sa (1,j) = minPath0 sa (1,j-1) + sa ! (1,j)
minPath0 sa (i,1) = minPath0 sa (i-1,1) + sa ! (i,1)
minPath0 sa (i,j) = minPath0 sa (i-1,j) `min` minPath0 sa (i,j-1) + sa ! (i,j)
--                  ^^^^^^^^^^^^^^^^^^^       ^^^^^^^^^^^^^^^^^^^ 

このminPath0の定義は仕様としては正しい定義ですが、実用には不向きです。

ただ改善の余地はあります。4つめの等式の右辺にあらわれた分岐の枝をさらに一段展開すると、

minPath0 sa (i,j) = (minPath0 sa (i-2,j) `min` minPath0 sa (i-1,j-1) + sa ! (i-1,j)) 
                    `min` --                  ^^^^^^^^^^^^^^^^^^^^^
                    (minPath0 sa (i,j-2) `min` minPath0 sa (i-1,j-1) + sa ! (i,j-1))
                    +     --                  ^^^^^^^^^^^^^^^^^^^^^
                    sa ! (i,j)

下線部の関数適用式が重なっています。なんども同じ計算が行われることになるので、これを１回ですませられるようにします。それには元の配列saのすべての点までのminPathの値で構成された、もとのsaと同じディメンションの配列taを構成します。これが作成できれば、(i,j)におけるminPathの値は、ta ! (i,j)になります。

minPath :: Array (Int,Int) -> (Int,Int) -> Int
minPath sa (i,j) = ta ! (i,j) where
    ta = listArray (bounds sa) (phi <$> assocs sa)
    phi = \ case
        ((1,1), s) -> s
        ((1,j), s) -> ta ! (1,j-1) + s
        ((i,1), s) -> ta ! (i-1,1) + s
        ((i,j), s) -> (ta ! (i,j-1)) `min` (ta ! (i-1,j)) + s

minPath0の再帰適用minPath0 sa (i,j-1)が、taのインデックスによる参照ta ! (i,j-1)に、minPath0 sa (i,j-1)が、ta ! (i,j-1)に、代っています。早速ためしてみましょう。

% for i in {10..15}
for> time stack exec -- ghc -e "flip minPath (${i},${i}) <$> getSample" src/MinPath.hs

832
--
Program : stack exec -- ghc -e "flip minPath (${i},${i}) <$> getSample" src/MinPath.hs
CPU     : 98%
user    : 0.363s
system  : 0.093s
total   : 0.462s
--
922
--
Program : stack exec -- ghc -e "flip minPath (${i},${i}) <$> getSample" src/MinPath.hs
CPU     : 96%
user    : 0.361s
system  : 0.087s
total   : 0.463s
--
989
--
Program : stack exec -- ghc -e "flip minPath (${i},${i}) <$> getSample" src/MinPath.hs
CPU     : 98%
user    : 0.350s
system  : 0.094s
total   : 0.453s
--
1152
--
Program : stack exec -- ghc -e "flip minPath (${i},${i}) <$> getSample" src/MinPath.hs
CPU     : 96%
user    : 0.336s
system  : 0.097s
total   : 0.447s
--
1149
--
Program : stack exec -- ghc -e "flip minPath (${i},${i}) <$> getSample" src/MinPath.hs
CPU     : 97%
user    : 0.345s
system  : 0.099s
total   : 0.455s
--
1282
--
Program : stack exec -- ghc -e "flip minPath (${i},${i}) <$> getSample" src/MinPath.hs
CPU     : 98%
user    : 0.353s
system  : 0.091s
total   : 0.453s
--

ほとんど、変化がありませんね。点$(M,N)$におけるminPathの値を計算量オーダーが指数オーダーから、$O(MN)$にまで低減されたからですね。

minPath は、指定の1点のみを計算するより、指定の点より前にある点もすべて求めてから、指定の点をインデックスアクセスするほうが効率がよいということになります。
このような計算構造は、1次元の例ですが、フィボナッチ数を求めるのにも（よりシンプルですが）現れます。

fib :: Array Int Integer -> Int -> Integer
fib sa n = ta ! n 
    where
        ta = listArray (bounds sa) (phi <$> assocs sa)
        phi = \ case
            (0,j) -> j
            (1,j) -> j
            (i,_) -> ta ! (i-2) + ta ! (i-1)

nats :: Array Int Integer
nats = listArray (0,1000) [0 .. 1000]

% time stack exec -- ghc -e 'fib nats 1000' src/Fib.hs
43466557686937456435688527675040625802564660517371780402481729089536555417949051890403879840079255169295922593080322634775209689623239873322471161642996440906533187938298969649928516003704476137795166849228875
--
Program : stack exec -- ghc -e 'fib nats 1000' src/Fib.hs
CPU     : 95%
user    : 0.274s
system  : 0.063s
total   : 0.352s
--

一般化： scanArray

minPath は配列sa :: Array (Int, Int) Intを別の配列ta :: Array (Int,Int) Intに変換して、taにインデックスアクセスするというものでした。この配列変換は蓄積演算をパラメータとして一般化できそうです。
$(i,j)$の直前点をすこし一般化して、$(i,j-1)$、$(i-1,j)$の2点に加えて、$(i-1,j-1)$も加えておくと応用が効きそうです。

scanArray :: (Ix i, Enum i)
          => (a -> b)                 -- 左上の点（直前点なし）
          -> (b -> a -> b)            -- 左辺または上辺上の点（直前点は(i,j-1)または(i-1,j)のどちらか1つ）
          -> (b -> b -> b -> a -> b)  -- 内部の点（直前点は(i-1,j-1)、(i,j-1）、(i-1,j)の3つ
          -> Array (i,i) a
          -> Array (i,i) b
scanArray f g h sa = ta where
    ta  = listArray (bounds sa) (phi <$> assocs sa)
    phi = \ case
        (ij@(i,j),a)
            | ij == ij0 -> f a
            | i  == i0  -> g (ta ! second pred ij) a
            | j  == j0  -> g (ta ! first  pred ij) a
            | otherwise -> h (ta ! (pred *** pred) ij) (ta ! first  pred ij) (ta ! second pred ij) a
            where
                ij0 = fst (bounds sa)
                i0  = fst ij0
                j0  = snd ij0

scanArray を使うと minPath は以下のように定義できます。

minPath2 :: Array (Int,Int) Int -> (Int,Int) -> Int
minPath2 = (scanArray f g h !) where
    f = id
    g = (+)
    h _ a b c = a `min` b + c

% time stack exec -- ghc -e 'flip minPath2 (100,100) <$> getSample' src/MinPath.hs

5544
--
Program : stack exec -- ghc -e 'flip minPath2 (100,100) <$> getSample' src/MinPath.hs
CPU     : 98%
user    : 0.521s
system  : 0.148s
total   : 0.679s
--

minPathが経路最小和でしたから、minの部分をmaxに置きかえれば、経路最大和maxPathができあがります。

maxPath :: Array (Int,Int) Int -> (Int,Int) -> Int
maxPath = (scanArray f g h !) where
    f = id
    g = (+)
    h _ a b c = a `max` b + c

scanArrayは、二次元の累積和の一般化とみなすことができるので、累積和行列を構成することもできます。

cumulativeSum2D :: Num a => Array (Int,Int) a -> Array (Int,Int) a
cumulativeSum2D = scanArray f g h where
    f = id
    g = (+)
    h a b c d = subtract a (b + c + d)

まとめ

minPath をプログラムするというところからはじめて、つれづれにプログラムして、汎用になりそうな、scanArray に到達しました。

scanArrayは結構便利そうなので、どこかのパッケージですでに、似たもの、あるいは、より一般化されたものが、採録されていそうですね。

筆者の他愛のないプログラミングの楽しみの断片を紹介しました。計測、分析はおおざっぱです。「ふう〜ん、そうするんだ」なんて思っていただければ、幸いです。

今回、書きちらしたコードは、筆者の書き散らし箱（今回つくった^^;） においてあります。御笑覧ください。