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