TSP（巡回セールスマン問題）：その２

前回のTSP（巡回セールスマン問題）の続きです。前回のKaggleコンペでは都市が197769もあり、計算するにもかなり時間がかかったので、都市数を減らしたサンプルで最適化アルゴリズムについて確かめてみました。

総当たりで計算すると都市数が10を超えるあたりからきつくなっていくので、真のスコアと比較検証するには少なめの都市数にする必要があります。


ライブラリ：

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from numba import jit, prange
from itertools import permutations
from tqdm import tqdm_notebook as tqdm
np.random.seed(seed=42)

まずは必要そうなライブラリ。numbaでどれだけ高速になるかはわかりませんが一応入れておきます。組み合わせのパターンを求めるにはitertools.permutationsを用いました。


都市数や座標の設定：

num = 10                                      # the number of cities
path = np.concatenate([np.arange(num), [0]])  # path starts from 0 and ends with 0
X = np.random.random(num)                     # X-coordinates of cities
Y = np.random.random(num)                     # Y-coordinates of cities
XY = X + Y * 1j                               # complex numbers of X and Y

都市は合計10箇所、0〜9までの番号を割り当てておきます。最後に0に戻るので経路としては
[0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 0]
の11要素になります。
それぞれにランダムでXとYの座標を与え、距離を求めるための複素数XYも用意しておきます。


スコア関数：

def score(path, a, b):
    c = XY[path[a: b+1]]
    sc = np.sum(np.abs(np.diff(c)))
    return sc

都市aからbに至るスコア（道のり）を求めるための関数を用意しておきます。複素数XYの差分をnp.abs()を使って求めます。


経路のプロット：

def plot_path(path):
    PX = [X[i] for i in path]
    PY = [Y[i] for i in path]
    plt.figure(figsize=(5,5))
    plt.axis('equal')
    plt.plot(PX, PY)
    plt.scatter(PX, PY)
    plt.scatter(PX[0], PY[0], s=80, c='r', marker='o')
    for i in range(num):
        plt.text(PX[i], PY[i]+0.03, s=i, fontsize=12, color='gray')
    print('Total Score: ', score(path, 0, len(path)))

開始点を赤点表示。巡回順のナンバリングもしておきます。

初期経路はランダムなので表示させるとこのような感じになります。
4と5がほぼ重なっていて見にくいですがこのままで。


真のスコアを求める：
総当たりで計算する場合、12都市以上になると何時間もかかるため、都市数を指定できるような関数にしました。10都市の場合は巡回すると11都市になり、開始点と到着点の二点を除いた9都市の組み合わせすべてについて計算します。9都市なら3分、10都市なら30分以上かかります。

def opt_n(path, n):
    for i in tqdm(range(1, len(path)-n)):
        perm = np.array(path[i:i+n])
        combs = np.array(list(permutations(perm)))
        sc1 = score(path, i-1, i+n+1)
        sc3 = np.inf
        for comb in combs:
            copy = np.concatenate([path[:i], comb, path[i+n:]])
            sc2 = score(copy, i-1, i+n+1)
            if sc1 > sc2:
                if sc3 > sc2:
                    sc3 = sc2
                    path = copy.copy()
                    print('Improved to: ', score(path, 0, num))
    return path

二方向の欲張り法：
ランダムな初期経路に対して最適化していってもいいかもしれませんが、欲張り法でもう少しまともな経路をつくっておきます。
スタート地点から最も近い都市をつなぎ続けていきます。スタート地点から二方向に都市をつなげていき、その二つの経路を中点でつなぎます。

@jit('i8[:](i8[:], i8, i8)')
def closest_num(path, num, index):
    copy = path.copy()
    d = np.array([np.abs(XY[index] - XY[i]) for i in copy])
    arg = d.argsort()
    closest_index = np.array([copy[i] for i in arg])[:num]
    return closest_index

@jit('i8[:](i8[:])', parallel=True)
def find_path(path):
    copy = path[1:-1].copy()
    half = len(copy)//2 + 1
    path1 = np.array([0])
    path2 = np.array([0])
    for i in tqdm(prange(half)):
        if len(copy) > 0:
            c1 = closest_num(copy, 1, path1[-1])
            path1 = np.append(path1, c1)
            copy = np.delete(copy, np.where(copy==c1))
        if len(copy) > 1:  
            c2 = closest_num(copy, 1, path2[-1])
            path2 = np.append(path2, c2)
            copy = np.delete(copy, np.where(copy==c2))
        
    print('Path_1:', path1, '\nPath_2:', path2[::-1])
    return np.concatenate([path1, path2[::-1]])

どのくらい効果的かはわかりませんが、一応使ってみます。

ランダムよりはまともな経路になりますが、まだまだ改良の余地があります。


二点間反転：
欲張り法で得た経路を最適化するアルゴリズムです。任意の二点間を選び、その区間を反転して繋ぎ直します。

@jit('i4[:](i4[:], i4)')
def opt_n_flip(path, n):
    for i in range(0, len(path)-n-1):
        n1 = i
        n2 = n1 + n
        sc1 = score(path, n1-1, n2+1)
        rev = path[n1: n2+1][::-1]
        copy = np.concatenate([path[:n1], rev, path[n2+1:]])
        sc2 = score(copy, n1-1, n2+1)
        if sc1 > sc2:
            path = copy.copy()
            print('Improved to: ', score(path, 0, len(path)-1))
    return path

@jit(parallel=True)
def opt_flip_loop(path, n, num):
    for i in tqdm(prange(n, n+num)): 
        path = opt_n_flip(path, i)
    return path

まずは隣り合わせの二都市、次は一つ離れた二都市、そして二つ離れた二都市という具合にN個分離れた二都市までを評価して、向上が見られれば随時更新していきます。

最適化の結果、都市数が少ないためか、ほぼ正解にたどり着いたようです。

一点移動挿入：
or-optというアルゴリズムに近い方法です。経路から一つだけずれている都市を最適な場所に移動してくれます。

def longest(path):
    sc_arr = - np.array([np.sum(np.abs(np.diff(np.array(XY[path[i-1: i+1]])))) for i in tqdm(range(1, len(path)-1))])
    arg = sc_arr.argsort()    
    return arg

def closest_n(path, n, index):
    copy = path[1:-1].copy()
    if path[index] in copy:
        copy = np.delete(copy, np.where(copy==path[index]))
    d = np.array([np.abs(XY[index] - XY[i]) for i in copy])
    arg = d.argsort()
    closest_index = np.array([copy[i] for i in arg])[:n]
    return closest_index

@jit('i8[:](i8[:])', parallel=True)
def or_opt(path):
    lon = longest(path)
    print('Initial Score: ', score(path, 0, num))
    for i in tqdm(lon):
        sc3 = np.inf
        c5 = closest_n(path, 10, path[i])
        for c in c5:
            j = np.where(path==c)[0][0]
            copy = np.delete(path, i)
            copy = np.insert(copy, j, path[i])
            if i >= j:
                sc1 = score(path, j-1, i+1)
                sc2 = score(copy, j-1, i+1)
            else:
                sc1 = score(path, i-1, j+1)
                sc2 = score(copy, i-1, j+1)
            if sc1 > sc2:
                if sc3 > sc2:
                    sc3 = sc2
                    path = copy.copy()
                    print('Improved to: ', score(path, 0, num))
    return path

まず二点間距離の降順リストをつくり、その都市の最も近くにある都市（上位５都市）の次の位置へ移動しスコアが向上すれば更新します。

最適化しても前回と結果は同じです。おそらくこれが正解の経路。

正解の経路：
最後に真のスコアを求めるためにつかったopt_nで再度最適化させます（都市数が多い場合はせいぜい10都市までの組み合わせによる最適化）。今回の場合は都市数が少ないのでopt_nに9を代入すれば正解が得られます（以下）。組み合わせは9!=362880通り（約3分）。

全体の流れ：
・ランダムに経路を生成
・欲張り法で経路を生成（find_path）
・二点間反転（opt_n_flip）
・一点移動挿入（or_opt）
・組み合わせ最適化（opt_n）

都市数が少ない場合は二点間反転だけでも真の経路になる場合がありました。ひとつのアルゴリズムを一度実行しただけでは最適化不十分なので、スコアが向上しなくなるまで複数回繰り返します。
都市数を２０、３０と増やしていくと、真の経路を求めるには時間がかかりすぎるので比較評価することができません。よって最適化によってスコアを徐々に向上していくしかありません。
複数のアルゴリズムをスコア向上しなくなるまで実行することで、ぱっと見は最適解になっているように見えますが、ランダムで生成した経路を最適化した場合と欲張り法のあとで最適化した場合とで結果が異なるので初期経路次第という感じです。今回は欲張り法でしか初期経路を生成していないので、もう少し他の方法も試してみる必要がありそうです。焼きなまし法というのもあるらしく、最初から最短都市や最少スコアを優先するアルゴリズムにするのではなく、許容値を設定して多少スコアオーバーしても採用するなどしたほうがいいかもしれません（この部分に関しては次回）。


以下は都市数100の場合（正解の経路とは限らない）：

初期ランダム経路
Total Score:  52.8760351727298

欲張り法（二方向）による経路
Total Score:  9.979255125288242
まだ交差している経路があります。

二点間反転による最適化
Total Score:  8.764946875917067
交差していた経路がなくなって、ほぼ最適解。

一点移動挿入による最適化
Total Score:  8.463720661190624
細部の修正と最適化（わずかながらスコアは向上しています）。
ここまでは数秒で終了。

組み合わせ数８による最適化（8-opt）
Total Score:  8.435705379605803
さらなる最適化（右下の部分が改善）。
経路の組み合わせが8!=40320通りあるので、この工程だけ時間がかかります（1分43秒）。
よく見ると、赤点近くの98と99を83と82につなげたほうがスコアが向上しそうです。まだ赤点付近が上手く最適化できないので改善の余地あり。

対策：
局所的には最適化されスコアは向上しましたが、初期経路によってその後の最適化が左右されてしまうという弱点があります。初期経路生成の工夫をするか、後からでも大規模な経路変更が可能になる方法を探さないといけないかもしれません。最初に多少スコアオーバーしても許容する大規模な経路変換を行ったのちに局所的な最適化を行って、結果的にスコアが向上するような二段階の最適化をするといいのかもしれません。

続き：convex hull：凸包/ギフト包装法/Graham scanアルゴリズム
関連：Traveling Salesman Problem:巡回セールスマン問題について（まとめ）