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