カテゴリー別アーカイブ: 制約論理プログラミング

CHRをGPUで実行させる論文を読みました。

Parallel Execution of Constraint Handling Rules on a Graphical Processing Unit
という論文を読みました。
リンク
リンク切れている場合は論文名で検索してください。

Constraint Handling Rulesの開発者のThom Fr¨uhwirth さん達が書いた論文で、CHRをGPU上で実行させる方法がCUDAのコード込みで書いてあります。

Constraint Handling Rulesに関して自分がわかっているレベル(≒あまりわかっていない)で簡単に説明すると、CHRは基本的には記述された一連の制約を論理的に同等かつよりシンプルなものにしたり、何かに変換(自作の制約で書いた記述 → 組み込み述語のみ使用した記述に変換など)したりします。

例えばECLiPSeではPropiaとは別にCHRを使用しても自作の制約が定義出来ます。

まず以下の3種類のいずれかの構文でルールを記述しておきます。(これが自作制約となります)
CHRには以下の3つしかルールがありません。
・simplification
  Head1,Head2 <=> Guard | Body
・propagation
  Head1,Head2 => Guard | Body
・simpagation
  Head1¥Head2 <=> Guard | Body
(HeadもGuardもBodyも普通のPrologのCallable Termで書きます。
GuardとBodyの部分は複数の述語を書いて良い。Headの部分は最大2つの述語のみ(ECLiPSeの場合。3つ以上が不可なのは以下に述べる述語のピックアップの際の効率の問題らしいです)→すみません訂正します。CHRの昔のライブラリでは2つのみでしたが、新しいECHライブラリでは3つ以上可能でした)

動作の仕組みなのですが、上記で定義した制約を含む述語群を連言(=AND。カンマ区切りの普通のPrologの表現)で記述すると、CHRの機構はそれらの述語群の中から2つピックアップし、
・上記のHead1,Head2にパターンマッチ
・GuardがTrueとなる
を満たすものを見つけると、以下のルールに従ってピックアップした2つの述語を「書き替えて」しまいます(これをfire[発火]と呼ぶ)そしてその際にBodyの実行可能な述語は実行してしまいます。
また、Guardの部分は常時チェックされる番犬のような動作をしていて、例えばguardにnonvar(A)みたいな記述をしていた際、別の処理でAに何か値を設定した場合、その瞬間にGuardのチェックが通り対応するBodyが発火されます(ECLiPSeのsuspend機構を使って実装されている)
simplificationの場合はBodyに書き換え
propagationの場合はHead1,Head2そのままでBodyの述語を追加
simpagationはHead2のみBodyに書き換え(Head1は残す)

上記の変換を適用できるルールがなくなるまで繰り返す(2つ述語ピックアップ→Head1,Head2にマッチングするかチェック→GuardがTrueかチェック→ルールに従いHead1,Head2を書き換え)
どうも書き換え中のルール群をconstraint storeと呼んでいるようです。

ちなみに上記一連の処理中にはバックトラックは発生しません。
CHRの動作を「コミッテドチョイス」と言うらしいのですが、CHRの機構がどの述語をピックアップするか勝手に決めるのでプログラマはコントロールできない(選択を委託する)という意味かもしれません。
また、Headの部分は2つじゃなくて1つでも良いです。

CHRは基本的には上記の書き換えを行うのみなので、用途は特に制約の表現のみにとどまらず、サンプルを見ると
 最大公約数の計算
 ガウスの消去法
 高速フーリエ変換
 マージソート
など、色々な応用があります。

めちゃめちゃ端折ってますが上記がCHRの基本動作です。

記述した制約群をどのような組み合わせでピックアップ・変換したとしても必ず最終的に生成される述語群が同じになるようなとき、Confluence(合流)性を持つと言うそうです。

Confluence性はどのようにすれば担保されるかなのですが、以下の論文に書いてありますが自分はまだ理解できていない(論理学の用語がたくさん使用されてて難しいため)どうも冗長なルールを追加することにより合流性が担保されるようです。

Theory and practice of constraint handling rules
リンク切れしている場合は論文名で検索してください。

Confluence性を満たす場合CHRが並列処理に非常に向いているとのことで、理屈では確かにそうだなと自分も納得しました。

前置きが長くなりましたが、上記を踏まえたうえでGPUでCHRを使った上記の論文を読んだわけです(GPUは数千コア積んでいるものもあり単純な計算を同時に行うのに非常に有力)
論文の内容はCUDAコードが実際に記載されていて実装イメージがしやすいものでした。
しかし、ルールの述語自体をCUDAコンパイル用のコードで書いているため、あまり汎用性がないのではないという気がしました。
使う側は当然任意の述語を使いたいのですが、論文のやり方ではそのたびに新たなプログラムをCUDA用に書いてコンパイルしなくてはならないのでは?と感じました。
あまりGPUのことわかってないですが、if文を多用するプログラムは遅くなるなど、性能を100%出すためにいろいろ条件があるようです。
現状ではCHRをGPUの速度を100%出し切って実行するのは難しいのかな、と感じました。
ただ理論的にCHRのいろいろな性質が証明されているので、ハードウェア側でCHRに向いた実装が実現されたら、かなり望みがあるような気がします。

余談ですが「CHRには3つしかルールがない」と書きましたが、
・simpagation
  Head1¥Head2 <=> Guard | Body
の、
・simplification はこれの Head1 がないバージョン(Head2がBodyに書き換えられる)
・propagation はこれの Head2 がないバージョン(Head1そのままでBodyが追加される)
となるので実質はsimpagationだけで良いみたいな文章をどこかで読みました。

labelingに関して

制約論理プログラミングで必ず出てくるlabelingに関して説明しようと思います。
だだーっとマニュアル見ないでいい加減に書くので多分嘘も書いてますが大体で読んでください。なるべくプログラム自体も書かない記述にします(めんどくさいので)

概要的にはECLiPSe CLP のtutorialに書いてあることを抜粋して書くので英語読める人はそちらを読んでもらう方が絶対に良いです。とても分かりやすいです。
ECLiPSeのtutorial

制約論理プログラミングでは基本的なプログラムの構造は以下のような形になってます。
1.制約の記述
2.解のサーチ
3.解の表示

labelingというのはこれの2.解のサーチに相当します。

labeling述語には制約変数のリストを渡します。
例えば、(便宜的に整数のみ扱います。というか整数しか自分は扱い方が良くわかりません。)
変数Aが1から3
変数Bが1から3
変数Cが1から3
といった制約をつけた後、labeling([A,B,C])というような呼び出しをすると、A=1,B=1,C=1という解が設定されます。
Prologをやっている人はピンとくると思いますが、labelingの呼び出しでバックトラックが発生し、failすると次のA=1,B=1,C=2という別解が設定されます。バックトラックするごとに制約を満たす解が順番に設定されます。

ちょっと余談で、ECLiPSe CLPの資料を読んでるとことあるごとに口酸っぱく書いてあるのですが
「1.の制約の記述の部分は、バックトラックを発生させてはいけない」とのことです。制約の部分でバックトラックを発生させる作りにすると正確な解が得られないとのことです。
バックトラックの有無がPrologの一般的な述語と制約論理プログラミングの制約の大きな違いだと思います。
ですので、ECLiPSe CLPにもSWI-PrologのCLPFDにもたくさん制約用の述語が定義されてますが、そのどれもがバックトラックをしない述語になってるはずです。

これもまた余談ですがECLiPSeには自分が任意に定義したPrologの述語を制約に変換することが出来るPropiaという仕組みがあり、とても便利です。これを使うとバックトラックを生成させてしまう術語がバックトラックを生成しない制約に変わります。

話は戻って、labelingなのですが、既定では与えられた変数群の解を生成する際、
・リストの先頭の変数から順番に値を確定
・値はドメイン(取りうる範囲)の小さいほうから大きいほうへ
決定しています。

ただ、これが一番良いやり方かどうかはわかりません。
SWI-PrologのclpfdでもECLiPSe CLPでも「どの変数の値を優先的に決定するか」「値をどのように設定するか」の2つの観点でいろいろなオプションが用意されてます。
例えば以下の観点のオプションなどあります。
・一番ドメインの範囲(取りうる値の範囲)が小さい変数から決定する
・値を中央から外側に向かうように設定する
・値をランダムに設定する

SWI-PrologのCLPFDではここまでの話でlabelingに関しての説明は完了なのですが、ECLiPSe CLPはこのlabelingの部分を完全にプログラム制御できます。tutorialのTree Search Methodの項を見てください。

まずlabelingを使用せずに、1つの制約変数を指定して、
indomain(変数)という書き方ができます。
これは変数の取りうる値を小さいほうから順に設定してfailするたびバックトラックをさせるという、labelingの1変数バージョンのような述語です。そしてindomainの2番目の引数にオプションを渡すと、いろいろな順番で値を決めることが出来ます(大きい順、ランダム、真ん中から、etc…)

また、deleteという紛らわしい名前の述語があるのですが、これにfirst_failなどのオプションを指定して制約変数のリストを渡すと、指定したオプションに応じたやり方で1つ変数を取り出し、それ以外の変数をRestに返してくれます。
first_failの場合はリスト内の一番ドメインが少ない変数を選択し、残りをRestで返します。
(「次に値を決める変数の選択」を行ってくれる述語)

ECLiPSeではこのdeleteとindomainを駆使して、細やかに変数のラベリングをコントロールできます。ECLiPSeのtutorialには、これらを使用してN-Queenの問題をいろいろなラベリング方法で解き、かかった時間を計測するという興味深い記事が出てますので是非読んでみてください。

あとsearch/6という似たようなことする術語もあるのですがまだよくわかってません。
あとECLiPSeでは整数以外の値が扱えるのですが、そのときにindomainの実数バージョンのようなlocateという術語があるのですが、これもまだ詳細に使いこなせてません。

これ読んでる方もぜひ一緒に勉強しましょう(そして教えて下さい)

DOMINO PUZZLEソルバーの解説

ここ最近、勉強で ECLiPSe CLPのExampleページのドミノタイル問題ソルバーのソースを解析していました。大体理解したので解説したいと思います。

問題の解説:
0-0、0-1、0-2…..6-6 の全部で28枚のドミノタイルがある。
tiles

このタイルが横8×縦7の長方形の中にランダムに敷き詰められている。
縦、横、左右逆、上下逆のどの置き方でも良い。

どのように置いたかは不明で、並べた結果全体の数字が以下のようになっていた場合、タイルがどのように配置されていたかを解け。
3 1 2 6 6 1 2 2
3 4 1 5 3 0 3 6
5 6 6 1 2 4 5 0
5 6 4 1 3 3 0 0
6 1 0 6 3 2 4 0
4 1 5 2 4 3 5 5
4 1 0 2 4 5 2 0

プログラムの解説:

プログラムの概要を説明します。
①述語one_by_two_tiling/4でLeftとUpというマス同士の接続状態を表す配列を作成、制約を設定します。

Leftの配列は以下の色のついた部分の接続状態を表す2次元配列です。
1つのドミノタイルとしてつながっていれば1、つながっていなければ0になります。
lefts

Upは以下の色のついた部分の接続状態を表す2次元配列です。
Ups

制約のつけ方なのですが、あるマスに注目し上下左右の接続状態を考えたとき、必ず上下左右のいずれか1つの方向のみ接続されていることになるので、そのような制約を設定します(上下左右の変数の和が1)

②2つのdo文で、隣接した2マスに注目しながら、「2つの数字-対応する接続変数」の組み合わせをどんどんリストに集めます。このとき、タイルは数字が小-大の順番になるようにして保存します。
その部分にタイルが配置されていた場合は、該当する接続状態の変数が1になるというわけです。

scan
ヨコ方向の隣接マスに注目した後は、タテ方向で同様の操作を行います。どちらも同じ配列にまとめます。
このようにまとめると当然 (1,4)-変数 というような組み合わせが複数見つかりますが、いくつか見つかった中のどれか1つの組み合わせだけが実際に配置されているタイルの情報だということになります。

③まとめた配列をkeysortでソートします。タイルごとにまとまります。

④group_same_key_values/2で、ソート後の
[...(1,4)-接続状態変数1,(1,4)-接続状態変数2,(1,4)-接続状態変数3,(1,5)-接続状態変数4,(1,5)-接続状態変数5...]
のような状態のリストを
[…(1,4)-[接続状態変数1,接続状態変数2,接続状態変数3],(1,5)-[接続状態変数4,接続状態変数5]…]
のような形式にまとめます。


[接続状態変数1,接続状態変数2,接続状態変数3] の和が1
[接続状態変数4,接続状態変数5] の和が1
というような制約を設定します(foreach)

⑥term_variableで全変数を1つのリストにまとめた後labelingですべての接続状態変数の解を探します。


prety_print出力結果:
pretty_print

Project Euler Problem 59

Problem 59

全体の文字数 – 「何文字がA-Z a-zの範囲だったか」 の数値をコストとして計算し、一番コストが小さくなるものをbb_minで調べました。
あんまり効率よくないプログラムで、1900秒くらいかかりました…遅いな…
あとテキストの文字数が3の倍数というところは前もって調べてて、それをあてにしたプログラムになってます(文字数が3の倍数でないとdecryptでこける)

ソース

:-lib(ic).
:-lib(ic_sets).
:-lib(propia).
:-lib(branch_and_bound).
:-lib(csv).

go(Str,Cost,Sum) :-
csv:(csv_read("p059_cipher.txt", [InLst], Option)),
char_code('a',AscA),
char_code('z',AscZ),
[Key1,Key2,Key3]#::AscA..AscZ,
decrypt(InLst,OutLst,[Key1,Key2,Key3]),
OutLst#::0..127,
OutLst#::0..127,
letter_count(OutLst,Cnt),
length(OutLst,Len),
Cost #= Len - Cnt,
flatten([Key1,Key2,Key3,OutLst],AllVals),
labeling(AllVals),
iso_light:(atom_codes(Str,OutLst)),
mysum(OutLst,Sum).

decrypt([],[],_).
decrypt([N1,N2,N3|NRest],[D1,D2,D3|DRest],[K1,K2,K3]):-
xor(N1,K1,D1) infers most,
xor(N2,K2,D2) infers most,
xor(N3,K3,D3) infers most,
decrypt(NRest,DRest,[K1,K2,K3]).

letter_count([],0).
letter_count([First|Rest],Cnt):-
letter_count(Rest,Cnt1),
#::(First, [65..90, 97..122],Hit),
Cnt #= Cnt1+Hit.

mysum([],0).
mysum([First|Rest],Sum):-
mysum(Rest,Sum1),
Sum is Sum1 + First.

実行結果

?- bb_min(go(A, Cost, Sum), Cost, _223).
A = 復号化した文字列(伏せます)
Cost = 伏せます
Sum = 伏せます
Yes (1893.05s cpu)

Project Euler Problem 52

Problem52

ECLiPSeで解きました。
処理時間は自分のマシン(Intel(R) Core(TM) i5-3340M CPU @ 270GHz)で6.5秒くらい。
見どころ?はECLiPSeのlogical-loopで書いた順列permutationと、propiaで数値を入れ替えた制約を作っているところ
logical loop はループ用のサブ述語みたいなアホっぽい定義がなくなってコードがエレガントになるのでSWI-Prologでもぜひ正式採用してほしい。
桁のところは一般化してなく、1桁、2桁、3桁… と手入力で増やしながら確認しました
(bb_minでそれぞれの桁での最小の数値が取れるようにしてます)。

プログラム:

:-lib(ic).
:-lib(propia).
:-lib(branch_and_bound).

go(Digits,[X1num,X2num,X3num,X4num,X5num,X6num]) :-
length(X1Lst,Digits),
X1Lst#::[0..9],
lst2num(X1Lst,X1num),
X1num #> 10 ^ (Digits-1),
X2num #= X1num * 2,
X3num #= X1num * 3,
X4num #= X1num * 4,
X5num #= X1num * 5,
X6num #= X1num * 6,
permutation(X1Lst,X2Lst) infers most,
lst2num(X2Lst,X2num),
permutation(X1Lst,X3Lst) infers most,
lst2num(X3Lst,X3num),
permutation(X1Lst,X4Lst) infers most,
lst2num(X4Lst,X4num),
permutation(X1Lst,X5Lst) infers most,
lst2num(X5Lst,X5num),
permutation(X1Lst,X6Lst) infers most,
lst2num(X6Lst,X6num),
labeling([X1num,X2num,X3num,X4num,X5num,X6num]).

lst2num(L,N):-
(foreach(Num,L),fromto(0,In,Out,N),count(I,0,_)
do
Out #= In + 10^I * Num
).

permutation(A,B):-
(fromto(A,AIn,AOut,[]),fromto(B,[BFirst|BRest],BRest,[])
do
select(BFirst,AIn,AOut)
).

実行結果(解答伏せます):

?- bb_min(go(XXXXXXX, [N1, N2, N3, N4, N5, N6]), N1, Option).
N1 = XXXXXXXX
N2 = XXXXXXXX
N3 = XXXXXXXX
N4 = XXXXXXXX
N5 = XXXXXXXX
N6 = XXXXXXXX
Option = bb_options(continue, -1.0Inf, 1.0Inf, 1, 1, 0, 0, one, _466, _467, _468)
Yes (6.59s cpu)

[Prolog]最近の勉強

最近あんまりPrologのブログ書いてませんが、ずっと研究してます(具体的にブログに書けるような結果が出てない)現在はECLiPse(IDEのほうじゃなくて制約プログラミングのほう)の勉強してます。

書くネタないので今までのPrologの勉強の流れを説明すると以下のような感じで進んできました(スパンは数年単位)

広井さんのページでPrologの基礎を把握しいくつかのパズルのソルバーを作る(20年近く前に大学の授業で受けていたがあまり深いところまで理解していなかった)

SWI-Prologのマニュアル読み込むうちCLPFDを発見し衝撃を受け、勉強開始

ProjectEulerでPrologで問題を解くうち、正答者スレでECLiPseを用いてSWIの自分のプログラムよりはるかに速いスピード(数十倍)で結果を出している人を発見して衝撃を受ける。プログラムの内容自体はCLPFDを使用した自分のとそんなに変わらなかった。

ECLiPseの勉強を始める。だんだんSWI-PrologのCLPFDがSicstusやECLiPseの制約ソルバーに比べると「オモチャ」レベルだったということがわかってくる。また、SWI-PrologのCLPFD開発者のmarkus triskaさんは開発の継続をsicstusで行うみたいな記述を最新のSWIのマニュアルで見つけた。

また、ECLiPse は node.jsを使用してjavascriptで呼び出せるらしく、ECLiPseの勉強はWeb系の仕事でのちのち面白い威力を発揮するかもしれないと期待を抱く。

ECLiPseには、自分がSWIのCLPFDでどうしてもできなかった自作制約の作成がPropiaとCHRの2種類の方法で可能なことを知り、そちらの勉強を開始。SWIのCLPFDの自作制約の作成方法はマニュアルが意味不明だった(まだ完成してないとも書いてある)

Propiaのほうは大体の概要を把握。任意の述語をそのまま制約にできて、Propagationの厳密さをinferという述語の引数でコントロールするという仕組みらしい

CHRの勉強に移る。ルール自体は少ないが結構難しくてどのように作ればよいかなかなか理解できない。CHRはHaskellやJavaなどでも使用できるようだ。

CHRが並列計算に向いていると知り、ハードウェアで実行したりGPUを使用したりして解くといった論文をネットで見つけて読み始め、少しECLiPseの勉強から浮気する。

「CHRをGPUを用いて計算する」方法を記載した論文を見つけるが実際にやった時の速度比較が載ったものを見つけることができなかったので、そういう論文が出ないかを定期的にウォッチすることに決め、とりあえず論文リサーチはやめる。CHRをハードウェアでやった論文では数万倍とかのスピードで解いてるの見つけたが良く読むとCHRの機構だけではなくプログラムのほうもハードウェア化していてちょっとちがうだろという感じだった(論文で使用している最大公約数求めるみたいな単純なロジックしかハードウェア化できなくて現実的ではないし、そもそもそれってCHR自体の速さではない)

またECLiPseのCHRの勉強に戻る

△いまここです。

まあ日本語の情報がないです。
ずーーっと英語読んでます。ECLiPse、日本で使ってる人100人もいないんじゃなかろうか。日本のECLiPse界では望まずして第一人者になりそうです。もしそういう人たちいれば仲間に入れてほしい。誰か一緒に勉強しませんか?情報交換して一緒に研究の成果を教えあいましょう

CHRはルール自体は3種類しかなくて単純なんだけど以下が良くわかってなくてまだ自作の制約を作るまで至ってない

・labelingの振る舞い
・制約の述語の変数に「値を設定した時?」と「新しい制約がセットされたとき?」に再度bodyがfireされるらしいがそこらへんの振る舞い

多分CHRの制約の形が「自分がCLPFDを勉強してきた上で持っている制約のイメージ」とかなりずれているのが理解が滞ってる理由だと思います。

自分の数年の目標であるTantrixというパズルのソルバーをECLiPseの自作制約でいつか作れれば良いなあ。手続き型だと現在のスキルで多分作れるが制約プログラミングで作るという目標がありわざとそういう手法では作ってない

…………………………うーん。

この記事書いてるうちになんで自分がCHRを必死に勉強してるのか分かんなくなってきた…
Propiaで間に合ってるような気がしてきたので、そっちでやるかな…。CHRは「つかうとすごいスピードが出るケースが見つかる(並列処理とかで)」とか「CHRで書かれためちゃくちゃエレガントなコードを見つける」みたいな状況になってから改めて勉強しようかな。
CHRの良さをご存じの方いらっしゃればぜひ koyahata*koyahatataku.com(*を@に変更) にメール下さい。そうでなくてもProlog関連で勉強中の人はぜひ連絡取り合って一緒に勉強しましょう

ブログまとめるのも自分の行動を見直す良い機会になりますね

よし方向転換