• 数独を解くプログラム

    昔いくつか書いたPrologのパズルソルバーのコードがリンク切れになってたので再掲します(これから4回分続きます)。

    まずは数独のソルバーです。

    これ書いたときはまだ制約論理プログラミングのことを知らなかったのでPure Prologで書いてます(SWI-Prolog)。

    その後SWI-Prologのマニュアルのclpfdライブラリの説明の箇所の数独ソルバーのプログラムを発見して衝撃を受け、制約論理プログラミングの勉強を開始したのであった。

    プログラム
    解説と実行結果

    以下は僕が衝撃を受けたSWI-Prologのマニュアルに載っている数独ソルバー(僕のプログラムよりエレガントかつ速い)

    :- use_module(library(clpfd)).
    sudoku(Rows) :-
    	length(Rows, 9), maplist(same_length(Rows), Rows),
    	append(Rows, Vs), Vs ins 1..9,
    	maplist(all_distinct, Rows),
    	transpose(Rows, Columns),
    	maplist(all_distinct, Columns),
    	Rows = [As,Bs,Cs,Ds,Es,Fs,Gs,Hs,Is],
    	blocks(As, Bs, Cs),
    	blocks(Ds, Es, Fs),
    	blocks(Gs, Hs, Is).
    
    blocks([], [], []).
    blocks([N1,N2,N3|Ns1], [N4,N5,N6|Ns2], [N7,N8,N9|Ns3]) :-
    	all_distinct([N1,N2,N3,N4,N5,N6,N7,N8,N9]),
    	blocks(Ns1, Ns2, Ns3).
    problem(1, [[_,_,_,_,_,_,_,_,_],
    	[_,_,_,_,_,3,_,8,5],
    	[_,_,1,_,2,_,_,_,_],
    	[_,_,_,5,_,7,_,_,_],
    	[_,_,4,_,_,_,1,_,_],
    	[_,9,_,_,_,_,_,_,_],
    	[5,_,_,_,_,_,_,7,3],
    	[_,_,2,_,1,_,_,_,_],
    	[_,_,_,_,4,_,_,_,9]]).

    実行結果:

    ?- problem(1, Rows), sudoku(Rows), maplist(portray_clause, Rows).
    [9, 8, 7, 6, 5, 4, 3, 2, 1].
    [2, 4, 6, 1, 7, 3, 9, 8, 5].
    [3, 5, 1, 9, 2, 8, 7, 4, 6].
    [1, 2, 8, 5, 3, 7, 6, 9, 4].
    [6, 3, 4, 8, 9, 2, 1, 5, 7].
    [7, 9, 5, 4, 6, 1, 8, 3, 2].
    [5, 1, 9, 2, 8, 6, 4, 7, 3].
    [4, 7, 2, 3, 1, 9, 5, 6, 8].
    [8, 6, 3, 7, 4, 5, 2, 1, 9].
    Rows = [[9, 8, 7, 6, 5, 4, 3, 2|...], ... , [...|...]].
  • Project Euler Problem 62

    Problem 62

    ECLiPSe CLP で解きました。
    ほぼPure Prolog で、制約論理の述語は使用してません。
    自分のマシンで約9秒。

    プログラム:

    :-lib(util).
    :-dynamic(cubic_sorted/2).
    
    go:-
    	retractall(cubic_sorted(_,_)),
    	between(1,1000000000000,I),
    		CubicNum is I * I * I,
    		digit_sort(CubicNum,Sorted),
    		assert(cubic_sorted(Sorted,I)),
    		findall(Num,cubic_sorted(Sorted,Num),NumLst),
    		length(NumLst,NumLen),
    		(
    			NumLen>=5->
    				sort(NumLst,NumLstSorted),
    				[First|_]=NumLstSorted,
    				CubeFirst is First * First * First,
    				printf(output, "CubicNum %d ,Lst %w", [CubeFirst, NumLst]),nl,
    				true
    				;
    				fail
    		)
    	.
    
    digit_sort(Num,Sorted):-
    	num2lst(Num,Lst),
    	msort(Lst,Sorted).
    
    % convert number to list(reversed)
    num2lst(0,[]):-!.
    num2lst(Num,[Dig|DigRst]):-
    	Dig is Num mod 10,
    	RestNum is Num // 10,
    	num2lst(RestNum,DigRst).

    コマンドライン:

    ?- time(go).
    Yes (9.22s cpu)

    標準出力:

    <解答伏せます。>

    Success, times: 9.218750s thread, 9.219+0.000s process, 9.224s real

  • 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)
  • Project Euler Problem79

    久しぶりにProject Eulerやってみました。
    制約論理プログラミングで解きました。ほとんど頭使ってない…
    SWI-Prologです。

    Problem79

    :-use_module(library(clpfd)).
    
    solve(A,Len):-
    	length(A,Len),
    	A ins 0..9,
    	choice(A,[3,1,9],IdxLst1,Len),
    	choice(A,[6,8,0],IdxLst2,Len),
    	choice(A,[1,8,0],IdxLst3,Len),
    	choice(A,[6,9,0],IdxLst4,Len),
    	choice(A,[1,2,9],IdxLst5,Len),
    	choice(A,[6,2,0],IdxLst6,Len),
    	choice(A,[7,6,2],IdxLst7,Len),
    	choice(A,[6,8,9],IdxLst8,Len),
    	choice(A,[7,6,2],IdxLst9,Len),
    	choice(A,[3,1,8],IdxLst10,Len),
    	choice(A,[3,6,8],IdxLst11,Len),
    	choice(A,[7,1,0],IdxLst12,Len),
    	choice(A,[7,2,0],IdxLst13,Len),
    	choice(A,[7,1,0],IdxLst14,Len),
    	choice(A,[6,2,9],IdxLst15,Len),
    	choice(A,[1,6,8],IdxLst16,Len),
    	choice(A,[1,6,0],IdxLst17,Len),
    	choice(A,[6,8,9],IdxLst18,Len),
    	choice(A,[7,1,6],IdxLst19,Len),
    	choice(A,[7,3,1],IdxLst20,Len),
    	choice(A,[7,3,6],IdxLst21,Len),
    	choice(A,[7,2,9],IdxLst22,Len),
    	choice(A,[3,1,6],IdxLst23,Len),
    	choice(A,[7,2,9],IdxLst24,Len),
    	choice(A,[7,2,9],IdxLst25,Len),
    	choice(A,[7,1,0],IdxLst26,Len),
    	choice(A,[7,6,9],IdxLst27,Len),
    	choice(A,[2,9,0],IdxLst28,Len),
    	choice(A,[7,1,9],IdxLst29,Len),
    	choice(A,[6,8,0],IdxLst30,Len),
    	choice(A,[3,1,8],IdxLst31,Len),
    	choice(A,[3,8,9],IdxLst32,Len),
    	choice(A,[1,6,2],IdxLst33,Len),
    	choice(A,[2,8,9],IdxLst34,Len),
    	choice(A,[1,6,2],IdxLst35,Len),
    	choice(A,[7,1,8],IdxLst36,Len),
    	choice(A,[7,2,9],IdxLst37,Len),
    	choice(A,[3,1,9],IdxLst38,Len),
    	choice(A,[7,9,0],IdxLst39,Len),
    	choice(A,[6,8,0],IdxLst40,Len),
    	choice(A,[8,9,0],IdxLst41,Len),
    	choice(A,[3,6,2],IdxLst42,Len),
    	choice(A,[3,1,9],IdxLst43,Len),
    	choice(A,[7,6,0],IdxLst44,Len),
    	choice(A,[3,1,6],IdxLst45,Len),
    	choice(A,[7,2,9],IdxLst46,Len),
    	choice(A,[3,8,0],IdxLst47,Len),
    	choice(A,[3,1,9],IdxLst48,Len),
    	choice(A,[7,2,8],IdxLst49,Len),
    	choice(A,[7,1,6],IdxLst50,Len),
    	label(A).
    
    
    choice(A,[P1,P2,P3],[N1,N2,N3],Cnt):-
    	[N1,N2,N3] ins 1..Cnt,
    	element(N1,A,P1),
    	element(N2,A,P2),
    	element(N3,A,P3),
    	N1 #< N2,
    	N2 #< N3.

    回答伏せます

  • ドローン

    なんか流れでドローンの制御をするソフトに携わることになってドローンの事情など何も知らない状態から苦労しながらどうにかこうにか設計にチャレンジしています。

    今のドローンはGPSと画像処理などを駆使してムチャクチャ頭良い処理ができるようになっていて、例えばコントローラーのタブレットで走る車などにターゲット指定したら移動するターゲットを常時距離を保ったまま追いかけ続け、カメラ(ジンバル)もターゲットに常時固定しながら動画を撮れる。

    素早く動くバイクなどが常時カメラに納まり続けているの動画を見てると、例えばもしこれが人を追いかけている状況だとしたらと考えると恐ろしい技術で、ジンバルに拳銃を固定したとして、遠隔で引き金を引くなんてプログラムは簡単にできるのだから、簡単にターミネーターのような殺人機械ができそうです。

    POIとかいうシステムで仏像とかのターゲットを指定したら、その周りをドローンがくるくる撮影しながらながら回り続けて、リアルな3dのモデリングデータ作成も出来る。

    ドローンは姿勢や位置の制御のために各種のセンサが搭載されているのだが、どうもこれらが安価に手に入るようになったことが安くドローンを作ることができるようになった理由らしくて、これらのセンサが安価につくれるようになった理由はスマホが大量生産されるようになったかららしい。

    確かにドローンに搭載されているGPSセンサ、気圧センサ、加速度センサ、磁気センサなどはどれもスマホに搭載されている。

    シヴィライゼーションというゲームで技術の順番を記述したテクノロジーツリーという考え方があり、例えば「法律」という技術は「アルファベット(文字)」という技術を取得していないと開発できないように順番があるのですが、さしずめ「ドローン」という技術の前提技術として「スマートフォン」が来る感じであろう。ちなみに「アルファベット」の前提技術は「筆記」となっている。

    20年前には自分がドローン制御のプログラムを書くことになるとは思いもしなかった。
    時代の進み方の速さにめまいがしそうだし、自分も末端の人間としてそのような進歩にほんの少しだけ関わっている

    勉強していてもオープンソース系でそのような技術のソフトウェア部分が構成されているのが多い印象で、技術進歩のスピードが無茶苦茶はやいのにはオープンソースコミュニティが1役も2役も買っていると思う

  • [Prolog] DCGで双方向性のある述語を作成

    最近stackoverflowでPrologの質問に良く回答しているのですが、面白そうな問題があり、自分でもプログラムを作ってみました。

    質問

    内容は、[1,2,3,4,1,4,7,5,7,8,9] のようなリストが与えられた際に「続いている数字は[最初の数字,最後の数字]の形で表し、そうでないものは単体の数字で表したようなリストを生成する という問題です。
    [1,2,3,4,1,4,7,5,7,8,9]は[[1,4],1,4,7,5,[7,9]] に変換されます。

    作っているうちにどんどんアイデアがわいてきて、はじめは長いプログラムだったのがDCG使ってエレガントな方法で実装する方法を思いつき、最後には入出力逆転してもOKなような非常に汎用的な形で作成できました。

    ECLiPSeで作成(多分他の環境でも問題なく動きます):

    chunk_list([])-->[],!.
    chunk_list([Chunk|Rest])-->chunk(Chunk),!,chunk_list(Rest).
    
    % sequence
    chunk([First,Last])-->sequence([ First , Last ] ),{ First < Last }. 
    
    % isolated number
    chunk(Val)-->sequence([Val,Val]).               
    
    sequence([First,Last]) --> 
                     { ( number ( First ) , number( Last ) ) -> First < Last ; true } ,
                     [ First ] ,
                     { succ ( First , Second ) } , sequence( [ Second , Last ] ) .
    
    sequence([Val,Val])-->[Val].

    実行結果:

    [eclipse 4]: phrase(chunk_list([2,[2,6],[3,5],3,7,[1,3]]),Ret).
    
    Ret = [2, 2, 3, 4, 5, 6, 3, 4, 5, 3, 7, 1, 2, 3]
    Yes (0.00s cpu)
    
    [eclipse 5]: phrase(chunk_list(Ret),[3,4,5,2,1,6,7,88,9,4,5,6,7,2,1,2,3]).
    
    Ret = [[3, 5], 2, 1, [6, 7], 88, 9, [4, 7], 2, [1, 3]]
    Yes (0.00s cpu)
    
    [eclipse 6]: phrase(chunk_list([[2,4],A,[2,8],3]),[2,B,4,6,2,C,4,5,6,7,D,E]).
    
    A = 6
    B = 3
    C = 3
    D = 8
    E = 3
    Yes (0.00s cpu)

    最後のマッチング例は「やってみたら、できた」という実験で、自分の予想を上回る汎用性が実現できていてびっくりしました。

    入出力がどちらもリストの述語を作成するときは、工夫すると双方向性が実現できることがあると思います。

    Prologの真髄に触れたような気がしました。

  • [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関連で勉強中の人はぜひ連絡取り合って一緒に勉強しましょう

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

    よし方向転換

  • C#のWPD(MTP)操作

    C#使ってMTPデバイスに読み書きするプログラム作ってるんですがこれがまあめんどい。

    WPDというAPIを使うんですがWeb上に全然動くサンプルがないし、Microsoftのページもかなり不親切な構成になっていて
    簡単なことやろうとしても非常に苦労します。

    極めつけはGetDevicesという接続デバイスを列挙する関数にバグがあって1デバイスしか取得できず、それを解決するためには
    Dllを逆アセンブル→修正→再アセンブルが必要という始末…
    参照

    フォルダ移動・アップロード・ダウンロード・削除といった基本機能ができるようになったプロジェクトをアップしますので、同じことやってる人はこれ使って時間無駄にしないようにしてください…いやーしょうもないことで時間使ってしまった。
    サンプルプロジェクト
    この程度のソースが全然手に入んないんです…

    Visual Studio 2012プロジェクトです。

  • PrologでForループ

    SicstusやECLiPSeではイテレータとかいうのを使ってDo-Loopが記述できるのですが、もともとECLiPSeで実現されていた論理プログラミングのループ(Logical Loop)を紹介する以下の論文を2002年にSchimpfさんという人が書いていて、これでLogical Loopが広まったようです?

    Logical Loop

    まだ深く読み込んでませんが、SWIにもなんかの方法で実装できるような気もする… lambdaライブラリとかを使う?