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/07/11

関手, 自然変換とHaskell

$\require{AMScd}$ 最近, Mathjaxで遊んでいるので, それを使って, 関手と自然変換について調べたことのメモ.

関手とHaskellのFunctor


 関手$F$は圏$X$から$Y$へ対応を表し, 圏$X$の各関数$f : A \rightarrow B$を $F(f) : F(A) \rightarrow F(B)$ として移すものです. 関手は元の圏の関数について合成則を保存し, $f : A \rightarrow B, g : B \rightarrow C$,を関手$F$によって移すと, $F(f) \circ F(g) = F(f \circ g)$という性質を持ちます. また, 元の圏の各対象Xについて, 移した先の関手$F(id_X)$の射もその先の対象の$id$射となります.

 Haskellにおける関手$F$は, 名前も性質もmap関数に似ています.

class Functor a where
  mapf : (a -> b) -> f a -> f b

 mapfの実装として, リスト上のmap関数を考えると, 関手の合成則の保存が成立し, (map f) . (map g) = map (f . g)という規則が成立するはずですが, これは, プログラムの融合変換では一般にmap分配則と呼ばれるもので, 一般に成立します.

 つまり, ある型からなる圏〜リストへの関手とは, ただのリストのmap関数だということになります. 上記の定義では, 対象についての対応関係を考えていないように見えますが, 射$F(f)$が定義されれば, 自動的に対象同士の対応関係も定義されるということでしょう. というわけで, 任意の型Aのリスト/Map関数は, AからAのリストへの関手だということになります(太字部分, 2016/05/02追記修正).

 一般にmap関数は, リストだけのものでなく, 木構造へ応用できます. 例えば, 2分木について, 木構造の各要素に関数fを適用するmap関数は容易に想像できます.

data BTree a = Leaf a | Fork (BTree a) (BTree a) deriving Show

を考えることにして,

instance Functor BTree where
  fmap f (Leaf a) = Leaf (f a)
  fmap f (Fork a b) = Fork (fmap f a) (fmap f b)

 タプルのリストで, 二番目の要素のみへの関数適用を考えるなら,

data TList a b = Nil | Cons a b (TList a b)

instance Functor (TList a) where
  fmap f Nil = Nil
  fmap f (Cons a b rest) = Cons a (f b) (fmap f rest)

とすれば, Functorのインスタンスとして具体的な関数が実装でき, 関手が定義できる流れになります. 代数的データ型なら, (fmap f) $= F(f)$の対応について, 自由に考えることが可能です. 関手というと身構えてしまいますが, 関数型言語上では, 単なる高階関数として表されることがわかります.

自然変換とHaskell上の関数


 関手$F, G$における自然変換$\mu$は, 関手$F, G : X \rightarrow Y$について, $\mu \circ F(f) = G(f) \circ \mu$を成り立たせるもの. というわけで, 可換図式を考えると,
$$\begin{CD} F(A) @>{F(f)}>> F(B)\\ @V{\mu_A}VV @V{\mu_B}VV \\ G(A) @>{G(f)}>> G(B) \end{CD}$$
 となります(上図 : 2016/05/18修正). 前述の条件を満たしていることがわかります. 自然変換の概念は, 2つの関手$F, G : X \rightarrow Y$の移した先の対応$\mu : F \dot{\rightarrow} G$を定義するものですが, $\mu \circ F(f)$のルートと, $G(f) \circ \mu$のルートの結果(計算)を同一にするものと考えると, わかりやすいですね. ちなみに, 上記の可換図式の個々の$\mu_A , \mu_B$は, 別の射(関数)である場合もあり, 個々の射は, 自然変換のコンポーネントと呼ばれるそうです.

 問題は, この自然変換を表すプログラム上の関数とは何かということです. 調べたら, 以下に例がありました.

 Category theory/Natural transformation - Haskell Wiki

 前述のデータ型BTreeとリスト([])で書くとすれば, 自然変換 $\mu :$ (BTree a) $\rightarrow$ [a]の自然変換とは, 例えば,

flatten2List :: (BTree a) -> [a]
flatten2List (Leaf a) = [a]
flatten2List (Fork a b) = (flatten2List a) ++ (flatten2List b)

 になりそうです.
 可換図式の主張する合成規則$\mu \circ F(f) = G(f) \circ \mu$は,

((flatten2List . fmap f) = (map f . flatten2List)

 が成立するはず. あまりにシンプルですね. ただのキャストに相当する(?)何かになりました.

 さらに自然変換$\mu$には, 別の自然変換$\eta$との垂直合成$\eta \circ \mu$を考えることも可能です. 前述の可換図式に下にもうひとつ自然変換をくっつけた形で, 以下のようになります.
$$\begin{CD} F(A) @>{F(f)}>> F(B)\\ @V{\mu_A}VV @V{\mu_B}VV \\ G(A) @>{G(f)}>> G(B) \\ @V{\eta_A}VV @V{\eta_B}VV \\ H(A) @>{H(f)}>> H(B) \end{CD} $$

 関手fmapの実装として, Maybeについて考えると,

data Maybe x = Nothing | Just x deriving Show

instance Functor Maybe where
  fmap f Nothing = Nothing
  fmap f (Just x) = Just (f x)

と定義することができて, この時の自然変換$\eta : $[a]$ \rightarrow$ (Maybe a)を考えれば,

list2Maybe :: [a] -> (Maybe a)
list2Maybe [] = Nothing
list2Maybe (x:xs) = (Just x)

となると思います. この時の, $\eta \circ \mu$による自然変換の垂直合成は,

(list2Maybe . flatten2List)

となり, 可換図式を満たすなら,

((list2Maybe . flatten2List . fmap f)
= (list2Maybe . map f . flatten2List)
= (fmap f . list2Maybe . flatten2List)

となるはずです.

参考文献

2015/05/29

Lispの新しい方言をいくつか

 Lispといえば方言ですが, 新しい(2000年以降に登場した)方言がたくさんあって, 覚えきれないのでメモ.


 によれば, ErlangVM上で動くLisp, LFEという方言があるそうです. 最初に登場したのは, 2008年なのだそうで, 2007年に登場したClojureとは一年違いです.

 英語版のWikipediaのLispの記事を見ると,
There are several new dialects of Lisp: ArcHyNuClojureLiskellLFE (Lisp Flavored Erlang) and Racket.
とあり,  Common Lisp, Emacs Lisp, Scheme, Clojure以外にも様々な方言が存在しています.

 Arcは, 特徴はいまいちよくわかりませんが, Paul Grahamによる方言とその実装. Hyは, PythonVM上で動作するLispで前に記事を書きました. Nuは, Objective-Cで書かれていて, MacOSXのアプリケーション向けのスクリプト言語で, Semantics的にはLispよりもRubyに近いそう. Liskellは, 名前からなんとなく想像がつきますが, HaskellのSyntaxをLisp風にしたものです. しかし, 実装は見当たらず, ネット上には論文しか落ちていません. (追記 : Rudolph Miller さんにGitHubに実装があると教えていただきました.) Racketは, もともとSchemeの処理系の一種でしたが, RnRS(Schemeの仕様書)から独立してRacketという言語になったものです.

 他にも調べると色々出てきました. どうやら, VM+Lispという組み合わせ(Clojure, Hy, LFE)が多いようなので, Luaで調べてみると, Moonlispというのがでてきましたが, これは, Luaのコード(VMのバイトコードなどではなく, 言語としてのLuaのコード)へコンパイルされるもののようです. CLR(.NET)では, L#, JavaScriptでは多すぎて調べきれませんが, ここ(JavaScript Lisp Implementations)に, そのリストがあります.

 また, Egisonというパターンマッチに重点をおいた方言がありますが. このあたりは, パターンマッチに否定的なClojureとは対照的かもしれません.

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/18

継承は必ずしもis-a関係ではない(OOP) - Inheritance is not subtyping問題

Inheritance is not subtyping」は, Cookらによる1989年の論文で, そのままですが, オブジェクト指向における継承≠派生型(is-aの関係)ということのようです.

 派生型(Subtyping)については, is-a - Wikipediaにその説明が載っていますが, 確かにここでも, 継承が必ずis-a関係を構築するとは書かれてはいません.

 数理科学的バグ撲滅方法論のすすめ - 第4回 関数型言語とオブジェクト指向,およびOCamlの"O"について - ITpro に, その反例を示したソースコードがありますが, OCamlなのであまり実感わきません. オブジェクト指向言語は, JavaかPythonくらいしか使ったことがないので, Pythonでとりあえず書き直してみました.

 Pointクラスを継承したColoredPointクラス.
class Point:
    def __init__(self, x, y):
        self.x = x
        self.y = y
    def equals(self, another):
        return self.x == another.x and self.y == another.y

class ColoredPoint(Point):
    def __init__(self, x, y, c):
        self.x = x
        self.y = y
        self.c = c
    def equals(self, another):
        xeq = self.x == another.x
        yeq = self.y == another.y
        ceq = self.c == another.c
        return xeq and yeq and ceq
 この例では, 確かにPointクラスが継承されていますが, equalsがオーバーライドされています. 一方で, Pointクラス間ではequals関数で比較できるように実装されています.

 PointクラスとColoredPointクラスのインスタンスがis-a関係ならば, 当然, Pointクラスのインスタンス⇔ColoredPointクラスのインスタンス間で, インスタンスのequalityの比較が可能なはずです.
>>> ColoredPoint(1,2,3).equals(ColoredPoint(1,2,3))
True
>>> ColoredPoint(1,2,3).equals(Point(1,2))
Traceback (most recent call last):
  File "", line 1, in 
  File "notisa.py", line 16, in equals
    ceq = self.c == another.c
AttributeError: Point instance has no attribute 'c'
 で, 比較してみると, エラーが発生する. という現象が反例としてあげられるため, 継承≠is-aということのようです.

 面白いのは, 継承しなくても, is-aの関係が成立してしまうケースが存在することで, 次のFruitとOarngeクラスの例が考えられます.
class Fruit:
    def __init__(self, price):
        self.price = price
    def equals(self, another):
        return self.price == another.price

class Orange:
    def __init__(self, price):
        self.price = price
    def equals(self, another):
        return self.price == another.price
厳密な意味(?)でis-a関係が成立しているわけではないですが, equalsによる比較ができてしまい, 一種の部分型(is-a関係)のような振る舞いが可能になります. OCamlなどでは, この辺を型推論により自動的に型付けしてくれるようで, このような型を構造的部分型(structual subtyping)と呼ぶようです(前述のITProの記事より).

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)))))
で, これが次のようにプロットされます.

2015/01/31

アンダースタンディング・コンピューテーション, メモと感想

 いつの間にか, 大学の図書館に『アンダースタンディング・コンピューテーション 単純な機会から不可能なプログラムまで』(Tom Stuart(著), 笹田耕一(監訳), 笹井崇司(訳), オライリー・ジャパン)があったので, 読みながらのメモと感想.

 読んでいて面白かったのは, 2章, 7章, 8章, 9章でした.

※以下のメモは, 各章に関する内容をまとめたものではなく, 各章を読んでいて, 気づいたこと/思ったこと/思い出したことをメモした内容です.

1章. Ruby

 触ったことなかったけど, Rubyの無名関数は,

 ->arguments { body }

 と書くらしい. 中括弧で関数本体を囲むのはともかくとして, なんで -> arguments なんだろう. Haskellでは, \arguments->bodyとかくけど関係があるのか......?

 Rubyの*演算子はapply関数相当 ... 準クォートのような文字列の式展開.

2章. プログラムの意味

 あやふやになりがちな, Big Step Semantics ↔ Small Step Semantics, 操作的意味論 ↔ 表示的意味論についての説明. 操作的意味論は, インタプリタ的な解釈で, 表示的意味論はコンパイラ的な解釈という感じでしょうか.

Small Step Semantics

 = プログラムの解釈の過程(計算過程)まで定義した反復的(iterative)な意味の定義.
 SchemeのR6RSの意味の定義は, Small Step Semantics. Schemeでも, R5RSやR7RSの定義は, Denotational Semanticsで, R7RSになって, 元に戻された感じ(戻されたといっても内容は違いますが).

Big Step Semantics

 = プログラムの意味を一度に表すための. 再帰的な定義.
以下は, Big Step Semanticsを使った例. Standard MLの定義なのだそうです.

3章 & 4章, 5章. 最も単純なコンピュータ & 能力を高める, 究極の機械

 いわゆるオートマトン(DFA, NFA, 正規表現, PDA, 構文解析など)の章. ここ最近, 形式言語関連はうんざりする程付き合ってきたので, ほとんど読み飛ばしてしまいました. 「形式言語とオートマトン」の教科書の内容が噛み砕いて解説されている印象. 学部生の頃に出会いたかったかな. 

 当然のように, (非)決定性オートマトン→正規表現→プッシュダウンオートマトン→構文解析→チューリングマシンという流れでした. 5章は, 究極の機械というタイトルですが, 内容はただのチューリングマシンです.

6章. 無からのプログラミング

 ここで遂に型なしラムダ計算が登場します. Rubyの無名関数, ->arguments { body }を使って, ラムダ計算をエミュレートします. pp.188-191にかけてのFizz Buzzは圧巻です.

 ちょっと特徴的だったのが, リストの扱いでした. ラムダ計算では, λxyf. f x yに相当する(Lispでいうところの)consは, PAIRと, car/cdrに相当する関数がLEFT/RIGHTと名付けられていました.

 これだけだと, 単に名前を変えただけなのですが, リストとconsセルを区別するために, 独自のリスト構造を定義しています. 初めて見る定義だったのですが, UNSHIFTというオペレータでリストを構築し, FIRST/RESTで要素を取り出します.

 UNSHIFTは, ほぼconsと機能的には同等のもので, UNSHIFTは, Rubyのunshiftメソッドに由来しているようです. Rubyのプログラムをそのままラムダ計算でエミュレートする流れであったため, リストのconsingは, UNSHIFTで行うということなのでしょう.

 例えば, (1 2)のようなリストは素朴なconsセル(この場合PAIR)を使った場合, データ構造としては, こんな感じのはずですが,
 本書では, 同じ (1 2)のリストは, 次のような構造になります. 素朴なconsセルのケースと同様に, PAIRにより構成されますが, 少し構造が違います.

  それで, UNSHIFTというオペレータで, 次のように構築しています.

 リストの各要素の手前にあるtrue/falseは, リストが空かどうかを判定するために用いられていました. この辺りがちょっと独特です.

7章. 至るところにある万能性

 このセクションは, チューリング完全な(計算)システムのコレクションです. 前半には, ラムダ計算, SKIコンビネータ, ι(イオタ)コンビネータ, 部分再帰関数のチューリング完全性についての説明があります.

 ラムダ計算で, チューリングマシンのエミュレータを書いて, SKIコンビネータで, ラムダ計算をエミュレートして, SKIコンビネータをιコンビネータでエミュレートして, という作業を繰り返すと, この3つがチューリング完全だという話でした.

 ところが, 中盤に, 普通の「形式言語」の教科書ではあまり扱われないタグシステムというのが登場します. 私にとっては, チューリングマシンやラムダ計算に比べると地味な存在でした. しかし, タグシステムに修正を加えた(同じようにチューリング完全な)循環タグシステムが, 「コンウェイのライフゲーム」や, 「ウルフラムの2, 3のチューリングマシン」がチューリング完全であることの証明に使われているという話があり, 面白いと感じました.

 余談ですが, このセクションで挙げられたもの以外では, 「マルコフアルゴリズム」なんかもチューリング完全です.

8章. 不可能なプログラム

 いわゆる停止性問題.

9章. おもちゃの国のプログラミング

 割りとふざけたタイトルに思えますが, 9章は, 抽象解釈(Abstract Interpretation)についての非形式的な説明が載っています. 「おもちゃの国」というのは, 要するに「抽象的なプログラムの世界」という意味なのでしょう.

 抽象解釈とは, プログラムの抽象的/近似的な実行(または解釈)のことです. あるプログラム(関数)が正しい結果を返すのかについて調べるために, そのプログラムへの入力をすべて試すというのは無理なので, よりシンプルな領域へ抽象化して, 抽象化された世界での結果を調べるという手法です.

 本書では, 掛け算における整数の符号の例の説明がありますが, 例えば, ある整数の元(z ∈ Z)同士の積は, 表で表すと, いわゆる無限に続く掛け算の九九のような形になります. それを正の数 : Pos, 負の数 : Neg, ゼロ : 0 からなる領域S = {0, Pos, Neg}で抽象的な表すと次のような表に収まりました. これが抽象解釈における抽象化です.


 で, 例えば, 1×-2 = -2のような計算が抽象化されて, Pos×Neg = Negのような計算(解釈)の世界が得られます.

 この時, 

f(x, y) = 1÷(if x * y < 0 then x else 2)

 というプログラムがあった時, この関数がゼロ除算のエラーを出さないことを調べたいとします.

 x, yがInteger型なら, Integer型の変数のペアの範囲すべてについて網羅的に計算すれば(とりあえず)OKです. つまり, こんな簡単なプログラムに, Integer.MAX_VALUE * Integer.MAX_VALUE回も実行する必要があり, 現実的ではありません.

 しかし, 抽象解釈なら, 上の領域S = {0, Pos, Neg}について, 3 * 3の9通りについて調べればよく, (f(0, 0) = Posとか, f(Neg, Pos) = Negといったふうに計算します.) 上記のプログラムは, ゼロ除算が起きないことを, (抽象解釈の世界にいながらにして), 知ることができます. 同じように網羅的に調べているだけですが, 抽象化おかげで計算量が格段に減っています.

 もっとも, 現実の問題がここまでうまくいく筈もありませんが, 少なくともこのようなケースが存在することは, 抽象解釈のアプローチが, ある種の問題に対して, 有効であることを示唆していると思います.

2015/01/05

Hy(lang), Lisp風味のPythonを使ってみた

 hy(lang)は, Lispです.

 公式にあるように, pipから簡単にトライできます.


 以下の記事は, hyのversion 0.10.1についての内容です.

 シンタックスは, 見た目, ほぼClojureにそっくりですが, cond, let式などをはじめとして, 微妙に異なる部分があるのと(どちらかといえばSchemeやRacketに近い), Python風の構文がいくつか含まれているなど, yet-another Clojureだと思ってプログラムを書こうとすると, 結構嵌まりました.

 ClojureからNumpyとか使ってみたい気もしますが, ベースがPythonであることもあり, hylangは, Clojureとは別の道を歩もうとしているようです.

 何回も使いまわしていますが, 8Queen問題の解法のプログラムを, Clojureで,書くと以下のような感じになるところを,
;; slove the 8 queen problem
(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 put-a-queen [x y others]
  (cond
   (< 7 x) false
   (not (conflict? x y others)) [x y]
   :else (put-a-queen (inc x) y others)))

(defn solve [x1 y1 answer1]
  (loop [x x1 y y1 answer answer1]
    (let [rests (rest answer)]
      (cond
       (and (< 7 x) (= 0 y))
       nil
       (< 7 y)
       (cons answer (solve (inc (ffirst answer)) (dec y) rests))
       :else
       (if-let [xy (put-a-queen x y answer)]
         (recur 0 (inc y) (cons xy answer))
         (recur (inc (ffirst answer)) (dec y) rests))))))

;; make the list of solutions
(def solutions (solve 0 0 []))
 hyで書くと, 次のようになりました.
;; slove the 8 queen problem
(require hy.contrib.loop)
(require hy.contrib.anaphoric)

(defn conflict? [x y others]
  (and (not (empty? others))
       (let [[x1 (car (car others))] [y1 (car (cdr (car others)))]]
          (or (= x x1) (= y y1)
              (= (- x x1) (- y y1)) (= (- x x1) (- y1 y))
              (conflict? x y (cdr others))))))

(defn put-a-queen [x y others]
  (cond
   [(< 7 x) false]
   [(not (conflict? x y others)) [x y]]
   [:else (put-a-queen (inc x) y others)]))

(defn solve [x1 y1 answer1]
  (loop [[x x1] [y y1] [answer answer1]]
    (let [[rests (cdr answer)]]
      (cond
        [(and (< 7 x) (= 0 y))
         nil]
        [(< 7 y)
         (cons answer (solve (inc (car (car answer))) (dec y) rests))]
        [:else
         (ap-if (put-a-queen x y answer)
           (recur 0 (inc y) (cons it answer))
           (recur (inc (car (car answer))) (dec y) rests))]))))

;; make the list of solutions
(def solutions (solve 0 0 []))
 ぱっと見は, Clojureとそっくりに書けることがわかります. もちろん, この書き方がhyの標準的なスタイルかどうかは分かりませんが. def/defnの記法は, Clojureと同じで, 大括弧([])を多様する書き方が, 全体をClojureっぽく見せているのだと思われます.

 condとlet式の部分は, 要素を偶数個並べる書き方ではありません. condは, 条件式とそれに対応する節のペアを括弧で囲います. letも, 変数を表すシンボルと, 割り当てる値を計算する式をリストで組にしたものを並べます(この辺は, Clojure感覚で書いていると, エラーが続出して, 結構嵌まりました).
 true/falseも, #t/#fではなく, true/falseのキーワードで表します. (Schemeと混同してました. 2015/02/02)
(defn put-a-queen [x y others]
  (cond
   [(< 7 x) false]
   [(not (conflict? x y others)) [x y]]
   [:else (put-a-queen (inc x) y others)]))
 さらに, car/cdrが復活しています(ただし, first/restも使えます).  loop/recurスペシャルフォームが使えますが, デフォルトでは使えません. contribからロードする必要があります.
(require hy.contrib.loop)
 こんな感じで, 割と表記上の問題は解決出来たのですが, リスト操作が不可解で, 結構難しかった印象です.

 一つは, empty?の特殊さで, これはリストの長さを計るのですが, (= (len ls) 0)の結果を返します. リストにrestをかけるとitertools.isliceというオブジェクトが生成されますが, このitertools.isliceオブジェクト, countableではないようで, (empty? (rest [1 2]))などと書くとエラーが出ます. そもそもempty?は, リストの長さを数える必要がないのですが, この辺はまだ, 実装途中といった感じなのかもしれません.
=> (rest [1 2])
<itertools.islice object at 0x271d208>
=> (empty? (rest [1 2]))
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/home/uks/Dropbox/langs/lisps/hy/example/env/local/lib/python2.7/site-packages/hy/core/language.hy", line 124, in is_empty
    (= 0 (len coll)))
TypeError: object of type 'itertools.islice' has no len()
=> (cdr [1 2])
[2L]
 ちなみにcdrが返すのはitertools.isliceではないようです.

 さらに, (憶測ですが)nil = Noneだったり, car ≠ first, cdr ≠ restだったり, 各関数がどのような内部クラスやオブジェクトの操作に対応しているのか知らないと, 嵌まるポイントがいくつかありました.
=> (first (rest (cons 1 None)))
=> (car (cdr (cons 1 None)))
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/home/uks/Dropbox/langs/lisps/hy/example/env/local/lib/python2.7/site-packages/hy/models/list.py", line 43, in __getitem__
    ret = super(HyList, self).__getitem__(item)
IndexError: list index out of range
 ネガティブなコメントが多くなってしまいましたが, 個人的には, Syntaxの観点から見ると, Clojureの書きにくかった部分が書きやすくなってるなーという印象です. Clojureは括弧の数を減らそうとしている箇所(cond, letなど)が多数見受けられますが, 逆にそこが仇となっている感じでしたがhylangでは元に戻っていました. リスト内包記法なんかも, Clojureのforとは比べ物にならないほど読みやすいです. リーダマクロも使い放題ですし.

 同じ機能を持つ関数やスペシャルフォームが, 複数名前があることが多いです. 例えば, 前述したcar/firstやdefn/defun, do/prognなどがあります.

 アナフォリックマクロが大量に定義されていて, contribをロードすると, 使えるようになります. anaphoricなのでapが頭につくみたいですね. ap-ifで, itに条件式の結果をbindします.
=> (require hy.contrib.anaphoric)
=> (ap-if true it false)
True