Moonmile Solutions Blog

数独をF#で解く – Moonmile Solutions Blog
http://www.moonmile.net/blog/archives/5304

の続き。二次元トラスは別途「Visual Basicでわかる やさしい有限要素法の基礎」が届いたので、これを見てから。この本、VB6 なのかと思ったら、VB2008 で書かれた新しい本でした。これだったら、結構流用がきくかも。元ネタは Fortran らしいのですが、グラフィックな部分とか .NET で書かれているので助かる。XAMLのPATHあたりに直すと結構いいかもという感じですね。

さて、

数独の問題を作るほうは、

1. 矛盾なくNxNを敷き詰める。
2. ランダムに1個空白にする。
3. 解が1つしかないことを確認する。再び2へ

を地道に実装したのが次です。やっぱり、再帰のところがうまく作れなくて、excepiton で抜けていますが、これはいずれ修正する…かも。

open System

let bsize = 3
let size = bsize * bsize
let size2 = size * size
let pazzle = Array2D.zeroCreate<int> size size

// ルールにマッチする数を取り出す
let rec rule (m:int[,]) x y v (lst:int list) =
    // 横一列をチェック
    let rule_x (m:int[,]) x y v =
        [|for i in 0..size-1 -> m.[y,i]|] |> Array.exists (fun t -> t = v)
    // 縦一列をチェック
    let rule_y (m:int[,]) x y v =
        [|for i in 0..size-1 -> m.[i,x]|] |> Array.exists (fun t -> t = v)
    // boxをチェック
    let rule_box (m:int[,]) x y v = 
        let x' = x/bsize*bsize
        let y' = y/bsize*bsize
        [| for i in 0..bsize-1 do
            for j in 0..bsize-1 -> m.[y'+j,x'+i] |] |> Array.exists (fun t -> t = v)

    if v = 0 then 
        []
    elif rule_x m x y v = false &&
        rule_y m x y v = false &&
        rule_box m x y v = false then
        v::rule m x y (v-1) lst
    else
        rule m x y (v-1) lst

// ランダムに1つ選ぶ
let rnd = new Random()

let rec remove i (l:'a list) =
     match i,l with
     | 0, x::xs -> xs
     | i, x::xs -> x::remove (i-1) xs
     | i, [] -> failwith "index out of range"

let rec shaffle (l:'a list) =
    match l with
    | x::[] -> [x]
    | x::xs -> 
        let i = rnd.Next(l.Length)
        let l' = remove i l
        l.[i]::shaffle l'
    | [] -> []

let rec quest (m:int[,]) x y =
    // printfn "%d %d" x y
    // printfn "%A" m
    if x = 9 && y = 8 then
        // 完成
        printfn "success."
        printfn "%A" m
        raise (new Exception("ok"))
    elif x = 9 then
        quest m 0 (y+1) 
    else
        let lst = rule m x y size []
        if lst.Length > 0 then
            let lst' = shaffle lst
            for i in 0..lst'.Length-1 do
                let v = lst'.[i]
                let m' = Array2D.copy m
                m'.[y,x] <- v
                quest m' (x+1) y 
        
let Quest (m:int&#91;,&#93;) =
    try 
        quest m 0 0
    with
        | _ -> printfn "ok."

// 数字の敷き詰めを作成
Quest pazzle

ごちゃごちゃしていますが、1番の数の敷き詰めを作っているところです。左上から順番に番号をいれて、ルールに沿って（縦/横/ボックス）候補の数を配置しています。配置するときにランダムに数値を選ばせるために、shaffle を使っています。単純に数え上げのところは、この部分はいらないのですが、それだと同じ問題になってしまうし。
3×3の場合には、結構なスピードで敷き詰めができあがります…が、4×4にすると途端に遅くなりますね。1分程度で終わる場合もあれば、30分やってやっと終わる場合もあります。このあたり、最適化の山を登っている感じ。達成しない準最適化の山に登ってしまうと一度降りるのに時間がかかります。将棋の枝切りとかこのあたりのロジックなのかも。そうそう、やねうらおさんの「探索の深さが強さのイコールとはならない」というのは、このあたりの話で、完全に探索ができない場合には「探索自体の深さ」は「強さ」＝正解そのものとは違っていて、さらに将棋の指し手が極めて少ない（２，３手とか）のであれば、その探索の深さ（指し手の多さ）は、強さに比例しない…だろう、ってことだと思います。まあ、実装できるのがすごいんですが。なんとなく想像で。

一旦敷き詰めができた配列をコピーして、今度は問題を作っていきます。

let ans' = 
    [[7; 3; 1; 2; 6; 4; 9; 8; 5]
     [9; 4; 2; 5; 3; 8; 7; 6; 1]
     [8; 6; 5; 7; 1; 9; 3; 2; 4]
     [2; 7; 9; 8; 5; 3; 4; 1; 6]
     [6; 5; 4; 9; 7; 1; 8; 3; 2]
     [1; 8; 3; 6; 4; 2; 5; 7; 9]
     [5; 2; 8; 3; 9; 6; 1; 4; 7]
     [4; 9; 6; 1; 8; 7; 2; 5; 3]
     [3; 1; 7; 4; 2; 5; 6; 9; 8]]
let A' = Array2D.init 9 9 (fun i j -> ans'.[i].[j])

let mutable scnt = 0

let rec solve (m:int[,]) x y = 
    // printfn "%d %d" x y
    if x = 9 && y = 8 then
        // 完成
        scnt <- scnt + 1
    elif x = 9 then
        solve m 0 (y+1)
    elif m.&#91;x,y&#93; <> 0 then
        solve m (x+1) y 
    else 
        // 縦/横/boxをチェック
        let mutable v = [|0..9|]
        for i in 0..8 do v.[m.[i,y]] <- 0
        for i in 0..8 do v.&#91;m.&#91;x,i&#93;&#93; <- 0
        let x0 = x/3*3
        let y0 = y/3*3
        for i in 0..2 do
            for j in 0..2 do
                v.&#91;m.&#91;x0+i,y0+j&#93;&#93; <- 0

        let m' = Array2D.copy m
        for i in 0..9 do
            if v.&#91;i&#93; <> 0 then
                m'.[x,y] <- v.&#91;i&#93;
                solve m' (x+1) y

let make_pazzle (A0:int&#91;,&#93;) =
    let A = Array2D.copy A0 
    
    // (x,y)をランダムに決める
    let slst = 
        shaffle &#91;for y in 0..size-1 do
                    for x in 0..size-1 -> (x,y) ]
    // 一つずつ取り出して、解法が2以上になった時にやめる
    for i in 0..size2-1 do
        printfn "%d" i
        printfn "%A" A'
        let A' = Array2D.copy A
        let (x,y) = slst.[i]

        A.[x,y] <- 0
        scnt <- 0
        solve A 0 0 
        if scnt > 1 then
            printfn "%A" A'
            raise (new Exception("ok"))

make_pazzle A'

問題の手順も簡単で、ランダムに1個空白（0にする）して、それをいちいち解いていきます。このとき、解法が2つ以上あれば、そこでストップしているだけです。これも非常に（コンピュータにとって）手間がかかる方法なのですが、なんか他に思いつかなかったので、これで。
空白にするコマをランダムに決めてしまうので、いわゆる数独の難しさを考慮していません。なので、これも低い最適化の山に登ってしまうと、途端につまらない≒空白の数の少ない問題を作成してしまいます。このあたりは、一定量の空白になるまで問題作成を繰り返せばいいんでしょうが、それだといつまでやっても終わらないって感じになりそうです。これも適当な足切りが必要でしょうね。

同じものをC#で書いたらどうなるんだろう？という疑問もありますが、このぐらいだと多分似た感じで書けそうです。
F#で書いてみてわかったんですが、配列を

F(x+1) = F(x) + a

の問題で書き直すとうまくいきます…つーか、逆な話で、この形の数式が出てくる場合は F# の再帰関数を使うと便利なわけです。実はこの形って、フィードバックそのもなので、F(1)の結果がF(2)に影響して、さらにF(2)の結果がF(3)に影響して、最終的にF(n)に影響する≒結果が出る、というパターンですね。F(F(F(F(…)))) な感じ。これを再帰の終端が、N から始まって 0 にして終わるのが F# の再帰関数のコツ（F#に限らないけど）、通常の for ループの場合は、0 から始まって max で終わるけど、再帰関数の場合は max を渡してやって減算していって 0 の時に終了、ってやるとうまくいきます…ってのは何処かにあるだろうか。

たとえば、フィボナッチ数列は

let rec fibo n =
    match n with
    | 0 -> 1
    | 1 -> 1
    | x -> fibo(x-1) + fibo(x-2) 

let fibonacci n =
    for i in 0..n do
        let ans = fibo i
        printf "%d," ans

fibonacci 20

とできるわけで、計算量はさておき、数式のシンプルさをそのまま表せるのが便利ですよね。

ふと、数独（sudoku）をF#で解く…で検索すると、

Sudoku solver in F#
http://www.ffconsultancy.com/dotnet/fsharp/sudoku/

があって、コードがええとよくわからないので、自作しました。
最初は素朴な方法を使っていて、左上から順番に数を入れて計算して、数を入れて計算してという方法をやっていたのですがいっこうに終わらないで、ちょっと考え直し。元ソースをよく見ると、左上からチェックして入れた後は、右に進むという方法をとっているので、ああ、これで十分かということで。

アルゴリズムとして、

1. 左上（0,0）から始める。
2. 0以外だったら、右へ進む。
3. 1 右端に達したら次の行のはじめへ
4. 0だったら、
5. 横一列を調べる
6. 縦一列を調べる
7. 3×3のボックスを調べる
8. 数が埋まれば、埋めて右に進んで再帰

な感じで進めればOK。そういえば、行列計算をちまちま自作していたときに思ったのですが、このあたり「アルゴリズム」≒計算機が得意とする手順は、がぜん生きていますね。普段、リスト構造とかツリー構造とかオブジェクトの構造とか、既存の.NETライブラリで賄えるものが多かったのですが、行列計算、線形代数の分野になると、高速化とか誤差も含めて数学的な解法とは違った古式ゆかしき「アルゴリズム」が生きている分野なのだな、と再確認しています。

さて、このアルゴリズムを直に実装したのが以下です。

let rec solve (m:int[,]) x y = 
    printfn "%d %d" x y
    if x = 9 && y = 8 then
        // 完成
        printfn "%A" m
        exit 0
    elif x = 9 then
        solve m 0 (y+1)
    elif m.[x,y] <> 0 then
        solve m (x+1) y 
    else 
        // 縦/横/boxをチェック
        let mutable v = [|0..9|]
        for i in 0..8 do v.[m.[i,y]] <- 0
        for i in 0..8 do v.[m.[x,i]] <- 0
        let x0 = x/3*3
        let y0 = y/3*3
        for i in 0..2 do
            for j in 0..2 do
                v.[m.[x0+i,y0+j]] <- 0

        let m' = Array2D.copy m
        for i in 0..9 do
            if v.[i] <> 0 then
                m'.[x,y] <- v.[i]
                solve m' (x+1) y

for ループの部分が「えー」とか、match を使ったほうがいいよとか、fold の使い方を知らないとか色々あるわけですが、素直に書いたのでこんな感じ。
ただし、この再帰関数 solve は元に帰ってこなくて、完成したとき（9,8 のとき）には、無理矢理 exit してプログラムを止めています。これが再帰関数を戻すと再び別の解を探し出すので無駄なんですよね。元ネタをみると、例外を発生させているので、それでもいいのですが。末尾再帰をしているわけではないので、スタックを食うのは仕方がなく（バックトラック法だし）、ただし、ツリー構造を探索していく中で各ノードはパラレルに動かせるので探索自体は深くなりそうなときは適当にスレッド化することも可能なのかなと。このあたりは、巨大な行列計算と似た感じになると思っています。

まあ、9×9の数独ならば1秒も経たないで解けてしまうのでスレッド化する意味はありません。が、10倍ぐらいの、100×100（数字は1から100まで）の数独ならばどうなんだろうというのがあります。ウィキペディアを見ると、

Mathematics of Sudoku – Wikipedia, the free encyclopedia
http://en.wikipedia.org/wiki/Mathematics_of_Sudoku

ボックス（bands）が5×5の場合の計算はあるのですが、それ以上はどうなるかわからんという状態みたいです。まあ、パズルとして人間が解ける限界を超えているし、パズルとしての面白味もないので解く意味はないのですが、コンピュータに解かせるときの最適化具合を見るのにはいいかも、と。

数独を解かせるためには、配列の配列を渡してやって、Array2D に直します。

// Input puzzle
let example = [|[|0;0;0;0;6;0;4;0;0|];
                [|0;5;0;0;0;3;6;0;0|];
                [|1;0;0;0;0;5;0;0;0|];
                [|0;4;1;0;0;0;0;0;0|];
                [|0;9;0;0;0;0;0;2;0|];
                [|5;0;2;0;0;0;3;4;0|];
                [|3;0;0;7;0;0;0;0;0|];
                [|0;0;6;5;0;0;0;9;0|];
                [|0;0;4;0;1;0;0;0;0|]|]

let pazzle = Array2D.init 9 9 (fun i j -> example.[i].[j])

solve pazzle 0 0 

答えはこんな感じ

[[2; 3; 7; 1; 6; 8; 4; 5; 9]
 [4; 5; 8; 9; 7; 3; 6; 1; 2]
 [1; 6; 9; 2; 4; 5; 7; 8; 3]
 [6; 4; 1; 3; 5; 2; 9; 7; 8]
 [7; 9; 3; 4; 8; 1; 5; 2; 6]
 [5; 8; 2; 6; 9; 7; 3; 4; 1]
 [3; 1; 5; 7; 2; 9; 8; 6; 4]
 [8; 2; 6; 5; 3; 4; 1; 9; 7]
 [9; 7; 4; 8; 1; 6; 2; 3; 5]]

実は、元ネタの数独は複数回答あって、問題としてよろしくないんですよね。
それはさておき、検算をするコードが以下

let ans = 
    [[2; 3; 7; 1; 6; 8; 4; 5; 9]
     [4; 5; 8; 9; 7; 3; 6; 1; 2]
     [1; 6; 9; 2; 4; 5; 7; 8; 3]
     [6; 4; 1; 3; 5; 2; 9; 7; 8]
     [7; 9; 3; 4; 8; 1; 5; 2; 6]
     [5; 8; 2; 6; 9; 7; 3; 4; 1]
     [3; 1; 5; 7; 2; 9; 8; 6; 4]
     [8; 2; 6; 5; 3; 4; 1; 9; 7]
     [9; 7; 4; 8; 1; 6; 2; 3; 5]]

let A = Array2D.init 9 9 (fun i j -> ans.[i].[j])

let check (m:int[,]) =
    
    for y in 0..8 do
        let mutable v = [|0..9|]
        for i in 0..8 do
            v.[m.[i,y]] <- 0
        if Array.sum v <> 0 then
            printfn "error y=%d" y

    for x in 0..8 do
        let mutable v = [|0..9|]
        for i in 0..8 do
            v.[m.[x,i]] <- 0
        if Array.sum v <> 0 then
            printfn "error x=%d" x

    for x in 0..3..8 do
        for y in 0..3..8 do
            let mutable v = [|0..9|]
            for i in 0..2 do
                for j in 0..2 do
                    v.[m.[x+i,y+j]] <- 0
            if Array.sum v <> 0 then
                printfn "error x=%d y=%d" x y
            
    printfn "OK"

# ああ、そうか。MSTest を使ってもいいんだっけ。どっかにサンプルページがあったはずなので後ほど。

で、bounds が 10×10の数独（全体が100×100）は、どうなるのかなと考えてみたものの…そんな数独の問題はどこにもないわけで、ということは自分で問題を作らないとダメなのですよ。どうせならば、解がひとつしかないパターンの自動作成をしたいわけで、ここは思案のしどころ。ざっと、問題作成の手順を書き下すと、

1. 矛盾なくNxNを敷き詰める。
2. ランダムに1個空白にする。
3. 解が1つしかないことを確認する。再び2へ

な感じかな。3の部分を最適化したいところですが、基本は先の解法コードを拡張するだけなので、後は計算時間の問題かな。

LU分解を F# らしく書き直してみます。
せっかくの F# なのに for ループを使っているのがもったいないし（？）、数式のΣから外れています。ここを sum 関数を使って書き直したいのですが… matrix に対しての Seq.sum とかは要素全体に掛け算してしまうので無理（だと思う）なので、

C#er のためのやさしい再帰入門: いげ太のブログ
http://igeta.cocolog-nifty.com/blog/2011/02/tailcall.html#more

を参考にして、再帰関数化しています。

/// LU分解
let LU2 ( A : matrix ) =
let n = A.NumCols
let L = Matrix.identity n
let U = Matrix.zero n n

// sum(k=0,n) L(i,k)*U(k,j) の計算を関数化
// 余計わかりづらいか？
let sum i j k =
let rec f i j k acc =
if k < 0 then 0.0 elif k > 0 then f i j (k-1) (acc + L.[i,k]*U.[k,j])
else L.[i,k]*U.[k,j]
f i j k 0.0

for j in 0..n-1 do
// Uを計算
for i in 0..j do
U.[i,j] <- A.[i,j] - sum i j (i-1) // Lを計算 for i in j+1..n-1 do L.[i,j] <- (A.[i,j] - sum i j (j-1))/U.[j,j] // 結果を返す (L,U) [/code]

LとUを計算するところで、ΣL.[i,k]*U.[k.j] が同じなので、括りだしてみたのですが、余計わかりづらいかも。最後の for ループに適用するとよいのかもしれません。まあ、この程度の簡単なものならば再帰を使うとかえってややこしいという感じですかね。

ならば、for の二重ループの中で、sum を定義して、ちょっと簡単にしてみたのが次のコードです。

/// LU分解
let LU3 ( A : matrix ) =
let n = A.NumCols
let L = Matrix.identity n
let U = Matrix.zero n n

// sum(k=0,n) L(i,k)*U(k,j) の計算を関数化
// 普通にforループで合計を計算する

for j in 0..n-1 do
// Uを計算
for i in 0..n-1 do
let sum max =
let mutable s = 0.0
for k in 0..max do
s <- s + L.[i,k]*U.[k,j] s if i <= j then // Uを計算 U.[i,j] <- A.[i,j] - sum (i-1) else // Lを計算 L.[i,j] <- (A.[i,j] - sum (j-1))/U.[j,j] // 結果を返す (L,U) [/code]

ループの中で関数が定義できるから、便利ってのがありますね。ループ変数のi,jもそのまま使えるので、sum 関数の引数が少なくて済みます。

お次は、このLU分解を使って、二次元トラスを具体的に解いてみます。

旧タイプなので、連立一次方程式を解くときに逆行列を使って解けばOK…と頭から思い込んでいたのですが、Gaussの消去法を使うと逆行列は出てこなくて「？」となるわけです。上三角化を使うのだから、LU分解でいいんですよね。

LU分解で行列式と逆行列の計算　その３: メモブログ
http://sssiii.seesaa.net/article/379560261.html

のところから、

連立１次方程式 I
http://nalab.mind.meiji.ac.jp/~mk/labo/text/linear-eq-1.pdf

のPDFを読んで、なるほど、LU分解で十分ですね。ってことで自前で実装してみます。

4 LU分解
http://akita-nct.jp/yamamoto/lecture/2004/5E/linear_equations/text/html/node4.html#eq:LU____________

を参考にして式を流用します。最初、手順がうまくわからなくて、自分で 4×4 の行列を作って手計算をした後に、数式に直して、元の式と一致することを確認。なるほど、行単位じゃなくて「列単位」で計算するので、計算する順番は、β11,α21,α31,α41 で計算した後に、β12,β22,α32,α42 と計算していくんですね。もう一度読み直したらそう書いてあるし…

出来上がったコードがこんな感じ。

/// LU分解
let LU ( A : matrix ) =
    let n = A.NumCols
    let L = Matrix.identity n
    let U = Matrix.zero n n
    
    for j in 0..n-1 do
        // Uを計算
        for i in 0..j do
            U.[i,j] <- A.&#91;i,j&#93;
            for k in 0..i-1 do
                U.&#91;i,j&#93; <- U.&#91;i,j&#93; - L.&#91;i,k&#93;*U.&#91;k,j&#93;
        // Lを計算
        for i in j+1..n-1 do
            L.&#91;i,j&#93; <- A.&#91;i,j&#93;
            for k in 0..j-1 do
                L.&#91;i,j&#93; <- L.&#91;i,j&#93; - L.&#91;i,k&#93;*U.&#91;k,j&#93;
            L.&#91;i,j&#93; <- L.&#91;i,j&#93; / U.&#91;j,j&#93;
    // 結果を返す
    (L,U)
&#91;/code&#93;
<p>
次のように行列を設定しておくと
</p>
[code lang="csharp"]
let A = matrix [[1.0;2.0;1.0];
                [2.0;1.0;0.0];
                [1.0;1.0;2.0]]

LU A

こんな感じで、L と U にして返してくれます。

val it : matrix * matrix =
  (matrix [[1.0; 0.0; 0.0]
           [2.0; 1.0; 0.0]
           [1.0; 0.3333333333; 1.0]], matrix [[1.0; 2.0; 1.0]
                                              [0.0; -3.0; -2.0]
                                              [0.0; 0.0; 1.666666667]])

先のページにも書いてあるのですが、Ujj で割ってしまうので、ピボット選択は必須…なのか？（ちょっとよくわからず）。浮動小数点だったら大丈夫な気もするのですが、これはあとで確認してみます。まあ、割り算のところにチェックを入れるのは数値計算の基本なので、なるべく割り算で誤差が出ないようにする方向でよいかと。

ちなみに F# MathProvider を使うと

let (p,L,U) = MathProvider.LinearAlgebra.lu A

として、次なる結果を得られます。

val p : (int -> int)
val U : matrix = matrix [[2.0; 1.0; 0.0]
                         [0.0; 1.5; 1.0]
                         [0.0; 0.0; 1.666666667]]
val L : matrix = matrix [[1.0; 0.0; 0.0]
                         [0.5; 1.0; 0.0]
                         [0.5; 0.3333333333; 1.0]]

MathProvider の LU 分解のコードはもっとシンプルで、こんな感じになっています。うまく交互に計算すると L のほうはほとんど計算せずに済みみたいです。行列数が多くなると swapR のパフォーマンスが問題になりますが、有限要素法以外だったら大丈夫かと。

for i=0 to nA-2 do
let mutable maxi = i // Find the biggest pivot element.
for k=i+1 to nA-1 do
if abs(U.[maxi,i]) < abs(U.[k,i]) then maxi <- k done //let maxi = maxIndex (i+1) (nA-1) (fun k -> abs(U.[k,i]))

if maxi <> i then
swapR U i maxi // Swap rows if needed.
swapR L i maxi // REVIEW can be made more performant.
let t = P.[i]
P.[i] <- P.[maxi] P.[maxi] <- t for j=i+1 to nA-1 do L.[j,i] <- U.[j,i] / U.[i,i] for k=i+1 to nA-1 do U.[j,k] <- U.[j,k] - L.[j,i] * U.[i,k] done U.[j,i] <- 0.0 done done [/code]

sum のところは、もっと F# らしく書き換えられないかな、と思案中です。コードが短くなるしΣを使ったほうが、数式を同じになるので。