ビネクラ杯について

こちらのリンクでご確認ください(ビネクラ社の皆様。楽しませていただきました！ありがとうございました)。
https://vigne-cla.com/category/optimization/vigne_cla_cup/
私自身はxtxtxという名前で出し、517という4位タイの結果を出しました。530以上のスコアを出した方々の戦略はわかりませんが、参考にしてもらえたら嬉しいです。
また、計算時間の短縮のためにCythonを用いたので、そこら辺も参考にしてもらえたらと感じています。

戦略

1. まずはマンハッタン距離順に並び替え

リンク先の18-5, 18-6の記事にあったように、for文のbreakをcontinueに変えたり、マンハッタン距離の短い順に並び替えたりすることによりスコアを470程度まで伸ばせることがわかります。ここからどのようにスコアを伸ばしていくかが、個性が出る部分であるように感じました。

マンハッタン距離順に並べるのは、Python組み込みのsorted関数を使うと18-6の記事より簡単に、

diff = generators - equipments
l = sorted(range(4000), key=lambda i: np.linalg.norm(diff[i], ord=1))

とすることでできます。
慣れてない方のために簡単に説明すると、まず、一行目で装置から発電機へのベクトル(相対的な座標)を並べた配列を作っています。
次に出てくるnp.linalg.normという関数は、与えられたベクトルのノルム(長さのようなもの)を計算してくれます。特にそれにord=1と引数を与えるとベクトルの根元から先までのマンハッタン距離を計算できます。
また、sortedは第1引数に受け取ったiterable(listのようなもの)の各要素をkeyに与えられた関数に入れた結果で並び替えます。
この2行だけで並び替えは完了です。

2. 探索の長さを変えてみる

次に私が考えてみたことは、いくらマンハッタン距離が近い順に並べたかと言って、他の道のせいで遠回りすることもあるだろうということです。つまり、マンハッタン距離の近い順に探索はするものの、他の道のせいで遠回りするようなコースを通るときには、(無理やりそんな道をつなげてもあとあと邪魔になってしまうだろうということで)後回しにしたほうがいいだろうと思ったのです。
確実に道が短い順につなげるためには、どうしたらいいでしょうか。少し考えてみてください。

正解の一例は、find_route関数をいじって探索する長さの上限を決められるようにしておいて、その上限を変化させて探索を繰り返すということです。
さきほど紹介した記事では999まで探索していましたが、それをあえて、max_rangeという変数でおいてあげてそれを変化させるということです。

また、今回は、この関数を繰り返し用いるのでその高速化のためにCythonを使ってコンパイルしておきます。

%load_ext cython

でCythonを読み込んだ後に、

%%cython -a
import numpy as np
import matplotlib.pyplot as plt
from scipy.ndimage.filters import minimum_filter, maximum_filter
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.cm as cm
import random
from tqdm import tqdm
# from tqdm import tqdm_notebook as tqdm
cimport numpy as np
#============================
# ひとつのルートを見つける関数
#============================
def find_route(np.ndarray[np.uint8_t, cast=True, ndim=3] barrier,
               np.ndarray[np.int_t, ndim=1] start,
               np.ndarray[np.int_t, ndim=1] goal, int max_range=999):
    #########################
    # 前処理
    #########################
    #たてとよこ
    h, w, l = len(barrier), len(barrier[0]), len(barrier[0][0])
    #コスト
    cdef np.ndarray[np.int_t, ndim=3] cost = np.zeros((h, w, l), dtype=int) + max_range
    #コストが書き込まれて探索が終了したマス（bool）
    cdef np.ndarray[np.uint8_t, cast=True, ndim=3] done = np.zeros((h, w, l), dtype=bool)
    
    #プーリング用のフィルタ
    cdef np.ndarray[np.int_t, ndim=3] g
    g = np.array([[[0, 0, 0],
                   [0, 1, 0],
                   [0, 0, 0]],
                  [[0, 1, 0],
                   [1, 1, 1],
                   [0, 1, 0]],
                  [[0, 0, 0],
                   [0, 1, 0],
                   [0, 0, 0]]], dtype=int)
    
    # 3の実装のために、startとgoalのbarrierを変えておく。
    barrier[start[0], start[1], start[2]] = False
    barrier[goal[0],goal[1], goal[2]] = False
    cost[start[0], start[1], start[2]] = 0
    done[start[0], start[1], start[2]] = True
    
    #########################
    # 幅優先探索
    #########################
    for i in range(1, max_range):
    
        #次に進出するマスのbool
        done_next = maximum_filter(done, footprint=g) * ~done
        #print('done_next\n{}'.format(done_next))
        
        #次に進出するマスのcost
        cost_next = minimum_filter(cost, footprint=g) * done_next
        cost_next[done_next] += 1
        #print('cost_next\n{}'.format(cost_next))
        
        #costを更新
        cost[done_next] = cost_next[done_next]
        #ただし障害物のコストは999とする
        cost[barrier] = max_range
        #print('cost\n{}'.format(cost))
        
        #探索終了マスを更新
        done[done_next] = done_next[done_next]
        #ただし障害物は探索終了としない
        done[barrier] = False
        #print('done\n{}'.format(done))
        
        #表示
        #show(i)
        
        #終了判定
        if done[goal[0], goal[1], goal[2]] == True:
            break
    
    #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    # 【追加】 ここで一旦終了判定、ゴールのコストが999であれば、ルートが見つかっていないので強制終了
    #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    if cost[goal[0], goal[1], goal[2]] == max_range:
        return None
    
    #########################
    # ゴールから逆順でルート計算
    #########################
    point_now = tuple(goal)
    cost_now = cost[goal[0], goal[1], goal[2]]
    route = [goal]
    #print('route\n{}'.format(route))
    while cost_now > 0:
        #for i in random.sample(list(range(6)), 6):
        for i in [2,3,0,1,4,5]: # この並び替えは気分です。ここを変えればすぐにランダムにできます(でもランダムの方がスコア下がりがちだったので…)。
            if i == 0:
                #x-から来た場合
                try:
                    if point_now[0] <= 0 : raise IndexError
                    if cost[point_now[0] - 1, point_now[1], point_now[2]] == cost_now - 1:
                        #更新
                        point_now = (point_now[0] - 1, point_now[1], point_now[2])
                        cost_now = cost_now - 1
                        #記録
                        route.append(point_now)
                except: pass
            elif i == 1:
                #x+から来た場合
                try:
                    if cost[point_now[0] + 1, point_now[1], point_now[2]] == cost_now - 1:
                        #更新
                        point_now = (point_now[0] + 1, point_now[1], point_now[2])
                        cost_now = cost_now - 1
                        #記録
                        route.append(point_now)
                except: pass
            elif i == 2:
                #y-から来た場合    
                try:
                    if point_now[1] <= 0 : raise IndexError
                    if cost[point_now[0], point_now[1] - 1, point_now[2]] == cost_now - 1:
                        #更新
                        point_now = (point_now[0], point_now[1] - 1, point_now[2])
                        cost_now = cost_now - 1
                        #記録
                        route.append(point_now)
                except: pass
            elif i == 3:
                #y+から来た場合
                try:
                    if cost[point_now[0], point_now[1] + 1, point_now[2]] == cost_now - 1:
                        #更新
                        point_now = (point_now[0], point_now[1] + 1, point_now[2])
                        cost_now = cost_now - 1
                        #記録
                        route.append(point_now)
                except: pass
            elif i == 4:
                #z-から来た場合    
                try:
                    if point_now[2] <= 0: raise IndexError
                    if cost[point_now[0], point_now[1], point_now[2] - 1] == cost_now - 1:
                        #更新
                        point_now = (point_now[0], point_now[1], point_now[2] - 1)
                        cost_now = cost_now - 1
                        #記録
                        route.append(point_now)
                except: pass
            else:
                #z+から来た場合
                try:
                    if cost[point_now[0], point_now[1], point_now[2] + 1] == cost_now - 1:
                        #更新
                        point_now = (point_now[0], point_now[1], point_now[2] + 1)
                        cost_now = cost_now - 1
                        #記録
                        route.append(point_now)
                except: pass
    
    #ルートを逆順にする
    route = route[::-1]
    #print('route\n{}'.format(route))
    
    return route

とします。NumPy配列や、関数自体の引数は型をつけておきました(全部につけなかったのは、自動である程度は型推論してくれますし、経験的にその手間に合う効果が得られないためです)¹。

これで、探索範囲を変化させることができるようになりました。ここでMain処理を
for max_range in range(80)などとしながら中でlのループをさせると最大で500弱程度のスコアが出ました。
また、そもそも1つ1つの探索を短くできるので、計算時間も短縮できます！(Cythonとは関係ない効果)

3. 発電機と装置をあえて避けるようなコースにする

ここからスコアを伸ばすにはどうしたらいいか、と考えたときに、初期の段階では発電機と装置も通らないような道を通るようにしてあげたらいいのではと考えました。なぜなら、発電機か装置を道で埋めてしまうと、もうその装置を繋ぐ可能性は無くなってしまうからです。ただ、これをやりすぎると本来障害物としないはずの発電機と装置を障害物とみなすようなものなので、かえってスコアが下がります。だから、l[:300](つまり近い順に300個)のみについて発電機と装置を通らないようにしました。300という数字は実験して決めました。
ここで、番号l[:300]に含まれる発電機と装置の位置のbarrierの要素について考えてみます。あとでそこにある発電機、装置を置く分には大丈夫なのでTrueのなかでも、通常の道によるTrueと区別する必要があります。そこで前者の値を2、後者を1、Falseを0とおくことにより、これらを区別しました。find_route関数内では1か2かの区別はないので、barrier.astype(bool)を引数に渡しました。

ただ、このままでは狙い通りに動かないので、(2の値になっていることもある)始点と終点について、find_route関数内で初めにbarrierの値をFalseとしました(この記事では2の時点でもう反映してあります)。

#============================
# メイン処理
#============================
l = sorted(range(4000), key=lambda i: np.linalg.norm(diff[i], ord=1))
#通った場所（Falseで初期化）
barrier = np.zeros((20, 20, 20), dtype=np.int)

for i in l[:300]:
    barrier[generators[i, 0], generators[i, 1], generators[i, 2]] = 2
    barrier[equipments[i, 0], equipments[i, 1], equipments[i, 2]] = 2

#記録用の配列
score = 0
routes = []
results = []
# l = l[:1000] + l[1000:][::-1] いろいろ実験した形跡。少し意味を考えてみてください
# random.seed(41)
# ll = random.sample(l[1000:4000], 3000)
for max_range in tqdm(list(range(5, 40)) + [80]): # 39までは1きざみで、そのあと80。適宜変更してみてください。
    
    for num, i in enumerate(l):

        #これからつなぐスタートやゴールが、すでに障害物となっていたら終わり
        if barrier[generators[i, 0], generators[i, 1], generators[i, 2]] == 1:
    #         print('double-g {}'.format(num))
            continue
        if barrier[equipments[i, 0], equipments[i, 1], equipments[i, 2]] == 1:
    #         print('double-e {}'.format(num))
            continue

        #ルートを見つける
        route = find_route(barrier.astype(bool), generators[i], equipments[i],
                          max_range)


        #ルートが見つからなかったら終わり
        if route is None:
            continue

        #ルートが見つかったら、スコアを記録する
        score += 1
        routes.append(route)
        results.append([max_range, len(route), num, i]) # あとで結果を見るため。

        #通った場所を障害物とする
        for [x, y, z] in route:
            barrier[x, y, z] = 1
            
#     if max_range > 45: 
#         random.shuffle(l) もしかしてshuffleするといいかな、と思った
    
print('\nscore:\n{}'.format(score))

これでスコア517が出ます！
今のところl[:300]というシンプルな方法で選んでいますが、具体的に何番とか指定してあげるとよりいいスコアが出せる気がします(時間がなくてできませんでしたが)。試してみてください！

今回は入門ということで深くは立ち入りませんでした。さらに速度を求めるのであればMain処理まで1つのセルに入れたり、listをC++のvectorで書き換えたりいろいろできます。
これで体感としては**30〜50%**くらい早くなったような感じです(このくらいの書き換えで済むならやって損はない！)。 ↩