はじめに
研究の過程でメッシュクラスを自作する必要に迫られたため、BGL(Boost Graph Library)やCGAL(Computational Geometry Algorithms Library)のドキュメントやソースコードを読んでメッシュクラスの実装を勉強しました。その過程でグラフ理論に興味を持ったので、自分でダイクストラ法をC++で書いてみました。
実装
ダイクストラ法の実装は、Wikipediaのダイクストラ法の擬似コードをそのまま使いました。
単体テストに[アルゴリズム] ダイクストラ法をやってみるの記事の問題を使わせていただき、実装が正しいことを確認しています。
ここでの実装方法は、[アルゴリズム] ダイクストラ法をやってみるの記事やその実装の元となるブログとは大きく異なっています。
これらの記事では、頂点の接続やスタートからの最短経路の距離などの情報を全て保持するnodeクラス
struct node {
std::vector<int> edges_to; // この頂点に接続している辺番号リスト
std::vector<double> edges_dist; // 接続されている辺の長さ
bool done; // 確定ノードか否か
double dist; // スタートからの最短経路の距離
};
を作成しています。
私の経験上、このようなオブジェクト指向の実装の仕方だと
- スケーラブルでない
- 変更に対して弱い
- パフォーマンスが低い
ので、ここではデータ指向かつBGL的な実装をしています。
注：以下では読みやすさを優先して、境界チェックを行うvector::at()やgsl::at()を使っていません。
インデックス
頂点と辺の番号を区別するため、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）をエイリアスとして定義します。
struct vertex_tag {};
struct edge_tag {};
using vertex_index = index_base<vertex_tag>;
using edge_index = index_base<edge_tag>;
頂点
今回は頂点＝頂点番号ということにして、頂点は定義しませんでした。
辺
（無向）辺は２つの頂点のペアとして定義しました。
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で使われるリストです。
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を使うとキャッシュミスが減り最適化されます。
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の実装は以下のとおりです。
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;
}
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の実装は次のようになります。
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{};
}
ダイクストラ法
以上でダイクストラ法を実装する準備ができました。グラフ、始点、終点、辺の長さを与えた時に最短経路の長さを返す関数を定義します。
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のダイクストラ法の擬似コード通りになっているため、理解しやすいと思います。
そのうちstd::priority_queueを使った実装もしてみようと思います。
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に相当するものを作ることにします。
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>;
最低限これだけ定義すれば十分です。ダイクストラ法の実装で
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>;
// ...
とエイリアスを入れ替えるだけで使用可能です。