tree_augmentation.hpp

Go to the documentation of this file.

//=======================================================================
// Copyright (c) 2013 Attila Bernath, Piotr Wygocki
//               2014 Piotr Godlewski, Piotr Smulewicz
// Distributed under the Boost Software License, Version 1.0. (See
// accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
//=======================================================================
#ifndef PAAL_TREE_AUGMENTATION_HPP
#define PAAL_TREE_AUGMENTATION_HPP


#include "paal/iterative_rounding/ir_components.hpp"
#include "paal/iterative_rounding/iterative_rounding.hpp"
#include "paal/utils/hash.hpp"

#include <boost/bimap.hpp>
#include <boost/graph/breadth_first_search.hpp>
#include <boost/graph/connected_components.hpp>
#include <boost/graph/filtered_graph.hpp>
#include <boost/graph/graph_traits.hpp>
#include <boost/graph/graph_utility.hpp>
#include <boost/graph/named_function_params.hpp>
#include <boost/graph/stoer_wagner_min_cut.hpp>
#include <boost/property_map/property_map.hpp>
#include <boost/range/as_array.hpp>
#include <boost/range/distance.hpp>

#include <algorithm>
#include <utility>

namespace paal {
namespace ir {

namespace detail {

template <class EdgeListGraph>
int filtered_num_edges(const EdgeListGraph & g) {
    return boost::distance(edges(g));
}

template <typename EdgeBoolMap> struct bool_map_to_tree_filter {
    bool_map_to_tree_filter() {}

    bool_map_to_tree_filter(EdgeBoolMap m) : ebmap(m) {}

    template <typename Edge> bool operator()(const Edge &e) const {
        return get(ebmap, e);
    }

    EdgeBoolMap ebmap;
};

template <typename EdgeBoolMap> struct bool_map_to_non_tree_filter {
    bool_map_to_non_tree_filter() {}

    bool_map_to_non_tree_filter(EdgeBoolMap m) : ebmap(m) {}

    template <typename Edge> bool operator()(const Edge &e) const {
        return !get(ebmap, e);
    }

    EdgeBoolMap ebmap;
};

template <typename KeyType, int num>
class const_int_map
    : public boost::put_get_helper<int, const_int_map<KeyType, num>> {
  public:
    using category = boost::readable_property_map_tag;
    using value_type = int;
    using reference = int;
    using key_type = KeyType;
    reference operator[](KeyType e) const { return num; }
};

} // namespace detail

namespace {
struct ta_compare_traits {
    static const double EPSILON;
};

const double ta_compare_traits::EPSILON = 1e-10;
}

struct ta_round_condition {
    ta_round_condition(double epsilon = ta_compare_traits::EPSILON)
        : m_round_half(epsilon) {}

    template <typename Problem, typename LP>
    boost::optional<double> operator()(Problem &problem, LP &lp,
                                       lp::col_id col) {
        auto res = m_round_half(problem, lp, col);
        if (res) {
            problem.add_to_solution(col);
        }
        return res;
    }

  private:
    round_condition_greater_than_half m_round_half;
};

struct ta_set_solution {
    ta_set_solution(double epsilon = ta_compare_traits::EPSILON)
        : m_compare(epsilon) {}

    template <typename Problem, typename GetSolution>
    void operator()(Problem &problem, const GetSolution &solution) {
        for (auto e :
             boost::as_array(edges(problem.get_links_graph()))) {
            if (!problem.is_in_solution(e)) {
                auto col = problem.edge_to_col(e);
                if (m_compare.e(solution(col), 1)) {
                    problem.add_to_solution(col);
                }
            }
        }
    }

private:
    const utils::compare<double> m_compare;
};

class ta_init {
  public:
    template <typename Problem, typename LP>
    void operator()(Problem &problem, LP &lp) {
        problem.init();
        lp.set_lp_name("Tree augmentation");
        lp.set_optimization_type(lp::MINIMIZE);

        add_variables(problem, lp);
        add_cut_constraints(problem, lp);
    }

  private:
    template <typename Problem, typename LP>
    void add_variables(Problem & problem, LP & lp) {
        for (auto e : boost::as_array(edges(problem.get_links_graph()))) {
            lp::col_id col_idx = lp.add_column(problem.get_cost(e), 0);
            problem.bind_edge_to_col(e, col_idx);
        }
    }

    template <typename Problem, typename LP>
    void add_cut_constraints(Problem &problem, LP &lp) {
        for (auto e :
             boost::as_array(edges(problem.get_tree_graph()))) {
            lp::linear_expression expr;
            for (auto pe : problem.get_covered_by(e)) {
                expr += problem.edge_to_col(pe);
            }

            lp::row_id row_idx = lp.add_row(std::move(expr) >= 1);
            problem.bind_edge_to_row(e, row_idx);
        }
    }
};

template <
    typename Init = ta_init,
    typename RoundCondition = ta_round_condition,
    typename RelaxContition = utils::always_false,
    typename SetSolution = ta_set_solution>
        using tree_augmentation_ir_components = IRcomponents<Init, RoundCondition,
                        RelaxContition, SetSolution>;

template <typename Graph, typename TreeMap, typename CostMap,
          typename VertexIndex, typename EdgeSetOutputIterator>
class tree_aug {
  public:
    using Edge = typename boost::graph_traits<Graph>::edge_descriptor;
    using Vertex = typename boost::graph_traits<Graph>::vertex_descriptor;
    using CostValue = double;

    using TreeGraph =
        boost::filtered_graph<Graph, detail::bool_map_to_tree_filter<TreeMap>>;
    using NonTreeGraph = boost::filtered_graph<
        Graph, detail::bool_map_to_non_tree_filter<TreeMap>>;

    using EdgeList = std::vector<Edge>;
    using CoverMap = std::unordered_map<Edge, EdgeList, edge_hash<Graph>>;

    // cross reference between links and columns
    using EdgeToColId = boost::bimap<Edge, lp::col_id>;
    using RowIdToEdge = std::unordered_map<lp::row_id, Edge>;

    using ErrorMessage = boost::optional<std::string>;

    tree_aug(const Graph & g, TreeMap tree_map, CostMap cost_map,
            VertexIndex vertex_index, EdgeSetOutputIterator solution) :
        m_g(g), m_tree_map(tree_map), m_cost_map(cost_map), m_index(vertex_index),
        m_solution(solution),
        m_tree(m_g, detail::bool_map_to_tree_filter<TreeMap>(m_tree_map)),
        m_ntree(m_g, detail::bool_map_to_non_tree_filter<TreeMap>(m_tree_map)),
        m_sol_cost(0)
    {}

    ErrorMessage check_input_validity() {
        // Num of edges == num of nodes-1 in the tree?
        int n_v = num_vertices(m_g);
        int n_e = filtered_num_edges(m_tree);

        if (n_e != n_v - 1) {
            return "Incorrect number of edges in the spanning tree. "
                        + std::string("Should be ") + std::to_string(n_v - 1)
                        + ", but it is " + std::to_string(n_e) + ".";
        }

        // Is the tree connected?
        std::vector<int> component(num_vertices(m_g));
        int num = boost::connected_components(m_tree, &component[0]);
        if (num > 1) {
            return ErrorMessage{ "The spanning tree is not connected." };
        }

        // Is the graph 2-edge-connected?
        detail::const_int_map<Edge, 1> const_1_edge_map;
        // TODO This stoer-wagner algorithm is unnecessarily slow for some reason
        int min_cut = boost::stoer_wagner_min_cut(m_g, const_1_edge_map);
        if (min_cut < 2) {
            return ErrorMessage{"The graph is not 2-edge-connected."};
        }

        return ErrorMessage{};
    }

    const NonTreeGraph &get_links_graph() const { return m_ntree; }

    const TreeGraph &get_tree_graph() const { return m_tree; }

    auto get_cost(Edge e)->decltype(get(std::declval<CostMap>(), e)) {
        return get(m_cost_map, e);
    }

    void add_to_solution(lp::col_id col) {
        *m_solution = m_edge_to_col_id.right.at(col);
        ++m_solution;
        m_sol_cost += m_cost_map[m_edge_to_col_id.right.at(col)];
        m_edge_to_col_id.right.erase(col);
    }

    void bind_edge_to_col(Edge e, lp::col_id col) {
        auto tmp =
            m_edge_to_col_id.insert(typename EdgeToColId::value_type(e, col));
        assert(tmp.second);
    }

    void bind_edge_to_row(Edge e, lp::row_id row) {
        auto tmp =
            m_row_id_to_edge.insert(typename RowIdToEdge::value_type(row, e));
        assert(tmp.second);
    }

    void init() {
        // We need to fill very useful auxiliary data structures:
        //\c m_covered_by - containing lists of
        // edges. For a tree edge \c t the list \c m_covered_by[t]
        // contains the list of links covering \c t.

        std::vector<Edge> pred(num_vertices(m_g));
        auto pred_map =
            boost::make_iterator_property_map(pred.begin(), m_index);
        std::set<Vertex> seen;
        for (auto u : boost::as_array(vertices(m_g))) {
            auto tmp = seen.insert(u);
            assert(tmp.second);
            boost::breadth_first_search(
                m_tree, u, boost::visitor(boost::make_bfs_visitor(
                               boost::record_edge_predecessors(
                                   pred_map, boost::on_tree_edge()))));

            for (auto e : boost::as_array(out_edges(u, m_ntree))) {
                auto node = target(e, m_ntree);
                if (!seen.count(node)) {
                    while (node != u) {
                        m_covered_by[get(pred_map, node)].push_back(e);
                        node = source(get(pred_map, node), m_tree);
                    }
                }
            }
        }
    }

    Edge row_to_edge(lp::row_id row) const { return m_row_id_to_edge.at(row); }

    lp::col_id edge_to_col(Edge e) const { return m_edge_to_col_id.left.at(e); }

    EdgeList &get_covered_by(Edge e) { return m_covered_by[e]; }

    bool is_in_solution(Edge e) const {
        return m_edge_to_col_id.left.find(e) == m_edge_to_col_id.left.end();
    }

    CostValue get_solution_cost() const { return m_sol_cost; }

  private:

    const Graph &m_g;
    TreeMap m_tree_map;
    CostMap m_cost_map;
    VertexIndex m_index;

    EdgeSetOutputIterator m_solution;

    EdgeToColId m_edge_to_col_id;

    TreeGraph m_tree;
    NonTreeGraph m_ntree;
    CostValue m_sol_cost;
    CoverMap m_covered_by;
    RowIdToEdge m_row_id_to_edge;
};

namespace detail {
template <typename Graph, typename TreeMap, typename CostMap,
          typename VertexIndex, typename EdgeSetOutputIterator>
tree_aug<Graph, TreeMap, CostMap, VertexIndex, EdgeSetOutputIterator>
make_tree_aug(const Graph & g, TreeMap tree_map, CostMap cost_map,
                VertexIndex vertex_index, EdgeSetOutputIterator solution) {
    return paal::ir::tree_aug<Graph, TreeMap, CostMap,
            VertexIndex, EdgeSetOutputIterator>(g, tree_map, cost_map, vertex_index, solution);
}

template <typename Graph, typename TreeMap, typename CostMap,
          typename VertexIndex, typename EdgeSetOutputIterator,
          typename IRcomponents = tree_augmentation_ir_components<>,
          typename Visitor = trivial_visitor>
IRResult tree_augmentation_iterative_rounding(
        const Graph & g,
        TreeMap tree_map,
        CostMap cost_map,
        VertexIndex vertex_index,
        EdgeSetOutputIterator solution,
        IRcomponents components = IRcomponents(),
        Visitor visitor = Visitor()) {
    auto treeaug = make_tree_aug(g, tree_map, cost_map, vertex_index, solution);
    return solve_iterative_rounding(treeaug, std::move(components), std::move(visitor));
}
} // detail

template <typename Graph, typename EdgeSetOutputIterator, typename P,
          typename T, typename R>
auto make_tree_aug(const Graph &g,
                   const boost::bgl_named_params<P, T, R> &params,
                   EdgeSetOutputIterator solution)
    ->tree_aug<
          Graph,
          decltype(choose_const_pmap(get_param(params, boost::edge_color), g,
                                     boost::edge_color)),
          decltype(choose_const_pmap(get_param(params, boost::edge_weight), g,
                                     boost::edge_weight)),
          decltype(choose_const_pmap(get_param(params, boost::vertex_index), g,
                                     boost::vertex_index)),
          EdgeSetOutputIterator> {
    return detail::make_tree_aug(
        g, choose_const_pmap(get_param(params, boost::edge_color), g,
                             boost::edge_color),
        choose_const_pmap(get_param(params, boost::edge_weight), g,
                          boost::edge_weight),
        choose_const_pmap(get_param(params, boost::vertex_index), g,
                          boost::vertex_index),
        solution);
}

template <typename Graph, typename EdgeSetOutputIterator>
auto make_tree_aug(const Graph &g, EdgeSetOutputIterator solution)
    ->decltype(make_tree_aug(g, boost::no_named_parameters(), solution)) {
    return make_tree_aug(g, boost::no_named_parameters(), solution);
}

template <typename Graph, typename EdgeSetOutputIterator,
          typename IRcomponents = tree_augmentation_ir_components<>,
          typename Visitor = trivial_visitor, typename P, typename T,
          typename R>
IRResult tree_augmentation_iterative_rounding(
    const Graph &g, const boost::bgl_named_params<P, T, R> &params,
    EdgeSetOutputIterator solution, IRcomponents components = IRcomponents(),
    Visitor visitor = Visitor()) {
    return detail::tree_augmentation_iterative_rounding(
        g, choose_const_pmap(get_param(params, boost::edge_color), g,
                             boost::edge_color),
        choose_const_pmap(get_param(params, boost::edge_weight), g,
                          boost::edge_weight),
        choose_const_pmap(get_param(params, boost::vertex_index), g,
                          boost::vertex_index),
        std::move(solution), std::move(components), std::move(visitor));
}

template <typename Graph, typename EdgeSetOutputIterator,
          typename IRcomponents = tree_augmentation_ir_components<>,
          typename Visitor = trivial_visitor>
IRResult tree_augmentation_iterative_rounding(const Graph &g,
                                              EdgeSetOutputIterator solution,
                                              IRcomponents components =
                                                  IRcomponents(),
                                              Visitor visitor = Visitor()) {
    return tree_augmentation_iterative_rounding(
        g, boost::no_named_parameters(), std::move(solution),
        std::move(components), std::move(visitor));
}

} // ir
} // paal

#endif // PAAL_TREE_AUGMENTATION_HPP