Mona Lisa TSPでモナリザの顔を表示させる -その１-

python TSP

少し前、研究室の先輩からこんな問題を教えてもらった。

http://www.tsp.gatech.edu/data/ml/monalisa.html

巡回セールスマン問題(TSP)の最適解を解くと、モナリザの顔が浮かび上がってくるらしい。段々モナリザの顔が…みたいな事をやりたくなったので、最適解は無理でもモナリザの顔が見えるくらいは頑張ってみることに。
とりあえず、簡単な最適化アルゴリズムで様子をみることにした。とりあえず簡単そうなアルゴリズムとして、二つ選んだ。最初に訪問順序をランダムに初期化して、経路順序{x_1,x_2,...,x_n}（nは頂点数）が得られたとすると、

    1. ランダムに選択した2頂点x_i、x_jの訪問順序を入れ替え、回路長が短くなったら入れ替えを保持。これを連続max_fail回失敗するまで続ける。
    2. ランダムに選択した2頂点x_iからx_jまでの経路をひっくり返す。回路長が短くなったら入れ替えを保持。これを連続max_fail回失敗するまで続ける。

実行速度が心配だったけどpythonで書いた。コードはこんな感じ

import random, math

def solve(max_fail, max_iteration):
    data = read_from("mona-lisa100K.tsp")
    min_answer = data
    for i in xrange(max_iteration):
        local_answer = opt_swap(list(data), max_fail)
        if dist_all(min_answer) > dist_all(local_answer):
            min_answer = local_answer
    print dist_all(min_answer)
    write_to("answer.txt", min_answer)

def choose_index_pair(size):
    index1 = random.randint(0, size - 1)
    index2 = random.randint(0, size - 1)
    while index1 == index2:
        index1 = random.randint(0, size - 1)
        index2 = random.randint(0, size - 1)
    return (index1, index2)

def dist_all(points):
    return sum(dist(points[i], points[i - 1]) for i in xrange(len(points) - 1, -1, -1))

def dist(p1, p2):
    return math.sqrt((p1[0] - p2[0]) ** 2 + (p1[1] - p2[1]) ** 2)

def opt_inverse(data, max_fail):
    random.shuffle(data)
    iteration = 0
    continuous_fail = 0
    while continuous_fail < max_fail:
        iteration += 1
        index1, index2 = choose_index_pair(len(data))
        if can_improve_inv(data, index1, index2):
            max_i, min_i = max(index1, index2), min(index1, index2)
            data[min_i:max_i] = data[max_i - 1:min_i - 1:-1] if min_i > 0 else data[max_i - 1::-1]
            continuous_fail = 0
        else:
            continuous_fail += 1
        if iteration % 50000 == 0:
            print iteration
    return data

def can_improve_inv(data, index1, index2):
    l = len(data)
    nIndex1 = index1 + 1 if index1 < l - 1 else 0
    nIndex2 = index2 + 1 if index2 < l - 1 else 0
    before_dist = dist(data[index1 - 1], data[index1]) + dist(data[index2], data[nIndex2])
    after_dist  = dist(data[index1 - 1], data[index2]) + dist(data[index1], data[nIndex2])
    return before_dist > after_dist

def opt_swap(data, max_fail):
    random.shuffle(data)
    iteration = 0
    continuous_fail = 0
    while continuous_fail < max_fail:
        iteration += 1
        index1, index2 = choose_index_pair(len(data))
        if can_improve_opt_swap(data, index1, index2):
            data[index1], data[index2] = data[index2], data[index1]
            continuous_fail = 0
        else:
            continuous_fail += 1
        if iteration % 5000000 == 0:
            print iteration
    return data

def can_improve_opt_swap(data, index1, index2):
    l = len(data)
    nIndex1 = index1 + 1 if index1 < l - 1 else 0
    nIndex2 = index2 + 1 if index2 < l - 1 else 0
    before_dist = dist(data[index1 - 1], data[index1]) + dist(data[index1], data[nIndex1]) + \
        dist(data[index2 - 1], data[index2]) + dist(data[index2], data[nIndex2])
    after_dist = dist(data[index1 - 1], data[index2]) + dist(data[index2], data[nIndex1]) + \
        dist(data[index2 - 1], data[index1]) + dist(data[index1], data[nIndex2])
    return before_dist > after_dist

def write_to(filename, data):
    fout = open(filename, 'w')
    fout.writelines([" ".join(map(str, item)) + "\n" for item in data])
    fout.close()    

def read_from(filename):
    lst = []
    for line in open(filename).readlines():
        coordinate = line.split()
        lst.append((int(coordinate[1]), int(coordinate[2])))
    return lst

if __name__ == "__main__":
    max_fail = 20000
    max_iteration = 1
    solve(max_fail, max_iteration)

標準の機能のみだったから、素のpythonじゃなくてずっと高速なpypyで実行。どうも数分程度で同一の実行時間ならswapの方が優秀らしい。とりあえず、最適化後の経路長はおおよそ1.1〜1.2億程度になる。最短経路はどんなもんかなーと見てみると何と570万...

嫌な予感しかしないが一応プロットしてみると、

何がなんだかwちなみに点だけを打った本物はこんな感じ

さすがに乱択のみでは厳しかったか。次は貪欲法っぽく探索してない点を結んで行くとか
を実装する予定。素で計算すると計算量がちょっとヤバそうなので、KD-Treeあたりを実装する予定。