pgRoutingのセットアップ
DockerかAmazon RDSを使うのが楽です。pgRoutingをはじめて使ってみたを参照。
グラフデータの準備
DIMACSというコンテストの問題データを使います。権利はパブリックドメイン。
http://www.dis.uniroma1.it/challenge9/download.shtml
から
http://www.dis.uniroma1.it/challenge9/data/USA-road-d/USA-road-d.NY.gr.gz
http://www.dis.uniroma1.it/challenge9/data/USA-road-d/USA-road-d.NY.co.gz
をダウンロードします。
grファイルは1行が1辺を表しており、2つの頂点IDとコストが書かれています。
coファイルは1行が1頂点を表しており、座標が書かれています。
v 1 -73530767 41085396
v 2 -73530767 41086098
v 3 -73519366 41048796
座標は1e-6倍すると経度、緯度になると思うのだけど（つまり-73530767 41085396が -73.530767 41.085396になる）、intのままでユークリッド平面座標として扱うことを想定してあるっぽい。実際、coファイルの方で頂点1と頂点2の距離は803になっています（Math.hypot(-73530767 - -73530767, 41085396 - 41086098) == 702で803とオーダーが合っている）。
これをRubyスクリプトでpgRouting用のSQLファイルに変換します：
#!/usr/bin/env ruby
COST_INFINITY = 99999999
table_name = "edges_astar"
co_filename = "USA-road-d.NY.co"
gr_filename = "USA-road-d.NY.gr"
def print_create_table(table_name)
puts <<~EOF
DROP TABLE IF EXISTS #{table_name};
CREATE TABLE #{table_name} (
id serial,
source int,
target int,
cost float,
reverse_cost float,
x1 float,
y1 float,
x2 float,
y2 float
);
EOF
end
def load_co_file(filename)
@vertices = {}
IO.foreach(filename, chomp: true) do |line|
a = line.split
if a[0] == "v"
@vertices[a[1]] = [a[2].to_i, a[3].to_i]
end
end
end
def load_gr_file(filename)
@edges = {}
IO.foreach(filename, chomp: true) do |line|
a = line.split
if a[0] == "a"
v1 = a[1].to_i
v2 = a[2].to_i
if v1 < v2
dir = :forward
key = "#{v1},#{v2}"
else
dir = :reverse
key = "#{v2},#{v1}"
end
@edges[key] ||= {}
@edges[key][dir] = a[3]
end
end
end
load_co_file(co_filename)
load_gr_file(gr_filename)
print_create_table(table_name)
puts <<EOF
INSERT INTO #{table_name} (
source,
target,
cost,
reverse_cost,
x1,
y1,
x2,
y2
) VALUES
EOF
i = 0
@edges.each do |key, value|
v1, v2 = key.split(",")
v1_co = @vertices[v1] || []
v2_co = @vertices[v2] || []
if i != 0
print ","
end
puts <<EOF
(#{v1}, #{v2}, #{value[:forward] || COST_INFINITY}, #{value[:reverse] || COST_INFINITY}, #{v1_co[0]}, #{v1_co[1]}, #{v2_co[0]}, #{v2_co[1]})
EOF
i += 1
end
puts <<EOF
;
EOF
$ ruby gr2sql_astar.rb > USA-road-d.NY.gr.sql
あとはpsqlでDBに接続して \i USA-road-d.NY.gr.sqlでインポート完了。
A*を試す
pgr_astarをやってみます。
# SELECT * FROM pgr_astar('SELECT id, source, target, cost, reverse_cost, x1, y1, x2, y2 FROM edges_astar', 91, 41, directed := true);
seq | path_seq | node | edge | cost | agg_cost
-----+----------+--------+--------+------+----------
1 | 1 | 91 | 56 | 456 | 0
2 | 2 | 90 | 363714 | 894 | 456
3 | 3 | 260469 | 358484 | 670 | 1350
4 | 4 | 260466 | 358483 | 1981 | 2020
5 | 5 | 5119 | 6030 | 1393 | 4001
6 | 6 | 5118 | 6031 | 201 | 5394
7 | 7 | 5116 | 6028 | 924 | 5595
8 | 8 | 5114 | 6026 | 1031 | 6519
9 | 9 | 5073 | 5977 | 3299 | 7550
10 | 10 | 5068 | 5972 | 3222 | 10849
11 | 11 | 5010 | 5973 | 805 | 14071
12 | 12 | 5011 | 5917 | 2458 | 14876
13 | 13 | 5005 | 5910 | 1966 | 17334
14 | 14 | 5003 | 5908 | 336 | 19300
15 | 15 | 5004 | 5911 | 473 | 19636
16 | 16 | 4997 | 5902 | 1027 | 20109
17 | 17 | 4995 | 5899 | 2879 | 21136
18 | 18 | 4985 | 5889 | 983 | 24015
19 | 19 | 4951 | 5888 | 2779 | 24998
20 | 20 | 4946 | 5882 | 889 | 27777
21 | 21 | 4947 | 5881 | 889 | 28666
22 | 22 | 4977 | 5878 | 980 | 29555
23 | 23 | 4978 | 5883 | 805 | 30535
24 | 24 | 4920 | 5810 | 1530 | 31340
25 | 25 | 4919 | 5809 | 1262 | 32870
26 | 26 | 4903 | 5794 | 1556 | 34132
27 | 27 | 4902 | 5805 | 7334 | 35688
28 | 28 | 4916 | 5804 | 1454 | 43022
29 | 29 | 4890 | 5778 | 5497 | 44476
30 | 30 | 4856 | 5741 | 940 | 49973
31 | 31 | 4857 | 7760 | 2843 | 50913
32 | 32 | 40 | 23 | 392 | 53756
33 | 33 | 41 | -1 | 0 | 54148
(33 rows)
頂点90から41までの経路が求まったようです。
pgr_astarではA*のヒューリスティック関数を数種類から選べます。
Heuristic number. Current valid values 0~5. Default 5
0: h(v) = 0 (Use this value to compare with pgr_dijkstra)
1: h(v) abs(max(dx, dy))
2: h(v) abs(min(dx, dy))
3: h(v) = dx * dx + dy * dy
4: h(v) = sqrt(dx * dx + dy * dy)
5: h(v) = abs(dx) + abs(dy)
座標が緯度経度でコストが距離（メートル）の場合は下記のようにfactor引数を指定するとうまい具合になるようです。
日本なら緯度が35くらいだからfactorは91288でしょうか。
Working with cost/reverse_cost as length in meters, x/y in lat/lon: Factor = would depend on the location of the points:
latitude conversion Factor
45 1 longitude degree is 78846.81 m 78846
0 1 longitude degree is 111319.46 m 111319
とりあえず今回は座標をintのまま入れているのでそこは関係ないとして、heuristicパラメータを変えて速度を見てみます。
| クエリ | 所要時間 |
|---|---|
| SELECT * FROM pgr_astar('SELECT id, source, target, cost, reverse_cost, x1, y1, x2, y2 FROM edges_astar', 91, 41, directed := true, heuristic := 0); | 861.512 ms |
| SELECT * FROM pgr_astar('SELECT id, source, target, cost, reverse_cost, x1, y1, x2, y2 FROM edges_astar', 91, 41, directed := true, heuristic := 1); | 1007.333 ms |
| SELECT * FROM pgr_astar('SELECT id, source, target, cost, reverse_cost, x1, y1, x2, y2 FROM edges_astar', 91, 41, directed := true, heuristic := 2); | 910.353 ms |
| SELECT * FROM pgr_astar('SELECT id, source, target, cost, reverse_cost, x1, y1, x2, y2 FROM edges_astar', 91, 41, directed := true, heuristic := 3); | 965.934 ms |
| SELECT * FROM pgr_astar('SELECT id, source, target, cost, reverse_cost, x1, y1, x2, y2 FROM edges_astar', 91, 41, directed := true, heuristic := 4); | 965.934 ms |
| SELECT * FROM pgr_astar('SELECT id, source, target, cost, reverse_cost, x1, y1, x2, y2 FROM edges_astar', 91, 41, directed := true, heuristic := 5); | 1001.475 ms |
あれ？heuristic := 0の場合はダイクストラ法と同じになるはずなのに、これが一番速いですね 🤔
なにか間違っているのでしょうか…誰か教えてください（教えてほしいから記事書いた）。
下記ページでもpgRoutingのAがダイクストラとほぼ同じか遅いと言われている。pgRoutingのAとダイクストラはBoost Graphを呼んでいるだけらしい。確かにそうだった。
https://gis.stackexchange.com/questions/84719/pgrouting-2-0-slow-execution-of-a-star-algorithm
https://gis.stackexchange.com/questions/69722/why-is-any-pgr-routing-function-taking-forever-based-on-osm-data-in-an-pgrouti/194328
boostを使ってダイクストラとAを比較するコードを書いて実験してみたところ、最短経路が短めの場合はダイクストラの方が早いが、長めの場合はAの方が早いという結果になった。