はじめに
この記事はSpeedcubing Advent Calendar 2020の2日目の記事です。
の続きをやっていきます。
後編では、index化による探索の効率化とTwo-Phase-Algorithmを実装して、現実的な時間で任意の状態のルービックキューブを解けるようにしていきます。
環境
環境については、前編中編と同様、Python3.6以上です。
前編,後編に引き続き、以下で出てくるコードをColaboratoryノートブック化したもの
も用意したので、[Play Groundで開く] をクリックし、順に実行しながら進めれば、環境構築など不要で楽だと思います。
Two-Phase-Algorithm
現実的な時間でルービックキューブを準最適に解くためのアルゴリズムとして、Two-Phase-Algorithm が広く使われています。
ルービックキューブの解法を2つのフェーズに分解し、それぞれをIDA*探索で解くというやり方です。
Phase 1では、ルービックキューブを自由に動かして、 <U,D,R2,L2,F2,B2> の動きだけで完成できる状態に持っていきます。
Phase 2では、そこから <U,D,R2,L2,F2,B2> の動きだけを使って完成まで持っていきます。
Phase 1は最悪で12手、だいたい10手くらい、
Phase 2は最悪で18手、だいたい13〜14手くらいで完成できます (cf. http://kociemba.org/math/symmetric.htm )。
フェーズを分けずに総当りするとだいたい平均18手ほど掛かります。
探索深さ (手順長) が深くなればなればなるほど指数関数的に探索範囲も増えていくので、10手の探索と13手の探索に分ければ速くなります。
更に、フェーズを2つに分けたことにより各フェーズの状態数が減り、状態を数値としてindex化し、予め状態遷移表と、枝刈り用の表も事前計算しておくことが可能になり、高速な探索を実装できます。
(補足: フェーズを分けない場合、そのままだとEPの状態が12!通りあり大きすぎてindex化が難しいですが、EPの状態を分割する工夫により実はindex化できます)
Phase 1の完成状態
Phase1のゴールは、<U,D,R2,L2,F2,B2> の動きだけで完成できる状態に持っていくことです。
つまり、完成状態から上面と下面は90度回してよいが、側面は180度しか回さないでたどり着ける状態に絞り込むということです。
完成状態のルービックキューブを、<U,D,R2,L2,F2,B2>の動きだけで崩してみてください。
例えばこんな状態になります。
Phase 1の完成状態の特徴としては、以下の3つがあります。
- COが全て揃っている
- EOが全て揃っている
- U面とD面に挟まれた中層 (以下E列と呼ぶことにします) のエッジパーツ (エッジパーツ0,1,2,3) はE列内でしか移動していない
逆に、この3つが満たされていれば、Phase 1の完成状態です。 CPやE列以外のEPの状態は無関係です。
Phase 1で考えなければいけない要素
Phase 1で考えなければいけない要素は、
- COの状態 2187 (3の7乗) 通り
- EOの状態 2048 (2の11乗) 通り
- E列の4つエッジパーツがどこにいるか 495 (12C4) 通り
です。 Phase 1の完成状態はこれらが全て合っていることです。
COは、ねじれの総和が3の倍数になると決まっているので、7個決まれば残りの1つは決まるので3の7乗となります。
EOも同様で、11個決まれば残りの1個は決まるので、2の11乗となります。
Phase 2で考えなければいけない要素
COが揃っている、EOが揃っている、E列のパーツがE列内にあることはPhase 1で達成済みで、Phase 2の操作では崩れないことが担保されているため、Phase 2で考えなければいけない要素は、
- CPの状態 40320 (8!) 通り
- UD面の8つのエッジパーツのEPの状態 40320 (8!) 通り
- E列の4つのエッジパーツEPの状態 24 (4!) 通り
となります。
定数の定義
この辺の定数は後で使うので定義しておきます
NUM_CORNERS = 8
NUM_EDGES = 12
NUM_CO = 2187 # 3 ** 7
NUM_EO = 2048 # 2 ** 11
NUM_E_COMBINATIONS = 495 # 12C4
NUM_CP = 40320 # 8!
# NUM_EP = 479001600 # 12! # これは使わない
NUM_UD_EP = 40320 # 8!
NUM_E_EP = 24 # 4!
状態のindex化
前編中編では、Stateクラスのインスタンスを操作しながら探索を進めていましたが、インスタンスの操作のコストが高いので、後編ではここも効率化します。
そのために、キューブの状態にindexを割り当てて、数値で状態を扱うようにします。
考える各要素を表すベクトルとindexの相互変換を実装していきます。
長くなりすぎるので、本記事では変換の理屈の説明は省略します。
より詳細な説明は Jaap's Puzzle Page の Computer Puzzling をご参照ください。
Phase 1 の状態のindex化
COのindex化
COの状態は、2187通りあるので、0〜2186の数値に各状態を割り当てていきます。これは、COを3進数として解釈してindex化すればできます。
先述の通り、最後の要素は勝手に決まるので無視します。
indexとcoの配列を相互に変換できる関数を書いておきます
def co_to_index(co):
index = 0
for co_i in co[:-1]:
index *= 3
index += co_i
return index
def index_to_co(index):
co = [0] * 8
sum_co = 0
for i in range(6, -1, -1):
co[i] = index % 3
index //= 3
sum_co += co[i]
co[-1] = (3 - sum_co % 3) % 3
return co
チェックしてみましょう。
print(co_to_index([0, 0, 0, 0, 0, 0, 0, 0]))
print(co_to_index([2, 1, 0, 0, 1, 0, 1, 1]))
print(co_to_index([2, 2, 2, 2, 2, 2, 2, 1]))
# 出力
0
1711
2186
print(index_to_co(0))
print(index_to_co(1711))
print(index_to_co(2186))
# 出力
[0, 0, 0, 0, 0, 0, 0, 0]
[2, 1, 0, 0, 1, 0, 1, 1]
[2, 2, 2, 2, 2, 2, 2, 1]
前編中編で扱ってきたcoの配列とindexの相互変換ができました。
index 0 がCOが揃った状態に対応します。
EOのindex化
COと同様です。 冗長ですが、分けて書いておきます。
def eo_to_index(eo):
index = 0
for eo_i in eo[:-1]:
index *= 2
index += eo_i
return index
def index_to_eo(index):
eo = [0] * 12
sum_eo = 0
for i in range(10, -1, -1):
eo[i] = index % 2
index //= 2
sum_eo += eo[i]
eo[-1] = (2 - sum_eo % 2) % 2
return eo
E列のエッジパーツの位置の組み合わせのindex化
E列のエッジパーツの位置の組み合わせを配列で [0, 0, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0] のように表すことにします。
1が立っているエッジの位置に、E列のエッジパーツのいずれかがいるということです。
def calc_combination(n, r):
"""nCrの計算"""
ret = 1
for i in range(r):
ret *= n - i
for i in range(r):
ret //= r - i
return ret
def e_combination_to_index(comb):
index = 0
r = 4
for i in range(12 - 1, -1, -1):
if comb[i]:
index += calc_combination(i, r)
r -= 1
return index
def index_to_e_combination(index):
combination = [0] * 12
r = 4
for i in range(12 - 1, -1, -1):
if index >= calc_combination(i, r):
combination[i] = 1
index -= calc_combination(i, r)
r -= 1
return combination
チェック
print(e_combination_to_index([1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0]))
print(e_combination_to_index([1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0]))
print(e_combination_to_index([0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1]))
# 出力
0
48
494
print(index_to_e_combination(0))
print(index_to_e_combination(48))
print(index_to_e_combination(494))
# 出力
[1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0]
[1, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 0]
[0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1]
組み合わせとindexの相互変換ができました。
前編で、E列のエッジパーツの位置が0〜3番目になるようにパーツ番号を定義しておいたので、index=0, [1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0] がE列のエッジパーツが全てE列にある、Phase 1の完成状態になります。
Phase 2 の状態のindex化
CPのindex化
CPの配列をindex化します。
def cp_to_index(cp):
index = 0
for i, cp_i in enumerate(cp):
index *= 8 - i
for j in range(i + 1, 8):
if cp[i] > cp[j]:
index += 1
return index
def index_to_cp(index):
cp = [0] * 8
for i in range(6, -1, -1):
cp[i] = index % (8 - i)
index //= 8 - i
for j in range(i + 1, 8):
if cp[j] >= cp[i]:
cp[j] += 1
return cp
チェック
print(cp_to_index([0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]))
print(cp_to_index([1, 7, 6, 2, 4, 5, 0, 3]))
print(cp_to_index([7, 6, 5, 4, 3, 2, 1, 0]))
出力
0
10000
40319
print(index_to_cp(0))
print(index_to_cp(10000))
print(index_to_cp(40319))
[0, 1, 2, 3, 4, 5, 6, 7]
[1, 7, 6, 2, 4, 5, 0, 3]
[7, 6, 5, 4, 3, 2, 1, 0]
CPとindexの相互変換ができました。
UD面のEPのindex化
パーツ数が8なのでCPと全く同じなのですが、一応別に書いておきます。
def ud_ep_to_index(ep):
index = 0
for i, ep_i in enumerate(ep):
index *= 8 - i
for j in range(i + 1, 8):
if ep[i] > ep[j]:
index += 1
return index
def index_to_ud_ep(index):
ep = [0] * 8
for i in range(6, -1, -1):
ep[i] = index % (8 - i)
index //= 8 - i
for j in range(i + 1, 8):
if ep[j] >= ep[i]:
ep[j] += 1
return ep
E列のEPのindex化
同様です。
def e_ep_to_index(eep):
index = 0
for i, eep_i in enumerate(eep):
index *= 4 - i
for j in range(i + 1, 4):
if eep[i] > eep[j]:
index += 1
return index
def index_to_e_ep(index):
eep = [0] * 4
for i in range(4 - 2, -1, -1):
eep[i] = index % (4 - i)
index //= 4 - i
for j in range(i + 1, 4):
if eep[j] >= eep[i]:
eep[j] += 1
return eep
遷移表の構築
後編の探索では状態をindexで扱っていきます。
各indexからある操作をしたときどのindexに遷移するかを予め計算して、それを引くだけで操作後のindexがわかるようにしておきます。
予め計算した遷移表は保存しておくことができるので、初回の初期化以降はそれを使って効率的に計算できます。
Phase 1 の遷移表
COの遷移表
NUM_CO × NUM_MOVES のサイズの二次元配列を用意し、その状態からその操作を行ったときにどの状態に遷移するかのindexを格納しておきます。
前編の実装で実際に全通り動かしてみることで、これを事前計算します。
COの遷移だけわかればいいので、cp, eo, cpには適当にゼロ配列を入れて計算しています。
import time
print("Computing co_move_table")
start = time.time()
co_move_table = [[0] * len(move_names) for _ in range(NUM_CO)]
for i in range(NUM_CO):
state = State(
[0] * 8,
index_to_co(i),
[0] * 12,
[0] * 12
)
for i_move, move_name in enumerate(move_names):
new_state = state.apply_move(moves[move_name])
co_move_table[i][i_move] = co_to_index(new_state.co)
print(f"Finished! ({time.time() - start:.5f} sec.)") # Finished! (0.32975 sec.)
EOの遷移表
同様に作ります。
print("Computing eo_move_table")
start = time.time()
eo_move_table = [[0] * len(move_names) for _ in range(NUM_EO)]
for i in range(NUM_EO):
state = State(
[0] * 8,
[0] * 8,
[0] * 12,
index_to_eo(i)
)
for i_move, move_name in enumerate(move_names):
new_state = state.apply_move(moves[move_name])
eo_move_table[i][i_move] = eo_to_index(new_state.eo)
print(f"Finished! ({time.time() - start:.5f} sec.)") # Finished! (0.30919 sec.)
E列エッジの組合せの遷移表
print("Computing e_combination_table")
start = time.time()
e_combination_table = [[0] * len(move_names) for _ in range(NUM_E_COMBINATIONS)]
for i in range(NUM_E_COMBINATIONS):
state = State(
[0] * 8,
[0] * 8,
index_to_e_combination(i),
[0] * 12,
)
for i_move, move_name in enumerate(move_names):
new_state = state.apply_move(moves[move_name])
e_combination_table[i][i_move] = e_combination_to_index(new_state.ep)
print(f"Finished! ({time.time() - start:.5f} sec.)") # Finished! (0.12020 sec.)
Phase 2の遷移表
Phase 2は使える操作が減っているので、これを定義しておきます。
move_names_ph2 = ["U", "U2", "U'", "D", "D2", "D'", "L2", "R2", "F2", "B2"]
CPの遷移表
print("Computing cp_move_table")
cp_move_table = [[0] * len(move_names_ph2) for _ in range(NUM_CP)]
start = time.time()
for i in range(NUM_CP):
state = State(
index_to_cp(i),
[0] * 8,
[0] * 12,
[0] * 12
)
for i_move, move_name in enumerate(move_names_ph2):
new_state = state.apply_move(moves[move_name])
cp_move_table[i][i_move] = cp_to_index(new_state.cp)
print(f"Finished! ({time.time() - start:.5f} sec.)") # Finished! (5.10473 sec.)
UD面エッジのEPの遷移表
print("Computing ud_ep_move_table")
ud_ep_move_table = [[0] * len(move_names_ph2) for _ in range(NUM_UD_EP)]
start = time.time()
for i in range(NUM_UD_EP):
state = State(
[0] * 8,
[0] * 8,
[0] * 4 + index_to_ud_ep(i),
[0] * 12
)
for i_move, move_name in enumerate(move_names_ph2):
new_state = state.apply_move(moves[move_name])
ud_ep_move_table[i][i_move] = ud_ep_to_index(new_state.ep[4:])
print(f"Finished! ({time.time() - start:.5f} sec.)") # Finished! (5.15420 sec.)
E列エッジのEPの遷移表
print("Computing e_edge_permutation_move_table")
e_ep_move_table = [[0] * len(move_names_ph2) for _ in range(NUM_E_EP)]
start = time.time()
for i in range(NUM_E_EP):
state = State(
[0] * 8,
[0] * 8,
index_to_e_ep(i) + [0] * 8,
[0] * 12,
)
for i_move, move_name in enumerate(move_names_ph2):
new_state = state.apply_move(moves[move_name])
e_ep_move_table[i][i_move] = e_ep_to_index(new_state.ep[:4])
print(f"Finished! ({time.time() - start:.5f} sec.)") # Finished! (0.00246 sec.)
枝刈り表の構築
Phase 1, Phase 2ともにまだ状態数が多すぎて、全ての状態について最短手順長を計算しておくことはできませんが、
考慮するパーツを限定して一部のパーツだけを揃えることを考えると、考慮している一部のパーツだけを揃えるための最短手順長を予め計算しておくことができます。
これを予め計算しておき、限定された一部のパーツですら残りの深さで揃わないので、全体も揃わないという知識を組み込めるようにして、枝刈りに使います。
Phase 1 枝刈り表
EOを無視して、COとE列だけ考えたときの最短手数表
NUM_CO × NUM_E_COMBINATIONS の二次元配列を用意して、各要素に、その状態からCOとE列をPhase1の完成状態に持っていくため最短手数を入れておきます。
print("Computing co_eec_prune_table")
start = time.time()
co_eec_prune_table = [[-1] * NUM_E_COMBINATIONS for _ in range(NUM_CO)]
co_eec_prune_table[0][0] = 0
distance = 0
num_filled = 1
while num_filled != NUM_CO * NUM_E_COMBINATIONS:
print(f"distance = {distance}")
print(f"num_filled = {num_filled}")
for i_co in range(NUM_CO):
for i_eec in range(NUM_E_COMBINATIONS):
if co_eec_prune_table[i_co][i_eec] == distance:
for i_move in range(len(move_names)):
next_co = co_move_table[i_co][i_move]
next_eec = e_combination_table[i_eec][i_move]
if co_eec_prune_table[next_co][next_eec] == -1:
co_eec_prune_table[next_co][next_eec] = distance + 1
num_filled += 1
distance += 1
print(f"Finished! ({time.time() - start:.5f} sec.)") # Finished! (8.46280 sec.)
COを無視して、EOとE列だけ考えたときの最短手数表
print("Computing eo_eec_prune_table")
start = time.time()
eo_eec_prune_table = [[-1] * NUM_E_COMBINATIONS for _ in range(NUM_EO)]
eo_eec_prune_table[0][0] = 0
distance = 0
num_filled = 1
while num_filled != NUM_EO * NUM_E_COMBINATIONS:
print(f"distance = {distance}")
print(f"num_filled = {num_filled}")
for i_eo in range(NUM_EO):
for i_eec in range(NUM_E_COMBINATIONS):
if eo_eec_prune_table[i_eo][i_eec] == distance:
for i_move in range(len(move_names)):
next_eo = eo_move_table[i_eo][i_move]
next_eec = e_combination_table[i_eec][i_move]
if eo_eec_prune_table[next_eo][next_eec] == -1:
eo_eec_prune_table[next_eo][next_eec] = distance + 1
num_filled += 1
distance += 1
print(f"Finished! ({time.time() - start:.5f} sec.)") # Finished! (8.01519 sec.)
Phase 2 枝刈り表
UD面のエッジを無視して、CPとE列エッジだけ揃えるときの最短手数表
print("Computing cp_eep_prune_table")
start = time.time()
cp_eep_prune_table = [[-1] * NUM_E_EP for _ in range(NUM_CP)]
cp_eep_prune_table[0][0] = 0
distance = 0
num_filled = 1
while num_filled != NUM_CP * NUM_E_EP:
print(f"distance = {distance}")
print(f"num_filled = {num_filled}")
for i_cp in range(NUM_CP):
for i_eep in range(NUM_E_EP):
if cp_eep_prune_table[i_cp][i_eep] == distance:
for i_move in range(len(move_names_ph2)):
next_cp = cp_move_table[i_cp][i_move]
next_eep = e_ep_move_table[i_eep][i_move]
if cp_eep_prune_table[next_cp][next_eep] == -1:
cp_eep_prune_table[next_cp][next_eep] = distance + 1
num_filled += 1
distance += 1
print(f"Finished! ({time.time() - start:.5f} sec.)") # Finished! (6.75380 sec.)
CPを無視して、UD面のエッジとE列エッジだけ揃えるときの最短手数表
print("Computing udep_eep_prune_table")
start = time.time()
udep_eep_prune_table = [[-1] * NUM_E_EP for _ in range(NUM_UD_EP)]
udep_eep_prune_table[0][0] = 0
distance = 0
num_filled = 1
while num_filled != NUM_UD_EP * NUM_E_EP:
print(f"distance = {distance}")
print(f"num_filled = {num_filled}")
for i_udep in range(NUM_UD_EP):
for i_eep in range(NUM_E_EP):
if udep_eep_prune_table[i_udep][i_eep] == distance:
for i_move in range(len(move_names_ph2)):
next_udep = ud_ep_move_table[i_udep][i_move]
next_eep = e_ep_move_table[i_eep][i_move]
if udep_eep_prune_table[next_udep][next_eep] == -1:
udep_eep_prune_table[next_udep][next_eep] = distance + 1
num_filled += 1
distance += 1
print(f"Finished! ({time.time() - start:.5f} sec.)") # Finished! (6.09583 sec.)
Phase 1の探索
まず Phase 1 の探索だけ書いてみましょう。
中編のコードを、状態遷移表と枝刈り表を使いつつ、Phase 1用に書き直してみます。
操作は相変わらず、操作名で扱おうと思うので、操作名から操作indexを引く辞書を作っておきます。
この辺はもうちょっと効率化ができると思いますが、それをやる余力がありませんでした。
move_names_to_index = {move_name: i for i, move_name in enumerate(move_names)}
class Search:
def __init__(self, state):
self.initial_state = state
self.current_solution_ph1 = []
def depth_limited_search_ph1(self, co_index, eo_index, e_comb_index, depth):
if depth == 0 and co_index == 0 and eo_index == 0 and e_comb_index == 0:
return True
if depth == 0:
return False
# 枝刈り
if max(co_eec_prune_table[co_index][e_comb_index], eo_eec_prune_table[eo_index][e_comb_index]) > depth:
return False
prev_move = self.current_solution_ph1[-1] if self.current_solution_ph1 else None
for move_name in move_names:
if not is_move_available(prev_move, move_name):
continue
self.current_solution_ph1.append(move_name)
move_index = move_names_to_index[move_name]
next_co_index = co_move_table[co_index][move_index]
next_eo_index = eo_move_table[eo_index][move_index]
next_e_comb_index = e_combination_table[e_comb_index][move_index]
if self.depth_limited_search_ph1(next_co_index, next_eo_index, next_e_comb_index, depth - 1):
return True
self.current_solution_ph1.pop()
def start_search(self, max_length=20):
co_index = co_to_index(self.initial_state.co)
eo_index = eo_to_index(self.initial_state.eo)
e_combination = [1 if e in (0, 1, 2, 3) else 0 for e in self.initial_state.ep]
e_comb_index = e_combination_to_index(e_combination)
depth = 0
while depth <= max_length:
print(f"# Start searching phase 1 length {depth}")
if self.depth_limited_search_ph1(co_index, eo_index, e_comb_index, depth):
return " ".join(self.current_solution_ph1)
depth += 1
return None
状態をindexで扱うようにして、状態遷移表と枝刈り表を使うようにした以外は、中編とほとんど同じです。
枝刈りでは、COとE列エッジを揃える最短手数と、EOとE列エッジを揃える最短手数のいずれかでも、残りのdepthで揃わないことがわかったらそこで打ち切るようにしています。
Phase1探索プログラムの動作確認
適当なランダムスクランブル R' U' F R' B' F2 L2 D' U' L2 F2 D' L2 D' R B D2 L D2 F2 U2 L R' U' F を用意しました。
scramble = "R' U' F R' B' F2 L2 D' U' L2 F2 D' L2 D' R B D2 L D2 F2 U2 L R' U' F"
scrambled_state = scamble2state(scramble)
search = Search(scrambled_state)
start = time.time()
solution = search.start_search()
print(f"Finished! ({time.time() - start:.2f} sec.)")
if solution:
print(f"Solution: {solution}.")
else:
print("Solution not found.")
# 出力
Finished! (0.0402 sec.)
Solution: U2 F2 D' U' L U' R' D2 R' F.
10手のPhase 1の解が0.04秒で見つかりました。中編では、6手の解を探すのに0.2秒かかっていたのでPhase1のみとはいえだいぶ高速化できていそうです。
スクランブルに続けて、見つけた解 U2 F2 D' U' L U' R' D2 R' F を回すと、
確かに、CO, EOが全部揃っており、E列エッジはE列エッジ内でしか移動していない状態、Phase 1の完成状態になっていることがわかります。
Phase 1を見つけたら、引き続き Phase 2を探すようにする
Phase 2についても、操作名から操作indexを引く辞書を作っておきます。
move_names_to_index_ph2 = {move_name: i for i, move_name in enumerate(move_names_ph2)}
Phase 1 の解を発見したら、Phase 1と同じように実装した Phase 2の探索を呼び出します。
class Search:
def __init__(self, state):
self.initial_state = state
self.current_solution_ph1 = []
self.current_solution_ph2 = []
self.max_solution_length = 9999
self.start = 0
def depth_limited_search_ph1(self, co_index, eo_index, e_comb_index, depth):
if depth == 0 and co_index == 0 and eo_index == 0 and e_comb_index == 0:
state = self.initial_state
for move_name in self.current_solution_ph1:
state = state.apply_move(moves[move_name])
return self.start_phase2(state)
if depth == 0:
return False
# 枝刈り
if max(co_eec_prune_table[co_index][e_comb_index], eo_eec_prune_table[eo_index][e_comb_index]) > depth:
return False
prev_move = self.current_solution_ph1[-1] if self.current_solution_ph1 else None
for move_name in move_names:
if not is_move_available(prev_move, move_name):
continue
self.current_solution_ph1.append(move_name)
move_index = move_names_to_index[move_name]
next_co_index = co_move_table[co_index][move_index]
next_eo_index = eo_move_table[eo_index][move_index]
next_e_comb_index = e_combination_table[e_comb_index][move_index]
if self.depth_limited_search_ph1(next_co_index, next_eo_index, next_e_comb_index, depth - 1):
return True
self.current_solution_ph1.pop()
def depth_limited_search_ph2(self, cp_index, udep_index, eep_index, depth):
if depth == 0 and cp_index == 0 and udep_index == 0 and eep_index == 0:
return True
if depth == 0:
return False
# 枝刈り
if max(cp_eep_prune_table[cp_index][eep_index], udep_eep_prune_table[udep_index][eep_index]) > depth:
return False
if self.current_solution_ph2:
prev_move = self.current_solution_ph2[-1]
elif self.current_solution_ph1:
prev_move = self.current_solution_ph1[-1]
else:
prev_move = None
for move_name in move_names_ph2:
if not is_move_available(prev_move, move_name):
continue
self.current_solution_ph2.append(move_name)
move_index = move_names_to_index_ph2[move_name]
next_cp_index = cp_move_table[cp_index][move_index]
next_udep_index = ud_ep_move_table[udep_index][move_index]
next_eep_index = e_ep_move_table[eep_index][move_index]
if self.depth_limited_search_ph2(next_cp_index, next_udep_index, next_eep_index, depth - 1):
return True
self.current_solution_ph2.pop()
def start_search(self, max_length=30):
self.start = time.time()
self.max_solution_length = max_length
co_index = co_to_index(self.initial_state.co)
eo_index = eo_to_index(self.initial_state.eo)
e_combination = [1 if e in (0, 1, 2, 3) else 0 for e in self.initial_state.ep]
e_comb_index = e_combination_to_index(e_combination)
depth = 0
while depth <= self.max_solution_length:
if self.depth_limited_search_ph1(co_index, eo_index, e_comb_index, depth):
return " ".join(self.current_solution_ph1) + " . " + " ".join(self.current_solution_ph2)
depth += 1
return None
def start_phase2(self, state):
cp_index = cp_to_index(state.cp)
udep_index = ud_ep_to_index(state.ep[4:])
eep_index = e_ep_to_index(state.ep[:4])
depth = 0
while depth <= self.max_solution_length - len(self.current_solution_ph1):
if self.depth_limited_search_ph2(cp_index, udep_index, eep_index, depth):
return True
depth += 1
動作確認
適当なスクランブルを与えて検証してみます。
scramble = "R' U' F R' B' F2 L2 D' U' L2 F2 D' L2 D' R B D2 L D2 F2 U2 L R' U' F"
scrambled_state = scamble2state(scramble)
search = Search(scrambled_state)
start = time.time()
solution = search.start_search()
print(f"Finished! ({time.time() - start:.4f} sec.)")
if solution:
print(f"Solution: {solution}.")
else:
print("Solution not found.")
# 出力
Finished! (0.0650 sec.)
Solution: U2 F2 D' U' L U' R' D2 R' F . U' B2 L2 D2 L2 U' B2 U B2 U' B2 U' L2.
0.065秒で、Phase 1が10手、Phase2が13手で、合計23手の解が見つかりました。
ここまでの実装の課題
0.065秒で適当なルービックキューブが解けたので、目標達成という感じもするのですが、イマイチな点があります。
こういうケースを考えてみます。
scramble = "U F2 D R' U2 R"
scrambled_state = scamble2state(scramble)
search = Search(scrambled_state)
start = time.time()
solution = search.start_search()
print(f"Finished! ({time.time() - start:.4f} sec.)")
if solution:
print(f"Solution: {solution}.")
else:
print("Solution not found.")
# 出力
Finished! (0.0114 sec.)
Solution: L' B2 L . R2 D L2 U B2 U L2 D2 U' R2.
6手のスクランブルに対して、13手の解を見つけてきました。
中編の総当りコードだったら6手のスクランブルに6手の解法をすぐ出せるのに、後編のコードがこれではイマイチ感が否めません。
これは、Phase 1の解を一番最初に見つけた解に固定してしまって、そこから続くPhase 2の解しか探していないのが原因です。
このプログラムは、Phase 1の解として最初に L' B2 L を見つけ、それを採用してPhase2の探索に移行し10手のPhase 2の解を発見し13手の解を出力しました。
しかし、実は R' U2 R もPhase 1の解となります。
後者を採用していれば、Phase 2は D' F2 U' となり、3手+3手＝6手の解が見つかっていました。
このように、Phase 1を最初に見つけた最適解に固定してしまうと、イマイチな場合があります。
Phase 1の最適解が他にもありそちらを採用したほうが良いというケース以外にも、Phase 1が最適より1手長くても、Phase 2が2手以上短ければお得になります。
なので、時間の許す範囲で余分に探索を続けることで、より良い解を見つけられる可能性があります。
それを実装してみます。
より良い解を求めて探索を続ける実装 (最終成果物)
次のように改修し、指定した時間内で、Phase 1が最適でなくとも、Phase 2がより短くなる解を探すようにしました。
すでに見つかっている解より長くなるとわかった場合は、Phase 2の探索を打ち切ることで効率化しています。
また、許容する最大解法長を予め渡しておくことによって、1つ目の解がすごく長くなってしまうようなケースを回避し効率化しています。23手くらいの解はすぐ見つかることがわかっているので、23手を許容する最大解法長としておけば、例えば最初に見つかったPhase 1が11手、それに続くPhase 2が12手まで探しても見つからなければとっとと別のPhase 1を探しに行くということができます。
import timeout_decorator
class Search:
def __init__(self, state):
self.initial_state = state
self.current_solution_ph1 = []
self.current_solution_ph2 = []
self.max_solution_length = 9999
self.best_solution = None
self.start = 0
def depth_limited_search_ph1(self, co_index, eo_index, e_comb_index, depth):
if depth == 0 and co_index == 0 and eo_index == 0 and e_comb_index == 0:
last_move = self.current_solution_ph1[-1] if self.current_solution_ph1 else None
if last_move is None or last_move in ("R", "L", "F", "B", "R'", "L'", "F'", "B'"):
# print(f"# Found phase 1 solution {" ".join(self.current_solution_ph1)} (length={len(self.current_solution_ph1)})")
state = self.initial_state
for move_name in self.current_solution_ph1:
state = state.apply_move(moves[move_name])
return self.start_phase2(state)
if depth == 0:
return False
# 枝刈り
if max(co_eec_prune_table[co_index][e_comb_index], eo_eec_prune_table[eo_index][e_comb_index]) > depth:
return False
prev_move = self.current_solution_ph1[-1] if self.current_solution_ph1 else None
for move_name in move_names:
if not is_move_available(prev_move, move_name):
continue
self.current_solution_ph1.append(move_name)
move_index = move_names_to_index[move_name]
next_co_index = co_move_table[co_index][move_index]
next_eo_index = eo_move_table[eo_index][move_index]
next_e_comb_index = e_combination_table[e_comb_index][move_index]
found = self.depth_limited_search_ph1(next_co_index, next_eo_index, next_e_comb_index, depth - 1)
self.current_solution_ph1.pop()
def depth_limited_search_ph2(self, cp_index, udep_index, eep_index, depth):
if depth == 0 and cp_index == 0 and udep_index == 0 and eep_index == 0:
# print(f"# Found phase 2 solution {" ".join(self.current_solution_ph2)} (length={len(self.current_solution_ph2)})")
solution = " ".join(self.current_solution_ph1) + " . " + " ".join(self.current_solution_ph2)
print(
f"Solution: {solution} "
f"({len(self.current_solution_ph1) + len(self.current_solution_ph2)} moves) "
f"in {time.time() - self.start:.5f} sec.")
self.max_solution_length = len(self.current_solution_ph1) + len(self.current_solution_ph2) - 1
self.best_solution = solution
return True
if depth == 0:
return False
# 枝刈り
if max(cp_eep_prune_table[cp_index][eep_index], udep_eep_prune_table[udep_index][eep_index]) > depth:
return False
if self.current_solution_ph2:
prev_move = self.current_solution_ph2[-1]
elif self.current_solution_ph1:
prev_move = self.current_solution_ph1[-1]
else:
prev_move = None
for move_name in move_names_ph2:
if not is_move_available(prev_move, move_name):
continue
self.current_solution_ph2.append(move_name)
move_index = move_names_to_index_ph2[move_name]
next_cp_index = cp_move_table[cp_index][move_index]
next_udep_index = ud_ep_move_table[udep_index][move_index]
next_eep_index = e_ep_move_table[eep_index][move_index]
found = self.depth_limited_search_ph2(next_cp_index, next_udep_index, next_eep_index, depth - 1)
self.current_solution_ph2.pop()
if found:
return True
def start_search(self, max_length=30, timeout=3):
@timeout_decorator.timeout(timeout)
def target():
self.start_search1(max_length)
try:
target()
return self.best_solution
except timeout_decorator.TimeoutError:
return self.best_solution
def start_search1(self, max_length=30):
self.start = time.time()
self.max_solution_length = max_length
co_index = co_to_index(self.initial_state.co)
eo_index = eo_to_index(self.initial_state.eo)
e_combination = [1 if e in (0, 1, 2, 3) else 0 for e in self.initial_state.ep]
e_comb_index = e_combination_to_index(e_combination)
depth = 0
while depth <= self.max_solution_length:
# print(f"# Start searching phase 1 length {depth}")
if self.depth_limited_search_ph1(co_index, eo_index, e_comb_index, depth):
# print(f"Found solution ({len(self.current_solution_ph1) + len(self.current_solution_ph2)} length)")
pass
depth += 1
return None
def start_phase2(self, state):
cp_index = cp_to_index(state.cp)
udep_index = ud_ep_to_index(state.ep[4:])
eep_index = e_ep_to_index(state.ep[:4])
depth = 0
while depth <= self.max_solution_length - len(self.current_solution_ph1):
if self.depth_limited_search_ph2(cp_index, udep_index, eep_index, depth):
return True
depth += 1
先程の例を試してみます。
scramble = "U F2 D R' U2 R"
scrambled_state = scamble2state(scramble)
search = Search(scrambled_state)
start = time.time()
solution = search.start_search()
print(f"Finished! ({time.time() - start:.4f} sec.)")
if solution:
print(f"Solution: {solution}.")
else:
print("Solution not found.")
# 出力
Solution: L' B2 L . R2 D L2 U B2 U L2 D2 U' R2 (13 moves) in 0.01765 sec.
Solution: R' U2 R . D' F2 U' (6 moves) in 0.01891 sec.
Finished! (0.0251 sec.)
Solution: R' U2 R . D' F2 U'.
Phase 1 の解として、 L' B2 L を見つけ、それに対する Phase 2の探索を行い13手の解を発見した後、引き続きPhase 1の他の解をとそれに続くPhase 2の解を探すことで、6手の解を見つけることができました。
より長いPhase 1の解が、全体として短い解を生む例
より長いPhase 1の解が全体として短い解を生む例ですが、ほとんどの場合はそのようになります。
適当なスクランブルに対して、90秒探索を続けさせてみます。
scramble = "R' U' F R' B' F2 L2 D' U' L2 F2 D' L2 D' R B D2 L D2 F2 U2 L R' U' F"
scrambled_state = scamble2state(scramble)
search = Search(scrambled_state)
start = time.time()
solution = search.start_search(max_length=22, timeout=90)
print(f"Finished! ({time.time() - start:.5f} sec.)")
if solution:
print(f"Solution: {solution}.")
else:
print("Solution not found.")
# 出力
Solution: U2 F2 D' U' F2 L' U B2 D' B . F2 U L2 U' B2 U' F2 L2 D2 F2 L2 D (22 moves) in 0.06983 sec.
Solution: D F' D2 U L' B2 D' B' F2 D' R' . F2 D L2 U R2 B2 D U R2 F2 (21 moves) in 1.21992 sec.
Solution: L2 B2 D B2 R2 U L' U' L2 D F' . L2 B2 D R2 F2 U' F2 U2 B2 (20 moves) in 2.37836 sec.
Solution: U2 R2 F2 D B D2 L2 B' L D2 F2 D L . U2 L2 D2 F2 D2 (18 moves) in 80.23398 sec.
Finished! (90.00034 sec.)
Solution: U2 R2 F2 D B D2 L2 B' L D2 F2 D L . U2 L2 D2 F2 D2.
最初は、Phase 1が10手、それに続くPhase 2が12手の計22手の解を見つけした。
その後、より長いPhase 1の解とそれに対応するPhase 2の解を探し続け、Phase 1は13手かかるが、Phase 2が5手で終わる、18手の解を見つけることができました。
このように、22,23手ほどの解は1秒未満で見つけることができ、時間をかければかけるほど短い解が見つかるような実装ができました。
この実装は、非常に長い時間をかければ必ず最適解を見つけることができます。
何故ならば、Phase 1を最大20手まで伸ばしていけば、Phase 2が0手となるようなPhase 1の解がいつかは見つかり、それは全体としての最適解だということだからです。
さらなる効率化
さらなる効率化として、以下のようなものがあるようです。私はまだ良くわかっていません。
- 対称性の考慮
- 対称性を考慮することで、探索を効率化できるようです
- 今回Phase 1は、U/D面のみ90度許容となる状態を完成状態としましたが、L/Rとしても同じ話ができるのでそちらも同時に考えた方が効率的みたいな話でしょうか?
- cf. http://kociemba.org/cube.htm
- pre-move
- 与えられたスクランブル冒頭になにか別の回転記号Xを付けたようなスクランブル状態を考え、解法を見つけた後、解法の最後にXをつければ、元のスクランブルの解法になるという技です
公式大会のスクランブルアルゴリズムは、上記の両方の工夫が入っているそうです。他の工夫もあると思います。
そのあたりも記事にできればと思っていますが、時間取れそうもないのでいつになるか、実現されるかはわかりません (誰か代わりにお願いします)。
おわりに
前編で目標として掲げた、ルービックキューブを大体1秒位 & 20手強くらいで解くプログラムが無事完成しました!
完結までアドベントカレンダー3年分となってしまいましたが、最後までお付き合いいただきありがとうございました。
ルービックキューブや類似のパズルを解くプログラムを書こうとしているどなたかのお役に立てれば幸いです。
コードはCC0ライセンスで扱っていただいて大丈夫です。
後編は息切れしてかなり説明を端折ってしまいましたが、中級編で飽き足らず後編まで見てくださった方ならば、一応考え方と動くプログラムは示してあるので解読して理解し応用してくださるのではないかなと思っています。
わからないところや、バグなどあれば、コメントをお寄せいただければと思います。