巡回セールスマン問題実践集③~非ランダムマップ編~

巡回セールスマン問題実践集①~全探索編~
巡回セールスマン問題実践集②~最近傍探索法編~

前回は近傍探索法により計算量を抑えつつTSPの近似解を求めた。実際の地図上で経路を求めるようなことをした方が楽しいだろう。

#非ランダムマップ
都道府県庁所在地について緯度経度から地図におこす。
##Coordinate_map
csv.readerを用いつつ、緯度経度をピックアップしてCityを生成していく関数Coordinate_mapを定義する。

def lines(text): return text.strip().splitlines()


def Coordinate_map(lines, delimiter=',', lat_col=0, long_col=1, lat_scale=91, long_scale=80):
    """文章から都市のセットを作る。
    経度・緯度をlong_scale,lat_scaleからx・y直交座標として扱う。
    ソースはファイルでも文章のリストでも良い。"""
    return frozenset(City(long_scale * float(row[long_col]), 
                          lat_scale  * float(row[lat_col]))
                     for row in csv.reader(lines, delimiter=delimiter, skipinitialspace=True))

lat_col=0やlong_col=1というのは与えたリストのどの部分を緯度、経度として取り込むかを選ぶパラメーターである。また、lat_scaleやlong_scaleは、緯度、経度から擬似的にxy直交座標系に変換する際の離れる度合いのパラメータである。例えば赤道付近での経度の1度の差と極付近における経度の1度の差は異なる。また、北緯は増えると北に移動するが、南緯は増えると南に移動する。西経の時はマイナスとかつけて微調整してね。

##都道府県庁所在地マップ
前述の通り都道府県庁所在地の情報からmapを設定する。JPN_mapを設定します。

JPN_map

JPN_map = Coordinate_map(lines("""
北海道,札幌市,43.06417,141.34694
青森県,青森市,40.82444,140.74
岩手県,盛岡市,39.70361,141.1525
宮城県,仙台市,38.26889,140.87194
秋田県,秋田市,39.71861,140.1025
山形県,山形市,38.24056,140.36333
福島県,福島市,37.75,140.46778
茨城県,水戸市,36.34139,140.44667
栃木県,宇都宮市,36.56583,139.88361
群馬県,前橋市,36.39111,139.06083
埼玉県,さいたま市,35.85694,139.64889
千葉県,千葉市,35.60472,140.12333
東京都,新宿区,35.68944,139.69167
神奈川県,横浜市,35.44778,139.6425
新潟県,新潟市,37.90222,139.02361
富山県,富山市,36.69528,137.21139
石川県,金沢市,36.59444,136.62556
福井県,福井市,36.06528,136.22194
山梨県,甲府市,35.66389,138.56833
長野県,長野市,36.65139,138.18111
岐阜県,岐阜市,35.39111,136.72222
静岡県,静岡市,34.97694,138.38306
愛知県,名古屋市,35.18028,136.90667
三重県,津市,34.73028,136.50861
滋賀県,大津市,35.00444,135.86833
京都府,京都市,35.02139,135.75556
大阪府,大阪市,34.68639,135.52
兵庫県,神戸市,34.69139,135.18306
奈良県,奈良市,34.68528,135.83278
和歌山県,和歌山市,34.22611,135.1675
鳥取県,鳥取市,35.50361,134.23833
島根県,松江市,35.47222,133.05056
岡山県,岡山市,34.66167,133.935
広島県,広島市,34.39639,132.45944
山口県,山口市,34.18583,131.47139
徳島県,徳島市,34.06583,134.55944
香川県,高松市,34.34028,134.04333
愛媛県,松山市,33.84167,132.76611
高知県,高知市,33.55972,133.53111
福岡県,福岡市,33.60639,130.41806
佐賀県,佐賀市,33.24944,130.29889
長崎県,長崎市,32.74472,129.87361
熊本県,熊本市,32.78972,130.74167
大分県,大分市,33.23806,131.6125
宮崎県,宮崎市,31.91111,131.42389
鹿児島県,鹿児島市,31.56028,130.55806
沖縄県,那覇市,26.2125,127.68111
"""), delimiter=',', lat_col=2, long_col=3,lat_scale=91, long_scale=80)

どうなっているか見てみましょう

plot_lines(JPN_map, 'bo')

いい感じですね。解いてみると

plot_tsp(repeated_altered_nn_tsp, JPN_map)

47 city tour with length 4830.0 in 0.125 secs for repeated_altered_nn_tsp

まぁいい感じ。ただ、日本が細長いのでパッとしませんね。

##福岡県の市区町村の役所・役場
福岡県の市区町村の役所・役場から福岡県の市区町村の役所・役場を全て巡る最短経路を探りたい。

FO

FO = lines("""
福岡県庁 130.25.05 33.36.23 
北九州市役所 130.52.31 33.53.00 
門司区役所 130.57.35 33.56.28 
若松区役所 130.48.40 33.54.20 
戸畑区役所 130.49.47 33.53.36 
小倉北区役所 130.52.25 33.52.51 
小倉南区役所 130.53.05 33.50.47 
八幡東区役所 130.48.43 33.51.49 
八幡西区役所 130.45.51 33.51.59 
福岡市役所 130.24.06 33.35.24 
東区役所 130.25.03 33.37.04 
博多区役所 130.24.54 33.35.29 
中央区役所 130.23.35 33.35.21 
南区役所 130.25.36 33.33.42 
西区役所 130.19.23 33.34.58 
城南区役所 130.22.12 33.34.33 
早良区役所 130.20.54 33.34.55 
大牟田市役所 130.26.46 33.01.49 
久留米市役所 130.30.30 33.19.10 
直方市役所 130.43.47 33.44.38 
飯塚市役所 130.41.28 33.38.47 
田川市役所 130.48.22 33.38.20 
柳川市役所 130.24.22 33.09.47 
八女市役所 130.33.28 33.12.43 
筑後市役所 130.30.08 33.12.44 
大川市役所 130.23.02 33.12.24 
行橋市役所 130.58.59 33.43.43 
豊前市役所 131.07.49 33.36.42 
中間市役所 130.42.33 33.49.00 
小郡市役所 130.33.20 33.23.47 
筑紫野市役所 130.31.34 33.29.15 
春日市役所 130.28.13 33.31.58 
大野城市役所 130.28.44 33.32.11 
宗像市役所 130.32.26 33.48.20 
太宰府市役所 130.31.26 33.30.46 
古賀市役所 130.28.12 33.43.44 
福津市役所 130.29.28 33.46.01 
うきは市役所 130.45.18 33.20.50 
宮若市役所 130.40.00 33.43.25 
嘉麻市役所 130.42.42 33.33.48 
朝倉市役所 130.39.56 33.25.24 
みやま市役所 130.28.29 33.09.09 
糸島市役所 130.11.44 33.33.26 
那珂川市役所 130.25.20 33.29.58 
宇美町役場 130.30.40 33.34.04 
篠栗町役場 130.31.35 33.37.26 
志免町役場 130.28.47 33.35.29 
須恵町役場 130.30.26 33.35.14 
新宮町役場 130.26.48 33.42.55 
久山町役場 130.30.00 33.38.48 
粕屋町役場 130.28.50 33.36.39 
芦屋町役場 130.39.50 33.53.38 
水巻町役場 130.41.41 33.51.17 
岡垣町役場 130.36.41 33.51.13 
遠賀町役場 130.40.06 33.50.53 
小竹町役場 130.42.46 33.41.33 
鞍手町役場 130.40.27 33.47.31 
桂川町役場 130.40.41 33.34.44 
筑前町役場 130.35.43 33.27.25 
東峰村役場 130.52.12 33.23.50 
大刀洗町役場 130.37.21 33.22.21 
大木町役場 130.26.23 33.12.38 
広川町役場 130.33.05 33.14.29 
香春町役場 130.50.50 33.40.05 
添田町役場 130.51.15 33.34.18 
糸田町役場 130.46.45 33.39.10
川崎町役場 130.48.55 33.36.00 
大任町役場 130.51.13 33.36.44 
赤村役場 130.52.15 33.37.00 
福智町役場 130.46.48 33.41.00 
苅田町役場 130.58.50 33.46.34 
みやこ町役場 130.55.14 33.41.57 
吉富町役場 131.10.34 33.36.10 
上毛町役場 131.09.52 33.34.42 
築上町役場 131.03.22 33.39.22 
""")

しかし、これは60進数で表されている。

##Coordinate_map2
まずは60進数を10進数になおす関数を準備。

def sixty2ten(x):
    x_list = x.split('.')
    t = float(x_list[0]) + float(x_list[1]) / 60 + float(x_list[2]) / 60 / 60
    return t

Coordinate_map2を定義。

def Coordinate_map2(lines, delimiter=' ', lat_col=2, long_col=1, lat_scale=91000, long_scale=81000):
    return frozenset(City(long_scale * sixty2ten(row[long_col]), 
                          lat_scale  * sixty2ten(row[lat_col]))
                     for row in csv.reader(lines, delimiter=delimiter, skipinitialspace=True))

FO_map = Coordinate_map2(FO)

plot_tsp(repeated_altered_nn_tsp,FO_map)

75 city tour with length 425110.1 in 0.262 secs for repeated_altered_nn_tsp

うん、いい感じ

##どの都市を巡ればいいのか
だが、しかし実際にどういう順番で巡れば良いのかわからない。
これまでは座標について並び替えていたが、座標の地点についての番号を振りどの番号を巡るのかを決めることで、巡る都市名を呼び出せるようにする。

そのためにはこれまで設定した関数を定義しなおす必要がある。
元のデータに都市名と都市の緯度軽度が連続して入力されているとする。
mat_lat_longに都市の緯度軽度の情報を抜き出し、
###distance2

def distance2(A, B): 
    global map_lat_long
    "2点間の距離"
    return abs(map_lat_long[A] - map_lat_long[B])

###tour_length2

def tour_length2(tour):
    "経路に沿った2点間の距離を足し合わせる。"
    return sum(distance2(tour[i], tour[i-1]) 
               for i in range(len(tour)))

今回はaltered_nn_tspが使えるところまで定義する。

###nn_tsp2

def nn_tsp2(cities):
    """ある都市から経路を始め、その都市から一番近い都市へと経路を進め、
　　　　その次の都市からさらに次の、まだ訪れていない都市へと経路を進める。"""
    start = 0
    tour = [start]
    unvisited = set(range(1,len(cities)))
    while unvisited:
        C = nearest_neighbor2(tour[-1], unvisited)
        tour.append(C)
        unvisited.remove(C)
    return tour

def nearest_neighbor2(A, cities):
    "citiesのうち、Aに一番近いものを見つける。"
    return min(cities, key=lambda c: distance2(c, A))

###altered_nn_tsp2

def reverse_segment_if_better2(tour, i, j):
    "もし逆転させたセグメント[i:j]によって経路が短くなるなら、逆転させる。"
    #元の経路[...A-B...C-D...]に対してB...Cを逆転させた[...A-C...B-D...]を考える。
    A, B, C, D = tour[i-1], tour[i], tour[j-1], tour[j % len(tour)]

    if distance2(A, B) + distance2(C, D) > distance2(A, C) + distance2(B, D):
        tour[i:j] = reversed(tour[i:j])

def alter_tour2(tour):
    "セグメントを逆転させることで経路をより短くしようと試す。"
    original_length = tour_length2(tour)
    for (start, end) in all_segments2(len(tour)):
        reverse_segment_if_better2(tour, start, end)
    # If we made an improvement, then try again; else stop and return tour.
    #改善があったならもう一度繰り返し、なければ終了しtourを出力する。
    if tour_length2(tour) < original_length:
        return alter_tour2(tour)
    return tour

def all_segments2(N):
    "長さNの経路についての全セグメントを表す(始点、終点)の組を出力"
    return [(start, start + length)
            for length in range(N, 2-1, -1)
            for start in range(N - length + 1)]

def altered_nn_tsp2(cities):
    "Run nearest neighbor TSP algorithm, and alter the results by reversing segments."
    return alter_tour2(nn_tsp2(cities))

###描画

#描画
def plot_lines2(points, style='bo-'):
    "Plot lines to connect a series of points."
    global map_lat_long
    plt.plot([map_lat_long[i].x for i in points], [map_lat_long[i].y for i in points], style)
    plt.axis('scaled'); plt.axis('off')

def plot_tour2(tour): 
    "Plot the cities as circles and the tour as lines between them."
    plot_lines2(list(tour) + [tour[0]])

###実行

JPN_map2

JPN_map2 = lines("""
北海道,札幌市,43.06417,141.34694
青森県,青森市,40.82444,140.74
岩手県,盛岡市,39.70361,141.1525
宮城県,仙台市,38.26889,140.87194
秋田県,秋田市,39.71861,140.1025
山形県,山形市,38.24056,140.36333
福島県,福島市,37.75,140.46778
茨城県,水戸市,36.34139,140.44667
栃木県,宇都宮市,36.56583,139.88361
群馬県,前橋市,36.39111,139.06083
埼玉県,さいたま市,35.85694,139.64889
千葉県,千葉市,35.60472,140.12333
東京都,新宿区,35.68944,139.69167
神奈川県,横浜市,35.44778,139.6425
新潟県,新潟市,37.90222,139.02361
富山県,富山市,36.69528,137.21139
石川県,金沢市,36.59444,136.62556
福井県,福井市,36.06528,136.22194
山梨県,甲府市,35.66389,138.56833
長野県,長野市,36.65139,138.18111
岐阜県,岐阜市,35.39111,136.72222
静岡県,静岡市,34.97694,138.38306
愛知県,名古屋市,35.18028,136.90667
三重県,津市,34.73028,136.50861
滋賀県,大津市,35.00444,135.86833
京都府,京都市,35.02139,135.75556
大阪府,大阪市,34.68639,135.52
兵庫県,神戸市,34.69139,135.18306
奈良県,奈良市,34.68528,135.83278
和歌山県,和歌山市,34.22611,135.1675
鳥取県,鳥取市,35.50361,134.23833
島根県,松江市,35.47222,133.05056
岡山県,岡山市,34.66167,133.935
広島県,広島市,34.39639,132.45944
山口県,山口市,34.18583,131.47139
徳島県,徳島市,34.06583,134.55944
香川県,高松市,34.34028,134.04333
愛媛県,松山市,33.84167,132.76611
高知県,高知市,33.55972,133.53111
福岡県,福岡市,33.60639,130.41806
佐賀県,佐賀市,33.24944,130.29889
長崎県,長崎市,32.74472,129.87361
熊本県,熊本市,32.78972,130.74167
大分県,大分市,33.23806,131.6125
宮崎県,宮崎市,31.91111,131.42389
鹿児島県,鹿児島市,31.56028,130.55806
沖縄県,那覇市,26.2125,127.68111
""")

map_lat_long = Coordinate_map2(JPN_map2, delimiter=',', lat_col=2, long_col=3,lat_scale=91, long_scale=80)
map = list((row[1])
    for row in csv.reader(JPN_map2, delimiter=','))
ant = altered_nn_tsp2(map_lat_long)
for i in range(len(map_lat_long)):
    print(map[ant[i]])
plot_tour2(ant)

###結果
以下のような都道府県庁所在地のリストが得られる。

都道府県庁所在地巡回順リスト

金沢市 富山市 長野市 新潟市 秋田市 青森市 札幌市 盛岡市 仙台市 山形市 福島市 宇都宮市 水戸市 千葉市 横浜市 新宿区 さいたま市 前橋市 甲府市 静岡市 名古屋市 岐阜市 津市 奈良市 大津市 京都市 大阪市 神戸市 和歌山市 徳島市 高知市 宮崎市 鹿児島市 那覇市 長崎市 佐賀市 福岡市 熊本市 大分市 山口市 広島市 松山市 高松市 岡山市 松江市 鳥取市 福井市

また今回の順番で巡った時の経路が 

#まとめ
非ランダム地図における巡回セールスマン問題を解けた。
都市の座標と都市名の紐付けをした上で結果を出力する方法は改善の余地があるだろう。
また、現実において経路を単純な座標的距離だ実際の遠さを表すわけではない。2点間の距離の表現をどう反映するのか(時間距離やコストの大小としての距離など)が課題である。