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頂点を表しており、座標が書かれています。

USA-road-d.NY.coの抜粋

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ファイルに変換します：

gr2sql_astar.rb

#!/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の方が早いという結果になった。