Commit fec365ec authored by Sébastien Piat's avatar Sébastien Piat
Browse files

eppstein: new algorithm

K Shortest Path algorithm implemented following the Eppstein implementation
described in the article: `Finding the k shortest paths`, David Eppstein (1997).

* tests/python/lightest.py: Test it here.
* vcsn/algos/dijkstra-node.hh: Here.
* vcsn/algos/eppstein.hh: Here.
* vcsn/algos/implicit-path.hh: Here.
* vcsn/algos/lightest.hh: Here.
* vcsn/algos/path.hh: Here.
* vcsn/algos/shortest-path-tree.hh: Here.
* vcsn/local.mk: Here.
parent a386d76b
......@@ -446,7 +446,7 @@ bench('a.shortest(5000)',
# lightest
ctx = 'lal_char(a-e), nmin'
r = "[a-e]?{150}"
for algo, number in [('auto', 10), ('yen', 500)]:
for algo, number in [('auto', 10), ('yen', 500), ('eppstein', 10)]:
bench('a.lightest(5, "{}")'.format(algo),
'a = std({}), c = {}'.format(r, ctx_signature(ctx)),
setup=['c = vcsn.context(ctx)',
......
......@@ -302,4 +302,8 @@ def PLAN():
print('PASS:', npass)
print('FAIL:', nfail)
def weightset_is(ctx, expected):
s = str(ctx).lower().split('-> ')
return 0 < len(s) and s[-1] == expected.lower()
atexit.register(PLAN)
......@@ -4,21 +4,22 @@ import vcsn
from test import *
algos = ['auto', 'a-star', 'bellman-ford', 'breadth-first', 'dijkstra', 'yen']
k_algos = ['auto', 'breadth-first', 'yen']
k_algos = ['auto', 'breadth-first', 'yen', 'eppstein']
def check(re, num, exp, tests = []):
a = ctx.expression(re, 'none').standard()
for algo in tests if tests else algos if num == 1 else k_algos:
print(algo + ':' + re)
p = a.lightest(num=num, algo=algo)
CHECK_EQ(exp, p)
if algo != 'eppstein' or weightset_is(ctx, 'nmin'):
print(algo + ':' + re)
p = a.lightest(num=num, algo=algo)
CHECK_EQ(exp, p)
ctx = vcsn.context('lal_char, nmin')
check(r'\z', 3, r'\z')
check(r'\e', 3, r'<0>\e')
check('a+b', 2, '<0>a + <0>b')
check('ababab', 10, '<0>ababab')
check('(<1>a+<1>b)*', 7, r'<0>\e + <1>a + <1>b + <2>aa + <2>ab + <2>ba + <2>bb', ['auto'])
check('<4>a+(<1>a<1>b)+<1>c+<2>d', 1, '<1>c')
check(r'\z', 3, r'\z', k_algos)
check(r'\e', 3, r'<0>\e', k_algos)
check('a+b', 2, '<0>a + <0>b', k_algos)
check('ababab', 10, '<0>ababab', k_algos)
check('(<1>a+<1>b)*', 7, r'<0>\e + <1>a + <1>b + <2>aa + <2>ab + <2>ba + <2>bb', ['auto', 'eppstein'])
check('<4>a+(<1>a<1>b)+<1>c+<2>d', 1, '<1>c', k_algos)
ctx = vcsn.context('lal_char(abcd), z')
check('[a-d]?{5}', 5, r'\e + aa + ab + ac + ad', ['auto'])
......
#pragma once
#include <limits>
#include <boost/optional.hpp>
#include <vcsn/ctx/traits.hh>
namespace vcsn
{
template <Automaton Aut>
class dijkstra_node
{
using automaton_t = Aut;
using weight_t = weight_t_of<automaton_t>;
using state_t = state_t_of<automaton_t>;
public:
dijkstra_node() = default;
dijkstra_node(state_t state, const boost::optional<weight_t>& weight,
state_t parent, unsigned depth = std::numeric_limits<unsigned>::max())
: state_{state}
, weight_{weight}
, parent_{parent}
, depth_{depth}
{}
bool
operator<(const dijkstra_node& other) const
{
if (!weight_)
return false;
else if (!other.weight_)
return true;
else
return *weight_ < *other.weight_;
}
weight_t
get_weight() const
{
// FIXME: max
return weight_ ? *weight_ : std::numeric_limits<weight_t>::max();
}
void
set_weight(weight_t weight)
{
weight_ = weight;
}
// FIXME: private
unsigned depth_;
state_t state_;
private:
boost::optional<weight_t> weight_;
public:
state_t parent_;
};
}
#pragma once
#include <vector>
#include <boost/heap/fibonacci_heap.hpp>
#include <vcsn/misc/unordered_map.hh>
#include <vcsn/algos/path.hh>
#include <vcsn/algos/shortest-path-tree.hh>
#include <vcsn/algos/implicit-path.hh>
namespace vcsn
{
/// Eppstein implementation of the K shortest path problem adapted to Vcsn.
///
/// Based on `Finding the k shortest paths`, David Eppstein (1997).
template <Automaton Aut>
class eppstein
{
using automaton_t = Aut;
using state_t = state_t_of<automaton_t>;
using transition_t = transition_t_of<automaton_t>;
using weight_t = weight_t_of<automaton_t>;
using implicit_path_t = implicit_path<automaton_t>;
using sidetrack_costs_t = std::unordered_map<transition_t, weight_t>;
public:
eppstein(const automaton_t& aut)
: aut_{aut}
{}
/// Elements of the heap
///
/// This is needed in order to unsure a max heap (operator< is actually
/// computing other < this).
struct profile
{
profile(const implicit_path_t& path)
: path_{path}
{}
profile(implicit_path_t&& path)
: path_{path}
{}
bool operator<(const profile& rhs) const
{
// We want a min-heap.
return rhs.path_ < path_;
}
implicit_path_t path_;
};
using queue_t = boost::heap::fibonacci_heap<profile>;
std::vector<path<automaton_t>>
k_shortest_path(state_t src, state_t dst, int K)
{
auto tree = compute_shortest_path_tree(aut_, dst);
auto sidetrack_edge_costs_map = sidetrack_edge_costs(aut_, tree);
auto res = std::vector<path<automaton_t>>{};
auto queue = queue_t{};
queue.emplace(implicit_path_t(aut_, aut_->null_transition(),
implicit_path_t::null_parent_path,
tree.states_[src].get_weight()));
for (int k = 0; k < K && !queue.empty(); k++)
{
auto k_path_implicit = std::move(queue.top().path_);
queue.pop();
auto k_path = k_path_implicit.explicit_path(res, tree, src);
if (k_path.path_.empty())
return res;
res.emplace_back(std::move(k_path));
add_children_to_queue_(sidetrack_edge_costs_map, src, k_path_implicit, k, queue);
}
return res;
}
sidetrack_costs_t
sidetrack_edge_costs(const automaton_t& aut, shortest_path_tree<automaton_t>& spt)
{
auto res = sidetrack_costs_t{};
for (auto tr : aut->all_transitions())
{
state_t parent = spt.get_parent_of(aut->src_of(tr));
if (parent == aut->null_state() || parent != aut->dst_of(tr))
res[tr] = (aut->weight_of(tr)
+ spt.states_[aut->dst_of(tr)].get_weight()
- spt.states_[aut->src_of(tr)].get_weight());
}
return res;
}
private:
void
add_children_to_queue_(sidetrack_costs_t& sidetracks, state_t src,
const implicit_path_t& k_path_implicit, int k, queue_t& queue)
{
auto k_path_cost = k_path_implicit.weight_;
auto transition_stack = std::vector<transition_t>{};
auto s = k == 0 ? src : aut_->dst_of(k_path_implicit.sidetrack_);
for (auto tr : all_out(aut_, s))
transition_stack.push_back(tr);
while (!transition_stack.empty())
{
transition_t curr = transition_stack.back();
transition_stack.pop_back();
if (has(sidetracks, curr))
queue.emplace(implicit_path_t(aut_, curr, k, k_path_cost + sidetracks[curr]));
else
for (auto tr : all_out(aut_, aut_->dst_of(curr)))
transition_stack.push_back(tr);
}
}
const automaton_t& aut_;
};
template <Automaton Aut>
std::enable_if_t<std::is_same<weightset_t_of<Aut>, vcsn::nmin>::value, std::vector<path<Aut>>>
compute_eppstein(const Aut& aut, state_t_of<Aut> src, state_t_of<Aut> dst, unsigned num)
{
auto ksp = eppstein<Aut>(aut);
return ksp.k_shortest_path(aut->pre(), aut->post(), num);
}
template <Automaton Aut>
std::enable_if_t<!std::is_same<weightset_t_of<Aut>, vcsn::nmin>::value, std::vector<path<Aut>>>
compute_eppstein(const Aut& aut, state_t_of<Aut> src, state_t_of<Aut> dst, unsigned num)
{
raise("eppstein: invalid algorithm with specified weightset: ",
*aut->weightset(), " is not nmin");
}
}
#pragma once
#include <vcsn/core/automaton.hh>
#include <vcsn/algos/shortest-path-tree.hh>
#include <vcsn/ctx/traits.hh>
namespace vcsn
{
template <Automaton Aut>
class implicit_path
{
using automaton_t = Aut;
using state_t = state_t_of<automaton_t>;
using transition_t = transition_t_of<automaton_t>;
using weight_t = weight_t_of<automaton_t>;
public:
static constexpr int null_parent_path = -1;
implicit_path(const automaton_t& aut, transition_t sidetrack,
int parent_path, weight_t weight)
: aut_{aut}
, sidetrack_{sidetrack}
, parent_path_{parent_path}
, weight_{weight}
{}
path<automaton_t>
explicit_path(const std::vector<path<automaton_t>>& ksp,
shortest_path_tree<automaton_t>& tree,
state_t src)
{
auto res = path<automaton_t>(aut_);
if (parent_path_ != null_parent_path)
{
auto& explicit_pref_path = ksp[parent_path_];
auto& path = explicit_pref_path.path_;
// The index of the last transition of the path before the sidetrack.
int last_transition_index = -1;
for (int i = path.size() - 1; 0 <= i; i--)
{
auto curr = path[i];
if (aut_->dst_of(curr) == aut_->src_of(sidetrack_))
{
last_transition_index = i;
break;
}
}
for (unsigned i = 0; i <= last_transition_index; i++)
res.emplace_back(aut_->weight_of(path[i]), path[i]);
res.emplace_back(aut_->weight_of(sidetrack_), sidetrack_);
}
auto s = parent_path_ == null_parent_path ? src : aut_->dst_of(sidetrack_);
while (s != tree.root_)
{
auto next = tree.get_parent_of(s);
if (s == aut_->null_state() || next == aut_->null_state())
return path<automaton_t>(aut_);
auto weight = tree.states_[s].get_weight() - tree.states_[next].get_weight();
auto t = find_transition(s, next);
assert(t != aut_->null_transition());
res.emplace_back(weight, t);
s = next;
}
return res;
}
transition_t
find_transition(state_t src, state_t dst) const
{
auto min_tr = aut_->null_transition();
for (auto tr : detail::outin(aut_, src, dst))
if (min_tr == aut_->null_transition()
|| aut_->weight_of(tr) < aut_->weight_of(min_tr))
min_tr = tr;
return min_tr;
}
bool operator<(const implicit_path& other) const
{
return weight_ < other.weight_;
}
// FIXME: private
const automaton_t& aut_;
transition_t sidetrack_;
int parent_path_;
weight_t weight_;
};
}
......@@ -7,6 +7,7 @@
#include <boost/heap/binomial_heap.hpp>
#include <vcsn/algos/lightest-path.hh>
#include <vcsn/algos/eppstein.hh>
#include <vcsn/algos/has-lightening-cycle.hh>
#include <vcsn/core/name-automaton.hh>
#include <vcsn/ctx/context.hh>
......@@ -200,6 +201,15 @@ namespace vcsn
ps.add_here(res, *m);
return res;
}
else if (algo == "eppstein" && std::is_same<weightset_t_of<Aut>, vcsn::nmin>::value)
{
const auto ps = make_word_polynomialset(aut->context());
auto res = ps.zero();
auto computed = compute_eppstein(aut, aut->pre(), aut->post(), num);
for (const auto& path : computed)
ps.add_here(res, path.make_monomial(ps));
return res;
}
else if ((algo == "auto" && num != 1) || algo == "breadth-first")
{
auto lightest = detail::lightest_impl<Aut>{aut};
......
#pragma once
#include <vector>
#include <vcsn/ctx/traits.hh>
namespace vcsn
{
template <Automaton Aut>
class path
{
using automaton_t = Aut;
using weight_t = weight_t_of<Aut>;
using transition_t = transition_t_of<automaton_t>;
using state_t = state_t_of<automaton_t>;
public:
weight_t
path_sum(const automaton_t& aut, const std::vector<transition_t>& elts)
{
weight_t res = weightset_t_of<automaton_t>::one();
for (auto tr : elts)
res += aut->weight_of(tr);
return res;
}
path(const automaton_t& aut, const std::vector<transition_t>& elts = {})
: path_{elts}
, weight_{path_sum(aut, elts)}
, aut_{aut}
{}
void
push_back(transition_t tr, weight_t weight)
{
path_.push_back(tr);
weight_ += weight;
}
template <typename... Args>
void emplace_back(weight_t weight, Args&&... args)
{
path_.emplace_back(std::forward<Args>(args)...);
weight_ += weight;
}
bool
operator<(const path& other) const
{
return weight_ < other.weight_;
}
template <typename PolynomialSet>
auto
make_monomial(const PolynomialSet& ps) const
-> typename PolynomialSet::monomial_t
{
const auto& pls = *ps.labelset();
const auto& pws = *ps.weightset();
const auto& ls = *aut_->labelset();
auto w = pws.one();
auto l = pls.one();
for (auto t : path_)
{
w = pws.mul(w, aut_->weight_of(t));
auto nl = aut_->label_of(t);
if (!ls.is_special(nl))
l = pls.mul(l, nl);
}
return {l, w};
}
// FIXME: private
std::vector<transition_t> path_;
weight_t weight_;
const automaton_t& aut_;
};
}
#pragma once
#include <vcsn/misc/map.hh>
#include <vcsn/algos/dijkstra-node.hh>
#include <vcsn/ctx/traits.hh>
namespace vcsn
{
template <typename Aut>
class shortest_path_tree
{
using automaton_t = Aut;
using state_t = state_t_of<automaton_t>;
using dijkstra_node_t = dijkstra_node<automaton_t>;
public:
shortest_path_tree(state_t root)
: root_{root}
{}
void
add(const dijkstra_node_t& n)
{
states_[n.state_] = n;
}
void
set_parent_of(state_t s, state_t parent)
{
if (has(states_, s))
states_[s].parent_ = parent;
else
states_[s] = dijkstra_node_t{s, {}, parent};
}
state_t
get_parent_of(state_t s)
{
auto it = states_.find(s);
if (it != states_.end())
return it->second.parent_;
else
return automaton_t::element_type::null_state();
}
state_t
get_parent_of(state_t s) const
{
auto it = states_.find(s);
if (it != states_.end())
return it->second.parent_;
else
return automaton_t::element_type::null_state();
}
// FIXME: private
std::map<state_t, dijkstra_node_t> states_;
state_t root_;
};
template <typename Aut>
shortest_path_tree<Aut>
compute_shortest_path_tree(const Aut& aut, state_t_of<Aut> src)
{
using automaton_t = Aut;
using dijkstra_node_t = dijkstra_node<automaton_t>;
struct profile
{
profile(const dijkstra_node_t& node)
: node_{node}
{}
profile(dijkstra_node_t&& node)
: node_{node}
{}
bool operator<(const profile& rhs) const
{
// We want a min-heap.
return rhs.node_ < node_;
}
dijkstra_node_t node_;
};
using queue_t = boost::heap::fibonacci_heap<profile>;
auto predecessor_tree = shortest_path_tree<automaton_t>(src);
auto queue = queue_t{};
for (auto s : aut->all_states())
predecessor_tree.add(dijkstra_node_t{s, {}, aut->null_state()});
auto src_node = predecessor_tree.states_[predecessor_tree.root_];
src_node.set_weight(weightset_t_of<automaton_t>::one());
src_node.depth_ = 0;
queue.emplace(src_node);
while (!queue.empty())
{
auto current = queue.top().node_;
queue.pop();
auto s = current.state_;
for (auto tr : all_in(aut, s))
{
auto& neighbor = predecessor_tree.states_[aut->src_of(tr)];
auto dist = neighbor.get_weight();
auto new_dist = current.get_weight() + aut->weight_of(tr);
if (new_dist < dist)
{
neighbor.set_weight(new_dist);
neighbor.depth_ = current.depth_ + 1;
neighbor.parent_ = s;
queue.emplace(neighbor);
}
}
}
return predecessor_tree;
}
}
......@@ -35,6 +35,7 @@ algo_headers = \
%D%/algos/detail/printer.hh \
%D%/algos/determinize.hh \
%D%/algos/dijkstra.hh \
%D%/algos/dijkstra-node.hh \
%D%/algos/distance.hh \
%D%/algos/divide.hh \
%D%/algos/divkbaseb.hh \
......@@ -42,6 +43,7 @@ algo_headers = \
%D%/algos/double-ring.hh \
%D%/algos/edit-automaton.hh \
%D%/algos/efsm.hh \
%D%/algos/eppstein.hh \
%D%/algos/epsilon-remover-distance.hh \
%D%/algos/epsilon-remover-lazy.hh \
%D%/algos/epsilon-remover-separate.hh \
......@@ -56,6 +58,7 @@ algo_headers = \
%D%/algos/has-lightening-cycle.hh \
%D%/algos/has-twins-property.hh \
%D%/algos/identities-of.hh \
%D%/algos/implicit-path.hh \
%D%/algos/inductive.hh \
%D%/algos/info.hh \
%D%/algos/insplit.hh \
......@@ -91,6 +94,7 @@ algo_headers = \
%D%/algos/pair.hh \
%D%/algos/partial-identity.hh \
%D%/algos/partial-identity-expression.hh \
%D%/algos/path.hh \
%D%/algos/prefix.hh \
%D%/algos/print.hh \
%D%/algos/project-automaton.hh \
......@@ -106,6 +110,7 @@ algo_headers = \
%D%/algos/reduce.hh \
%D%/algos/scc.hh \
%D%/algos/shortest.hh \