text: // Copyright (c) 2014  GeometryFactory (France).  All rights reserved.
  //
  // This file is part of CGAL (www.cgal.org)
  //
  // $URL$
  // $Id$
  // SPDX-License-Identifier: LGPL-3.0-or-later OR LicenseRef-Commercial
  //
  //
  // Author(s)     : Ilker O. Yaz, Simon Giraudot
  
  #ifndef CGAL_BOOST_GRAPH_ALPHA_EXPANSION_GRAPHCUT_H
  #define CGAL_BOOST_GRAPH_ALPHA_EXPANSION_GRAPHCUT_H
  
  #include <CGAL/Iterator_range.h>
  #include <CGAL/assertions.h>
  #include <CGAL/property_map.h>
  #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
  #include <CGAL/Timer.h>
  #endif
  #include <CGAL/IO/trace.h>
  
  #include <CGAL/boost/graph/Named_function_parameters.h>
  #include <CGAL/boost/graph/named_params_helper.h>
  
  #include <boost/version.hpp>
  
  #include <boost/graph/adjacency_list.hpp>
  #include <boost/graph/compressed_sparse_row_graph.hpp>
  #include <boost/graph/boykov_kolmogorov_max_flow.hpp>
  
  #include <vector>
  
  
  
  
  namespace CGAL
  {
  
  /// \cond SKIP_IN_MANUAL
  namespace internal
  {
  
  struct Alpha_expansion_old_API_wrapper_graph
  {
    typedef std::size_t vertex_descriptor;
    typedef std::size_t edge_descriptor;
    typedef boost::directed_tag directed_category;
    typedef boost::disallow_parallel_edge_tag edge_parallel_category;
    typedef boost::edge_list_graph_tag traversal_category;
  
    typedef boost::counting_iterator<std::size_t> counting_iterator;
    typedef CGAL::Iterator_range<counting_iterator> counting_range;
  
    typedef CGAL::Identity_property_map<std::size_t> Vertex_index_map;
    typedef CGAL::Pointer_property_map<std::size_t>::type Vertex_label_map;
  
    struct Vertex_label_cost_map
    {
      typedef std::size_t key_type;
      typedef std::vector<double> value_type;
      typedef value_type reference;
      typedef boost::readable_property_map_tag category;
  
      const std::vector<std::vector<double> >* cost_matrix;
  
      Vertex_label_cost_map (const std::vector<std::vector<double> >* cost_matrix)
        : cost_matrix (cost_matrix)
      { }
  
      friend reference get (const Vertex_label_cost_map& pmap, key_type idx)
      {
        std::vector<double> out;
        out.reserve (pmap.cost_matrix->size());
        for (std::size_t i = 0; i < pmap.cost_matrix->size(); ++ i)
          out.push_back ((*pmap.cost_matrix)[i][idx]);
        return out;
      }
  
    };
  
    typedef CGAL::Pointer_property_map<double>::const_type Edge_cost_map;
  
    const std::vector<std::pair<std::size_t, std::size_t> >& edges;
    const std::vector<double>& edge_costs;
    const std::vector<std::vector<double> >& cost_matrix;
    std::vector<std::size_t>& labels;
  
    Alpha_expansion_old_API_wrapper_graph (const std::vector<std::pair<std::size_t, std::size_t> >& edges,
                                           const std::vector<double>& edge_costs,
                                           const std::vector<std::vector<double> >& cost_matrix,
                                           std::vector<std::size_t>& labels)
      : edges (edges), edge_costs (edge_costs), cost_matrix (cost_matrix), labels (labels)
    { }
  
    friend counting_range vertices (const Alpha_expansion_old_API_wrapper_graph& graph)
    {
      return CGAL::make_range (boost::counting_iterator<std::size_t>(0),
                               boost::counting_iterator<std::size_t>(graph.labels.size()));
    }
  
    friend std::size_t num_vertices (const Alpha_expansion_old_API_wrapper_graph& graph) { return graph.labels.size(); }
  
    friend counting_range edges (const Alpha_expansion_old_API_wrapper_graph& graph)
    {
      return CGAL::make_range (boost::counting_iterator<std::size_t>(0),
                               boost::counting_iterator<std::size_t>(graph.edges.size()));
    }
  
    friend vertex_descriptor source (edge_descriptor ed, const Alpha_expansion_old_API_wrapper_graph& graph)
    { return graph.edges[ed].first; }
    friend vertex_descriptor target (edge_descriptor ed, const Alpha_expansion_old_API_wrapper_graph& graph)
    { return graph.edges[ed].second; }
  
    Vertex_index_map vertex_index_map() const { return Vertex_index_map(); }
    Vertex_label_map vertex_label_map() { return CGAL::make_property_map(labels); }
    Vertex_label_cost_map vertex_label_cost_map() const
    { return Vertex_label_cost_map(&cost_matrix); }
    Edge_cost_map edge_cost_map() const { return CGAL::make_property_map(edge_costs); }
  };
  
  ////////////////////////////////////////////////////////////////////////////////////////
  //   Comments about performance:
  //
  // 1) With BGL:
  //     * Using adjacency_list:
  //     ** Without pre-allocating vertex-list
  //       | OutEdgeList | VertexList | Performance |
  //       |    listS    |   listS    |     25.2    |
  //       |    vecS     |   listS    |     22.7    |
  //       |    listS    |   vecS     |     30.7    |
  //       |    vecS     |   vecS     |     26.1    |
  //
  //     ** With pre-allocating vertex-list with max-node size
  //        (Note: exact number of vertices are not certain at the beginning)
  //       | OutEdgeList | VertexList | Performance |
  //       |    listS    |   vecS     |     25.2    |
  //       |    vecS     |   vecS     |     23.4    |
  //
  //     * Didn't try adjacency_matrix since our graph is sparse
  //     ( Also one can check BGL book, performance section )
  //
  //    Decision:
  //     * Alpha_expansion_graph_cut_boost: use adjacency_list<vecS, listS> without
  //       pre-allocating vertex-list.
  //
  // 2) With Boykov-Kolmogorov MAXFLOW software:
  //   (http://pub.ist.ac.at/~vnk/software/maxflow-v2.21.src.tar.gz)
  //                                  | Performance |
  //                                  |     3.1     |
  //     * Alpha_expansion_graph_cut_boykov_kolmogorov provides an implementation.
  //       MAXFLOW does not provide any option for pre-allocation (It is possible with v_3.02 though).
  //
  // Typical Benchmark result provided by Ilker
  //                                 | construction of vertices  |  construction of edges    | graph cut  | Total
  //   -----------------------------------------------------------------------------------------------------------
  //   boost with an adjacency list  |         1.53              | 1.51                      |  3.00      | 6.04
  //   boost with CSR                | 0.11 (gather in a vector) | 0.15 (gather in a vector) |  2.67      | 2.93
  //   MaxFlow                       |       0.042               | 0.076                     |  1.043     | 1.161
  //
  // The main issue for now with CSR is the construction of the opposite edge map that is too costly,
  // since it is done by exploring all edges to find opposite
  ////////////////////////////////////////////////////////////////////////////////////////
  
  } // namespace internal
  
  /**
   * @brief Implements alpha-expansion graph cut algorithm.
   *
   * For representing graph, it uses adjacency_list with OutEdgeList = vecS, VertexList = listS.
   * Also no pre-allocation is made for vertex-list.
   */
  class Alpha_expansion_boost_adjacency_list_impl
  {
  private:
    typedef boost::adjacency_list_traits<boost::vecS, boost::listS, boost::directedS>
    Adjacency_list_traits;
  
    typedef boost::adjacency_list<boost::vecS, boost::listS, boost::directedS,
            // 4 vertex properties
            boost::property<boost::vertex_index_t, std::size_t,
            boost::property<boost::vertex_color_t, boost::default_color_type,
            boost::property<boost::vertex_distance_t, double,
            boost::property<boost::vertex_predecessor_t, Adjacency_list_traits::edge_descriptor >
            > > >,
            // 3 edge properties
            boost::property<boost::edge_capacity_t, double,
            boost::property<boost::edge_residual_capacity_t, double,
            boost::property<boost::edge_reverse_t, Adjacency_list_traits::edge_descriptor> >
            > > Graph;
  
    typedef boost::graph_traits<Graph> Traits;
    typedef boost::color_traits<boost::default_color_type> ColorTraits;
  
  public:
  
    typedef Traits::vertex_descriptor Vertex_descriptor;
    typedef Traits::vertex_iterator   Vertex_iterator;
    typedef Traits::edge_descriptor   Edge_descriptor;
    typedef Traits::edge_iterator     Edge_iterator;
  
  private:
  
    Graph graph;
    Vertex_descriptor cluster_source;
    Vertex_descriptor cluster_sink;
  
  public:
  
    void clear_graph()
    {
      graph.clear();
      cluster_source = boost::add_vertex(graph);
      cluster_sink = boost::add_vertex(graph);
    }
  
    Vertex_descriptor add_vertex()
    {
      return boost::add_vertex(graph);
    }
  
    void add_tweight (Vertex_descriptor& v, double w1, double w2)
    {
      add_edge (cluster_source, v, w1, 0);
      add_edge (v, cluster_sink, w2, 0);
    }
  
    void init_vertices()
    {
      // initialize vertex indices, it is necessary since we are using VertexList = listS
      Vertex_iterator v_begin, v_end;
      Traits::vertices_size_type index = 0;
      for(boost::tie(v_begin, v_end) = vertices(graph); v_begin != v_end; ++v_begin) {
        boost::put(boost::vertex_index, graph, *v_begin, index++);
      }
    }
  
    double max_flow()
    {
      return boost::boykov_kolmogorov_max_flow(graph, cluster_source,
                                                        cluster_sink);
    }
  
    template <typename VertexLabelMap, typename InputVertexDescriptor>
    void update(VertexLabelMap vertex_label_map,
                const std::vector<Vertex_descriptor>& inserted_vertices,
                InputVertexDescriptor vd,
                std::size_t vertex_i,
                std::size_t alpha)
    {
      boost::default_color_type color = boost::get(boost::vertex_color, graph,
                                                   inserted_vertices[vertex_i]);
      if(std::size_t(get (vertex_label_map, vd)) != alpha
         && color == ColorTraits::white()) //new comers (expansion occurs)
        put (vertex_label_map, vd,
             static_cast<typename boost::property_traits<VertexLabelMap>::value_type>(alpha));
    }
  
    void add_edge (Vertex_descriptor& v1, Vertex_descriptor& v2, double w1, double w2)
    {
      Edge_descriptor v1_v2, v2_v1;
      bool v1_v2_added, v2_v1_added;
  
      boost::tie(v1_v2, v1_v2_added) = boost::add_edge(v1, v2, graph);
      boost::tie(v2_v1, v2_v1_added) = boost::add_edge(v2, v1, graph);
  
      CGAL_assertion(v1_v2_added && v2_v1_added);
      //put edge capacities
      boost::put(boost::edge_reverse, graph, v1_v2, v2_v1);
      boost::put(boost::edge_reverse, graph, v2_v1, v1_v2);
  
      //map reverse edges
      boost::put(boost::edge_capacity, graph, v1_v2, w1);
      boost::put(boost::edge_capacity, graph, v2_v1, w2);
  
    }
  };
  
  // another implementation using compressed_sparse_row_graph
  // for now there is a performance problem while setting reverse edges
  // if that can be solved, it is faster than Alpha_expansion_graph_cut_boost
  class Alpha_expansion_boost_compressed_sparse_row_impl
  {
  private:
    // CSR only accepts bundled props
    struct VertexP {
      boost::default_color_type vertex_color;
      double vertex_distance_t;
      // ? do not now there is another way to take it, I think since edge_descriptor does not rely on properties
      // this should be fine...
      boost::compressed_sparse_row_graph<boost::directedS>::edge_descriptor
      vertex_predecessor;
    };
  
    struct EdgeP {
      double edge_capacity;
      double edge_residual_capacity;
      boost::compressed_sparse_row_graph<boost::directedS>::edge_descriptor
      edge_reverse;
    };
  
    typedef boost::compressed_sparse_row_graph<boost::directedS,
            VertexP, EdgeP> Graph;
  
    typedef boost::graph_traits<Graph> Traits;
    typedef boost::color_traits<boost::default_color_type> ColorTraits;
  
  public:
  
    typedef Traits::vertex_descriptor Vertex_descriptor;
    typedef Traits::vertex_iterator   Vertex_iterator;
    typedef Traits::edge_descriptor   Edge_descriptor;
    typedef Traits::edge_iterator     Edge_iterator;
  
  private:
  
    Graph graph;
    std::size_t nb_vertices;
    std::vector<std::pair<std::size_t, std::size_t> > edge_map;
    std::vector<EdgeP>                edge_map_weights;
  
  public:
    void clear_graph()
    {
      nb_vertices = 2;
      edge_map.clear();
      edge_map_weights.clear();
      // edge_map.reserve(labels.size() *
      //                  8); // there is no way to know exact edge count, it is a heuristic value
      // edge_map_weights.reserve(labels.size() * 8);
    }
  
    Vertex_descriptor add_vertex()
    {
      return (nb_vertices ++);
    }
  
    void add_tweight (Vertex_descriptor& v, double w1, double w2)
    {
      add_edge (0, v, w1, 0);
      add_edge (v, 1, w2, 0);
    }
  
    void init_vertices()
    {
      graph = Graph(boost::edges_are_unsorted, edge_map.begin(), edge_map.end(),
                    edge_map_weights.begin(), nb_vertices);
  
      // PERFORMANCE PROBLEM
      // need to set reverse edge map, I guess there is no way to do that before creating the graph
      // since we do not have edge_descs
      // however from our edge_map, we know that each (2i, 2i + 1) is reverse pairs, how to facilitate that ?
      // will look it back
      Graph::edge_iterator ei, ee;
      for(boost::tie(ei, ee) = boost::edges(graph); ei != ee; ++ei) {
        Graph::vertex_descriptor v1 = boost::source(*ei, graph);
        Graph::vertex_descriptor v2 = boost::target(*ei, graph);
        std::pair<Graph::edge_descriptor, bool> opp_edge = boost::edge(v2, v1, graph);
  
        CGAL_assertion(opp_edge.second);
        graph[opp_edge.first].edge_reverse =
          *ei; // and edge_reverse of *ei will be (or already have been) set by the opp_edge
      }
    }
  
    double max_flow()
    {
      // since properties are bundled, defaults does not work need to specify them
      return boost::boykov_kolmogorov_max_flow
        (graph,
         boost::get(&EdgeP::edge_capacity, graph),
         boost::get(&EdgeP::edge_residual_capacity, graph),
         boost::get(&EdgeP::edge_reverse, graph),
         boost::get(&VertexP::vertex_predecessor, graph),
         boost::get(&VertexP::vertex_color, graph),
         boost::get(&VertexP::vertex_distance_t, graph),
         boost::get(boost::vertex_index,
                    graph), // this is not bundled, get it from graph (CRS provides one)
         0, 1);
    }
  
    template <typename VertexLabelMap, typename InputVertexDescriptor>
    void update(VertexLabelMap vertex_label_map,
                const std::vector<Vertex_descriptor>&,
                InputVertexDescriptor vd,
                std::size_t vertex_i,
                std::size_t alpha)
    {
      boost::default_color_type color =  graph[vertex_i + 2].vertex_color;
      if(get(vertex_label_map, vd)!= alpha
         && color == ColorTraits::white()) //new comers (expansion occurs)
        put(vertex_label_map, vd, alpha);
    }
  
    void add_edge(Vertex_descriptor v1, Vertex_descriptor v2, double w1, double w2)
    {
      edge_map.push_back(std::make_pair(v1, v2));
      EdgeP p1;
      p1.edge_capacity = w1;
      edge_map_weights.push_back(p1);
  
      edge_map.push_back(std::make_pair(v2, v1));
      EdgeP p2;
      p2.edge_capacity = w2;
      edge_map_weights.push_back(p2);
    }
  
  
  };
  
  // tags
  struct Alpha_expansion_boost_adjacency_list_tag { };
  struct Alpha_expansion_boost_compressed_sparse_row_tag { };
  struct Alpha_expansion_MaxFlow_tag { };
  
  // forward declaration
  class Alpha_expansion_MaxFlow_impl;
  
  /// \endcond
  
  // NOTE: latest performances check (2019-07-22)
  //
  // Using a random graph with 50000 vertices, 100000 edges and 30 labels:
  //
  // METHOD                 TIMING           MEMORY
  // Boost Adjacency list   49s              122MiB
  // Boost CSR              187s             77MiB
  // MaxFlow                12s              717MiB
  
  /**
     \ingroup PkgBGLPartition
  
     regularizes a partition of a graph into `n` labels using the alpha
     expansion algorithm \cgalCite{Boykov2001FastApproximate}.
  
     For a graph \f$(V,E)\f$, this function computes a partition `f`
     that minimizes the following cost function:
  
     \f[
     \mathrm{C}(f) = \sum_{\{v0,v1\} \in E} C_E(v0,v1) + \sum_{v \in V} C_V(f_v)
     \f]
  
     where \f$C_E(v0,v1)\f$ is the edge cost of assigning a different
     label to \f$v0\f$ and \f$v1\f$, and \f$C_V(f_v)\f$ is the vertex
     cost of assigning the label \f$f\f$ to the vertex \f$v\f$.
  
     \tparam InputGraph a model of `VertexAndEdgeListGraph`
  
     \tparam EdgeCostMap a model of `ReadablePropertyMap` with
     `boost::graph_traits<InputGraph>::%edge_descriptor` as key and `double`
     as value
  
     \tparam VertexLabelCostMap a model of `ReadablePropertyMap`
     with `boost::graph_traits<InputGraph>::%vertex_descriptor` as key and
     `std::vector<double>` as value
  
     \tparam VertexLabelMap a model of `ReadWritePropertyMap` with
     `boost::graph_traits<InputGraph>::%vertex_descriptor` as key and
     `std::size_t` as value
  
     \tparam NamedParameters a sequence of named parameters
  
     \param input_graph the input graph.
  
     \param edge_cost_map a property map providing the weight of each
     edge.
  
     \param vertex_label_map a property map providing the label of each
     vertex. This map will be updated by the algorithm with the
     regularized version of the partition.
  
     \param vertex_label_cost_map a property map providing, for each
     vertex, an `std::vector` containing the cost of this vertex to
     belong to each label. Each `std::vector` should have the same size
     `n` (which is the number of labels), each label being indexed from
     `0` to `n-1`. For example, `get(vertex_label_cost_map,
     vd)[label_idx]` returns the cost of vertex `vd` to belong to the
     label `label_idx`.
  
     \param np an optional sequence of \ref bgl_namedparameters "Named Parameters" among the ones listed below
  
     \cgalNamedParamsBegin
       \cgalParamNBegin{vertex_index_map}
         \cgalParamDescription{a property map associating to each vertex of `input_graph` a unique index between `0` and `num_vertices(input_graph) - 1`}
         \cgalParamType{a class model of `ReadablePropertyMap` with `boost::graph_traits<InputGraph>::%vertex_descriptor`
                        as key type and `std::size_t` as value type}
         \cgalParamDefault{an automatically indexed internal map}
         \cgalParamExtra{If this parameter is not passed, internal machinery will create and initialize
                         a face index property map, either using the internal property map if it exists
                         or using an external map. The latter might result in  - slightly - worsened performance
                         in case of non-constant complexity for index access.}
       \cgalParamNEnd
  
       \cgalParamNBegin{implementation_tag}
         \cgalParamDescription{a tag used to select which implementation of the alpha expansion should be used.
                               Available implementation tags are:
                               - `CGAL::Alpha_expansion_boost_adjacency_list`
                               - `CGAL::Alpha_expansion_boost_compressed_sparse_row_tag`
                               - `CGAL::Alpha_expansion_MaxFlow_tag`}
         \cgalParamDefault{`CGAL::Alpha_expansion_boost_adjacency_list`}
       \cgalParamNEnd
     \cgalNamedParamsEnd
  
     \note The `MaxFlow` implementation is provided by the \ref PkgSurfaceMeshSegmentationRef
     under GPL license. The header `<CGAL/boost/graph/Alpha_expansion_MaxFlow_tag.h>`
     must be included if users want to use this implementation.
  */
  template <typename InputGraph,
            typename EdgeCostMap,
            typename VertexLabelCostMap,
            typename VertexLabelMap,
            typename NamedParameters>
  double alpha_expansion_graphcut (const InputGraph& input_graph,
                                   EdgeCostMap edge_cost_map,
                                   VertexLabelCostMap vertex_label_cost_map,
                                   VertexLabelMap vertex_label_map,
                                   const NamedParameters& np)
  {
    using parameters::choose_parameter;
    using parameters::get_parameter;
  
    typedef boost::graph_traits<InputGraph> GT;
    typedef typename GT::edge_descriptor input_edge_descriptor;
    typedef typename GT::vertex_descriptor input_vertex_descriptor;
  
    typedef typename GetInitializedVertexIndexMap<InputGraph, NamedParameters>::type VertexIndexMap;
    VertexIndexMap vertex_index_map = CGAL::get_initialized_vertex_index_map(input_graph, np);
  
    typedef typename GetImplementationTag<NamedParameters>::type Impl_tag;
  
    // select implementation
    typedef typename std::conditional
      <std::is_same<Impl_tag, Alpha_expansion_boost_adjacency_list_tag>::value,
       Alpha_expansion_boost_adjacency_list_impl,
       typename std::conditional
       <std::is_same<Impl_tag, Alpha_expansion_boost_compressed_sparse_row_tag>::value,
        Alpha_expansion_boost_compressed_sparse_row_impl,
        Alpha_expansion_MaxFlow_impl>::type>::type
      Alpha_expansion;
  
    typedef typename Alpha_expansion::Vertex_descriptor Vertex_descriptor;
  
    Alpha_expansion alpha_expansion;
  
    // TODO: check this hardcoded parameter
    const double tolerance = 1e-10;
  
    double min_cut = (std::numeric_limits<double>::max)();
  
  #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
    double vertex_creation_time, edge_creation_time, cut_time;
    vertex_creation_time = edge_creation_time = cut_time = 0.0;
  #endif
  
    std::vector<Vertex_descriptor> inserted_vertices;
    inserted_vertices.resize(num_vertices (input_graph));
  
    std::size_t number_of_labels = get(vertex_label_cost_map, *(vertices(input_graph).first)).size();
  
    bool success;
    do {
      success = false;
  
      for (std::size_t alpha = 0; alpha < number_of_labels; ++ alpha)
      {
        alpha_expansion.clear_graph();
  
  #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
        Timer timer;
        timer.start();
  #endif
  
        // For E-Data
        // add every input vertex as a vertex to the graph, put edges to source & sink vertices
        for (input_vertex_descriptor vd : CGAL::make_range(vertices(input_graph)))
        {
          std::size_t vertex_i = get(vertex_index_map, vd);
          Vertex_descriptor new_vertex = alpha_expansion.add_vertex();
          inserted_vertices[vertex_i] = new_vertex;
          double source_weight = get(vertex_label_cost_map, vd)[alpha];
          // since it is expansion move, current alpha labeled vertices will be assigned to alpha again,
          // making sink_weight 'infinity' guarantee this.
          double sink_weight = (std::size_t(get(vertex_label_map, vd)) == alpha ?
                                (std::numeric_limits<double>::max)()
                                : get(vertex_label_cost_map, vd)[get(vertex_label_map, vd)]);
  
          alpha_expansion.add_tweight(new_vertex, source_weight, sink_weight);
        }
  #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
        vertex_creation_time += timer.time();
        timer.reset();
  #endif
  
        // For E-Smooth
        // add edge between every vertex,
        for (input_edge_descriptor ed : CGAL::make_range(edges(input_graph)))
        {
          input_vertex_descriptor vd1 = source(ed, input_graph);
          input_vertex_descriptor vd2 = target(ed, input_graph);
          std::size_t idx1 = get (vertex_index_map, vd1);
          std::size_t idx2 = get (vertex_index_map, vd2);
  
          double weight = get (edge_cost_map, ed);
  
          Vertex_descriptor v1 = inserted_vertices[idx1],
            v2 = inserted_vertices[idx2];
  
          std::size_t label_1 = get (vertex_label_map, vd1);
          std::size_t label_2 = get (vertex_label_map, vd2);
          if(label_1 == label_2) {
            if(label_1 != alpha) {
              alpha_expansion.add_edge(v1, v2, weight, weight);
            }
          } else {
            Vertex_descriptor inbetween = alpha_expansion.add_vertex();
  
            double w1 = (label_1 == alpha) ? 0 : weight;
            double w2 = (label_2 == alpha) ? 0 : weight;
            alpha_expansion.add_edge(inbetween, v1, w1, w1);
            alpha_expansion.add_edge(inbetween, v2, w2, w2);
            alpha_expansion.add_tweight(inbetween, 0., weight);
          }
        }
  #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
        edge_creation_time += timer.time();
  #endif
  
        alpha_expansion.init_vertices();
  
  #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
        timer.reset();
  #endif
  
        double flow = alpha_expansion.max_flow();
  
  #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
        cut_time += timer.time();
  #endif
  
        if(min_cut - flow <= flow * tolerance) {
          continue;
        }
        min_cut = flow;
        success = true;
        //update labeling
        for (input_vertex_descriptor vd : CGAL::make_range(vertices (input_graph)))
        {
          std::size_t vertex_i = get (vertex_index_map, vd);
          alpha_expansion.update(vertex_label_map, inserted_vertices, vd, vertex_i, alpha);
        }
      }
    } while(success);
  
  #ifdef CGAL_SEGMENTATION_BENCH_GRAPHCUT
    CGAL_TRACE_STREAM << "vertex creation time: " << vertex_creation_time <<
      std::endl;
    CGAL_TRACE_STREAM << "edge creation time: " << edge_creation_time << std::endl;
    CGAL_TRACE_STREAM << "max flow algorithm time: " << cut_time << std::endl;
  #endif
  
    return min_cut;
  }
  
  
  /// \cond SKIP_IN_MANUAL
  // variant with default NP
  template <typename InputGraph,
            typename EdgeCostMap,
            typename VertexLabelCostMap,
            typename VertexLabelMap>
  double alpha_expansion_graphcut (const InputGraph& input_graph,
                                   EdgeCostMap edge_cost_map,
                                   VertexLabelCostMap vertex_label_cost_map,
                                   VertexLabelMap vertex_label_map)
  {
    return alpha_expansion_graphcut (input_graph, edge_cost_map,
                                     vertex_label_cost_map, vertex_label_map,
                                     CGAL::parameters::all_default());
  }
  
  // Old API
  inline double alpha_expansion_graphcut (const std::vector<std::pair<std::size_t, std::size_t> >& edges,
                                          const std::vector<double>& edge_costs,
                                          const std::vector<std::vector<double> >& cost_matrix,
                                          std::vector<std::size_t>& labels)
  {
    internal::Alpha_expansion_old_API_wrapper_graph graph (edges, edge_costs, cost_matrix, labels);
    return alpha_expansion_graphcut(graph,
                                    graph.edge_cost_map(),
                                    graph.vertex_label_cost_map(),
                                    graph.vertex_label_map(),
                                    CGAL::parameters::vertex_index_map (graph.vertex_index_map()));
  }
  
  template <typename AlphaExpansionImplementationTag>
  double alpha_expansion_graphcut (const std::vector<std::pair<std::size_t, std::size_t> >& edges,
                                   const std::vector<double>& edge_costs,
                                   const std::vector<std::vector<double> >& cost_matrix,
                                   std::vector<std::size_t>& labels,
                                   const AlphaExpansionImplementationTag&)
  {
    internal::Alpha_expansion_old_API_wrapper_graph graph (edges, edge_costs, cost_matrix, labels);
  
    return alpha_expansion_graphcut(graph,
                                    graph.edge_cost_map(),
                                    graph.vertex_label_cost_map(),
                                    graph.vertex_label_map(),
                                    CGAL::parameters::vertex_index_map (graph.vertex_index_map()).
                                    implementation_tag (AlphaExpansionImplementationTag()));
  }
  /// \endcond
  
  }//namespace CGAL
  
  namespace boost
  {
  
  template <>
  struct property_map<CGAL::internal::Alpha_expansion_old_API_wrapper_graph, boost::vertex_index_t>
  {
    typedef CGAL::internal::Alpha_expansion_old_API_wrapper_graph::Vertex_index_map type;
    typedef CGAL::internal::Alpha_expansion_old_API_wrapper_graph::Vertex_index_map const_type;
  };
  }
  
  #endif //CGAL_BOOST_GRAPH_ALPHA_EXPANSION_GRAPHCUT_H
  