ラベル Clojure の投稿を表示しています。 すべての投稿を表示
ラベル Clojure の投稿を表示しています。 すべての投稿を表示

2017/12/19

Binding: Implicit parameters with Clojure

この記事は、Clojureアドベントカレンダーの19日目の記事です。今日は、Clojureの動的スコープの話です。

ScalaのImplicit parameter

ここ最近、業務でScalaがメインなのでScalaしか書いてないのですが、たまに、仕事から帰ってきて家でClojureのコードを書いてると、ClojureにはImplicit parametersの機能が無い事に気づきます。というか、そもそもScala以外だと、Implicit parameterは、HaskellのGHC拡張くらいにしかない?(ちゃんと調べてないですが)と思われますが。。。Implicit parametersとは、関数の引数を暗黙的(implicit)に受け渡す機能です。

例えば、DBConnectionのような関数呼び出し毎に同じ変数を渡すようなシチュエーションが発生した時、Implicit parameterでは、関数間で受け渡す、共有したいパラメータ(ここだとConnectionのような変数)について、呼び出し元と呼び出し先のメソッドでそれぞれimplicitであることを宣言すると、暗黙的に変数の受け渡しをやってくれます。例えば、以下のコード。このコードではcallerという関数からcallee1を呼び出し、callee1からcallee2を呼び出しています。

def caller():Unit = {
  implicit val message = "sample message"
  callee1("from caller")
}

def callee1(text1: String)(implicit message: String):Unit ={
  println(text1)
  callee2(text1 + "2")
}

def callee2(text2: String)(implicit message: String):Unit ={
  println(text2)
  println(message)
}
この時、callee2で使用される変数messageは暗黙的に、caller->callee1->callee2へと受け渡されていきます。明示的に引数を宣言していませんが、callee1及びcallee2の関数呼び出しは、暗黙な引数となったmessageを補完されて呼び出されています。この時、補完される値は、implicitキーワードにより定義された暗黙のパラメータであり、caller内ではimplicitキーワードにより宣言された変数(implicitキーワードのないval messageは普通の変数宣言です。)で、callee1呼び出しで不足しているパラメータを補います。さらに、callee1内では、暗黙的に渡されたパラメータがさらに暗黙的にcallee2に渡されています。

implicit parameterのより正確で詳しい説明は、以下にあります。


実行時のコンテクストなど、種々のパラメータを明示的に記述しないことで、不要な(つまり、例えば、その変数に代入されたオブジェクトが直接使用されない)箇所での余計な記述が減り、その関数で記述すべきロジックが明確になるというメリットがあります。特に、関数の引数で変数引きずり回しプログラミングをやっていると、必須の機能になってくるのではないでしょうか(?)
Implicit parameterも一長一短の機能で手放しで喜べる機能ではないですが、その機能を意識してコントロールしている限りでは、便利な機能と言えます。

ClojureでImplicit parameter的なものを考える

そんな、機能がClojureにあったら良いなと度々思うようになりました。Clojureでは、Lexical Closureのように、環境を捕縛することは出来ても、一々引数として指定せずに別のコンテクストから、値を渡す、といったことはできません。マクロで無理矢理implicitを表す、implicit-letのような構文を作って代入する、みたいなことも出来そうではあります。例えば以下のような感じに。

(implicit-let
 [param1 (get-context1)]
 (function4
    (function1 arg1 arg2)
    (function2 arg1 arg2) ;; 本当はfunction3に3番目のパラメータとしてparam1が入る
    (function3 arg1 arg2))) ;; 本当はfunction3に3番目のパラメータとしてparam1が入る
implicit-letマクロ拡張時に、scopeにバインドされているfunction2とfunction3のスコープを調べ(調べるならdocかsourceマクロあたりでしょうか)、引数が不足していてかつ、パラメータ名(上記で言う所のparam1)が同一であれば、(function2 arg1 arg2)を(function2 arg1 arg2 aram1)に拡張するようなS式変換を行うことで、静的に暗黙のパラメータを受け渡すことが出来そうです。ここに型チェックでも入ればScalaが提供しているImplicit parameterとほぼ同様のものになると思われます。
ただ、これには、マルチメソッドはどうするのか、とか、マクロ内マクロはどのような扱いにするのか、などletのbody部の各名前や構造の静的解決が簡単ではありません。

Binding(Clojureの動的スコープ)

そんなことを考えながら仕事をしていたある日、Implicit Parametersって動的スコープとよく似ているなあ、と思いググったところ、まさにこれ、といった文献がヒットしました。


上記の文献によると、"Implicit Parameters"とは、"Static Types"で"Dynamic Scoping"をするためのものであるようです。タイトルそのまんまですが。だとすると、ClojureでImplicit parameter使いたい問題は、Clojureで動的スコープ使いたい問題と言い換える事ができ、なんとも、都合の良いことにClojureには動的スコープを操作するためのマクロが存在します。
^:dynamicキーワードとbindingです。以下のように動的に値を束縛します。

(def ^:dynamic x "α")
(def ^:dynamic y "β")

(println (format "point1: x = %s, y = %s" x y))

(defn call-xy []
  (println (format "in call-xy: x = %s, y = %s" x y)))

(binding [x "a" y "b"]
  (println (format "after binding: x = %s, y = %s" x y))
  (call-xy))

(println (format "point2: x = %s, y = %s" x y))

(call-xy)
そして、これに対して以下のような出力結果が得られます。

point1: x = α, y = β
after binding: x = a, y = b
in call-xy: x = a, y = b
point2: x = α, y = β
in call-xy: x = α, y = β
bindingにより束縛した、後と束縛中に呼び出した関数呼び出しの中のみ、束縛した値に変更されます。そして、bindingのスコープ外では束縛した値が束縛前のもとの値に戻っています。bindingは、bindingスコープの外側に出ると、もとの値、束縛前の値を再度束縛するという機能があります。これにより引数よる明示的な値渡し以外の方法で値を渡すことを実現します。bindingされた時の呼び出しから呼び出される関数内でのdynamicパラメータを持つ変数のみに動的に値が設定されます。そして、さらに都合の良いことに、スレッドセーフでもあるのです。。。このあたりが、ref(transactionalなclojureのグローバル変数)との使い分けのポイントになると思います。

そして、気になる静的スコープとの競合ですが、以下のようになります。

user=> (binding [x "2"] x)
"2"
user=> (binding [x "2"] (let [x "3"] x))
"3"
user=> (let [x "3"] (binding [x "2"] x))
"3"
基本的に、静的スコープにより遮蔽されるようです。しかし以下のようにvar-get関数により取得することは可能です。
(単にletを飛び越えてトップレベルのdefで定義されたvarを取ってきているだけですが。)

user=>  (let [x "3"] (binding [x "2"] (var-get #'x)))
"2"

Binding(動的スコープ)を使う

動的スコープは副作用ありきの機能というか参照透過性を破壊のような概念ではありますが、Clojureで使われているbindingによる動的スコープなら、例えば、DBコネクションを複数の関数に渡したい、とか、コードの実行コンテクストを複数の関数どうしで引き回したいとか、再帰で書く必要があるものの、一々引数指定したくない(けど、代入する必要がある)などに有用だと思われます。このような、bindingの使い方は、例えば、標準ライブラリのcl-formatや、IOのreadなどで使われています。

例えば、以下のようなコードがある時、関数f, g,でそれぞれg, hの呼び出し以外にparam1, param2が使用されない時、hの変数の値を解決しようとすると、呼び出した先の関数の先の関数で使う値を、一旦、呼び出す関数に渡し、さらにその先での呼び出しでもまた同様に値を渡す必要が出てきます。

(defn h [param1 param2]
  param1, param2を使用)

(defn g [b param1 param2]
  ...
  (h param1 param2)
  ...)

(defn f [a param1 param2]
  ...
  (g abc param1 param2)
  ...)

(f param0 param1 param2)
このようなコードだとかなり厳しさがありますが、例えば動的スコープなら、以下のように解決することが出来ます。

(defn h []
  param1, param2を使用)

(defn g [b]
  ...
  (h)
  ...)

(defn f [a]
  ...
  (g abc)
  ...)

(binding [param1 .... param2 ...]
  (f param0))
また、再帰呼出しを行う以下のような関数がある時、再帰呼出し中、自分自身を呼び出す時に、再帰呼出し中では値が変更されない引数を自分自身の中で共有するために再帰中に値を引き回す必要が出てきます。以下の例では、不変なリストpoints (ex. [[10 20] [30 40] ... ]のような値が入ります)に対して、find-minとweightで、不変なリストを複数の関数内で引き回して使っています。
(ちなみに、以下のコードはの前回の記事から。)

(defn weight [points i j]
  (st/euclidean-distance (points i) (points j)))

(defn +m [x y size]
  (mod (+ x y size) size))

(defn min-path [args]
   (first (sort-by second (filter not-empty args))))

(defn find-min [points i j]
  (if (<= (Math/abs (- i j)) 2)
    [[] 0]
    (let [size (count points)
          i+1 (+m i 1 size)
          j-1 (+m j -1 size)
          [lines w] (find-min points i+1 j)
          skip-i [(cons [i+1 j] lines) (+ w (weight points i+1 j))]
          [lines w] (find-min points i j-1)
          skip-j [(cons [i j-1] lines) (+ w (weight points i j-1))]]
      (letfn [(search-inter [k]
                (let [[lines1 w1] (find-min points i k)
                      [lines2 w2] (find-min points k j)]
                  [(concat [[i k] [k j]] lines1 lines2)
                   (+ w1 w2 (weight points i k) (weight points k j))]))]
        (->> (map search-inter (range (+m i 2 size) (+m j -2 size)))
             (cons skip-i)
             (cons skip-j)
             min-path)))))

(defn triangulate [points]
  (find-min points 0 (- (count points) 1)))
勿論、上記のようなコードは、一旦、Closureを作って解決するという方法も取れますが、そうなると、letfnなどを使い、関数中にさらに不変な変数を共有するためだけに、関数を再定義する必要がでてきます。このようなコードは動的スコープにより以下のように書き直せます。

(def ^:dynamic points [])

(defn weight [i j]
  (st/euclidean-distance (points i) (points j)))

(defn +m [x y size]
  (mod (+ x y size) size))

(defn min-path [args]
   (first (sort-by second (filter not-empty args))))

(defn find-min [i j]
  (if (<= (Math/abs (- i j)) 2)
    [[] 0]
    (let [size (count points)
          i+1 (+m i 1 size)
          j-1 (+m j -1 size)
          [lines w] (find-min i+1 j)
          skip-i [(cons [i+1 j] lines) (+ w (weight i+1 j))]
          [lines w] (find-min i j-1)
          skip-j [(cons [i j-1] lines) (+ w (weight i j-1))]]
      (letfn [(search-inter [k]
                (let [[lines1 w1] (find-min i k)
                      [lines2 w2] (find-min k j)]
                  [(concat [[i k] [k j]] lines1 lines2)
                   (+ w1 w2 (weight i k)) (weight k j)]))]
        (->> (map search-inter (range (+m i 2 size) (+m j -2 size)))
             (cons skip-i)
             (cons skip-j)
             min-path)))))

(defn triangulate [points]
  (binding [points points]
    (find-min 0 (- (count points) 1))))
単にグローバル変数に変数を束縛しているだけである事との違いは、呼び出し元(この場合だとtriangulate)の関数内で、pointsの値を書き換えてる点です。これにより、(当然ですが)triangulateに違うリストを渡した場合に違う結果が返ってくることとなり、あたかも副作用のない通常の関数呼び出しを行っているかのような記述が可能になります。

動的スコープのデメリットとしては、静的にはどこの値が束縛されるのか分からないという問題も発生します。上記のような書き方だと、(トップレベルには記述があるとは言え)どうしてもどこにも所属しない自由変数が関数内に登場してしまうことにはなります。しかし、名前の解決のためだけに値を引数で引き回したり、再帰的な関数をletfnなどで再定義するよりは、シンプルなコードが実現できると思います。

おわりに

前述の文献のように、Implicit parametersが静的型付言語でDynamic Scopingの役割を果たすものだとするならば、その逆もまたありえる、つまり、ScalaでImplicit paramtersで記述されている箇所は、Clojureの動的スコープを使用することで値渡しが可能になるのではないでしょうか。

濫用するとコードが読みにくく扱いづらくなりますが、(これに関してはScalaのimplicit Parameterも同様です)適度に使用することで、無用な引数渡しや、煩雑なデータ渡しが簡潔に記述できて、見通しがよいコードが書けそうです。特にScalaのコードにおいてImplicit parametersで値を渡しているような箇所は、Clojureではdynamic bindingで解決するという方法が取れるのかも知れません。

言語の名前こそClojure(Closure)ですが、動的スコープもまた利用可能で、Lexical Closureのみによって束縛されない変数を含む関数を定義できます。Scalaのimplicit Parameterが良いなと思った時、Clojureではbindingを試してみるのはいかがでしょうか。

参考文献

  1. Let vs. Binding in Clojure - Stack Over Flow
  2. Jeffrey R. Lewis, Mark B. Shields, Erik Meijer, John Launchbury, Iplicit Parameters: Dynamic Scoping with Static Types, 2000, [PDF]
  3. On the Perils of Dynamic Scope - Digital Digressions by Stuart Sierra
  4. Clojure勉強日記(その21 varを使用したスレッドローカルな状態管理 - 夢とガラクタの集積場
  5. Scala 言語仕様 Version 2.8 DRAFT July 13, 2010

2017/12/11

凸多角形の最適三角形分割

はいどうも! この記事は、数学とコンピュータ・アドベントカレンダーの11日目の記事です。 凸多角形の最適三角形分割の話をします。

凸多角形の定義とは、文献1によれば、「多角形の内部のどの2点を結んでも、その2点を両端点とする線分がその多角形の内部に含まれるとき、その多角形は凸である」ということだそうです。「多角形Pが凸多角形であるための必要十分条件は、Pの全ての対角線がPの内部にあるときである」と言い換えられます。

こんな感じに三角分割する

凸多角形を分割する

凸多角形(単調な多角形)を任意の三角形で分割するだけなら、特別なアルゴリズムを考える必要もなさそうです。凸多角形の定義から、多角形の点を一つ選び、その線から多角形の他の点に対して線を引いていくだけでも凸多角形の三角形分割は可能です。というわけで、単に凸多角形を三角形に分解するだけなら簡単なので、凸多角形を三角形分割した時に、分割に使った対角線の総和が最小となるような分割を求めてみます。

この問題を分割統治法で解きます(後述しますが、この方法は、動的計画法によって解くのが一般的です)。

凸多角形を対角線で分割するため、解は必ず多角形の頂点のペアの列として表せます。つまり、頂点のペアの列を求めるわけですが、より正確には、合計の長さが最小となる対角線を表す頂点のインデックスのリストつまり、i番目とj番目の頂点のペアのリストが求めたい解となります。

凸多角形がどんな形で分割されたとしても、必ず、四角形以上であれば、必ず対角線が一本引かれます。つまり、$n$個の頂点を持つ多角形の対角線のうち、最初の1つの対角線を一本、どこかに引く必要があります。
そのような直線が一本決まると、残った部分の多角形の三角分割は、元の多角形の部分多角形の三角分割の問題となります。(三角分割なので、当然、対角線どうしが交差することはない為、決まった対角線によって分割された2部分多角形の三角分割の問題になります。)

まず始めに、凸多角形の全体を考えるのではなく、凸多角形の部分多角形の分割について考えてみたいと思います。話を簡単にするために、頂点を$v_i$(凸多角形の$i$番目の頂点)の$i$番目〜$j$番目間の分割を表す関数を$L[i, j]$とします。$L[i, j] (|i - j| < 2)$となる場合については、分割出来ないため、分割なし、分割を表す対角線の長さは0とします。また、インデックスですが、頂点$n$の多角形に対し$i+1$や$j-1$などは$\pmod{n}$で循環するものとします。

こんな風に分割

部分多角形の分割を表す$L[i, j]$に対して、部分多角形内の頂点$v_i, v_j$を含むような三角形を考えます。今求めようとしている部分多角形の性質上、頂点$v_i$と$v_j$間が別の対角線によって遮られる事はあり得ません。したがって、$v_i$, $v_j$をともうひとつの頂点から成る三角形によって分割される筈です。そのような三角形の頂点は、$v_i+1 $cdots v_j-1$の中から必ず1つ選ばれる事になります。(選ばれなかった場合、三角分割として成り立たないため)
このような三角形の選び方は3通りあります。

  1. $(v_i, v_j, v_{j-1})$で三角形が構成されるパターン
  2. $(v_i, v_{i+1}, v_j)$で三角形が構成されるパターン
  3. $(v_i, v_k, v_j)$で三角形が構成されるパターン(ただし、$i+1 < k < j-1$)

3通りの分割

1, 2, 3は3の形式でひとまとめに出来そうですが、1, 2は、引かれる対角線が1つであるのに対して、3では2つ対角線が引かれるという点が異なります。そして、それぞれ、以下のように3通りの部分多角形の三角分割を行った結果と足し合わせる形で表せます。

  1. $L[i, j] = (i, j-1)$ 及び $L[i, j-1]$
  2. $L[i, j] = (i+1, j)$ 及び $L[i+1, j]$
  3. $L[i, j] = (i, k), (k, j)$ 及び $L[i, k]$ と $L[k, j]$

上記の条件をまとめると、三角分割した部分多角形$L[i, j]$の対角線の長さ$L'[i,j]$は、

$L'[i,j] = 0 (|i - j| < 2)$
$L'[i,j] = min \bigl\{$
    $L'[i + 1, j] + w[i + 1, j],$
    $L'[i, j - 1] + w[i, j - 1],$
    $\displaystyle \min_{k = (i+2) \dots (j-2)} \{ L'[i, k] + L'[k, j] + w[k, j] + w[k, j] \bigr\} \}$

であると言えます。$w[i, j]$は、頂点$v_i$, $v_j$間の重さ、$min { 〜 }$は、与えられた集合の中で最も小さい物を還す関数。添字の$k$も頂点を表すインデックスです。これにより、再帰的な関数によって線分$L'[i, j]$の長さが定義出来ることになりました。

ここでもう一度、元の問題、$n$個の頂点を持つ凸多角形の線分の長さが最小となる三角分割ですが、上記の部分多角形の三角分割$L[i, j]$の$L[0, n-1]$のケースであると見なせます。$L[0, n-1]$を解くことで、$n$個の頂点を持つ凸多角形の三角分割が可能になるのです。
ちなみに上記の$L'$ですが、コードに落とす時は、重さだけでなく、$L[i, j]$の要素のリストについても考慮する必要があります。

実装

というわけでClojureのコードです。
(実行結果を可視化するのが一番簡単だったため、ClojureとIncanterを使いました。)
(def weight st/euclidean-distance)

(defn +m [x y size]
  (mod (+ x y size) size))

(defn min-path [args]
   (first (sort-by second (filter not-empty args))))

(defn find-min [points i j]
  (if (<= (Math/abs (- i j)) 2)
    [[] 0]
    (let [size (count points)
          i+1 (+m i 1 size)
          j-1 (+m j -1 size)
          [lines w] (find-min points i+1 j)
          skip-i [(cons [i+1 j] lines) (+ w (weight (points i+1) (points j)))]
          [lines w] (find-min points i j-1)
          skip-j [(cons [i j-1] lines) (+ w (weight (points i) (points j-1)))]]
      (letfn [(search-inter [k]
                (let [[lines1 w1] (find-min points i k)
                      [lines2 w2] (find-min points k j)]
                  [(concat [[i k] [k j]] lines1 lines2)
                   (+ w1 w2
                      (weight (points i) (points k))
                      (weight (points k) (points j)))]))]
        (->> (map search-inter (range (+m i 2 size) (+m j -2 size)))
             (cons skip-i)
             (cons skip-j)
             min-path)))))

(defn triangulate [points]
  (find-min points 0 (- (count points) 1)))
weightは、重さ(長さ)を表す関数(incanterの関数)、前述の$w[i, j]$に相当するものです、+mは$\pmod{n}$を計算する関数、find-minが$L[i, j]$、triangulateは、$L[0, n-1]$を呼び出すだけの関数です。

find-minは、頂点のリストを引数に受け取り、頂点列のインデックスと対角線の長さの合計のペア[[$(i, j) \cdots$] $'L[i, j]$]を返します。skip-iとskip-jがそれぞれパターン1, 2の箇所で、(map search-inter 〜)で書かれている所がパターン3になります。最後にmin-pathで最も小さい対角線を求めています。ちなみに、紛らわしいかもしれませんが、スペースが入っていないi+1とj-1は変数です。

実行結果

こんな感じになります。一番上の画像と同じですが。

おわりに

さて、出来上がった関数ですが、計算の重複が多く、見事にメモ化再帰の実装対象になっています。文献1によれば、動的計画法でメモ化しながら、分割統治法的に問題を解くようです。(多角形の頂点の数を増やしていくと、めっちゃ重くなる。。。)つまり、メモ化再帰をテストするためのサンプルに使用できるということであります。メモ化再帰の実装サンプルに、凸多角形の最適三角分割を使用してみてはいかがでしょうか。
それでは!

全ソースコード

(require '[incanter.core :as in]
         '[incanter.stats :as st]
         '[incanter.charts :as ch])

;; ---- generate polygon: Jarvis's March to find convex hull

(defn left-most-point [[x1 y1 :as p1] [x2 y2 :as p2]]
  (if (< x1 x2) p1 p2))

(defn ex-product
  "exterior product"
  [[x1 y1 :as p1] [x2 y2 :as p2]]
  (- (* x1 y2) (* x2 y1)))

(defn -vec [[x1 y1 :as p1] [x2 y2 :as p2]]
  [(- x1 x2) (- y2 y1)])

(defn min-angle [current-p point-b point-c]
  (let [v (ex-product (-vec current-p point-b) (-vec current-p point-c))]
    (if-not (or (= point-b current-p)
                (< 0 v)
                (and (< -0.001 v 0.001)
                     (< (st/euclidean-distance current-p point-b)
                        (st/euclidean-distance current-p point-c))))
      point-b point-c)))

(defn jarvis-march
  "find convex hull, points = [[x1, y1], ...]"
  [points]
  (let [start-p (reduce left-most-point points)]
    (loop [current-p start-p hull-ps [start-p]]
      (let [next-p (reduce (partial min-angle current-p) points)]
        (if (and (= next-p start-p))
          hull-ps
          (recur next-p (cons next-p hull-ps)))))))

(defn generate-sample-set []
  (letfn [(generate-by [m]
            (st/sample-normal
             (+ 100 (rand-int 100)) :mean m :sd (+ 0.15 (rand 0.1))))]
    (map list (generate-by (dec (rand 2))) (generate-by (dec (rand 2))))))

(defn generate-sample []
  (let [ps (generate-sample-set)
        result (jarvis-march ps)]
    (vec (map (fn [[x y]] [(int (* x 100)) (int (* y 100))]) result))))

;; --- triangulation -----------------------
(def weight st/euclidean-distance)

(defn +m [x y size]
  (mod (+ x y size) size))

(defn min-path [args]
   (first (sort-by second (filter not-empty args))))

(defn find-min [points i j]
  (if (<= (Math/abs (- i j)) 2)
    [[] 0]
    (let [size (count points)
          i+1 (+m i 1 size)
          j-1 (+m j -1 size)
          [lines w] (find-min points i+1 j)
          skip-i [(cons [i+1 j] lines) (+ w (weight (points i+1) (points j)))]
          [lines w] (find-min points i j-1)
          skip-j [(cons [i j-1] lines) (+ w (weight (points i) (points j-1)))]]
      (letfn [(search-inter [k]
                (let [[lines1 w1] (find-min points i k)
                      [lines2 w2] (find-min points k j)]
                  [(concat [[i k] [k j]] lines1 lines2)
                   (+ w1 w2
                      (weight (points i) (points k))
                      (weight (points k) (points j)))]))]
        (->> (map search-inter (range (+m i 2 size) (+m j -2 size)))
             (cons skip-i)
             (cons skip-j)
             min-path)))))

(defn triangulate [points]
  (find-min points 0 (- (count points) 1)))

;; --- plot --------------------
(defn plot-triangles [points lines]
  (let [ys (concat points [(first points) (last points)])
        poly (-> (ch/scatter-plot (map first points) (map second points))
                 (ch/add-lines (map first ys) (map second ys) :auto-sort false))
        graph (reduce #(ch/add-lines %1 (map first %2) (map second %2) :auto-sort false)
                      poly
                      (partition 2 lines))]
    (in/view graph)))

(defn point-selector [points indices]
  (apply concat (map (fn [[start end]] [(points start) (points end)]) indices)))

(defn show [points]
  (let [[indices score] (triangulate points)
        xs (point-selector points indices)]
    (plot-triangles points xs)
    points))

参考文献

  1. 浅野 哲夫, 計算幾何―理論の基礎から実装まで (アルゴリズム・サイエンスシリーズ―数理技法編), 共立出版, 2007

2016/03/02

Incanterのグラフのフォントが文字化けする問題(Clojure)

ClojureのIncanter, 便利なライブラリですが, 日本語などを使った場合, 表示されるフォントが文字化けするという問題があります.

ちなみに, こんな感じ.

元は, 次のようなデータ.
(def top-family-names
  [["佐藤" 313079]
   ["鈴木" 268103]
   ["高橋" 226090]
   ["田中" 205560]
   ["伊藤" 171104]
   ["山本" 167767]
   ["渡辺" 166831]
   ["中村" 163789]
   ["小林" 162481]
   ["加藤" 133341]])
※データは, 同姓同名探しと名前ランキングより.

勿論, 日本語を使わないとか, ローマ字表記にするいう選択肢もありますが, 日本語の処理を行う場合, 各軸の値が日本語の場合は避けて通れません.

上記の場合, 日本の苗字トップ10を表示させたつもりなのですが, 凡例は勿論, 各データの内容もわからなくなってしまいます. 唯一正しく表示されているのは, 数値のみです.

調べると, Google Groupsに解決策がありました.
https://groups.google.com/forum/#!topic/incanter/s_5OR-OP9gw

Fontを変えれば良いようです. というわけでMSゴシックへフォントを変更しました.
(use '(incanter core charts))
(import 'java.awt.Font)

(defn view-incanter [xy-values]
  (let [xs (map first xy-values)
        ys (map second xy-values)
        chart (bar-chart xs ys :vertical false
                         :x-label "family names" :y-label "number of people")
        _ (-> chart .getCategoryPlot .getDomainAxis
              (.setTickLabelFont (new Font "MS ゴシック" Font/BOLD 10)))]
    (save chart "family-names.png")))

(view-incanter top-family-names)
という形です.

これを描写すると,

こんな感じで, 軸ラベルは変更できませんが, 縦軸の値は日本語化できました.
(-> chart .getCategoryPlot .getDomainAxis
    (.setTickLabelFont (new Font "MS ゴシック" Font/BOLD 10)))
の部分で縦軸のフォントを変えることができます. 横軸を変える場合は, .getDomainAxis→.getRangeAxisで変えられると思われます.

しかし, ここまで来ると, 軸ラベルもすべて一括して日本語で表示できるようにしたくなります. そして, 調べていて分かったのですが, Incanterのプロットは, JavaのJFreeChartを使ってプロットしているようです. そして, JFreeChartには, グラフが文字化けした場合のおまじないがありました.

JFreeChart – グラフ中の日本語が文字化けしないようにするには - TK FACTORY

そしてそれは,
ChartFactory.setChartTheme(StandardChartTheme.createLegacyTheme());
だそうです. というわけで,
(import 'org.jfree.chart.ChartFactory)
(import 'org.jfree.chart.StandardChartTheme)
(ChartFactory/setChartTheme (StandardChartTheme/createLegacyTheme))
を描画前に入れると, 文字化けが修正されます(ただし, メソッド名にあるように, テーマも変更されてしまいますが).
これは, 次のようにプロットされます.

見た目を気にしないなら, とりあえず, これでOKですね.

2015/08/11

Convex Hull(凸包)を求める(Jarvis's March, Quickhull, Clojure)

 凸包を求めるアルゴリズムを2つ(Jarvis's MarchとQuickhull)調べたので, そのメモ.

 どちらもアルゴリズム的には, シンプルですが, Quickhullの方は, 理解するのに少し時間がかかりました. Jarvis's Marchは, ギフト包装法という呼び方の方が一般的なようです. どちらも基本的には幾何学の問題で, 高校数学の基礎的な内容と外積を駆使します.

 画像としては以下のような感じ. 点の集合を囲む凹みの無いつつみが凸包です. 凸包は, 一部の線形計画法の問題などで使われますが, 抽象解釈などにも凸多面体などとして, 応用されているという点が面白いですね.

以下では, 2次元について考えます.
 Jarvis's Marchから.

 求めるのは, 凸包を構成する点列です. そしてその点列に含まれるのが自明な点は, 最左, 最右と, 最上, 最下の点でしょう. Jarvis's Marchのアルゴリズムは, この4つのどの点でもいいのですが, とりあえずとっかかりとして, 最左の点(x座標が一番小さい点)から始めます.

1. 現在の点から他の全ての点へのベクトルの角度を求め, 角度が最小の点を凸包の点列に追加し, これを次の点とする.
2. 始めた点へ戻ってくるまで1を続ける.

 というのがアルゴリズムで, 非常に簡単でした. 1で角度が同じ場合は, 現在の点から距離が近い方を選びます. 角度を求める手法は何でも良いらしいですが, Wikipediaに従い, 外積を用いることにしました. (外積の意味については, 参考 : 基礎の基礎その1 内積と外積の使い方 - ○×つくろーどっとコム)

Wikipediaの擬似コードを無理やり関数型言語へ押し込むと以下のようになりました. 多分, 手続き型言語よりも整理出来ていると思います. (incanterを使っています.)
(require '[incanter.core :as in]
         '[incanter.stats :as st]
         '[incanter.charts :as ch])

;; -- Jarvis's March to find convex hull

(defn left-most-point [[x1 y1 :as p1] [x2 y2 :as p2]]
  (if (< x1 x2) p1 p2))

(defn ex-product
  "exterior product"
  [[x1 y1 :as p1] [x2 y2 :as p2]]
  (- (* x1 y2) (* x2 y1)))

(defn -vec [[x1 y1 :as p1] [x2 y2 :as p2]]
  [(- x1 x2) (- y2 y1)])

(defn min-angle [current-p point-b point-c]
  (let [v (ex-product (-vec current-p point-b) (-vec current-p point-c))]
    (if-not (or (= point-b current-p)
                (< 0 v)
                (and (< -0.001 v 0.001)
                     (< (st/euclidean-distance current-p point-b)
                        (st/euclidean-distance current-p point-c))))
      point-b point-c)))

(defn jarvis-march
  "find convex hull, points = [[x1, y1], ...]"
  [points]
  (let [start-p (reduce left-most-point points)]
    (loop [current-p start-p hull-ps [start-p]]
      (let [next-p (reduce (partial min-angle current-p) points)]
        (if (and (= next-p start-p))
          hull-ps
          (recur next-p (cons next-p hull-ps)))))))

;; --- visualization --- (use incanter)
(defn plot-convex-hull [points hull-ps]
  (let [hull-ps (cons (last hull-ps) hull-ps)]
    (ch/add-lines
     (ch/scatter-plot (map first points) (map second points))
     (map first hull-ps)
     (map second hull-ps)
     :auto-sort false)))

;; --- test ---
(defn generate-sample-set []
  (letfn [(generate-by [m]
            (st/sample-normal
             (+ 100 (rand-int 100)) :mean m :sd (+ 0.15 (rand 0.1))))]
    (map list
         (generate-by (dec (rand 2)))
         (generate-by (dec (rand 2))))))

(defn trial []
  (let [ps (generate-sample-set)
        result (jarvis-march ps)]
    (in/view (plot-convex-hull ps result))))
 このアルゴリズムは, 計算量が大きそうですが, 自明な凸包の輪郭とならない点を削ることで, 計算量を減らすことは可能なようです. 例えば, 前述の上下左右4点の内側の点は全て凸包の輪郭には含まれないのでこれらの点は無視できますね.

 意外に気づかずにハマったのは, incanterのchart/add-linesが引数として渡されたxy各座標のリストを勝手にソートしてくれることです. :auto-sortをfalseに設定しています.

次に, Quickhullのアルゴリズム. これは, Quicksort風の分割統治法によるアルゴリズムで, やはり, このアルゴリズムも自明な点から始めます.


 最左と最右の点を初期値にします. 上記のWikipediaの記事とは書き方が少し異なりますが, 内容としては同じ(はず)です.

1. 最左/最右の点を求める.
2. 前述の二点(p-left-most, p-right-most)から線分を求め, これで点の集合を二分し, 上半分についてのみ考え, 3.に上記二点(p-left-most, p-right-most)と上半分の集合を渡す.

3. (関数. 引数は 二点(p1, p2)と点の集合)
3-1. 線分から最も遠い点fpを求める
3-2. これを凸包の点の一つとし, 線分とfpからなる三角形の内側の点を除外.
3-3. 三角形の外側に残る点のうち, fpより左側の点の集合と, fpより右側の点の集合に分ける.
3-4. 3.へp1, fpと左半分を渡し, p1〜pf間の凸包の部分列(輪郭の点列)を求める.
3-5. 3.へfp, p2と右半分を渡し, pf〜p2間の凸包の部分列を求める.
3-6. 3-4, pf, 3-5の結果を合わせて, p1〜p2間の凸包の部分列として返す.

4. 下半分も3.へ渡す.
5. 上半分, 下半分について求めた結果と, 最左最右の点を纏める => 凸包を構成する点列.

 Jarvis's Marchと比較して, 少し手間がかかりますが, こんなイメージ.

 全体として分割統治法で, 再帰的に求めていくと, 答えが求まるという流れになります. 線分からの距離は, y=ax+b(直線)からの距離を求めればOK. 三角形に含まれるかどうかも, 外積から求まります.
(require '[incanter.core :as in]
         '[incanter.stats :as st]
         '[incanter.charts :as ch])

(defn get-ab
  "calc a,b for y = ax + b from p1, p2"
  [[x1 y1] [x2 y2]]
  (let [a (/ (- y1 y2) (- x1 x2))]
    [a (- y2 (* a x2))]))

(defn ex-product
  "exterior product"
  [[x1 y1 :as p1] [x2 y2 :as p2]]
  (- (* x1 y2) (* x2 y1)))

(defn -vec
  "p1 - p2 (p1 & p2 is vector)"
  [[x1 y1 :as p1] [x2 y2 :as p2]]
  [(- x1 x2) (- y2 y1)])

(defn in-triangle?
  [[x1 y1 :as p1] [x2 y2 :as p2] [x3 y3 :as p3] [tx ty :as tp]]
  (let [from-a (ex-product (-vec p3 p1) (-vec p1 tp))
        from-b (ex-product (-vec p1 p2) (-vec p2 tp))
        from-c (ex-product (-vec p2 p3) (-vec p3 tp))
        abc [from-a from-b from-c]]
    (or (reduce #(and %1 %2) (map pos? abc))
        (reduce #(and %1 %2) (map neg? abc)))))

(defn furthest-point [[x1 y1 :as xy1] [x2 y2 :as xy2] points]
  (let [[a b] (get-ab xy1 xy2)
        sqa2b2 (. Math sqrt (+ (. Math pow a 2) (. Math pow b 2)))]
    (letfn [(get-dist [[x0 y0]] (/ (. Math abs (+ (* a x0) (* -1 y0) b)) sqa2b2))]
      (reduce #(if (< (get-dist %1) (get-dist %2)) %2 %1) points))))

(defn find-points [p1 p2 points] ;; find points between p1 & p2
  (if (empty? points)
    []
    (let [[xf yf :as pf] (furthest-point p1 p2 points) ;; furthest point
          rest-points (remove #(or (in-triangle? p1 p2 pf %) (= pf %)) points)]
      (concat (find-points p1 pf (filter #(< (first %) xf) rest-points))
              [pf]
              (find-points pf p2 (filter #(<= xf (first %)) rest-points))))))

(defn quick-hull [points]
  (let [p-left-most (reduce #(if (< (first %1) (first %2)) %1 %2) points)
        p-right-most (reduce #(if (< (first %1) (first %2)) %2 %1) points)
        [a b] (get-ab p-left-most p-right-most)]
    (letfn [(upper? [[x y]] (< (+ (* a x) b) y))]
      (concat [p-left-most]
              (find-points p-left-most p-right-most (filter upper? points))
              [p-right-most]
              (->> (filter #(not (upper? %)) points)
                   (find-points p-left-most p-right-most)
                   reverse)))))

;; --- visualization --- (use incanter)
(defn plot-convex-hull [points hull-ps]
  (let [hull-ps (cons (last hull-ps) hull-ps)]
    (ch/add-lines
     (ch/scatter-plot (map first points) (map second points))
     (map first hull-ps)
     (map second hull-ps)
     :auto-sort false)))

;; --- test ---
(defn generate-sample-set []
  (letfn [(generate-by [m]
            (st/sample-normal
             (+ 1000 (rand-int 500)) :mean m :sd (+ 0.15 (rand 0.1))))]
    (map list
         (generate-by (dec (rand 2)))
         (generate-by (dec (rand 2))))))

(defn trial []
  (let [ps (generate-sample-set)
        result (quick-hull ps)]
    (in/view (plot-convex-hull ps result))))
filterやconcatの使い方が, Haskell系の偽Quicksortにそっくりですね.

参考文献
Quickhullアルゴリズム Project ASURA

2015/07/28

loop/recurで相互再帰(Clojure)

 Clojureのtrampolineを使った相互再帰の例として, even/odd関数が有りますが, しかし,  簡単な相互再帰なら, trampolineよりloop-recur内に自分でディスパッチさせる処理を書いたほうがわかりやすいような気がします.

 例えば, even/oddならこんな感じ.
(defn even-odd? [even-or-odd n]
  (loop [status even-or-odd n n]
    (cond
     (= status :even)
     (if (= n 0)
       true
       (recur :odd (dec n)))
     (= status :odd)
     (if (= n 0)
       false
       (recur :even (dec n))))))
こう書いておいて,
user> (even-odd? :odd 9)
と呼び出すとか. cond節が冗長になりますが, declare, defnと書いて, 無名関数で囲って, trampolineで待ち受けさせるのと大差ないかなと....思ったのですが, 比べてみると, やっぱり読みにくい. trampoline版のほうがシンプルですね.
(declare my-odd?)

(defn my-even? [n]
  (if (= n 0)
    true
    #(my-odd? (dec n))))

(defn my-odd? [n]
  (if (= n 0)
    false
    #(my-even? (dec n))))
適当に測ってみると少し早くなった程度.
user> (time (trampoline my-even? 10000000))
"Elapsed time: 2127.23756 msecs"
true
user> (time (even-odd? :even 10000000))
"Elapsed time: 1781.678006 msecs"
true
しかし, 相互再帰を書く時の選択肢としてはありかなと思ったりしてます.

2015/07/15

dorothy(Graphviz)で木構造(構文解析木)を可視化

 dorothyは, GraphvizのClojure向けのラッパーです. シーケンスによる(Graphvizの)DOT言語を用いて, グラフ構造を記述すると, そこから, そのグラフを可視化します.
 dorothyの基本的な使い方については, dorothyのREADMEか, 以下のページの説明が詳しいです.


 木構造もグラフの一種なので, パーサ(opennlp)によって生成された構文解析木を表示させてみることにします. Clojureのopennlpについては, にも触れましたが, 今回は, これにdorothyを合わせて使ってみます. 長文を構文解析し, dorothyで可視化すると, 以下のようになります.

 opennlpからdorothyへの変換は, 基本的に, 木構造のリンク関係を[src dst]のリストへ変換すればOKで, 書いてみると, こんな感じになります.
(require '[dorothy.core :refer :all];; [dorothy "0.0.6"]
         '[opennlp.treebank :refer :all]) ;; [clojure-opennlp "0.3.2"]

(def treebank-parser
  (make-treebank-parser "en-parser-chunking.bin"))

(defn gensym- [tag]
  (keyword (gensym (str (name tag) "-"))))

(defn ungensym- [tag]
  (first (clojure.string/split (name tag) #"-")))

(defn tree->links [tree]
  (if (string? tree)
    [(gensym- tree)]
    (let [{tag :tag chunk :chunk} tree
          children (map tree->links chunk)
          top (keyword (gensym- tag))]
      (concat [top]
              (map (fn [x] [top (first x)]) children)
              (apply concat (map rest children))))))

(defn treeviz [en-text]
  (letfn [(remove-sym [node]
            (subgraph [[:node {:label (ungensym- node)}] node]))]
    (let [link-pairs
          (->> (treebank-parser [en-text]) first make-tree tree->links)
          subnodes
          (->> link-pairs flatten distinct (map remove-sym))
          params
          [(graph-attrs {:size "5, 5"}) (node-attrs {:fontsize 10})]]
      (->> link-pairs (concat params subnodes) digraph dot show!))))
 Graphviz側では, 親ノードから, 子ノードへのリンクの表現のリストで, 木構造を表すわけですが, ノードを表す識別子(POS(品詞)や名詞)が重複する可能性があるのでgensymなどを使ってノードのアイデンティティを確立します. 再帰的にリンクを表すペアのリストを作って上に投げるという処理を繰り返すだけで生成できるdorothyの記法は, かなりお手軽です.
 dorothyには, 以下のような[srs dst]で表されるペアのリストを渡します.
([:TOP-2782 :S-2783] [:S-2783 :NP-2784] [:S-2783 :VP-2791] [:S-2783 :.-2801] [:NP-2784 :DT-2785] [:NP-2784 :NN-2787] [:NP-2784 :NN-2789] [:DT-2785 :This-2786] [:NN-2787 :parse-2788] [:NN-2789 :tree-2790] [:VP-2791 :VBZ-2792] [:VP-2791 :NP-2794] [:VBZ-2792 :visualizes-2793] [:NP-2794 :DT-2795] [:NP-2794 :NN-2797] [:NP-2794 :NN-2799] [:DT-2795 :a-2796] [:NN-2797 :parse-2798] [:NN-2799 :tree-2800] [:.-2801 :.-2802])
 これだけでも, 木構造は可視化できますが, これだと, シノニム部分がくっついて見づらいため, 各ノードのラベルは数値の部分を抜いた名前で表すことにします. 各ノードにラベル(ノードとして表示される名前)をつけるには, subgraphを使います. その他, 表示用のオプションがいくつか設定できますが, 今回は, フォントサイズと, 表示される画像のサイズを設定しました.
user> (treeviz "This tree visualizes a parse tree .")
こんな感じで実行するとして,
 上図の木構造が出てきます. 自然変換言語処理の教科書に乗ってそうな図ですね.
 TOPの下のSは, Statement(文)のSです. SubjectのSではないです. 主語は, NPで表されている左側の部分木です.

 ちなみに, 要素の順番は重要で, 上記の例だと, 綺麗に元の文の単語列が左から順に表示されていますが, リストの順序が異なれば, 可視化の結果も異なり, 木構造としては成立していても変な形で表示されたりもします.

2015/05/27

Clojure.core.logicでレイトン教授の問題を解く

 Alloy Analyzerについて教えてもらった時, Alloyなら, なぞなぞを解けるという話になって, でもそれってPrologでもほとんど同じことができるんじゃないかと思ったのですが, とりあえず, 以下のページにAlloyのなぞなぞSolverの実装があります.


 Prologで同じ問題を解いてみようと思ったのですが, ほとんど書いたことのないPrologで一からやるより, 身近なClojureでやってみたほうが簡単だろうと思って, Clojure.core.logicを使って, 上記のページの問題を1問解いてみました.


 Clojure.core.logicは, Clojureの論理プログラミング用のライブラリで, 主な使い方は, 以下のページに載っています.


 解いたのは, レイトン教授のNo.50の問題で,

 牛が5頭(それぞれa~eとする)いて, うち, 3頭は, ウソしか言わない種であり, 残りの2頭は, ホントのことしか言わない種であるそうです. このとき, 各牛は, 次のように発言しています.

a : dはウソをついている.
b : cはウソをついている.
c : aはウソをついていない.
d : eはウソをついている.
e : bはウソをついている.

 さて, このとき, 嘘をついているのは誰でしょうという問題.

 というわけで, 小学生の頃よくやった推理算ですが, 元の記事では, この問題をAlloyを使ってコンピュータに考えさせようという話でした.

 core.logicで書くとAlloyで書いた場合よりも少し冗長になっている気がしました(マクロを駆使すれば解決できるのでしょうが/また, 書きなれていないことも原因かもしれません).

 core.logicのオペレータは非常にシンプルで, ==, conde, fresh, consoが主でした. それぞれ, 同値であるという述語, 条件分岐, 変数の束縛, リストの構築を表す述語(?)です.run*でUnificationを実行し, 各変数のmost general unifier(最も一般化された単一子, 要するに述語を満たすような最も抽象的な値)を求めます. 結果は, 束縛した変数に対応する値のリストとして返されます. 答えの候補が複数ある場合は, 一つの変数に対して, リストが返ってきます. run*の次のシンボルのベクタが, 求めたい変数のリストで, その後ろのS式が, 制約条件となります.

 問題の発言を述語にするのは, 簡単で, 「aがdはウソをついている」 という発言は, 「a = "正直"かつb = "嘘つき"」または, 「a = "嘘つき"かつb = "正直"」といった感じで, 発言者について場合分けすれば, きれいに分解できます. もっともAlloyのコードでは, 逆のパターンの記述はする必要がないようですが.

 5頭のうち, 2頭が正直であるという条件の制約の定義は手間取りました. もう少し, 上手い書き方があるような気がするのですが... 以下のソースコードのhonest-cows-inという関数がそれに相当します. しかし, このような書き方が適切なのかよくわかっていません.

以下がソースコード.
(use 'clojure.core.logic)

(defn honest-cows-in [ls n]
  (fresh [xs]
    (conde
     [(== ls []) (== n 0)]
     [(conso 'honest xs ls)
      (honest-cows-in xs (dec n))]
     [(conso 'liar xs ls)
      (honest-cows-in xs n)])))

(defn layton50 []
  (run* [a b c d e]
     (fresh [ls]
       (conso a [b c d e] ls)
       (honest-cows-in ls 2))
     ;; A said D is a liar. (A is honest and D is liar) or (A is liar and D is honest).
     (conde [(== a 'honest) (== d 'liar)]  [(== a 'liar) (== d 'honest)])
     ;; B said C is a liar.
     (conde [(== b 'honest) (== c 'liar)]  [(== b 'liar) (== c 'honest)])
     ;; C said A is honest.
     (conde [(== c 'honest) (== a 'honest)]  [(== c 'liar) (== a 'liar)])
     ;; D said E is a liar.
     (conde [(== d 'honest) (== e 'liar)]  [(== d 'liar) (== e 'honest)])
     ;; E said B is a liar.
     (conde [(== e 'honest) (== b 'liar)]  [(== e 'liar) (== b 'honest)])))
実行してみると,
user> (layton50)
([liar honest liar honest liar])
 それぞれ, run*で束縛した[a b c d e]という各牛のパラメータで, aとc, eが嘘つきだということが判明しました. というわけで, 答えは出ましたが, 正解かどうかは知りません. 少なくとも, clojure.core.logicでなぞなぞSolverは記述できるようです.

2015/04/20

Clojureで単純パーセプトロンを書いてみた.

 機械学習を学習しているので, とりあえず, 単純パーセプトロンを書いてみました. 理論と比べると, 実装は非常に簡単(というかシンプル)でした.

 機械学習系のアルゴリズムは, immutableなClojure向きではないような気もしますが, 単純パーセプトロンなら, 大量の配列のupdateの作業がないので, 大丈夫でした.

 参考文献は, 以下のブログ中のスライド.

単純パーセプトロンをPythonで組んでみる - 銀座で働くデータサイエンティストのブログ

 単純パーセプトロンは, 線形分離可能な分類問題において, 識別関数を求めるアルゴリズムです.

 求める分離関数を g(x) = w . + bとすると, アルゴリズムとしては(参考文献のスライドにもありますが),
R = maxi(||xi||)
repeat
  modified? ← false  for i = 0 to (サンプルの総数) - 1 do
    if (labeli * (w . xi) + b < 0) then
      modified? ← true
      ww + ((学習係数) * labeli) xi
      b ← b + (学習係数) * labeli * R
    end if
  end for
until (not modified? or ループ上限に達する.)

xi : i番目のサンプル(教師データ)を表すベクトル(2次元なら(x, y)など)
labeli : i番目のサンプルのラベル(どちらのクラスタに属するか, -1 or 1)
(学習係数) : 1未満の値

 で良いようです.
(2014/11/20 : modified?と終了判定が間違っていたので修正)

 ソースコードは以下のようになりました.(incanterを使ってます.)  というわけで, 実行してみると, 以下のようになりました.


2015/03/24

Clojureのdefun, パターンマッチ付きの関数定義

 defunというと, Common LispやEmacs Lispの関数定義を連想しますが, Clojureでもdefunを使って関数定義ができるようです. さっき, 偶然発見したのですが.

 このdefun, Clojureの普通の関数定義であるdefnとの違いは, パターンマッチ付きの関数定義ができることです. GitHub中の説明には,
A macro to define clojure functions with parameter pattern matching just like erlang or elixir.
 とあります.

 Clojureのパターンマッチングといえば, core.matchかdefmulti/defmethodのマルチディスパッチですが, 関数の引数に対して, そのままパターンを記述できるので, defn+matchの組み合わせよりもシンプルになります.

 .lein/profile.clj (or project.clj)の:pluginsに[defun "0.2.0-RC"]を追加してとりあえず, quicksortもどきを書いてみると,
 core.match版
(use '[clojure.core.match :refer [match]])

(defn quicksort1 [ls]
  (match ls
   ([] :seq) []
   ([x & xs] :seq)
   (concat (quicksort1 (filter #(< % x) xs))
           [x] (quicksort1 (filter #(>= % x) xs)))))
 と, defun版
(use '[defun :only [defun]])

(defun quicksort2
  ([([] :seq)] [])
  ([([x & xs] :seq)]
   (concat (quicksort2 (filter #(< % x) xs))
           [x] (quicksort2 (filter #(>= % x) xs)))))
 余分なシンボルが減り, 括弧だけになってすっきりした印象になりました. matchの一旦変数を定義して, それをそのままmatchマクロに渡すという処理がなくなり1行分減り, 変数名が一つ減っています. しかし, 括弧は, 条件一つにつき, 2ペア増えてしまいました. どうやら, この括弧は(大括弧, 小括弧ともに)外せないようです.
(defun sum
  ([0 m] m)
  ([n m] (recur (dec n) (+ n m))))
こんな感じで, recurが内臓させることも可能なようで, loop-recurの省略にも使えます.

 パターンやガードの記述は, core.matchと同じ(というかcore.matchそのもの)なので, core.matchを使う要領で定義できます.

2015/03/05

doallで遅くなる(Clojure)

 mapなどリスト処理の間に, doallを間に挟むと, その分遅くなっているような気がしたので, 調べてみると, 確かにdoallを挟むことで実行時間が伸びています.

 副作用を考えない場合, 意味的にはidentityと同じなので, 忘れがちですが, しっかり計算コストとメモリ確保のためのコストがかかっているようです.

 以下がその証拠. (Linux Mint 17 Cinnamon 64-bit, Intel Core i7-4790, メモリ : 31.4GB)
user> (defn dub [n] (* n 2))
#'user/dub
user> (time (reduce + (map inc (map dub (range 10000000)))))
"Elapsed time: 1661.92639 msecs"
100000000000000
user> (time (reduce + (map (comp inc dub) (range 10000000))))
"Elapsed time: 1388.861015 msecs"
100000000000000
user> (time (reduce + (map inc (doall (map dub (range 10000000))))))
"Elapsed time: 1892.659523 msecs"
100000000000000
user> (time (reduce + (map inc (map dub (range 1000000000)))))
user> (time (reduce + (map inc (map dub (range 100000000)))))
"Elapsed time: 15394.220163 msecs"
10000000000000000
user> (time (reduce + (map (comp inc dub) (range 100000000))))
"Elapsed time: 13912.628054 msecs"
10000000000000000
user> (time (reduce + (map inc (doall (map dub (range 100000000))))))
"Elapsed time: 54838.803745 msecs"
10000000000000000
 一回目のトライは, 普通にmapを2回リストに適用した場合で1.6秒程度.

 次にmap分配則を用いて2回分のmapを1回に縮めてみた場合. 若干, 早くなりました. Haskellだとこの最適化は, 効果的なのですが, Clojureではおそらくリストの実装上の違いによりあまり効果がないようです.

 三行目が, 問題のdoallを挟んだリスト処理です. 通常のリスト操作と比べて0.2秒ほど, 遅くなっています. 遅くはなっていますが, 1.1倍程度です.

 下三行の入力は, リストのサイズを大きくした(1ケタ上げた)場合. map分配則による変換の効果の割合はリストサイズが10分の1以下の時と変わらないようです.

 一方で, doallを挟んだ場合は, 3倍以上の時間がかかっています. 与えられたリストのサイズに対して実行時間が線形に増加していないのは, メモリ確保に時間がかかっていることが原因だと思われます.

 ちなみに, 上記のreplのプロセスが終了すると約5.5GiBのメモリ領域が開放されました. というわけで, メモリサイズがシビアなノートPCとかだと, (GCなどに)よりその影響が顕著に出てくるかもしれません.

2015/02/14

quilを使って巡回セールスマン問題の解の経路を表示する(Clojure) + 矢印の表示


 は, ProcessingのClojureラッパーです. グラフィック系のライブラリで, Clojureから, ProcessingのAPI(メソッド等)を簡単に触ることができます. もともと, Processing自体が非プログラマを対象としたヴィジュアルアート向きの, Javaベースのプログラミング言語なので, その機能をClojureで使えるようにしてくれます. Clojure/Swingでやるには少し煩雑なグラフィック処理を, Processing風の書き方で書くことができます. マウスやキーボードの情報の取得も可能なので, ゲームも作れそうです.  以下がAPI(関数)の一覧.


 少し前に, 巡回セールスマン問題を遺伝的アルゴリズムを使って解いたのですが, その時の解(となる経路)の可視化がいまいちだったので, quilを使ってやり直してみました. もともとは, Processing同様, quilも 次のquil.infoにあるようなヴィジュアルアートを作成することを目的としているようですが, 容易に絵が描かけるという点には変わりありません.

 ちなみに, quilのコードを記述する部分はほぼ100%(?)手続き型言語風の書き方になります. map等でループを書くことはできますが.

 他のClojureのプラグイン同様, プロジェクトを作成しているなら, project.cljの :dependenciesに[quil "2.2.5"]を追加(quilでは試していないため保証はできませんが, 一応これ) でOKのはずです.
 lein repl 単体ならば, "~/.lein/profiles.clj" (または, それに相当するファイル)の:pluginsに対応するリストへ  [quil "2.2.5"] を加えれば利用できるようになります.

 実際にプロットさせると次のようになりました. (ソースコードは一番下)

 上記の表示は,
(def sample-points
  {:a [20 70] :b [30 20] :c [50 50] :d [60 55] :e [40 30]
   :f [20 50] :g [80 60] :h [10 30] :i [65 10] :j [40 80]})
のようなXY上に点がばらついているときに, 次のような解が生成された時に, 表示されるルートです.
(def route [:a :f :j :d :g :i :c :e :b :h])
 単に線を引いて, 丸を描いて, 文字を打ち込んでいるだけです. なので非常に簡単でした. しかし, Swingではなく, quil独特(Processing特有?)の関数名がいくつかあり, 目的の関数を探すのに時間がかかりました.

 API一覧はあるのですが, 例えば, 丸を描く関数名がellipseだったり(SwingだとfillOvalとかですが), その丸を描く関数を使うときに円の内側の色を指定する関数名が, fillだったりします. 線の太さは, stroke-weightで設定します.

 一方で, Processingのラッパーだけあっって, ネット上にあるProcessingのサンプルをClojure/quilに読み替えれば, 様々なことができるので, その点は非常に良い点です. quilのAPI一覧を見てよくわからない場合は, ProcessingのAPIやサンプルを見れば良いわけですから.

 ウィンドウのみを表示させたければ, 次のように書けます.
user> (require '[quil.core :as q])
nil
user> (q/defsketch skt2)
#'user/skt2
という形で表示できます. このdefsketch関数のオプションで, 描画の設定等を行います.
(defn plot-route []
  (q/defsketch skt1
    :title "routes"
    :setup (fn [] (q/smooth)) ;; anti-aliased
    :draw (draw-route sample-points route)
    :size [400 400]))
のようにdefsketchを定義します. タイトル等はオプションです. :setupは, ウィンドウ生成時に, 実行される関数, :drawは, 定期的に実行される関数ということでしょう. smooth関数でアンチエイリアスが指定できます. 上記の例を見ればわかるかとは思いますが, 滑らかな表示に切り替わります. :sizeのパラメータは, ウィンドウのサイズ.

 こういう経路情報の表示で意外と面倒なのが, 矢印の表示です 点(x1, y1)と点(x2, y2)間に線を引いた時, 頂点の片方に矢印の三角形を付け加えて矢印にする作業です. そして, quilにも多分矢印を表示する関数は無かったと思います.

 仮想的なXY空間を用意し, そこに矢印の先となる三角形を想定した3つの点をおいて, 回転行列/アフィン変換等で適当な角度に回転させます. 次にその3つの点に必要なdrawPolygonメソッド等で描画するというステップが必要になります.

 色々ググって調べてみたのですが, Swing/Awtでは, やはり上記の手法が一般的なようでした. 回転の処理を直接書くのは面倒なのですが, Processingには, matrix stackというオブジェクトがあり, これを使うことで簡単に三角形のオブジェクトを回転させることができることができるそうです.

 これで回転の処理の問題が解決したので, あとは与えられた座標から矢印の先端を表す三角形を回転させる角度がほしいのですが, これはatan(アークタンジェント)で解決します. しかし, ゼロ除算の問題でatanをそのまま使うよりも, atan2という(x, y)の値を渡すとその角度(単位はラジアン)を返してくれる関数の方が有用です. 知らなかったのですが, atan2は, MayaのMELエクセルとか, C++にも標準でついてくるような一般的な関数のようです.
(defn draw-arrow-head [xy1 xy2 top-mergin bot-mergin]
      (let [x1 (first xy1) x2 (first xy2)
            y1 (second xy1) y2 (second xy2)
            size (+ 25 (- bot-mergin))]
        (set-white)
        (q/push-matrix)
        (q/translate x1 y1)
        (q/rotate (+ (q/radians 270) (q/atan2 (- y2 y1) (- x2 x1))))
        (q/triangle 0 (+ top-mergin 0) 5 size -5 size)
        (q/pop-matrix)))
上記のような実装で, 矢印の頭の部分が書けます. 矢印の本体の線は, line関数で書けばよいので, 特に問題はないかと思います.

 rotateは, 時計回りにオブジェクトを回転させるようです. GUIプログラミングではお馴染みですが, Y軸が反転しているので, その反転分の270度を足すことで, 全体として, 適切に傾いた矢印が描けます. 以下の図が参考.
 atan2は, x軸に対する(x, y)座標の角度θを返します. translateには, オブジェクトの原点に相当する点の位置を与えます. 矢印の三角形を原点から少し離しているのは, 原点のまわりに各点を表すノードを表す円(ellipse)を表示させるためです. こんな感じで, ノード間の矢印を表しました.

 quilを使っていて少し困ったのは, エラーメッセージがreplに出力されないことです, 一度defsketchでオブジェクトを生成されてしまうと, :setup/:drawに指定した関数内で実行時エラーが発生した場合, エラーメッセージが出力されず, エラーが発生しているのかすら, 分かりませんでした. これはおそらく, defsketchを実行した時点で別プロセス(スレッド)でその処理が実行されるため, replのあるスレッド上ではエラーが発生していないことに由来するようですが, デバッグが大変なので, 対策を考える必要がありそうです.

 全ソースコードは以下から.
uid0130 / plot-graph.clj - Gist

参考文献

2015/02/10

ClojureでリストのHylomorphismの変換マクロ

 前に書いた記事(Clojureでunfoldとhylomorphism) の続き.

 上の記事では, ClojureにおいてList版のHylomorphismを中間のデータ構造を生成しない純粋な再帰呼び出しで書きなおすことで高速化しました. しかし, 可読性が最悪で, はっきりいって読めない...
(defn correct-answer? [arrangement]
  (reduce andf
    (map not-conflict?
      (take-while first (iterate rest arrangement)))))

(defn correct-answer?-v2 [arrangement]
  (hylol andf true #(not (first %)) not-conflict? rest arrangement))
 上の関数が, リストを生成しながら計算するHylomorphism(に相当するもの), 下の関数が, リスト(中間のデータ構造)なしで同じ計算をするものです. (iterateで)リストを生成せずにプログラムが実行される分, 下の関数の方が高速であるという結果になったという話でした.

 hylolとはどいういった関数かというと, 次のような関数です.
(defn hylol
  [f e p g h b]
  (if (p b)
    e
    (f (g b) (hylol2 f e p g h (h b)))))
 ピカソの絵画ような抽象的な定義ですが, しかし抽象的であるがために, ある程度普遍性があります(様々なプログラムに出現する可能性が高いという意味です).

 できれば可読性を維持したまま(つまり, reduceとかmapとかが書かれた状態のまま) Hylomorphismに相当するプログラムを純粋な再帰に変換したいですね. というわけで, 関数でhylolを実装するのではなく, マクロで実装すれば良いのでは? という結論に至りました.
(hylol-expand
  (reduce + 0
       (take-while (partial not= 0)
                    (map dec (iterate dec 100)))))
とか
(hylol-expand
 (->> (iterate dec 100)
      (take-while (partial not= 0))
      (reduce + 0)))
みたいな書き方なら, 可読性に関しては問題ありません.

 というわけで上記のような書き方を維持したまま, Hylomorphismを再帰呼び出しに書き戻すマクロを作ってみました. 以下がそのマクロの実装. hylol-formで展開後の構造を定義し, hylol-expandで記法のパターンマッチを行い, それに対応するマクロフォームに展開します. 素朴な実装なので, 特に柔軟性などはありません.
(use '[clojure.core.match :only (match)])

(defn hylol-form
  "returns quoted hylol form"
  ([f-in-reduce init-in-reduce pred use-map-option f-in-map
   f-in-iterate init-in-iterate]
     `(letfn [(~'f [~'b]
                (if (~pred ~'b)
                  (~f-in-reduce
                   ~(if use-map-option (list f-in-map 'b) 'b)
                   (~'f (~f-in-iterate ~'b)))
                  ~init-in-reduce))]
        (~'f ~init-in-iterate))))

(defmacro hylol-expand
  "hylomorphism of lists without generating list"
  [hylol-exp]
  (match [hylol-exp]
    [(['reduce f-in-reduce init-in-reduce ;; no map & natural
       (['map f-in-map
         (['take-while pred
           (['iterate f-in-iterate init-in-iterate] :seq)]
            :seq)] :seq)] :seq)]
    (hylol-form f-in-reduce init-in-reduce pred true f-in-map
                f-in-iterate init-in-iterate)
    [(['reduce f-in-reduce init-in-reduce ;; no map & natural
       (['take-while pred
         (['iterate f-in-iterate init-in-iterate] :seq)]
            :seq)] :seq)]
    (hylol-form f-in-reduce init-in-reduce pred false nil
                f-in-iterate init-in-iterate)
    [(['->> (['iterate f-in-iterate init-in-iterate] :seq) ;; map & thread
       (['take-while pred] :seq)
       (['map f-in-map] :seq)
       (['reduce f-in-reduce init-in-reduce] :seq)] :seq)]
    (hylol-form f-in-reduce init-in-reduce pred true f-in-map
                f-in-iterate init-in-iterate)
    [(['->> (['iterate f-in-iterate init-in-iterate] :seq) ;; no map & thread
       (['take-while pred] :seq)
       (['reduce f-in-reduce init-in-reduce] :seq)] :seq)]
    (hylol-form f-in-reduce init-in-reduce pred false nil
                f-in-iterate init-in-iterate)
     [the-other-cases]
     (throw (Exception. (str "malformed hylol : " the-other-cases)))))
 上記のようなマクロ定義で, 次のような4通りの記述が可能になります.
(println (hylol-expand
          (reduce + 0
                  (map dec
                       (take-while (partial not= 0)
                                   (iterate dec 100))))))

(println (hylol-expand
          (reduce + 0
                  (take-while (partial not= 0) (iterate dec 100)))))

(println (hylol-expand
          (->> (iterate dec 100)
               (take-while (partial not= 0))
               (map inc)
               (reduce + 0))))

(println (hylol-expand
          (->> (iterate dec 100)
               (take-while (partial not= 0))
               (reduce + 0))))
 reduce, map,  take-while, iterateの位置は固定で, これ以外の書き方はできません. 引数となる「関数/(繰り返しの)初期値」のみが変更可能です. 上の二例は, 引数渡しで普通に書けるケース, 下の二例は, スレッディングマクロで書かれたものを再帰に書きなおすものです.

 簡単なベンチマークテストのため, 前の記事で作成した8queen問題でテストしてみました. 前回書いたのは, 次のような書き方(関数コール)でした.
(defn correct-answer?-v2 [arrangement]
  (hylol andf true #(not (first %)) not-conflict? rest arrangement))
 hylolの変換マクロの場合, 同じ意味のプログラムが,
(defn correct-answer?-mbt [arrangement]
  (hylol-expand ;; = identify
   (->> (iterate rest arrangement)
        (take-while first)
        (map not-conflict?)
        (reduce andf true))))
で書けます. 見た目は, 大分改善されました.

 以下が簡単なベンチマーク結果. correct-answers-with-hylol-macroがhylolの変換のマクロ実装, correct-answers-with-hylol-functionは, 前に書いた記事のcorrect-answers-v2, correct-answersは, ナイーブな(hylolの変換が入らない)実装でこれも前回書いた記事と同じ.
user> (time (count correct-answers-with-hylol-macro))
"Elapsed time: 590.56785 msecs"
92
user> (time (count correct-answers-with-hylol-function))
"Elapsed time: 594.688247 msecs"
92
user> (time (count correct-answers))
"Elapsed time: 739.69738 msecs"
92
 これで, 可読性を維持したまま, Hylomorphismのマクロ展開により高速化が可能になりました.

2015/02/07

Clojure, incanterを使ってクラスタリング, Single Linkage(単連結法) / Complete Linkage(完全連結法)

 Clojureのincanter(統計解析用のライブラリ)が使えるようになったので, とりあえず, クラスタリングをやってみました. Single Linkage (単連結法)および, Complete Linkage(完全連結法)は, 階層的クラスタリングの手法の一種で, 距離空間に散らばった多数の点をその位置からそれぞれグループ(クラスタ)ごとに分類する手法です.

 例えば, 下の図は, 5つのガウス分布に従ってプロットされた二次元上の点です (各ガウス分布ごとに色で分類されています).

 これをComplete Linkageによって, 5つのクラスタへ分類すると, 次のようになります.
 元のガウス分布にいくらか従いながら, 近い(近そうな)点どうしが, それぞれのクラスタとして, ひとまとめにされていることが観察できると思います.

 クラスタリングでは, 空間上の点の位置を元にして同じようなどうしを一つにまとめ, 元の分布を予測したり, グループごとに分類します.

データクラスタの可視化

 incanterでは, 二次元空間上のグラフのプロットや行列の処理が簡単に使えるので, このようなクラスタリングが, 簡単に実装できます.

(require '[incanter.core :as core]
         '[incanter.stats :as st]
         '[incanter.charts :as ch]
         '[clojure.math.numeric-tower :as tower])

;; ------------------------- ------------------------- -------------------------
;; show clusters
;; ------------------------- ------------------------- -------------------------

(defn plot-with-clusters [clusters] ;; clusters = [[[x1 y1] [x2 y2] ...] ...]
  (let [head (first clusters)
        init-plot (ch/scatter-plot (map first head) (map second head))]
    (reduce #(ch/add-points %1 (map first %2) (map second %2))
            init-plot (rest clusters))))
 名前の衝突を避けるためにrequireでincanterをロードします. incanterではscatter-plotで点を二次元座標上にプロットできます(一番目のクラスタをscatter-plotでプロット). add-pointsでchartオブジェクトに点を付加, reduceで各クラスタを一つにまとめます.

サンプルとなるデータ点のリスト

 これに, ランダムなパラメータを持つ複数のガウス分布に従うデータセットを作ります.
(defn generate-random-parameter []
  {:sd (+ 0.15 (rand 0.1)) :x (dec (rand 2)) :y (dec (rand 2))
   :size (+ 20 (rand-int 10))})

(defn generate-sample-set [parts]
  (let [generate-by #(st/sample-normal (:size %1) :mean (%2 %1) :sd (:sd %1))]
    (map (fn [p] (map list (generate-by p :x) (generate-by p :y)))
         (take parts (repeatedly generate-random-parameter)))))
 sample-normalは, ガウス分布に従うサンプルデータを生成します. 標準偏差等のパラメータの範囲は適当です. generate-sample-setでサンプルのデータ点(x, y)のリストを生成します.

 データ的には, [ [[0.1 0.2] [-0.3 0.1] ......] [[0.3 -0.02] ......] ......]のように, 各クラスタごとに(x, y)座標がひとつのリストに.

 (coreの)view関数で, (scatter)プロットされたオブジェクトを表示できます.
user> (core/view (plot-with-clusters (generate-sample-set 3)))
 このリストを(apply concat ls)で, 結合してサンプルのデータ点のリストを作ります. 何か面白いデータがあれば, それでクラスタリングができるのですが...... このデータ列が本来クラスタリングにより分析されるデータになります.

距離

 距離空間には, マンハッタン距離やマハラノビス距離などありますが, 今回は典型的なユークリッド距離を使います. 簡単にいえば, 視覚的な距離がそのまま, データ間の距離になります.
(defn distance [x1 y1 x2 y2]
  (double (tower/sqrt (+ (core/sq (- x1 x2)) (core/sq (- y1 y2))))))

(defn compare-pairs [xy1 xy2]
  (distance (first xy1) (second xy1) (first xy2) (second xy2)))

Single Linkage, Complete Linkageと階層的クラスタリング

 Single/Complete Linkageは, 階層的クラスタリングの一種です. 階層的クラスタリングとは, クラスタ群の中から2つクラスタを選びその2つを結合していき, 目標となるクラスタ数に達するまで結合を続けます. クラスタ群の初期状態では, データ点一つに対し, 一つクラスタを定義します. 階層的クラスタリングの欠点の一つは, 最終的なクラスタ数を指示する必要があることです. 最終的なクラスタ数の適切な値がわからない時にはこの設定が問題となります. しかし, 裏を返せば, 好きなようにクラスタ数を決められるということでもあります.

 Single LinkageとComplete Linkageの違いは, クラスタ間の距離の定義です. Single Linkageでは, 2つのクラスタ間どうしで最も近い要素間の距離がそのクラスタ間の距離になります. 一方, Complete Linkageは, 最も遠い要素間の距離がクラスタ間の距離とします.
(defn single-linkage "cluster a & cluster b -> distance" [a b]
  (reduce min (map #(reduce min (map (partial compare-pairs %) b)) a)))

(defn complete-linkage "cluster a & cluster b -> distance" [a b]
  (reduce max (map #(reduce max (map (partial compare-pairs %) b)) a)))
 プログラム上の違いはmin/maxだけですが, Single Linkageの考え方としては, 2つのクラスタが同じクラスタなら要素間に繋がりがあるはずだということで, Complete Linkageの方は, 同じクラスタに含まれるということは要素がそれなりに凝集しているはずだ, ということでしょう.

IncanterのMatricesと距離行列

 incanterを使っていて, せっかくなので距離行列(クラスタ間の距離のデータを保存している行列)もincanterの行列を使いました.

 Clojureのシーケンス(リスト)から, matrixオブジェクトを生成できます. 同じクラスタ群の距離を保存するための行列なので, 対象行列になり, 対格要素は0になります. ここでの計算において, 対格要素は使いませんが, わかりにくいので対格要素は無限大としました.
 こんな感じ.
user> (generate-init-matrix (map (fn [x] [x]) (apply concat (generate-sample-set 5))))
 A 118x118 matrix
 -----------------
 Infinity  1.79e-01  3.13e-01  .  1.96e+00  2.09e+00  2.26e+00 
 1.79e-01  Infinity  1.78e-01  .  1.85e+00  1.99e+00  2.14e+00 
 3.13e-01  1.78e-01  Infinity  .  1.93e+00  2.08e+00  2.21e+00 
 ... 
 1.96e+00  1.85e+00  1.93e+00  .  Infinity  2.07e-01  3.39e-01 
 2.09e+00  1.99e+00  2.08e+00  .  2.07e-01  Infinity  3.69e-01 
 2.26e+00  2.14e+00  2.21e+00  .  3.39e-01  3.69e-01  Infinity
 matrixオブジェクトは, デフォルトでtoStringメソッドが実装されているため, replでは, 間を省略した行列が表示されます. 上の例は, 対格要素をPOSITIVE_INFINITYにした118×118の行列. 巨大な行列でも間引かれているので表示されるときはコンパクトになります.

(def inf Double/POSITIVE_INFINITY)

(defn distances-xy-xys [xy xys]
  (cons inf (map (fn [xy1] (compare-pairs (first xy) (first xy1))) (rest xys))))

(defn generate-init-matrix [xys]
  (let [dist-vec (map distances-xy-xys (drop-last xys) (iterate rest xys))]
    (core/symmetric-matrix
     (apply concat (concat dist-vec [[inf]])) :lower false)))

(defn find-conj-pair [dist-matrix width]
  (loop [x (dec width) y (dec width) conj-pair [width width] min-value inf]
    (cond
     (< y 0) conj-pair
     (< x 0) (recur y (dec y) conj-pair min-value)
     :else (let [value (core/sel dist-matrix :cols x :rows y)]
             (if (< value min-value)
               (recur (dec x) y [x y] value) ;; update pair
               (recur (dec x) y conj-pair min-value)))))) ;; not update pair

(defn update-dists
  "updates distance matrix
   dist-matrix : distance matrix
   i, j : delete target
   new-column : newly generated distances from i,j"
  [dist-matrix i j new-column]
  (->> (core/sel dist-matrix :except-rows [i j] :except-cols [i j])
       (core/bind-rows new-column)
       (core/bind-columns (cons inf new-column))))
 上2つの関数は, 初期状態の距離行列を作るためのものです. find-conj-pairで距離行列から結合する2つのクラスタのペアをインデックスで返します. update-distsで距離行列を更新します. 新しく生成されたクラスタと他のクラスタ間の距離を行列の上の行(と左側の列)に付け加えて更新します.

階層的クラスタリング

 クラスタ群は, 簡単のため, i,jのインデックスで扱います.
(defn remove-clusters [i j clusters]
  (let [;; supposed to be c_0 c_1 ... c_n ... c_m ... c_max
        n (min i j) m (max i j)
        ;; [[c_0  ... c_n ... c_m-1] [c_m ... c_max]]
        first-m-last (split-at m clusters)
        ;; [[c_0 ... c_n-1] [c_n ... c_m-1]]
        first-n-m (split-at n (first first-m-last))]
    ;; (concat [c_0 ... c_n-1] [c_n+1 ... c_m-1] [c_m+1 ... c_max])
    (concat (first first-n-m) (rest (second first-n-m)) (rest (second first-m-last)))))
 remove-clustersでクラスタ群のi,j番目のクラスタを除去します.

 以下は, クラスタリングのループです. 一番近いクラスタを探し出し, 結合して更新という処理を指定されたクラスタ数になるまで続けます.
(defn clustering [calc-dists xys target-size]
  (loop [clusters (map (fn [x] [x]) xys)
         dist-matrix (generate-init-matrix clusters)
         cluster-size (count xys)]
    (if (<= cluster-size target-size)
      clusters
      (let [conj-pair (find-conj-pair dist-matrix cluster-size)
            i (first conj-pair) j (second conj-pair)
            new-cluster (concat (nth clusters i) (nth clusters j))
            other-clusters (remove-clusters i j clusters)]
       (recur (cons new-cluster other-clusters)
              (update-dists dist-matrix i j (calc-dists new-cluster other-clusters))
              (dec cluster-size))))))

(defn link-with [dist-f new-cluster other-clusters]
  (map (fn [cluster] (dist-f new-cluster cluster)) other-clusters))

(defn clustering-with-single-linkage [xys number-of-clusters]
  (clustering (partial link-with single-linkage) xys number-of-clusters))

(defn clustering-with-complete-linkage [xys number-of-clusters]
  (clustering (partial link-with complete-linkage) xys number-of-clusters))

実行結果

こんな感じでサンプルのデータ点を生成して,
user> (def A (generate-sample-set 5))
user> (core/view (plot-with-clusters A))

こんな感じのガウス分布になっているものに, 最終的なクラスタ数を5に設定してSingle Linkageをつかうと,
user> (core/view (plot-with-clusters (clustering-with-single-linkage (apply concat A) 5)))

となり, 同様に最終的なクラスタ数を5に設定したComplete Linkageでは,
user> (core/view (plot-with-clusters (clustering-with-complete-linkage (apply concat A) 5)))
になります.

 Complete Linkageの方が, うまく(?)クラスタリングされているように見えます. 元の分布を完全に予測しているわけではありませんが, 生成されたときの分布4に近いかたちで分類できています. Single Linkageは, 孤立した点が一つのクラスタに分類されてしまう傾向があるようです. 全体がほぼ2つのクラスタに分類されてしまっています.

参考文献
 全ソースコード.
(require '[incanter.core :as core]
         '[incanter.stats :as st]
         '[incanter.charts :as ch]
         '[clojure.math.numeric-tower :as tower])

;; ------------------------- ------------------------- -------------------------
;; show clusters
;; ------------------------- ------------------------- -------------------------

(defn plot-with-clusters [clusters] ;; clusters = [[[x1 y1] [x2 y2] ...] ...]
  (let [head (first clusters)
        init-plot (ch/scatter-plot (map first head) (map second head))]
    (reduce #(ch/add-points %1 (map first %2) (map second %2))
            init-plot (rest clusters))))

;; ------------------------- ------------------------- -------------------------
;; generate data set
;; ------------------------- ------------------------- -------------------------

(defn generate-random-parameter []
  {:sd (+ 0.15 (rand 0.1)) :x (dec (rand 2)) :y (dec (rand 2))
   :size (+ 20 (rand-int 10))})

(defn generate-sample-set [parts]
  (let [generate-by #(st/sample-normal (:size %1) :mean (%2 %1) :sd (:sd %1))]
    (map (fn [p] (map list (generate-by p :x) (generate-by p :y)))
         (take parts (repeatedly generate-random-parameter)))))

;; ------------------------- ------------------------- -------------------------
;; clustering
;; ------------------------- ------------------------- -------------------------

(defn distance [x1 y1 x2 y2]
  (double (tower/sqrt (+ (core/sq (- x1 x2)) (core/sq (- y1 y2))))))

(defn compare-pairs [xy1 xy2]
  (distance (first xy1) (second xy1) (first xy2) (second xy2)))

(defn single-linkage "cluster a & cluster b -> distance" [a b]
  (reduce min (map #(reduce min (map (partial compare-pairs %) b)) a)))

(defn complete-linkage "cluster a & cluster b -> distance" [a b]
  (reduce max (map #(reduce max (map (partial compare-pairs %) b)) a)))

(def inf Double/POSITIVE_INFINITY)

(defn distances-xy-xys [xy xys]
  (cons inf (map (fn [xy1] (compare-pairs (first xy) (first xy1))) (rest xys))))

(defn generate-init-matrix [xys]
  (let [dist-vec (map distances-xy-xys (drop-last xys) (iterate rest xys))]
    (core/symmetric-matrix
     (apply concat (concat dist-vec [[inf]])) :lower false)))

(defn find-conj-pair [dist-matrix width]
  (loop [x (dec width) y (dec width) conj-pair [width width] min-value inf]
    (cond
     (< y 0) conj-pair
     (< x 0) (recur y (dec y) conj-pair min-value)
     :else (let [value (core/sel dist-matrix :cols x :rows y)]
             (if (< value min-value)
               (recur (dec x) y [x y] value) ;; update pair
               (recur (dec x) y conj-pair min-value)))))) ;; not update pair

(defn update-dists
  "updates distance matrix
   dist-matrix : distance matrix
   i, j : delete target
   new-column : newly generated distances from i,j"
  [dist-matrix i j new-column]
  (->> (core/sel dist-matrix :except-rows [i j] :except-cols [i j])
       (core/bind-rows new-column)
       (core/bind-columns (cons inf new-column))))

(defn remove-clusters [i j clusters]
  (let [;; supposed to be c_0 c_1 ... c_n ... c_m ... c_max
        n (min i j) m (max i j)
        ;; [[c_0  ... c_n ... c_m-1] [c_m ... c_max]]
        first-m-last (split-at m clusters)
        ;; [[c_0 ... c_n-1] [c_n ... c_m-1]]
        first-n-m (split-at n (first first-m-last))]
    ;; (concat [c_0 ... c_n-1] [c_n+1 ... c_m-1] [c_m+1 ... c_max])
    (concat (first first-n-m) (rest (second first-n-m)) (rest (second first-m-last)))))

(defn clustering [calc-dists xys target-size]
  (loop [clusters (map (fn [x] [x]) xys)
         dist-matrix (generate-init-matrix clusters)
         cluster-size (count xys)]
    (if (<= cluster-size target-size)
      clusters
      (let [conj-pair (find-conj-pair dist-matrix cluster-size)
            i (first conj-pair) j (second conj-pair)
            new-cluster (concat (nth clusters i) (nth clusters j))
            other-clusters (remove-clusters i j clusters)]
        (recur (cons new-cluster other-clusters)
               (update-dists dist-matrix i j (calc-dists new-cluster other-clusters))
               (dec cluster-size))))))

(defn link-with [dist-f new-cluster other-clusters]
  (map (fn [cluster] (dist-f new-cluster cluster)) other-clusters))

(defn clustering-with-single-linkage [xys number-of-clusters]
  (clustering (partial link-with single-linkage) xys number-of-clusters))

(defn clustering-with-complete-linkage [xys number-of-clusters]
  (clustering (partial link-with complete-linkage) xys number-of-clusters))

2015/02/03

incanterでhistgramをプロットする

 Clojureのincanterについては, ずいぶん前から知っていたのですが, 使う機会があまりなかった為, 特に触っていませんでした. 最近, Clojureで生成したデータからグラフをプロットするという作業が必要になったこともあって, incanterを弄っているので, メモです.


 とりあえず, histgramを使うところまで.

 以下のコマンドは, REPLに
user> (use '(incanter core stats charts pdf))
を叩き込んだ前提で. (2015/02/05, pdf追加)
 名前的に他のライブラリ等の関数名と名前が重複する可能性がある(特にnumeric-towerなど, な)ので, prefixなしで直にロードするというのはあまりよくありませんが, 弄ってみるだけなら, これでも問題ない筈です.

 公式のチュートリアル(Github)にも載っていますが, ヒストグラムの表示は,
user> (view (histogram (sample-normal 10000)))
でいけます.

 ヒストグラムは, n個のデータサンプルをカテゴリ別に分類し, そのデータの分布を可視化するものなので, カテゴリ数を調整したくなります. カテゴリ数は, :nbinsオプションで指定できます. カテゴリ数を1000にする場合, 以下のように書きます.
user> (view (histogram (sample-normal 10000) :nbins 1000))
これを実行すると, 次のようになります.

 (sample-normal 10000)は, 意味的には,
user> (import 'java.util.Random)
user> (def r (Random .))
user> (map (fn [x] (.nextGaussian r)) (range 10000))
と一緒のようです. 平均0, 標準偏差1の正規分布に従い乱数を生成するものです.

 sample-normalが生成するのは, 単なる数値型(double)のリストです.
user> (sample-normal 3)
(1.2821074070154617 1.638783357428357 -0.2801943679292413)
 サンプルのリストをconcatでつなげれば, 別々のサンプルを合成できます. :sdは標準偏差(Standard Deviation), :meanは平均のオプション.
user> (view (histogram (concat (sample-normal 10000 :mean -5 :sd 4)
                               (sample-normal 10000 :mean 4)) :nbins 1000))

 ヒストグラムが, 正規分布に従っているかどうかを調べるために, 正規分布の曲線を重ねるのは, histgramで生成したオブジェクトに, add-line関数で重ねます. pdf-normalは, 正規分布の関数. histgramの:densityオプションは, 密度(Density)か頻度(Frequency)のオプション.
user> (view (doto (histogram (sample-normal 10000) :nbins 100 :density true)
                  (add-lines (range -4 4 0.01) (pdf-normal (range -4 4 0.01)))))
で, これが次のようにプロットされます.

2014/12/31

ClojureでunfoldとHylomorphism

unfold

 unfoldは, fold(reduce)の双対として定義される関数で, リストを生成する一種のジェネレータのような働きをするものです. 畳み込み関数に対応する, 要素を展開する関数とみなすことができます. Clojureだとiterateが一番近い関数でしょう.
user> (take 10 (iterate (partial + 2) 0))
(0 2 4 6 8 10 12 14 16 18)
 Haskellでは, foldが二種類あり, foldrとfoldlで, 右から畳み込みと左から畳み込みに対応しています. Clojureのreduceは, foldlに相当します. そして, foldr/foldlそれぞれにunfoldr/unfoldlが定義できますが, 以下の文章では, foldlを前提に話を進めます. (Haskellのfoldとunfoldrについてはここで書いています.)

 unfoldlをClojureで書くとしたら以下のような感じになると思います.
(defn unfoldl [pred g h start]
   (map g (take-while #(not (pred %)) (iterate h start))))
 計算を始めるオブジェクトと繰り返し適用するh, 生成したリストにmapをかけるためのgと, リストを終わらせる条件式predからなります. map gがtake-whileよりも後の処理になるのも重要な点です.

Hylomorphism

 Hylomorphismとは, (代数的データ型へ)一般化されたfoldとunfoldの合成のことで, リストにおけるfoldlとunfoldlについて,

hylol f e p g h = foldl f e . unfoldl p g h

と書ける関数です. (太字部分 2015/01/08追記) unfoldlで生成したリストをfoldlで畳み込むような計算になります. map/reduceならぬ, generate(リスト生成)/map/reduceの処理が一つにまとまりました. Clojureだと次のようになります.
(defn hylol1
  "hylol f e p g h = foldl f e . unfoldl p g h"
  [f init pred g h start]
  (reduce f init (unfoldl pred g h start)))
 これだけだと, 抽象的な関数を組み合わせたまた別の抽象的なだけの関数といった印象しか残りません. しかし, Hylomorphismは, これを別のリストを使わないコードへ書き直すことができるのです.
(defn hylol2
  [f e p g h b]
  (if (p b)
    e
    (f (g b) (hylol2 f e p g h (h b)))))
 これが, HylomorphismLの別の実装です.

 foldlとunfoldlを組み合わせると当然, unfoldlで生成されたリストが, foldlにおいて, 畳み込まれてしまうため, 結果的に, 最終的な計算結果(戻り値)に不要なリストを計算途中で生成する必要が出てきてしまいます. しかし, 上記の2番目の実装では, リストを生成する様子がありません.

 Hylomorphismの変換では, foldとunfoldの合成関数を, これと同等の意味を持つ「リストを生成しないプログラム」へ書き換えることでプログラムを高速化します. この辺はmap分配則と同じような感じですね.

 例えば, 1〜100の和を計算するようなプログラムについて, foldl/unfoldlを使って,
(defn sum-from1to [n]
  (reduce + 0 (unfoldl (partial = 0) identity dec n)))
 と書けるわけですが, これをHylomorphismで書きなおしてみると, hylol2の実装を使って,
user> (hylol2 + 0 (partial = 0) identity dec 10)
55
と書けます. これをhylol2の実装へ展開すると,
(defn hylol2-sum
  [b]
  (if ((partial = 0) b)
    0
    (+ (identity b) (hylol2-sum (dec b)))))
 のように書きなおすことができて, foldl/unfoldlを使った場合に比べて素直な書き方になっていることがわかるかと思います.

 ある意味では, Hylomorphismの変換とは, リストを使って抽象化した処理を, 素朴な再帰関数に書き戻す手法だとみなすことができます.

Hylomorphismの変換例

 さて, Hylomorphismの変換に関する一番の関心事は, 本当に高速化するのかという点と, 具体的にプログラムへ適用できるのかという点です. というわけで, 以下は, 8Queen問題のプログラムでHylomorphismを使ってみたケースです.

 前に書いた記事の8Queen問題の解法では, バックトラックにより計算していましたが, 今回は, 素朴な総当り法を使います.

 8Queen問題の解とは, 究極的には, 1〜8の数字の並び, だとみなすことができます. 例えば, x座標を固定した場合, 鳩の巣論法的に, 各x座標(= 1, 2, ... 8)には必ず一つQueenが入っている必要があるからです. というわけで, x座標の情報は不要になり, 順序を持ったy座標列が, 解となります.

 というわけで, このプログラムでは, y座標のpermutationを与え, そこから, xy座標の組を計算し, その座標列が正しいかどうかを判定し, 正しければ解を返します. 以下がコードです.
(use 'clojure.math.combinatorics)

;; unfoldl & hylomorphism
(defn unfoldl [pred g h start]
   (map g (take-while #(not (pred %)) (iterate h start))))

(defn hylol2
  [f e p g h b]
  (if (p b)
    e
    (f (g b) (hylol2 f e p g h (h b)))))

;; 8 queen
(defn conflict? [x y others]
  (and (not (empty? others))
       (let [x1 (ffirst others) y1 (second (first others))]
         (or (= x x1) (= y y1)
             (= (- x x1) (- y y1)) (= (- x x1) (- y1 y))
             (conflict? x y (rest others))))))

(defn not-conflict? [ls]
  (let [xy (first ls) others (rest ls)]
    (not (conflict? (first xy) (second xy) others))))

(defn andf
  ([a] a)
  ([a b] (and a b)))

(defn correct-answer? [arrangement]
  (reduce andf
    (map not-conflict?
      (take-while first (iterate rest arrangement)))))

(defn correct-answer?-v2 [arrangement]
  (hylol2 andf true #(not (first %)) not-conflict? rest arrangement))

(defn check-and-cons [checker yss]
  (let [arrangement (map list (range 8) yss)]
    (if (checker arrangement) arrangement [])))

(defn get-correct-answers [permutations checker]
  (filter first (map (partial check-and-cons checker) permutations)))

(def perm (doall (permutations (range 8))))

(def correct-answers
  (get-correct-answers perm correct-answer?))

(def correct-answers-v2
  (get-correct-answers perm correct-answer?-v2))
 実行時間を計測してみると, 確かに高速化されていることがわかります. correct-answersの方が, foldl/unfoldl(に相当するもの)を使った場合で, correct-answers-v2がhylol2によりプログラム変換を行ったケースです.
user> (time (count correct-answers))
"Elapsed time: 733.504269 msecs"
92
user> (time (count correct-answers-v2))
"Elapsed time: 592.678338 msecs"
92
 上記の結果を見た場合, それ以外の部分でボトルネックになっている部分もありそうですが, Hylomorphismの変換により, 一応, 高速化ができていることが観察できました.

 さて, 問題のプログラム変換は以下の部分で行っています. correct-answer?が, 素朴なfoldl/unfoldlで, correct-answer?-v2がリストを生成しないHylomorphismを使ったコードになります. correct-answer?の方にunfoldlというキーワードは見当たりませんが, map/take-while/iterateの組み合わせがそれに相当します.
(defn correct-answer? [arrangement]
  (reduce andf
    (map not-conflict?
      (take-while first (iterate rest arrangement)))))

(defn correct-answer?-v2 [arrangement]
  (hylol2 andf true #(not (first %)) not-conflict? rest arrangement))
 ここでは, [[1 2] [3 4] ...]のようなQueenの位置を表すリスト(arrangement)を受け取り, そこから, 部分リストを生成し(iterate/take-while), 各部分リストの「先頭の要素」と「他の要素」が衝突しないか(互いが取り合える位置にいないか)をチェックして(map not-conflict?), 矛盾しているケースがないかどうかをandf(andの関数版)で確認し(reduce), その答えの正解/不正解をboolで返すというものです.

 例えば, 座標のリスト[p1 p2 ... p8]があるときに互いに衝突していないかを判定するためには, p1と[p2 ... p8], p2と[p3 ... p8], p3と[p4 ... p8]と判定していくため, 上記のようなrestをiterateする実装になっています.

 reduce, map. take-while, iterateが一堂に会しているので, 図らずともHylomorphismの変換が簡単にできました.

 Hylomorphismの変換から得られる教訓(?)は, プログラム中に, ジェネレータによるリスト生成を行う部分があり, その先で, reduce(畳み込み関数)がそのリストを畳もうと待ち構えているとき, そのコードは, リストを使わないコードへ書き換えることができるかもしれない, ということです. そして, 書き換えたコードは, 余分な(本来の計算結果に不要な)データ構造の生成が不要になることで, 高速化できる筈だという結論です.

参考文献
  • 関数プログラミングの楽しみ, Jeremy Gibbons & Oege de Moor(編集), 山下 伸夫 (翻訳), オーム社, 2010/6/23
    • ifphの続編(?), この記事の内容は本書のhylomorphismのセクションを元にかかれています.
  • 構成的アルゴリズム論
    • hyloだけでなく, map分配則などプログラム運算の様々な規則が紹介されています. チュートリアル.
  • Anamorphic Adventure in Clojure - Matlux
    • Clojureで, Anamorphism(unfold相当)を定義する話. 

2014/12/26

Clojureでダイクストラ法


 ダイクストラ法, 与えられたグラフにおいて, ある頂点からそれ以外の各頂点への最短のパス(ルート)を求めるアルゴリズムの一つです.

  1. グラフの頂点を, 到達済み頂点のグループと, 未到達頂点のグループへ分割します. 到達済み頂点のグループの初期値は, 開始頂点1つからなるグループで, 未到達の頂点グループの初期値は, それ以外の頂点すべてを含んだグループになります.
  2. 「到達済み頂点のグループのうち, どれかの頂点」(=vn)との間に「最も短い辺を持つ未到達の頂点」(=vm)を到達済み頂点のグループへ追加する. また, この辺を新しいルートとする. (開始点から新しく追加されたvmまでのパスは, 「開始点~vnまでのパス + 新しく追加された辺」となります.)
  3. 未到達頂点のグループに要素がなければアルゴリズム終了. 未到達頂点のグループに要素があれば, 2へ.

 日本語で書くとこんな感じで, 未到達の頂点が1つづつ減っていって, すべての頂点へのルートが決定します.

 なんとなく階層的クラスタリングに似ている気もしますが. 以下のような感じでダイクストラ法の実行結果が得られます. a~jがグラフ中の各点で隣の数字が開始点からの最短距離で, 赤線が開始点からの各点へのルートを表します. 黒線は, グラフに含まれているが, 最短ルートには含まれない辺. 例えば, c→jのルートは, c → d → a → j となります.

 以下の図において, 各点は, (x,y)座標が与えられており, 各点間の距離(=辺の長さ)は二乗距離で定義されます. (つまり, 視覚的な距離と, 頂点間の距離(辺の長さ)がおおよそ対応しています.)
 ナイーブなダイクストラ法の実装が以下のコードです. ダイクストラ法本体の実装より, グラフと実行結果可視化のコードの方がなぜか, 分量が多いですが.
(use '[clojure.core.match :only (match)]
     '[clojure.math.numeric-tower :as math])

;;---------------------------------------------------------------
;; find path

(defn mean-square [x1 x2 y1 y2]
  (math/sqrt (+ (* (- x1 x2) (- x1 x2)) (* (- y1 y2) (- y1 y2)))))

(defn get-weights [routes points-with-coord]
  (defn get-distance
    [point-A point-B]
    (let [a (get points-with-coord point-A)
          b (get points-with-coord point-B)]
      (mean-square (first a) (first b) (second a) (second b))))
  (->> routes
       (map #(get-distance (first %) (second %)))
       (interleave routes)
       (apply hash-map)))

(defn get-neighbors [routes]
  (->> (concat (map reverse routes) routes)
       (map vec)
       sort
       (group-by first)
       (map (fn [x] (list (first x) (map second (second x)))))
       (reduce concat)
       (apply hash-map)))

(defn find-path [start routes vs weights]

  (defn length [u v]
    (weights (sort (list u v))))

  (defn get-dist [x dp]
    (get-in dp [x :dist]))

  (defn find-nearest [xs dp]
    (->> (map (fn [x] [(get-dist x dp) x]) xs)
         sort first second))

  (defn short-route? [u v dp]
    (< (+ (length u v) (get-dist u dp)) (get-dist v dp)))

  (defn set-shorter [u v dp]
    (assoc-in dp [v] {:dist (+ (get-dist u dp) (length u v)) :prev u}))

  (defn update-at [u dp v]
    (if (short-route? u v dp) (set-shorter u v dp) dp))

  (let [distance-previous-data
        (->> vs
             (map #(list % {:dist Double/POSITIVE_INFINITY :prev nil}))
             (apply concat) (apply hash-map)
             (#(assoc-in % [start] {:dist 0 :prev nil})))]
   (loop [q vs dp distance-previous-data]
      (if (empty? q)
        dp
        (let [u (find-nearest q dp)
              new (reduce (partial update-at u) dp (u (get-neighbors routes)))]
          (recur (remove #(= u %) q) new))))))

;;---------------------------------------------------------------
(def spwc ;;sample-points-with-coord
  {:a [20 80] :b [35 20] :c [80 40] :d [60 100] :e [70 10]
   :f [10 60] :g [120 80] :h [20 35] :i [110 20] :j [0 120]})

(def sample-routes
  (let [routes
        [[:j :f] [:j :a] [:f :a] [:f :h]
         [:a :d] [:d :g] [:g :c] [:c :d]
         [:g :i] [:c :i] [:c :e] [:e :i]
         [:a :h] [:h :b] [:b :e] [:b :c]]]
    (map sort routes)))

;; samples
;; (find-path :j sample-routes (keys spwc) (get-weights sample-routes spwc))
;; (find-path :c sample-routes (keys spwc) (get-weights sample-routes spwc))

;;---------------------------------------------------------------
;; graph plot

(import (javax.swing JFrame))
(import (java.awt Color))
(import (java.awt Graphics))
(import (java.awt Graphics2D))
(import (java.awt Font))

(def radius 8)

(defn x-extend [x]
  (+ (* (- x radius) 3) 50))

(defn y-extend [y]
  (+ (- (* (- y radius) 3)) 400))

(defn vertex->str [key]
  (apply str (rest (str key))))

;; (plot-graph spwc sample-routes)
(defn plot-graph [points route-pairs answer]
  (let [window-size 450
        font  (Font.  "Monospace" Font/BOLD 14)
        frame (JFrame. "Route")
        get-x #(+ radius (x-extend (first %1)))   ;; for draw-line
        get-y #(+ radius (y-extend (second %1)))] ;; for draw-line

    (doto frame
      (.setSize window-size window-size)
      (.setVisible true)
      (.setResizable false))
    (Thread/sleep 100) ;; wait for the generateion of Window

    (def graphics (.. frame (getGraphics)))

    (defn plot-point [name]
      (let [coord (name points)
            r (* radius 2)
            x-pos (x-extend (first  coord))
            y-pos (y-extend (second coord))]

        (doto (cast Graphics2D graphics)
          (.setColor (Color. 255 100 100))
          (.fillOval x-pos y-pos r r)
          (.setColor (Color. 100 190 190))
          (.drawOval x-pos y-pos r r)
          (.setColor Color/blue)
          (.setFont font)
          (.drawString
           (->> answer name :dist round (str (vertex->str name) "="))
           (+ x-pos 20) (+ y-pos 10)))))

    (defn draw-line [point-a point-b]
      (let [a (point-a points) b (point-b points)]
        (if (or (= point-b (:prev (point-a answer)))  ;; from b to a
                (= point-a (:prev (point-b answer)))) ;; from a to b
          (.setColor graphics Color/red)
          (.setColor graphics Color/black))
        (.drawLine graphics (get-x a) (get-y a) (get-x b) (get-y b))))

    ;; draw background
    (doto graphics
      (.setColor Color/white)
      (.fillRect 0 0 window-size window-size))

    ;; draw route
    (doall (map draw-line (map first route-pairs) (map second route-pairs)))

    ;; draw points
    (doall (map plot-point (keys points)))

    (println "done!")))

;;---------------------------------------------------------------
;; execution example
(def from-a (find-path :a sample-routes (keys spwc) (get-weights sample-routes spwc)))
(def from-c (find-path :c sample-routes (keys spwc) (get-weights sample-routes spwc)))
(def from-e (find-path :e sample-routes (keys spwc) (get-weights sample-routes spwc)))
 こんな感じになりました. find-path関数で, 最短ルートとその距離を求めます. find-pathの一番目の引数は, 開始点. get-weights関数は, 各辺の重み(=距離)を計算しています. plot-graphでSwingにより, グラフを表示します.

 以下のように実行します.
user> (load-file "find-path.clj")
#'user/from-e
user> (plot-graph spwc sample-routes from-e)
done!
nil
user> (plot-graph spwc sample-routes from-a)
done!
nil
というわけで実行結果. 開始点をaとすると,
となり, 開始点をeとすると,

こんな感じ.