2014/08/21
R > 数独を解くシミュレーション（4）
R > 数独を解くシミュレーション（4）
一昨日までは入力候補を list に入れていたが、データフレームを使う方法に変え、解き方も改良した。実行環境は 2014/08/13 の実機を参照。下が実際の様子。一昨日は数十秒かかっていた問題を数秒で解き、続いてもう一問試す。問題は一昨日と同様 ↓ から。スクリプトは末尾に。
■ Challenging Sudoku Puzzles 97-100 - Target Time: 18 mins.
http://puzzles.about.com/library/sudoku/blprsudokuh25.htm
動画 : https://kenpg.up.seesaa.net/image/20140821_sud_1.flv
一問目の結果は https://kenpg.up.seesaa.net/image/20140821_sud_1.gif
下は二問目。その結果は https://kenpg.up.seesaa.net/image/20140821_sud_4.gif
今回から、試行錯誤で入力する数字をグラフィックウィンドウへ描く・描かないを選べるようにした。動画には両方あり。下は描かない方。一回以上試行錯誤したセルには四角形が入り、完成したら数字に変わる。この方が描画負荷が少ないので速度とメモリ消費の点で有利。ただし試行錯誤が長いと画面が変わらず、止まっているように見える。
↓ 今回のデータフレーム。各セルを一行にして、行番号 i、列番号 j、左上からのブロック番号 b、確定した数字 fix、入力候補の数字 can をまとめる。NA になっている fix をなるべく効率的に入力していくのが基本の流れ。
まず can が一つだけの行を探す。↓ の場合、セル A の can は 6 のみ。そこでデータフレームの A の行で fix = 6 を入力し、全セルの can を更新する。そうするとセル B の can は 2 だけになるので fix = 2 を入力。この繰り返しを最初に行う。
can が一つだけのセルがなくなったら、行、列、ブロック別に情報をまとめたデータフレーム ↓ を作り、入力候補 can の数が最小の所に着目して、埋める組み合わせを列挙する。下の場合、j - 9 つまり最右列を対象にする。
↓ 最右列の X、Y、Z に 1、8、9 が入るので、その組み合わせを考える。
↓ 可能な組み合わせを列挙した様子。X = 9 と Z = 8 は除外でき (1, 8, 9)、(8, 1, 9)、(8, 9, 1) の三通りになる。それぞれにつき、とりあえず数字を入れて先ほどと同じことを繰り返す。順列の生成には e1071 パッケージの permutations 関数を使う。
このように試行錯誤し、データフレームの中で can が入らないセルが出たら「矛盾」ということで捨てる。可能な組み合わせを全て列挙できれば、いずれは解答にたどり着く。下が今回のコード。従前に比べれば改良できたが、本当に可能な組み合わせを全て列挙できているか分からないし、難しい問題だとやっぱり時間がかかる。というか、待ち切れずに中断した問題も結構ある。
■ Challenging Sudoku Puzzles 97-100 - Target Time: 18 mins.
http://puzzles.about.com/library/sudoku/blprsudokuh25.htm
動画 : https://kenpg.up.seesaa.net/image/20140821_sud_1.flv
一問目の結果は https://kenpg.up.seesaa.net/image/20140821_sud_1.gif
下は二問目。その結果は https://kenpg.up.seesaa.net/image/20140821_sud_4.gif
今回から、試行錯誤で入力する数字をグラフィックウィンドウへ描く・描かないを選べるようにした。動画には両方あり。下は描かない方。一回以上試行錯誤したセルには四角形が入り、完成したら数字に変わる。この方が描画負荷が少ないので速度とメモリ消費の点で有利。ただし試行錯誤が長いと画面が変わらず、止まっているように見える。
↓ 今回のデータフレーム。各セルを一行にして、行番号 i、列番号 j、左上からのブロック番号 b、確定した数字 fix、入力候補の数字 can をまとめる。NA になっている fix をなるべく効率的に入力していくのが基本の流れ。
まず can が一つだけの行を探す。↓ の場合、セル A の can は 6 のみ。そこでデータフレームの A の行で fix = 6 を入力し、全セルの can を更新する。そうするとセル B の can は 2 だけになるので fix = 2 を入力。この繰り返しを最初に行う。
can が一つだけのセルがなくなったら、行、列、ブロック別に情報をまとめたデータフレーム ↓ を作り、入力候補 can の数が最小の所に着目して、埋める組み合わせを列挙する。下の場合、j - 9 つまり最右列を対象にする。
↓ 最右列の X、Y、Z に 1、8、9 が入るので、その組み合わせを考える。
↓ 可能な組み合わせを列挙した様子。X = 9 と Z = 8 は除外でき (1, 8, 9)、(8, 1, 9)、(8, 9, 1) の三通りになる。それぞれにつき、とりあえず数字を入れて先ほどと同じことを繰り返す。順列の生成には e1071 パッケージの permutations 関数を使う。
このように試行錯誤し、データフレームの中で can が入らないセルが出たら「矛盾」ということで捨てる。可能な組み合わせを全て列挙できれば、いずれは解答にたどり着く。下が今回のコード。従前に比べれば改良できたが、本当に可能な組み合わせを全て列挙できているか分からないし、難しい問題だとやっぱり時間がかかる。というか、待ち切れずに中断した問題も結構ある。
rm(list=ls())
options(show.error.messages=T, showWarnCalls=T, stringsAsFactors=F)
library(e1071)
sud_ini =
' 2 8 56 1 6 9 7 4 2 136 3 7 2 4 1 3 84 5 3 '
f_main = function(){
assign('sud_ini_num', unlist(strsplit(sud_ini, '')), env=.GlobalEnv)
f_plt(sud_ini_num)
readline('ENTER to start')
dat = f_ini(sud_ini_num)
assign('dsp', 0, env=.GlobalEnv) # 1 : 解きながら数字を描画。メモリ消費大
assign('slp', 0, env=.GlobalEnv) # 探索ごとの待機秒
assign('sta', proc.time(), env=.GlobalEnv)
dat = f_slv(dat) }
f_plt = function(sud){
pxs = 400
lty = c(3,3,1)
if(length(dev.list()) == 0) windows(height=(pxs-45)/96, width=(pxs-7.5)/96)
par(plt=rep(0:1, 2))
plot(c(0.5, 9.5), c(0.5, 9.5), axes=F, type='n', xaxs='i', yaxs='i')
abline(h=seq(1.5, 8.5, by=1), lty=lty)
abline(v=seq(1.5, 8.5, by=1), lty=lty)
f_num(sud, 1) }
f_num = function(sud, col){
skp = which(sud_ini_num != ' ')
for(i in 1:81){
xys = f_xys(i)
n = sud[i]
if(is.na(n) | (col != 1 & is.element(i, skp))) next
points(xys, cex=4, pch=15, col='white')
text(xys, cex=1.5, col=col, fon=1, lab=n)
} }
f_xys = function(ord){
return(list(
x = (ord - 1) %% 9 + 1, y = 9 - floor((ord - 1) / 9)
)) }
f_ini = function(sud){
df1 = data.frame()
for(n in 1:length(sud)){
if(sud[n] == ' '){
fix = NA
can = f_can(sud, n)
if(length(can)==0) return()
can = list(can)
} else {
fix = sud[n]
can = NA
}
i = floor((n - 1) / 9) + 1
j = (n - 1) %% 9 + 1
b = floor((n - 1) / 27) * 3 + floor((n - 1) %% 9 / 3) + 1
df1 = rbind(df1, cbind(i, j, b, fix, can))
}
return(df1) }
f_slv = function(df1){
# まず、候補一つのセルを埋める
# 最初を除き、即座に全部埋めると矛盾する場合ある
while(1){
flg = F
for(n in 1:nrow(df1)){
if(!is.na(df1[n,]$fix)) next
if(length(unlist(df1[n,]$can)) > 1) next
f_poi(n)
df1[n,]$fix = df1[n,]$can
df1 = f_upd(df1)
if(is.null(df1)) return()
flg = T
}
if(!flg) break
}
# 空白セル最小の行・列・ブロック ijb を探索
ijb = f_ijb(df1)
tmp = subset(ijb, !is.na(can))
ord = order(apply(tmp, 1, function(x){ length(unlist(x$can)) }))[1]
tmp = tmp[ord,]
ijb = as.character(tmp$ijb)
val = as.numeric(as.character(tmp$val))
can = unlist(tmp$can)
# あり得る埋め方を列挙
df2 = subset(df1[df1[,ijb]==val,], is.na(fix))
mat = apply(e1071::permutations(nrow(df2)), 2, function(x){ can[x] })
for(i in 1:nrow(df2)){
del = which(!is.element(mat[,i], unlist(df2[i,]$can)))
if(length(del) > 0) mat = mat[-del,]
}
if(nrow(mat)==0) return()
# 可能な埋め方を一つずつ試行
for(r in 1:nrow(mat)){
df3 = df1
for(c in 1:ncol(mat)) df3[rownames(df2)[c],]$fix = mat[r,c]
df3 = f_upd(df3)
if(is.null(df3)) next
df3 = f_slv(df3) # 再帰
if(is.null(df3)) next
}
return() }
f_poi = function(n, col='pink'){
points(f_xys(n), cex=4, pch=15, col=col)
Sys.sleep(slp) }
f_can = function(sud, i){
row = floor((i - 1) / 9) * 9 + 1
row = row : (row + 8)
col = (i - 1) %% 9 + 1
col = col + (0 : 8) * 9
blk = floor((i - 1) / 27) * 27 + 1
blk = blk + floor((i - 1) %% 9 / 3) * 3
blk = blk + rep((0 : 2) * 9, each=3) + (0 : 2)
exs = sud[c(row, col, blk)]
res = 1:9
return(res[which(!is.element(res, exs))]) }
f_ijb = function(df1){
tmp = data.frame()
for(lab in paste(rep(c('i','j','b'), each=9), 1:9)){
ijb = substr(lab, 1, 1)
val = substr(lab, 3, 3)
whi = unname(which(df1[,ijb]==val))
if(length(whi) == 0) next
fix = unlist(df1[whi,]$fix)
fix = sort(as.numeric(fix[which(!is.na(fix))]))
if(length(fix)==0){
fix = NA
} else {
fix = list(fix)
}
can = unlist(df1[whi,]$can)
can = sort(as.numeric(unique(can[which(!is.na(can))])))
if(length(can)==0){
can = NA
} else {
can = list(can)
}
tmp = rbind(tmp, cbind(ijb, val, fix, can))
}
return(tmp) }
f_upd = function(df1){
sud = unname(unlist(ifelse(is.na(df1$fix), ' ', df1$fix)))
if(!any(is.na(df1$fix))){
f_num(sud, 2)
message('COMPLETED')
print(proc.time() - sta)
options(show.error.messages = F)
stop()
} else if(dsp > 0){
f_num(sud, 2)
}
return(f_ini(sud)) }
f_main()