convex hull：凸包/ギフト包装法/Graham scanアルゴリズム

引き続きTSP（巡回セールスマン問題）についてというか、下準備の派生アルゴリズムについてです。
前回、初期経路生成に「欲張り法（貪欲法）」を用いてスタート地点から最も近い都市を次々とつないでいきましたが、その他の方法として「挿入法」というものもあるらしい。
「挿入法」では、最初に簡単な初期経路つくっておき、コストが最小になるような地点を徐々に加えていくようです。しかし、初期経路を何にするかという問題があり、通常は「凸包」で都市全体を外側から包み込む経路を初期経路とするといいようです（参照）。

今回はTSPの下準備となる、この凸包アルゴリズムについて試してみました。
凸包には「ギフト包装法」や「グラハムスキャン」などあるようです。この手のアルゴリズムは計算コストの工夫も同時に考慮しなければいけないようですが、今回はとりあえず解を求めることを目標に自分なりにコーディングしてみました。

追記：
グラハムスキャン、ならびに挿入法（経路生成）も追加しました。

これはランダムに配置した20個（0〜19）の点の場合。
全体を包み込むような経路（黄色）を生成します。
TSPに利用する場合は、ここからコストが最小になるような点を加えていき徐々に複雑な経路にしていきます。
赤点「0」はTSPにおけるスタート地点であり、今回はあまり意味はありません。


ギフト包装法アルゴリズム：
コード：

def giftwrap():
    lowest_y = np.where(Y==min(Y))[0][0]
    plist = [lowest_y]
    prev_theta = np.inf
    
    for _ in range(NUM): 
        best_theta = -np.inf
        best_idx = None
        for i in range(NUM):
            if i not in  plist:
                theta = np.arctan2(Y[i]-Y[plist[-1]], X[i]-X[plist[-1]])
                if theta >= best_theta and theta <= prev_theta:
                    best_idx = i
                    best_theta = theta
                    
        if best_idx is not None:
            theta_to_start = np.arctan2(Y[plist[0]]-Y[best_idx], X[plist[0]]-X[best_idx])
            if best_theta >= theta_to_start:
                plist.append(best_idx)
                prev_theta = best_theta
            else:
                break
        else:
            break

    plist.append(plist[0])
    
    return plist

・まずは点群の中でY座標値が最小となる点を選びます（複数ある場合はXも最小値）。
　上記の場合であれば「14」。これを開始点にしリストに加えておきます。
・「14」を起点に、その他のすべての点に向けての角度を比較します。
・角度については（時計で例えるなら）最小値→最大値は半時計回りになっています。
　9時：-180度
　6時：-90度
　3時：0度
　12時：90度
　9時：180度
・「14」を起点に角度が最大となる点を選びます。この場合「8」を選択。
　次回この(14, 8)の辺の角度を利用するので記録しておきます。
・「8」をリストに加えておき、次の起点を「8」にします。
・先程と同様に起点「8」からその他の点に向けての角度を比較します。
　この際、リストに加えた点以外という条件設定にしておきます。
・基本的にはこれを繰り返していくだけです。

ただし、もう少し条件を加えないと循環経路として閉じてくれません。
リストに加えていく候補点の条件としては：
（１）リスト以外の点
（２）起点からの角度が最大の点
（３）前回の辺の角度以下の角度となる点
（４）候補点から開始点（「14」）に向けた角度以上の角度となる点
という感じでしょうか。
最後に開始点をリストに加えて終了します。
各辺の角度は囲い込んでいくにつれ減少していく特徴があるので、その部分を条件式として利用しています。


今回は自前でコーディングしたので、正当な「ギフト包装法」のコードとは違うかもしれません。しかしながらこのようなアルゴリズムはTSPの初期経路生成に利用するだけでなくいろんな分野でも活用できそうです。



グラハムスキャン（Graham Scan）：（追記）
もうひとつの凸包アルゴリズムであるグラハムスキャンも試してみました。計算量はグラハムスキャンの方が少ないらしいのですが、実際どうなのか？
アルゴリズムは、各点をスキャンして行く際に行ったり来たりする挙動なので最初わかりにくかったのですが、処理の順番を理解すればそれほど難しいものでもありませんでした。

コード：

def graham_scan():
    lowest_y = np.where(Y==min(Y))[0][0]
    plist = [lowest_y]
    
    theta_list = [] 
    for i in range(NUM):
        theta = np.arctan2(Y[i]-Y[plist[-1]], X[i]-X[plist[-1]])
        theta_list.append(theta)
        
    sorted_index = np.argsort(theta_list)
    
    for i in range(1,NUM):
        plist.append(sorted_index[i])
        if i > 1:
            v1 = [X[plist[-2]] - X[plist[-3]], Y[plist[-2]] - Y[plist[-3]]]
            v2 = [X[plist[-1]] - X[plist[-2]], Y[plist[-1]] - Y[plist[-2]]]
            while np.cross(v1, v2) <= 0:
                plist.pop(-2)
                v1 = [X[plist[-2]] - X[plist[-3]], Y[plist[-2]] - Y[plist[-3]]]
                v2 = [X[plist[-1]] - X[plist[-2]], Y[plist[-1]] - Y[plist[-2]]]
        
    plist.append(plist[0])
    return plist

・Y座標が最小の点を選択し起点とします。
・最終的な経路リストを用意して、起点を入れておきます。
・起点から全ての点に対する角度のリストをつくります。
・角度のリストをソートしなおして角度順のリストをつくります。
・角度順リストを元に角度の小さい方から各点を経路リストに加えていきます。
・その際に現在の辺と一つ前の辺の外積を求め、以下をチェックします。
　　外積=0：一つ手前の辺と現在の辺が平行（延長線上にある）
　　外積<0：一つ手前の辺からみて現在の辺が右曲がり
　　外積>0：一つ手前の辺からみて現在の辺が左曲がり
・外積<0の際は一つ手前の点を経路リストから削除します（外積>=0になるまで繰り返す）。
・次のノードへ移動し、上記を繰り返す。
・最後に経路を閉じるため経路リストの最初の値を加えて終了。

起点を中心とすると、起点の左側から反時計周りに角度が小→大となります。
循環経路なので追加される各辺は常に左曲がりとなります。つまり徐々に角度の値が大きくなっていくということです。
左曲/右曲の判定は外積（現在の辺と一つ手前の辺による）を使います。np.cross(v1, v2)を利用しましたが、x1*y2-x2*y1でも求められます。
もし途中で現在の辺が右曲がりなら、その手前の点を経路リストから削除します。そうすると、二つ手前の辺が繰り上がって一つ手前の辺になります。再度、現在の辺と一つ手前の辺（繰り上がった二つ手前の辺）の外積を調べて左曲/右曲の判定をします。現在の辺が右曲であるかぎり、この外積の左曲/右曲判定と点の削除をwhile文で繰り返します。現在の辺が左曲がりになれば次の点へ移動していきます。

以下の動画が分かりやすかったです。

ということで、ギフト包装法とグラハムスキャンの二つのコード（以下）。
100000点で比較してみると、
ギフト包装法：6.80秒
グラハムスキャン：4.76秒
という結果になりました。


挿入法によるTSP初期経路生成（追記）：
TSPに適用するには、凸包上のノードAとBからなる任意の辺、加えようとするノードをCとし、
L=AC+BC-AB
という長さを各辺において比較し最小値のノードを加えていきます。一点加えるごとに形が変化していくので、辺と点のすべての組み合わせの二重forループで一点追加していきます。その分計算量は多くなりますが、貪欲法でありながらも無駄な経路は少なくなります。
要は外側から徐々に内部の方に侵食していくように経路を形成していきます。

def hull_insert(path, hull):
    Path = list(path[:-1].copy())
    Hull = list(hull[:-1].copy())
    for h in Hull:
        Path.remove(h)
        
    def tri_length(p0, p1, pt): 
        d0t = abs(XY[p0]-XY[pt])
        d1t = abs(XY[p1]-XY[pt])
        d01 = abs(XY[p0]-XY[p1])
        L = d0t + d1t - d01
        return L
       
    while len(Path):
        best_len = np.inf
        best_node = None
        for i in range(len(Hull)):
            for p in Path:
                if i == len(Hull)-1:
                    L = tri_length(Hull[i], Hull[0], p)
                else:
                    L = tri_length(Hull[i], Hull[i+1], p)
                if L < best_len:
                    best_len = L
                    best_node = [Hull[i], p]
                        
        Hull.insert(Hull.index(best_node[0])+1, best_node[1])
        Path.remove(best_node[1])

    zero = Hull.index(0)
    tour = Hull[zero:] + Hull[:zero] + [0]
    return tour

最終的には2-OPTなどによる最適化も必要ですが、単純な貪欲法よりは整った経路ができあがります。しかしながら、その後の最適化によって改善しやすいかどうかは分かりません。

ノード20個の場合。

ノード100個の場合。

交差している箇所も少しありますが、最終型に近い経路になります。

追記：

以下は経路生成（ノード：50）のアニメーション。下のGistの最後にコードを追加しました。

外側から内側へ向かって一つずつ最短距離の点を加えていくので一見無駄がなさそうですが、まだ選択されていない点によって形成される最短経路を無視して繋いでいくため、必ずしも最終的な最短経路になるとは限りません。

関連：Traveling Salesman Problem:巡回セールスマン問題について（まとめ）

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:巡回セールスマン問題について（まとめ）