Boost1.55でMultiple Source Dijkstra

この記事はC++ (fork) Advent Calendar 2013の22日目の記事になります．
今回も軽めの記事です．ご了承ください．

Boost1.55ではBoost.Graphの関数boost::dijkstra_shortes_pathsがマルチソースに対応しましたので，本記事ではこれを紹介します．

dijkstra_shortest_paths

関数dijkstra_shortest_pathsはその名前からも分かる通り，ダイクストラ法を用いた最短経路木を探索するためのアルゴリズムです．通常，ダイクストラ法ではスタート地点となる節点を一つ選択し，この節点から他の節点までの最短経路を求めます．

マルチソースの場合はスタート地点(ルート)となる節点が複数になります．実際に図で見てみるとこのようになります(見た目上枝の長さは違いますが，論理的にはすべて同じにしてます)．
二重線で囲まれている節点(ラベルが0, 20, 40, 60, 80の節点)がルートになります．
マルチソースダイクストラでは，ルート以外の節点vと，各ルート間のすべての最短経路が計算されるのではなく，経路距離が最小となるルートとの経路のみが計算されます．そのため図中ではルート以外の節点は経路距離が最小となるようなルートと同じ色で塗り分けしています．
この図はグラフのボロノイ分割と呼ばれ，勢力図等を作成するのに役立ちます．例えば，ルートを消防署と考えれば，ある節点で火事が起きた際にどこの消防署から消防車を出せばよいかの指標になります．

マルチソース版の実装で，シングルソース版との違いは探索開始時に優先度付きキューの中身に複数のルートがあるという点だけですので，容易に実装できます(私も昔作りました)．

上記のような図をboost::dijkstra_shortes_pathsを用いて作成するコードは以下になります．

A example of multiple source dijkstra.

ルート節点の指定

注目すべき箇所は関数multisource_dijkstra内のboost::dijkstra_shortest_pathsの呼び出しです．
第一引数はシングルソースの場合と同じで探索対象のグラフですが，第二，第三引数がルートとなる節点の集合を示すイテレータペアになっています(Rangeで渡したいですね!!)．シングルソースの場合は，第二引数にルートとなる節点一つを指定していた場合とは異なります．それ以外の引数はシングルソースの場合と同じです．

各節点から最短のルートの記録

ここではparent_recorderという名前のEventVisitorを新たに作成しています(root_recorderという名前にすればよかった)．これは新たに探索された節点がどのルートから探索されたかを記憶します．親(predecessorであってparentではない)のルートを自節点のルートに設定するだけです．あらかじめルートとなる節点のルートをその節点に設定しておく必要はありますが，そこは我慢します．

void operator()(Edge const e, Graph const& g) const {
    auto value = get(parent_map_, source(e, g));
    put(parent_map_, target(e, g), value);
}

getとputを合わせて一行にせずにわざわざgetで取得した値をコピーしているのは，put時にparent_map_の内部で保持するvectorのresize()が呼ばれてgetで取得した値が破棄されるのを防ぐためです(というか一行で書いて嵌りました)．

引数の数

上記のdijkstra_shortes_pathsの呼び出しでは12個の引数を指定していますが，これらを省略することはできません．
通常Boost.Graphの関数はdjikstra_shortest_pathsも含め，Named Parameter Idiomを用いた名前付き引数によってデフォルト引数を持つパラメータに対して引数を指定する必要はありません．しかしながら，マルチソースの場合名前付き引数を持つオーバロードが一つもないので，すべての引数を指定する必要があります(さらにVertexColorMapを指定するオーバロードならありました．これで引数の数は13個!!)．
使い勝手はかなり悪いので，名前付き引数版をオーバロード可能かどうか調査して追加するか独自にラッパを作成する必要がありそうです．

グラフの視覚化

出力されたファイルはgraphvizを使えばいい感じにレイアウトしてくれます.

fdp -Nstyle=filled -Nshape=circle -Nperipheries=peripheries -Tpng -O graph.dot

まとめ

マルチソースダイクストラは現状undocumentedですが，シングルソース版の使い方が分かっていれば使用方法が分からないことはないと思います．また，Boost1.55ではダイクストラだけではなく，幅優先探索もマルチソースに対応しているので利用してみるといいかもしれません(今回の例では枝の長さがすべて等しいので本当はこっちを使うべきでした)．

アルゴリズムとか意味がなかった

https://paiza.jp/poh/ec-campaign に挑戦してみましたので，そのコードを載せておきます．

#include <iostream>
#include <vector>
#include <algorithm>

template <class Iterator>
inline int calc(Iterator const first, Iterator const last, int const value)
{
    if (*(last - 1) + *(last - 2) <= value) {
        return *(last - 1) + *(last - 2);
    }
    Iterator upper = std::lower_bound(first + 1, last, value / 2);
    if (upper + 1 != last && *upper + *(upper + 1) == value) {
        return value;
    }
    int best_price = 0;
    for (Iterator lower = std::upper_bound(first, upper, value - *upper);
         lower != first && upper != last;
         lower = std::upper_bound(first, lower - 1, value - *upper)) {
        upper = std::upper_bound(upper + 1, last, value - *(lower - 1));
        best_price = std::max(best_price, *(upper - 1) + *(lower - 1));
        if (best_price == value) {
            break;
        }
    }
    return best_price;
}

int main()
{
    std::cin.tie(0);
    std::ios::sync_with_stdio(false);
    int N, D;
    std::cin >> N >> D;
    std::vector<int> prices;
    prices.reserve(N);
    for (int i = 0; i < N; ++i) {
        int p;
        std::cin >> p;
        prices.push_back(p);
    }
    std::sort(prices.begin(), prices.end());
    for (int i = 0; i < D; ++i) {
        int m;
        std::cin >> m;
        std::cout << calc(prices.begin(), prices.end(), m) << '\n';
    }
}

結果は http://paiza.jp/poh/ec-campaign/result/16051cf4bd6d13df012bc008d1fdad13 です．case3は0.13secでした．

boost::adjacency_listのVertexListテンプレート引数とvertex_indexの関係

id:faith_and_brave:20131003:1380788496に関して少しばかし補足.

まず，知るべきことはboost::adjaceny_listのテンプレート引数を省略した場合である．この場合，以下の引数を指定したものとして扱われる．

boost::adjacency_list<
    boost::vecS,		// OutEdgeList
    boost::vecS,		// VertexList
    boost::directedS,	// Directed
    boost::no_property,	// VertexProperty
    boost::no_property,	// EdgeProperty
    boost::no_property,	// GraphProperty
    boost::listS		// EdgeList
>

ここでは各引数の詳細は説明しないが，ここで注目していただきたいのは2番目のVertexListとコメントしたテンプレート引数である．

このVertexListテンプレート引数にはグラフの節点（Vertex）を格納するコンテナを指定する．デフォルトではvecSとなっており，これはstd::vectorに対応している．ここで注意して頂きたいことはVertexListがvecSの場合，adjaceny_listは自動的に各節点にindex番号が割り当てられたIndexedGraph*1になるということである．

このときもう一つ注意すべきことは，このIndexedGraphではindexの値が自動で割り振られ，外部からは変更不可能ということである．このときの節点のindex値は[0, num_vertices(g))中の値に一対一に対応する（gはadjacency_listのオブジェクトである）．つまり，indexに任意の値を割り当てることができないことになる．

よって，VertexListがvecSでない場合はデフォルトではIndexedGraphにならない．以下のコードはエラーになる．

#include <iostream>
#include <boost/graph/adjacency_list.hpp>
#include <boost/range/algorithm/for_each.hpp>
#include <boost/range/counting_range.hpp>

int main()
{
    using Graph = boost::adjacency_list<
        boost::listS,
        boost::listS,		// VertexListはlistS
        boost::directedS,
        boost::no_property	// VertexPropertyはなし
	>;

    Graph g;
    for (auto i : boost::counting_range(0, 5)) {
        add_vertex(g);
    }

    using Vertex = boost::graph_traits<Graph>::vertex_descriptor;
    boost::for_each(vertices(g), [&](Vertex const v) {
        std::cout << get(boost::vertex_index, g, v) << std::endl;	// error gはIndexedGraphではない
    });
}

この場合，以下のように手動でグラフにindexを持たせる必要がある．

#include <iostream>
#include <boost/graph/adjacency_list.hpp>
#include <boost/range/algorithm/for_each.hpp>
#include <boost/range/counting_range.hpp>

int main()
{
    using IndexProp = boost::property<boost::vertex_index_t, std::size_t>;
    using Graph = boost::adjacency_list<
        boost::listS,
        boost::listS,	// VertexListはListS
        boost::directedS,
        IndexProp	// VertexPropertyはvertex_indexを持つproperty
    >;

    Graph g;
    for (auto i : boost::counting_range(0, 5)) {
        add_vertex(IndexProp(i * 2), g);	// ついでなのでindexを2の倍数とする
    }

    using Vertex = boost::graph_traits<Graph>::vertex_descriptor;
    boost::for_each(vertices(g), [&](Vertex const v) {
        std::cout << get(boost::vertex_index, g, v) << std::endl;
    });
}

ここではadjacency_listのVertexPropertyにvertex_indexを持つpropertyを指定した．これにより，get(boost::vertex_index, g, v)で節点vからindexを取得することが可能となる．

また，ここでは関数add_vertexの呼び出し時に各節点のindexの値を設定している．

実際はvertex_indexの値は上記のような不連続な値にせず，[0, num_vertices(g))に対応するように割り当てるべきである*2．こういったindex付けをしたい場合は，vertex_indexの代わりにvertex_index1やvertex_index2が用意されているのでこちらを用いるとよい．

*1:公式にはこのような呼び方はない

*2:この理由についてはいつか述べるかもしれない

Propety_iter_range

BGLといえばPropertyMap，PropertyMapといえばBGLなのだが，BGLのアダプタ関数はイテレータのペアであるRangeを返すなのでBGLとPropertyMap，Rangeをうまく組み合わせて使いたい．
で少し調べてみたら，BGL内にproperty_iter_rangeというのを発見．property_iteratorはイテレータが指し示す値に，その値をキーとするプロパティマップに適用するイテレータアダプタ．このペアを作るのがproperty_iter_range．
以下サンプル

#include <iostream>
#include <random>
#include <ctime>
#include <boost/graph/adjacency_list.hpp>
#include <boost/graph/plod_generator.hpp>
#include <boost/graph/property_iter_range.hpp>  // property_iter_range
#include <boost/graph/dijkstra_shortest_paths.hpp>
#include <boost/property_map/property_map.hpp>
#include <boost/range/adaptor/filtered.hpp>
#include <boost/range/algorithm/generate.hpp>
#include <boost/range/numeric.hpp>

int main()
{
	typedef boost::adjacency_list<
				boost::vecS, boost::vecS, boost::undirectedS,
				boost::property<boost::vertex_root_t, bool,
					boost::property<boost::vertex_distance_t, int>>,
				boost::property<boost::edge_weight_t, int>>
			Graph;

	typedef boost::plod_iterator<std::mt19937, Graph> PlodGen;

	std::size_t const num_v = 20;
	std::mt19937 gen{
		static_cast<std::mt19937::result_type>(time(nullptr))
	};
	// ランダムグラフの生成
	Graph graph{
		PlodGen{gen, num_v, 1.5, 200},
		PlodGen{},
		num_v
	};

	std::uniform_int_distribution<> bool_dist{0, 1};
	boost::generate(
			boost::get_property_iter_range(graph, boost::vertex_root),
			std::bind(bool_dist, gen));

	std::uniform_int_distribution<> dist{0, 15};
	boost::generate(
			boost::get_property_iter_range(graph, boost::edge_weight),
			std::bind(dist, gen));

	boost::dijkstra_shortest_paths(
			graph,
			vertex(std::uniform_int_distribution<>{0, num_v - 1}(gen), graph),
			boost::weight_map(get(boost::edge_weight, graph)).
			distance_map(get(boost::vertex_distance, graph)));

        // rootから各vertexまでの距離の合計の計算
	auto const is_root_func
		= boost::make_property_map_function(get(boost::vertex_root, graph));
	auto const vrange
		= vertices(graph) | boost::adaptors::filtered(is_root_func);
	int total = 0;
	for (auto const v : vrange) {
		total += get(boost::vertex_distance, graph, v);
	}
	std::cout << total << std::endl;

	return 0;
}

正直これはむむむ...propertyp_iteratorのレンジを得るにはget_property_iter_rangeという関数にグラフオブジェクトとPropertyTagのオブジェクトを渡す必要がある．property_iter_range.hppのファイルがPropertyMapではなくBGL内にあるからグラフに特化したインタフェースなのかもしれないが，これだと使いにくい．
例えば上のサンプルの最後の距離の合計の計算の場合，フィルタを掛けたvrangeにプロパティマップでマッピングできた方がいい．絶対にいい！
というわけでRangeで一般的なAdaptorを書いてみた．
https://github.com/amedama41/Canard/blob/master/property_map/property_map_adaptor.hpp
ついてでにBoostのproperty_iteratorはWritable PropertyMapには適用できなかったので，Writable PropertyMapにも対応させてみた．
これを使えばサンプルの計算も

boost::accumulate(
		vrange | Canard::adaptors::mapped(get(boost::vertex_distance, graph)),
		0.0);

みたいに書けていい感じなのではないだろうか

GitHubのリポジトリからのみファイルを削除する

ローカルリポ上のファイルは残しながら，GitHubのリポからファイルを削除する方法が分からなかったのでメモ．

git rm --cached filename

でOK．初めは-nかと思ってたのだが何も変わらなくて？だった．

Graph Layoutが使えない件

BGLのGraph Layoutは名前そのままグラフのレイアウトを行うアルゴリズム群．
Random LayoutやCircle Layoutなど基本的なものから，Kameda Kawaiアルゴリズムといったものまで揃ってる．
これらのアルゴリズムはTopologyと呼ばれるTraitsを持つことで，いろいろなグラフに適用可能となっている．

Topologyが満たすべき要件はたぶんこれ↓
http://www.boost.org/doc/libs/1_48_0/libs/graph/doc/topology.html
なのだが，これだけだとコンパイルエラーになるアルゴリズムがある．

graph/topology.hppのソースを見ると多分Valid Expressioは多分こんな感じになる

Expression	Type	Description
Topology::point_type	type	空間内の座標を表す型
space.random_point()	point_type	戻り値：空間内からランダムに選択した点
space.distance(p1, p2)	double	戻り値：座標間の距離
space.move_position_toward(p1, fraction, p2)	point_type	戻り値：p1,p2線分上の点p．p1,p:p:p2 = fraction : 1-fraction となる
point_type::dimension	int型に変換可能な型	空間の次数
Topology::point_difference_type	type	二つのpoint_type p1,p2の差を表す型
point_difference_type::dimension	int型に変換可能な型	空間の次数
space.difference(p1, p2)	point_difference_type	戻り値：p1とp2の差
space.adjust(p1, delta)	point_type	戻り値：difference(p, p1) == deltaとなるようなp
space.pointwise_min(p1, p2)	point_type	戻り値：成分iがmin(p1[i], p1[i])となるようなp
space.pointwise_max(p1, p2)	point_type	戻り値：成分iがmax(p1[i], p1[i])となるようなp
space.norm(delta)	double	戻り値：deltaのユークリッドノルム sqrt(delta[0]^2 + delta[1]^2 + ...)
space.volume(delta)	double	戻り値：deltaの各成分の積？

あってるかは知らん

追記

これだけじゃ動かん...boundとかextendがいる模様

コンテナとは違うのだよコンテナとは

はにゃー．修論提出・発表と論文提出が被って停滞してしまっていました．修論まだ提出できる段階じゃないけど徐々に復活でふ！

What is the difference between a std::map and a boost::associative_property_map ?
When a nonexisting key value is given (e.g. "khwgfruiwgp" in the example above), the boost::associative_property_map's internal asserts will find it. For the std::map you have to write these asserts yourself every time.

You cannot iterate through a boost::associative_property_map, where with a std::map you can.

boost::associative_property_mapはstd::mapをLvaluePropertyMapに変換するアダプタクラス．公式ドキュメントにもサンプルコードがあるのだけれども，これだけ見ちゃうとstd::mapにgetとputを適用可能にしただけに見えて有り難味を感じないのよね．
http://www.boost.org/doc/libs/1_48_0/libs/property_map/doc/associative_property_map.html
例えば，

template <class AddressMap, class Key, class Value>
void write(AddressMap& address_map, const Key& key, const Value& value)
{
    address_map[key] = value;
}

のようなoperator経由で値を書き込む関数を定義すると，std::mapにもassociative_property_mapにもoperatorが定義されているので，

int main()
{
    typedef std::map<std::string, std::string> Name2Address;
    Name2Address name2address = {
        {"Fred", "710 West 13th Street"},
        {"Joe",  "710 West 13th Street"}
    };

    typedef boost::associative_property_map<Name2Address> AddressMap;
    AddressMap address_map(name2address);    // PropertyMap

    write(name2address, "Fred", "325 Cushing Avenue");
    write(address_map,  "Joe",  "325 Cushing Avenue");

のように呼び出すことができます．
これだと上のように何が違うのっていう疑問がでても仕方がないかなと思います．

しかしstd::mapを直接書き換えるか，associative_property_mapを通して書き換えるかはカプセル化の話を抜きにしても実際は全く違います．

const性の違い

例えば例えば，std::mapのインスタンスのconst参照を先ほどの書き込み関数に渡した場合，

    const Name2Address& const_name2address = name2address;
    write(const_name2address, "Fred", "Dokoka To-ku");    // エラー！

になるのは当然ですね？（書き込み以前にoperator[]を呼び出せませんが）
一方，PropertyMapインスタンスのconst参照を渡すと，

    const AddressMap& const_address_map = address_map;
    write(const_address_map, "Joe", "Dokoka To-ku");     // エラー...じゃない！

ちゃんと書き換えれます！不思議！

値渡しの違い

今度は先ほどのwriteの値渡しバージョンを使います．

template <class AddressMap, class Key, class Value>
void write_by_value(AddressMap address_map, const Key& key, const Value& value)
{
    address_map[key] = value;
}

std::mapに対して呼び出しても，

    write_by_value(name2address, "Fred", "Arumikan No Ue");
    std::cout << name2address["Fred"] << std::endl;    // "710 West 13th Street"を出力

もちろん値は書き換わりません．
しかし，PropertyMapに対しては，

    write_by_value(address_map, "Joe", "Arumikan No Ue");
    std::cout << address_map["Joe"] << std::endl;      // "Arumikan No Ue"を出力！
    assert(address_map["Joe"] == name2address["Joe"]);

書き換わってしまいます！

まとめ

このようなconst渡し，値渡しによる値の変化はassociative_property_mapだけでなく，他のプロパティマップの場合でも同様に起こります．これはPropetyMapがstd::mapのような値の集合を持つコンテナではないことを意味しています．
つまり，↓このようにPropetyMapから直接値を取り出すような関係でなく，

↓このようにPropertyMapはkeyとContainer（もちろんContainerの部分は関数でもなんでもかまいません）とを仲介するような関係になっています．

PropetyMapはkeyと値をどのように対応づけるか，どのように取り出すかしか関係しません．そのためPropetyMap自身のconst性はContainer自身のconst性には関係せず，PropetyMapのコピーでも同じContainerを参照できます．

PropertyMapのこのconst性のため，const関数でPropertyMapを用いると意図しない値の書き換えが起きてしまいそうですが，この場合，ReadOnlyのPropertyMapの使用を強制することで問題を回避できます．PropertyMapは一つのContainerに対して複数定義することができるのです．これが上の図にPropertyMapを二つ書いた理由にもなるのですが，一方はReadWrite可能にし，もう一方はReadOnlyに限定するということができるのです．

PropertyMapについての利点，概念等はCryolite先生がわかりやすく纏めてくれているので，そちらは見ましょう！
http://www.slideshare.net/Cryolite/boostpropertymap-pptx