グラフ問題：ダイクストラ法

# はじめに

研究の過程でメッシュクラスを自作する必要に迫られたため、BGL(Boost Graph Library)やCGAL(Computational Geometry Algorithms Library)のドキュメントやソースコードを読んでメッシュクラスの実装を勉強しました。その過程でグラフ理論に興味を持ったので、自分でダイクストラ法をC++で書いてみました。

# 実装

ダイクストラ法の実装は、[Wikipediaのダイクストラ法](https://ja.wikipedia.org/wiki/%E3%83%80%E3%82%A4%E3%82%AF%E3%82%B9%E3%83%88%E3%83%A9%E6%B3%95)の擬似コードをそのまま使いました。
単体テストに[[アルゴリズム] ダイクストラ法をやってみる](https://qiita.com/edo_m18/items/0588d290a19f2afc0a84)の記事の問題を使わせていただき、実装が正しいことを確認しています。

ここでの実装方法は、[[アルゴリズム] ダイクストラ法をやってみる](https://qiita.com/edo_m18/items/0588d290a19f2afc0a84)の記事や[その実装の元となるブログ](http://www.deqnotes.net/acmicpc/dijkstra/)とは大きく異なっています。
これらの記事では、頂点の接続やスタートからの最短経路の距離などの情報を全て保持する`node`クラス

```cpp
struct node {
  std::vector<int> edges_to;       // この頂点に接続している辺番号リスト
  std::vector<double> edges_dist;  // 接続されている辺の長さ
  bool done;                       // 確定ノードか否か
  double dist;                     // スタートからの最短経路の距離
};
```

を作成しています。
私の経験上、このようなオブジェクト指向の実装の仕方だと

1. スケーラブルでない
2. 変更に対して弱い
3. パフォーマンスが低い

ので、ここではデータ指向かつBGL的な実装をしています。

注：以下では読みやすさを優先して、境界チェックを行う`vector::at()`や`gsl::at()`を使っていません。

## インデックス

頂点と辺の番号を区別するため、`index_base`というクラスを定義しました。

```cpp:index_base
template <typename Tag>
class index_base {
 public:
  index_base() = default;
  explicit index_base(std::size_t index) : index_{index} {}
  operator std::size_t() const noexcept { return index_; }
  bool is_valid() const noexcept {
    return index_ != std::numeric_limits<std::size_t>::max();
  }

 private:
  std::size_t index_ = std::numeric_limits<std::size_t>::max();
};
```

インデックスを使ってリストに簡単にアクセスできるように、`std::size_t()`演算子を定義し、暗黙的に`std::size_t`型へ変換可能にしています。
頂点および辺のタグ（`vertex_tag`/`edge_tag`）を作り、頂点番号と辺番号（`vertex_index`/`edge_index`）をエイリアスとして定義します。

```cpp
struct vertex_tag {};
struct edge_tag {};

using vertex_index = index_base<vertex_tag>;
using edge_index   = index_base<edge_tag>;
```

## 頂点

今回は頂点＝頂点番号ということにして、頂点は定義しませんでした。

## 辺

（無向）辺は２つの頂点のペアとして定義しました。

```cpp
struct edge {
  vertex_index first;
  vertex_index second;

  bool has_vertex(vertex_index v) const noexcept {
    return this->first == v || this->second == v;
  }
};
```

## グラフ

グラフ理論では、グラフ$G$は頂点の集合$V$および辺の集合$E$からなります。

$$ G = (V, E) $$

そこで、頂点のリスト`std::vector<vertex_index>`と（無向）辺のリスト`std::vector<edge>`をデータとして持つクラス`graph`を作ります。さらに、頂点間の接続を表す`adjacency_list`と、頂点につながっている辺を表す`vertex_edge_list`も加えます。`adjacency_list`はダイクストラ法のアルゴリズムに必要なリストで、`vertex_edge_list`は２つの頂点$u, v \in V$を与えられたときに対応する辺番号`edge_index`を探す関数`find_index`で使われるリストです。

```cpp
using vertex_list      = std::vector<vertex_index>;
using edge_list        = std::vector<edge>;
using adjacency_list   = std::vector<std::vector<vertex_index>>;
using vertex_edge_list = std::vector<std::vector<edge_index>>;

struct graph {
  vertex_list vertices;
  edge_list edges;
  adjacency_list adjacency;
  vertex_edge_list vertex_edges;
};

adjacency_list make_adjacency_list(const vertex_list& vs, const edge_list& es);
vertex_edge_list make_vertex_edge_list(const vertex_list& vs, const edge_list& es);

inline graph make_graph(const vertex_list& vs, const edge_list& es) {
  return {vs, es, make_adjacency_list(vs, es), make_vertex_edge_list(vs, es)};
}

edge_index find_index(const graph& g, const edge& e);
```

ここでは簡単のため、`adjacency_list`と`vertex_edge_list`を`std::vector`を使った二次元リストとして定義していますが、`boost::container::small_vector`を使うとキャッシュミスが減り最適化されます。

```cpp
using adjacency_list   = std::vector<boost::container::small_vector<vertex_index, 4>>;
using vertex_edge_list = std::vector<boost::container::small_vector<edge_index, 4>>;
```

`make_adjacency_list`と`make_vertex_edge_list`の実装は以下のとおりです。

```cpp
adjacency_list make_adjacency_list(const vertex_list& vs, const edge_list& es) {
  adjacency_list adj(vs.size());
  for (auto&& e : es) {
    adj[e.first].push_back(e.second);
    adj[e.second].push_back(e.first);
  }
  return adj;
}
```

```cpp
vertex_edge_list make_vertex_edge_list(const vertex_list& vs, const edge_list& es) {
  vertex_edge_list ves(vs.size());
  for (std::size_t i = 0; i < es.size(); ++i) {
    const auto ei = edge_index{i};
    const auto& e = es[ei];
    ves[e.first].push_back(ei);
    ves[e.second].push_back(ei);
  }
  return ves;
}
```

`find_index`の実装は次のようになります。

```cpp
edge_index find_index(const graph& g, const edge& e) {
  const auto& eis = g.vertex_edges[e.first];
  for (auto&& ei : eis)
    if (g.edges[ei].has_vertex(e.second) return ei;
  return edge_index{};
}
```

## ダイクストラ法

以上でダイクストラ法を実装する準備ができました。グラフ、始点、終点、辺の長さを与えた時に最短経路の長さを返す関数を定義します。

```cpp
using edge_length_list = std::vector<double>;

double min_distance(const graph& g,
                    const vertex_index start,
                    const vertex_index end,
                    const edge_length_list& lengths);
```

ダイクストラ法の実装は下記の通りになります。ほぼ[Wikipediaのダイクストラ法の擬似コード](https://ja.wikipedia.org/wiki/%E3%83%80%E3%82%A4%E3%82%AF%E3%82%B9%E3%83%88%E3%83%A9%E6%B3%95)通りになっているため、理解しやすいと思います。
そのうち`std::priority_queue`を使った実装もしてみようと思います。


```cpp
double min_distance(const graph& g,
                    const vertex_index start,
                    const vertex_index end,
                    const edge_length_list& lengths) {
  // 初期化
  const auto& vs = g.vertices;
  vertex_list q = vs;

  std::vector<double> dist(vs.size());
  dist.resize(vs.size(), std::numeric_limits<double>::max());
  dist[start] = 0.0;

  std::vector<vertex_index> prev(vs.size());
  prev.resize(vs.size());

  const auto& adj = g.adjacency;

  // 本計算
  while (!q.empty()) {
    // qからdist[u]が最小である頂点uを探し取り除く
    auto it = std::min_element(q.begin(), q.end(), [&dist](auto u, auto v) {
      return dist[u] < dist[v];
    });
    const auto u = *it;
    q.erase(it);

    // uに接続している各頂点に対しスタートからの経路の長さを計算し、最短経路ならば情報を更新
    const auto& vis_u = adj[u];
    for (auto&& v : vis_u) {
      const auto ei = find_index(g, {u, v});
      const auto d = dist[u] + lengths[ei];
      if (dist[v] > d) {
        dist[v] = d;
        prev[v] = u;
      }
    }
  }
  return dist[end];
}
```

## `property_map`

今回は簡単のため`std::vector`を使いましたが、各リストが頂点または辺のどちらに紐付けられたデータなのかわかりにくいという問題があります。__グラフの頂点または辺に紐付けられたデータである__ということを型を使って表現できると便利ですよね。そこでBoostライブラリのPropertyMapに相当するものを作ることにします。

- [Boost.PropertyMap](https://www.boost.org/doc/libs/1_49_0/libs/property_map/doc/property_map.html)
- [Boost.勉強会 #3 関西「Boost.PropertyMap」](https://www.slideshare.net/Cryolite/property-map)

```cpp
template <typename T, typename Key, typename Allocator = std::allocator<T>>
struct property_map : std::vector<T, Allocator> {
  using base_type               = typename std::vector<T, Allocator>;
  using value_type              = typename base_type::value_type;
  using allocator_type          = typename base_type::allocator_type;
  using size_type               = typename base_type::size_type;
  using difference_type         = typename base_type::difference_type;
  using reference               = typename base_type::reference;
  using const_reference         = typename base_type::const_reference;
  using pointer                 = typename base_type::pointer;
  using const_pointer           = typename base_type::const_pointer;
  using iterator                = typename base_type::iterator;
  using const_iterater          = typename base_type::const_iterator;
  using reverse_iterator        = typename base_type::reverse_iterator;
  using const_reverse_iterater  = typename base_type::const_reverse_iterator;

  // vectorのコンストラクタを使用可能にする
  using base_type::vector;

  // 要素へのアクセスはKeyを通してのみ可能
  // Keyは少なくともsize_typeへ明示的に変換可能であることを要求する
  // ついでにgslのExpectsを使って境界チェックを行う
  reference operator[](const Key& key) {
    const auto i = static_cast<size_type>(key);
    Expects(i < this->size());
    return base_type::operator[](i);
  }

  const_reference operator[](const Key& key) const {
    const auto i = static_cast<size_type>(key);
    Expects(i < this->size());
    return base_type::operator[](i);
  }
};

# vertex_indexまたはedge_indexをKeyに指定することで、頂点または辺のどちらに属するかを示す。
template <typename T>
using vertex_property = property_map<T, vertex_index>;
template <typename T>
using edge_property = property_map<T, edge_index>;
```

最低限これだけ定義すれば十分です。ダイクストラ法の実装で

```cpp
using edge_list        = edge_property<edge>;
using adjacency_list   = vertex_property<std::vector<vertex_index>>;
using vertex_edge_list = vertex_property<std::vector<edge_index>>;
using edge_length_list = edge_property<double>;
// ...
```

とエイリアスを入れ替えるだけで使用可能です。