(**************************************************************************)
(*                                                                        *)
(*  Copyright (C) 2012 Johannes 'josch' Schauer <j.schauer@email.de>      *)
(*  Copyright (C) 2012 Pietro Abate <pietro.abate@pps.jussieu.fr>         *)
(*                                                                        *)
(*  This library is free software: you can redistribute it and/or modify  *)
(*  it under the terms of the GNU Lesser General Public License as        *)
(*  published by the Free Software Foundation, either version 3 of the    *)
(*  License, or (at your option) any later version.  A special linking    *)
(*  exception to the GNU Lesser General Public License applies to this    *)
(*  library, see the COPYING file for more information.                   *)
(**************************************************************************)

module OCAMLHashtbl = Hashtbl

open! ExtLib
open Dose_common
open Dose_algo

#define __label __FILE__
let label =  __label ;;
include Util.Logging(struct let label = label end) ;;

let timer_johnson = Util.Timer.create "GraphUtils.FindCycles.johnson"

module GraphUtils (G : Graph.Sig.I) = struct
  module C = Graph.Components.Make(G)
  module VSet = Set.Make(G.V)
  module O = Defaultgraphs.GraphOper(G)

  let v_list_mem v l = List.exists (fun p -> (G.V.compare p v) = 0) l;;

  let scc_with_vertex g v =
    let _,scc = C.scc g in
    let n = scc v in
    let sg = G.create () in
    G.iter_vertex (fun v1 ->
      if (scc v1) = n then
        List.iter (fun e ->
          if scc (G.E.dst e) = n then
            G.add_edge_e sg e
        ) (G.succ_e g v1)
    ) g;
    sg
  ;;

  let copy_graph g =
    let g1 = G.create () in
    G.iter_edges_e (fun e -> G.add_edge_e g1 e) g;
    g1
  ;;

  let find_strong_articulation_points g =
    let scc = List.filter_map (function [] | [_] -> None | s -> Some(O.subgraph g s)) (C.scc_list g) in
    List.fold_left (fun acc scc ->
      (* test if this scc has a strong articulation point *)
      let l = G.fold_vertex (fun v acc ->
        (* test if the scc without this vertex is split or not *)
        let og = copy_graph g in
        G.remove_vertex og v;
        let num_scc = List.length (C.scc_list og) in
        if num_scc = 1 then acc else (v,num_scc)::acc
      ) scc [] in
      l@acc
    ) [] scc
  ;;

  let find_strong_bridges g =
    let scc = List.filter_map (function [] | [_] -> None | s -> Some(O.subgraph g s)) (C.scc_list g) in
    List.fold_left (fun acc scc ->
      (* test if this scc has a strong articulation point *)
      let l = G.fold_edges_e (fun e acc ->
        (* test if the scc without this vertex is split or not *)
        let og = copy_graph g in
        G.remove_edge_e og e;
        let num_scc = List.length (C.scc_list og) in
        if num_scc = 1 then acc else (e,num_scc)::acc
      ) scc [] in
      l@acc
    ) [] scc
  ;;

  let get_partial_order g =
    let module Dfs = Graph.Traverse.Dfs(G) in
    if Dfs.has_cycle g then
      failwith "need a DAG without cycles as input for partial order";
    let module Hashtbl = OCAMLHashtbl.Make(G.V) in
    (* find all vertices with no successors (those that have all dependencies
     * fulfilled) *)
    let init = G.fold_vertex (fun v acc ->
      match G.succ g v with
        | [] -> v::acc
        | _ -> acc
    ) g [] in
    let result = ref [init] in
    let processed = Hashtbl.create (G.nb_vertex g) in
    let tocheck = Hashtbl.create (G.nb_vertex g) in
    (* fill the two hashtables, starting from the initial result *)
    List.iter (fun v ->
      Hashtbl.replace processed v ();
      G.iter_pred (fun pred ->
        Hashtbl.replace tocheck pred ()
      ) g v
    ) init;
    while Hashtbl.length tocheck > 0 do
      let localprocessed = Hashtbl.create (G.nb_vertex g) in
      (* iterate over the to-be-checked vertices and check if all of their
       * dependencies are already in the result *)
      Hashtbl.iter (fun v _ ->
        let satisfied = G.fold_succ (fun succ acc ->
          if acc then Hashtbl.mem processed succ else acc
        ) g v true in
        (* if yes, remove this vertex from tocheck, add it to the result and 
         * add its predecessors to tocheck*)
        if satisfied then begin
          Hashtbl.remove tocheck v;
          Hashtbl.replace localprocessed v ();
          G.iter_pred (fun pred ->
            Hashtbl.replace tocheck pred ()
          ) g v;
        end;
      ) tocheck;
      (* add results of this round to global variables *)
      let l = Hashtbl.fold (fun v _ acc ->
        Hashtbl.replace processed v ();
        v::acc
      ) localprocessed [] in
      result := l::!result
    done;
    List.rev !result
  ;;

  exception Found of G.V.t
  let find_vertex_option prop g =
    try
      G.iter_vertex (fun v -> if prop v then raise (Found v)) g;
      None
    with Found v -> Some(v)
  ;;

  exception Empty
  let first_vertex g =
    try
      G.iter_vertex (fun v -> raise (Found v)) g;
      raise Empty
    with Found v -> v
  ;;
end

module FindCycles (G : Graph.Sig.I) = struct
  module C = Graph.Components.Make(G)
  module GraphUtilsG = GraphUtils(G)
  module O = Defaultgraphs.GraphOper(G)

  let edge_cycle_from_vertex_cycle g cycle =
    let fst, tail = match cycle with
      | h::t -> h,t
      | [] -> raise List.Empty_list
    in
    let rec get_edges vertices last acc =
      match vertices with
        | [] ->
            (* edge from last to fst to close the cycle *)
            let e = G.find_edge g last fst in
            e::acc
        | h::t ->
            (*edge from last to h *)
            let e = G.find_edge g last h in
            get_edges t h (e::acc)
    in
    let edges = get_edges tail fst [] in
    List.rev edges
  ;;

  let get_cycles_per_edge_ht g cycles =
    let hist = Hashtbl.create (G.nb_edges g) in
    List.iter (fun cycle ->
      let edges = edge_cycle_from_vertex_cycle g cycle in
      List.iter (fun edge ->
        Hashtbl.replace hist edge ((Hashtbl.find_default hist edge 0)+1);
      ) edges;
    ) cycles;
    hist
  ;;

  let to_edges g cycle =
  let rec aux (acc,first,last) = function
    |[] ->
        begin match (acc,first,last) with
        |acc,Some f,Some l -> List.rev ((G.find_edge g l f) :: acc)
        |[],None,None -> []
        |_,_,_ -> failwith "not a cycle (empty)" end
    |h1::h2::t ->
        begin match (acc,first,last) with
        |[],None,None ->
            let e = G.find_edge g h1 h2 in
            aux (e::acc,Some h1,Some h2) t
        |acc,Some f,Some l ->
            let e1 = G.find_edge g l h1 in
            let e2 = G.find_edge g h1 h2 in
            aux (e2::e1::acc,Some f,Some h2) t
        |_,_,_ -> failwith "not a cycle (full)" end
    |_ -> failwith "not a cycle (impossible)"
  in aux ([],None,None) cycle


  let johnson ?(maxlength=0) ?(maxamount=0) g =
    Util.Timer.start timer_johnson;

    let path = Stack.create () in
      (* vertex: blocked from search *)
    let blocked = Hashtbl.create 1023 in
      (* graph portions that yield no elementary circuit *)
    let b = Hashtbl.create 1023 in
      (* list to accumulate the circuits found  *)
    let result = ref [] in
    let resultlen = ref 0 in
    let i = ref 0 in

    let rec unblock n =
      if Hashtbl.find blocked n then begin
        Hashtbl.replace blocked n false;
        List.iter unblock (Hashtbl.find b n);
        Hashtbl.replace b n [];
      end
    in

    let stack_to_list s =
      let l = ref [] in
      Stack.iter (fun e -> l:= e::!l) s;
      !l
    in

    let vertex_set = G.fold_vertex (fun v l -> v::l) g [] in
    let vertex_set_len = List.length vertex_set in

    let rec circuit thisnode startnode component =
      let closed = ref false in
      if (maxlength = 0 || (Stack.length path) < maxlength) && (maxamount = 0 || !resultlen < maxamount) then begin
        Stack.push thisnode path;
        Hashtbl.replace blocked thisnode true;
        G.iter_succ (fun nextnode ->
          if G.V.equal nextnode startnode then begin
            result := ((stack_to_list path))::!result;
            incr resultlen;
            debug "[%d/%d] found cycles: %d\r%!" !i vertex_set_len !resultlen;
            closed := true;
          end else begin if not(Hashtbl.find blocked nextnode) then
            if circuit nextnode startnode component then begin
              closed := true;
            end
          end
        ) component thisnode;
        if !closed then begin
          unblock thisnode
        end
        else
          G.iter_succ (fun nextnode ->
            let l = Hashtbl.find b nextnode in
            if not(List.mem thisnode l) then
              Hashtbl.replace b nextnode (thisnode::l)
          ) component thisnode;
        ignore(Stack.pop path);
      end;
      !closed
    in
    debug "[%d/%d] found cycles: %d\r%!" 0 vertex_set_len 0;

    let non_degenerate_scc = List.filter (function [] | [_] -> incr i; false | _ -> true) (C.scc_list g) in

    let non_degenerate_subgraphs = List.map (fun scc -> O.subgraph g scc) non_degenerate_scc in

    List.iter (fun g ->
      let vertex_set = G.fold_vertex (fun v l -> v::l) g [] in
      List.iter (fun s ->
        incr i;
        debug "[%d/%d] found cycles: %d\r%!" !i vertex_set_len !resultlen;
        if maxamount = 0 || !resultlen < maxamount then begin
          let component = GraphUtilsG.scc_with_vertex g s in
          if G.nb_edges component > 0 then begin
            G.iter_vertex (fun node ->
              Hashtbl.replace blocked node false;
              Hashtbl.replace b node [];
            ) component;
            ignore(circuit s s component);
          end;
          G.remove_vertex g s;
        end
      ) (List.sort ~cmp:G.V.compare vertex_set);
    ) non_degenerate_subgraphs;

    debug "[%d/%d] found cycles: %d\n%!" !i vertex_set_len !resultlen;
    Util.Timer.stop timer_johnson (List.rev !result)
  ;;
end

module FAS (G : Graph.Sig.I) = struct
  module EdgeSet = Set.Make(G.E)
  module FindCyclesG = FindCycles(G)
  module Dfs = Graph.Traverse.Dfs(G)
  module GraphUtilsG = GraphUtils(G)
  module Int = struct type t = int let compare = Pervasives.compare end
  module IntSet = Set.Make(Int)
  module T = Graph.Topological.Make(G)

  (* decide whether a given feedback arc set is minimal for the supplied graph.
   * A feedback arc set is minimal if the reinsertion of any edge in the set
   * into the graph would make the graph cyclic again
   * this is FALSE but written as such in:
   *   Brandenburg, K. H. F., & Hanauer, K. (2011). Sorting Heuristics for the
   *   Feedback Arc Set Problem. *)
  let isminimal fas g =
    let g = GraphUtilsG.copy_graph g in
    EdgeSet.iter (G.remove_edge_e g) fas;
    if Dfs.has_cycle g then
      failwith "the input is not a feedback arc set for the supplied graph";
    EdgeSet.for_all (fun e ->
      G.add_edge_e g e;
      let res = Dfs.has_cycle g in
      G.remove_edge_e g e;
      res
    ) fas
  ;;

  (* turn a feedback arc set into a vertex ordering 
   * since there are many topological orderings for a given acyclic graph, take
   * care to choose the order which keeps the feedback arc set small 
   * for this reason, instead of using Graph.Topological,
   * GraphUtils.get_partial_order is used. All vertices within each group
   * returned by this function are then ordered such that the cardinality of the
   * given feedback arc set is reduced. *)
  let getorder fas g =
    (* fasverts is a hashtable which maps vertices which are the source of edges
     * in the feedback arc set to a list of vertices which are destinations of
     * edges in the feedback arc set. *)
    let fasverts = Hashtbl.create (EdgeSet.cardinal fas) in
    EdgeSet.iter (fun e ->
      let src = G.E.src e in let dst = G.E.dst e in
      if src <> dst then (* ignore selfcycles as they don't influence the order *)
        Hashtbl.add fasverts src (dst,e)
    ) fas;
    (* go through all vertex lists returned by GraphUtilsG.get_partial_order and
     * sort all vertices in this list which make edges in the feedback arc set,
     * such that those edges are removed *)
    List.fold_left (fun acc l ->
      (* get all the edges that are part of this list for lookup later *)
      let localverts = Hashtbl.create (List.length l) in
      List.iter (fun v -> Hashtbl.add localverts v ()) l;
      (* create a graph that only contains those feedback arcs whose source and
       * destination are in the current vertex list. The resulting graph might
       * be cyclic if the feedback arc set was really bad but at this point we
       * don't care*)
      let g = G.create () in
      List.iter (fun v1 ->
        List.iter (fun (v2,e) ->
          if Hashtbl.mem localverts v2 then G.add_edge_e g e
        ) (Hashtbl.find_all fasverts v1)
      ) l;
      if Dfs.has_cycle g then
        warning "fas has forward and backward edge (creating a cycle) we don't handle this yet";
      (* get the topological order of the vertices in the graph *)
      let vlist = T.fold (fun v acc -> v::acc) g [] in
      (* concatenate all vertices that are not part of the graph above *)
      let acc = List.fold_left (fun acc v ->
        if not (List.mem v vlist) then v::acc else acc
      ) acc l in
      List.rev_append vlist acc
    ) [] (GraphUtilsG.get_partial_order g)
  ;;

  let ordertofas order g =
    (* check if the order can match the graph *)
    if (List.length order) <> (G.nb_vertex g) then
      failwith "invalid vertex order (length differs)";
    let seen = Hashtbl.create (List.length order) in
    List.fold_left (fun acc v ->
      Hashtbl.add seen v (); (* because of edges in self cycles *)
      G.fold_succ_e (fun edge acc ->
        if Hashtbl.mem seen (G.E.dst edge) then
          EdgeSet.add edge acc (* this vertex has already been processed, so it is a backarc*)
        else
          acc
      ) g v acc
    ) EdgeSet.empty order
  ;;

  (* cost function for an order equals the size of the fas *)
  let costoforder order g =
    EdgeSet.cardinal (ordertofas order g)
  ;;

  let cyclefas ?(maxlength=4) g =
    let g = GraphUtilsG.copy_graph g in
    (* given a graph and a list of cycles in it, return a set of edges that
     * remove all those cycles by iteratively removing the edge that is shared
     * by most cycles. *)
    let calculate_partial_fas cycles =
      let hist = Hashtbl.create (G.nb_edges g) in
      (* create a hashtable mapping edges to a set of integers where each
       * integer maps to the cycle that this edge is part of *)
      List.iteri (fun i cycle ->
        let edges = FindCyclesG.edge_cycle_from_vertex_cycle g cycle in
        List.iter (fun edge ->
          Hashtbl.replace hist edge (IntSet.add i (Hashtbl.find_default hist edge IntSet.empty));
        ) edges;
      ) cycles;
      let rec remove_most_popular_edge acc =
        (* get the edge that is part of the most cycles *)
        match List.of_enum (Hashtbl.enum hist) with
          | [] -> acc (* it might be that no cycles of this length can be broken *)
          | hd::tl -> begin
              let max_edge,cids = List.fold_left (fun (k1,v1) (k2,v2) ->
                if (IntSet.cardinal v1) > (IntSet.cardinal v2) then k1,v1 else k2,v2
              ) hd tl in
              (* end if the edge with the most cycles has zero cycles *)
              if (IntSet.cardinal cids) = 0 then
                acc
              else begin
                (* remove those cycle ids from all sets *)
                Hashtbl.iter (fun edge set ->
                  Hashtbl.replace hist edge (IntSet.diff set cids)
                ) hist;
                (* add edge to feedback arc set *)
                remove_most_popular_edge (EdgeSet.add max_edge acc)
              end
            end
      in
      remove_most_popular_edge EdgeSet.empty
    in
    let remove_edgeset = EdgeSet.iter (G.remove_edge_e g) in
    (* find selfcycles *)
    let selfcycles = G.fold_edges_e (fun e acc ->
      if (G.E.src e) = (G.E.dst e) then EdgeSet.add e acc else acc
    ) g EdgeSet.empty in
    (* remove the found edges from the graph *)
    remove_edgeset selfcycles;
    (* apply calculate_partial_fas on the graph, remove the resulting edges and
     * increment the max cycle length each time until the graph is loop free *)
    let rec foo ml acc =
      if Dfs.has_cycle g then begin
        let cycles = FindCyclesG.johnson ~maxlength:ml g in
        match cycles with
          | [] -> foo (ml+2) acc
          | l ->
              let partial_fas = calculate_partial_fas l in
              remove_edgeset partial_fas;
              foo (ml+2) (EdgeSet.union partial_fas acc)
      end else
        acc
    in
    let fas = foo maxlength EdgeSet.empty in
    (* instead of calculating the union of the set of selfcycle edges and the
     * edges in the calculated feedback arc set, calculate the induced vertex
     * order. The back arcs in this order are of either equal or less amount
     * than the edges picked above. *)
    getorder fas g
  ;;

  let sortingfaswrapper algo g_orig =
    let g = GraphUtilsG.copy_graph g_orig in

    (* first, remove all selfcycles *)
    let selfcycles = G.fold_edges_e (fun e acc ->
      if (G.E.src e) = (G.E.dst e) then EdgeSet.add e acc else acc
    ) g EdgeSet.empty in
    (* remove the found edges from the graph *)
    EdgeSet.iter (G.remove_edge_e g) selfcycles;

    algo g
  ;;

  (* implementation of algorithm presented in
   *   Peter Eades, Xuemin Lin, W.F. Smyth, A fast and effective heuristic for
   *   the feedback arc set problem, Information Processing Letters, Volume 47,
   *   Issue 6, 18 October 1993, Pages 319-323, ISSN 0020-0190,
   *   10.1016/0020-0190(93)90079-O.
   *
   * by setting improved:true it becomes the implementation of an improvement
   * of the eades algo as presented in
   *   Tom Coleman and Anthony Wirth. 2010. Ranking tournaments: Local search
   *   and a new algorithm. J. Exp. Algorithmics 14, Article 6 (January 2010),
   *   .62 pages. DOI=10.1145/1498698.1537601
   *
   * the implementation is rather slow as focus was put on correctness and
   * readability rather than execution speed
   * *)
  let eadesfas ?(improved=false) g =
    let aux g =
      (* instead of appending to s1, we will prepend to s1 and reverse s1 in the
       * end*)
      let s1 = ref [] in
      let s2 = ref [] in
  
      let is_sink v = G.out_degree g v = 0 in
      let is_source v = G.in_degree g v = 0 in
      let delta = if improved then begin
        fun v -> abs ((G.out_degree g v) - (G.in_degree g v))
      end else begin
        fun v -> (G.out_degree g v) - (G.in_degree g v)
      end in
  
      while G.nb_vertex g > 0 do
        (* add all sinks to s2 *)
        try
          while true do
            match GraphUtilsG.find_vertex_option is_sink g with
              | Some v -> begin
                  s2 := v::(!s2);
                  G.remove_vertex g v;
                end
              | None -> raise Exit
          done
        with Exit -> ();
        (* add all sources to s1 *)
        try
          while true do
            match GraphUtilsG.find_vertex_option is_source g with
              | Some v -> begin
                  s1 := v::(!s1);
                  G.remove_vertex g v;
                end
              | None -> raise Exit
          done
        with Exit -> ();
        (* add vertex with maximum δ = d_{out} - d_{in} to s1 *)
        if G.nb_vertex g > 0 then begin
          let v = GraphUtilsG.first_vertex g in
          let v,_ = G.fold_vertex (fun v (vmax, d) ->
            let dn = delta v in
            if dn > d then (v,dn) else (vmax,d)
          ) g (v, delta v) in
          if improved && (G.in_degree g v) > (G.out_degree g v) then
            s2 := v::(!s2) (* treat as sink *)
          else
            s1 := v::(!s1); (* treat as source *)
          G.remove_vertex g v;
        end
      done;
      List.rev_append !s1 !s2
    in
    (* now figure out the feedback arc set from this ordering *)
    sortingfaswrapper aux g
  ;;

  (* cost function for an order equals the size of the fas *)
  let costoforder2l l e pos g =
    let seen = Hashtbl.create ((List.length l) + 1) in
    let res = ref EdgeSet.empty in (* could be speeded up by using a list instead of a set *)
    let gather = G.iter_succ_e (fun edge ->
      if Hashtbl.mem seen (G.E.dst edge) then
        res := EdgeSet.add edge !res (* this vertex has already been processed, so it is a backarc*)
    ) g in
    List.iteri (fun i v ->
      (* once we are at position i, add e in *)
      if i = pos then begin
        Hashtbl.add seen e ();
        gather e;
      end;
      Hashtbl.add seen v (); (* because of edges in self cycles *)
      gather v;
    ) l;
    (* if e was to be added at the end: do it now *)
    if pos = (List.length l) then begin
      Hashtbl.add seen e ();
      gather e;
    end;
    EdgeSet.cardinal !res
  ;;

  let handlevertexorder g = function
    | None -> G.fold_vertex (fun v acc -> v::acc) g []
    | Some l ->
        let l = List.filter (G.mem_vertex g) l in
        if (List.length l) <> (G.nb_vertex g) then
          failwith "invalid vertex order (length differs)"
        else l
  ;;


  let sortingfas ?(order=None) algo g =
    let aux g =
      let l = handlevertexorder g order in
      algo l g
    in
    sortingfaswrapper aux g
  ;;

  (* find the position to insert v into l with the least cost *)
  let insert g v = function
    | [] -> [v]
    | l -> begin
        let min = ref max_int in
        let res = ref 0 in
        for i=0 to (List.length l) do
          let m = costoforder2l l v i g in
          if m < !min then begin
            min := m;
            res := i;
          end
        done;
        match !res with
          | 0 -> v::l
          | _ -> begin
              let l1, l2 = List.split_nth !res l in
              l1@[v]@l2
            end
      end
  ;;


  (* TODO: move this further down to the insertsort section *)
  (* this version of insertion sort has an order as input and output
   * it does not care about selfcycles *)
  let insertorderfas order g =
    (* find the position to insert v into l with the least cost *)
    let rec sort l =
      match l with
        | h::tl -> insert g h (sort tl)
        | [] -> []
    in

    sort order
  ;;

  let rec convergefas g algo res min =
    let res = algo res g in
    let newmin = costoforder res g in
    if newmin = min then
      res (* convergence! *)
    else if newmin < min then begin
      convergefas g algo res newmin (* try harder! *)
    end else
      failwith "minimum increased" (* apparently algo is not monotone... *)
  ;;

  let revfas ?(order=None) algo g =
    let aux g =
      let res = handlevertexorder g order in
  
      (* then use algo until convergence is reached *)
      let c1 = convergefas g algo res max_int in
  
      (* now use this solution as an input to an algorithm which repeatedly
       * reverses that solution, applies one step of insertion sort and then
       * applies algo until convergence *)
      let rec converge2 res min =
        let res = List.rev res in
        (* after a reversal MUST come an insertionsort - otherwise an improvement
         * of the result is not guaranteed *)
        let res = insertorderfas res g in
        let newmin = costoforder res g in
        let res = convergefas g algo res newmin in
        let newmin = costoforder res g in
        if newmin = min then
          res (* convergence! *)
        else if newmin < min then begin
          debug "current best result: %d\n%!" newmin;
          converge2 res newmin
        end else
          failwith "minimum increased"
      in
      converge2 c1 max_int
    in
    sortingfaswrapper aux g
  ;;

  let multifas ?(order=None) algo g =
    let aux g =
      let res = handlevertexorder g order in
      convergefas g algo res max_int
    in
    sortingfaswrapper aux g
  ;;

  (************************ INSERT HEURISTIC ************************)

  (*
   * simplified version (no reversal, just iterative insertion) of
   *   Chanas, S., & Kobylański, P. (1996). A new heuristic algorithm solving
   *   the linear ordering problem. Computational optimization and applications,
   *   6(2), 191-205.
   *
   * everything is very imperative but this way of doing it outperformed any
   * recursive solution
   *
   * the complexity is O(n^2) as usual for insertion sort
   *
   * runs insertion sort on a default order
   * *)

  let insertfas g =
    sortingfas insertorderfas g
  ;;

  let insertmultifas g =
    multifas insertorderfas g
  ;;

  let insertrevfas g =
    revfas insertorderfas g
  ;;

  (************************ SIFTING HEURISTIC ************************)

  (* TODO: find a paper with the origin of sifting 
   * the term sifting comes from the paper
   *   Brandenburg, K. H. F., & Hanauer, K. (2011). Sorting Heuristics for the
   *   Feedback Arc Set Problem.
   * *)
  let siftorderfas order g =
    (* iterate over all vertices of the graph and try to find the best position
     * in the list for each vertex *)
    G.fold_vertex (fun v acc ->
      (* remove v from the list - TODO: work around this somehow...*)
      let acc = List.filter (fun e -> e <> v) acc in
      insert g v acc
    ) g order
  ;;

  let siftfas g =
    sortingfas siftorderfas g
  ;;

  let siftmultifas g =
    multifas siftorderfas g
  ;;

  (*
   * this combined algorithm was introduced in
   *   Brandenburg, K. H. F., & Hanauer, K. (2011). Sorting Heuristics for the
   *   Feedback Arc Set Problem.
   *
   * similar to insertrevfas with the difference of
   *  - using sifting instead of insertion sort
   *  - applying one step of insertion sort after reversal
   *      o this is necessary because sifting a reversed order does not
   *        necessarily improve the result but insertion sorting a reversed 
   *        order does
   *      o that reasoning can also be read about in
   *          Brandenburg, K. H. F., & Hanauer, K. (2011). Sorting Heuristics
   *          for the Feedback Arc Set Problem.
   *)
  let siftrevfas g =
    revfas siftorderfas g
  ;;

  (************************ MOVE HEURISTIC ************************)

  (*
   * the move heuristic comes from
   *   Coleman, T., & Wirth, A. (2009). Ranking tournaments: Local search and a
   *   new algorithm. Journal of Experimental Algorithmics (JEA), 14, 6.
   * and is there called chanas-both because of it's similarity to the insertion
   * sort algorithm above. It differs in that it allows each element to be moved
   * to either side. It was called move in
   *   Brandenburg, K. H. F., & Hanauer, K. (2011). Sorting Heuristics for the
   *   Feedback Arc Set Problem.
   * *)
  let moveorderfas order g =
    (* iterate over all positions in the order and try to find the best position
     * for the vertex in that postion *)
    let rec aux acc i =
      if List.length acc = i then acc else begin
        let l1, v, l2 = match List.split_nth i acc with
          | _, [] -> failwith "impossible"
          | l1, h::tl -> l1,h,tl
        in
        aux (insert g v (l1@l2)) (i+1)
      end
    in
    aux order 0
  ;;

  let movefas g =
    sortingfas moveorderfas g
  ;;

  let movemultifas g =
    multifas moveorderfas g
  ;;

  let moverevfas g =
    revfas moveorderfas g


  (************************ CYCLE HYBRIDS ************************)
  let cycleinsert ?(maxlength=4) g =
    let order = cyclefas ~maxlength g in
    sortingfas ~order:(Some order) insertorderfas g
  ;;
  let cycleinsertmulti ?(maxlength=4) g =
    let order = cyclefas ~maxlength g in
    multifas ~order:(Some order) insertorderfas g
  ;;
  let cycleinsertrev ?(maxlength=4) g =
    let order = cyclefas ~maxlength g in
    revfas ~order:(Some order) insertorderfas g
  ;;
  let cyclesift ?(maxlength=4) g =
    let order = cyclefas ~maxlength g in
    sortingfas ~order:(Some order) siftorderfas g
  ;;
  let cyclesiftmulti ?(maxlength=4) g =
    let order = cyclefas ~maxlength g in
    multifas ~order:(Some order) siftorderfas g
  ;;
  let cyclesiftrev ?(maxlength=4) g =
    let order = cyclefas ~maxlength g in
    revfas ~order:(Some order) siftorderfas g
  ;;
  let cyclemove ?(maxlength=4) g =
    let order = cyclefas ~maxlength g in
    sortingfas ~order:(Some order) moveorderfas g
  ;;
  let cyclemovemulti ?(maxlength=4) g =
    let order = cyclefas ~maxlength g in
    multifas ~order:(Some order) moveorderfas g
  ;;
  let cyclemoverev ?(maxlength=4) g =
    let order = cyclefas ~maxlength g in
    revfas ~order:(Some order) moveorderfas g
  ;;

  (************************ EADES HYBRIDS ************************)
  let eadesinsert g =
    let order = eadesfas g in
    sortingfas ~order:(Some order) insertorderfas g
  ;;
  let eadesinsertmulti g =
    let order = eadesfas g in
    multifas ~order:(Some order) insertorderfas g
  ;;
  let eadesinsertrev g =
    let order = eadesfas g in
    revfas ~order:(Some order) insertorderfas g
  ;;
  let eadessift g =
    let order = eadesfas g in
    sortingfas ~order:(Some order) siftorderfas g
  ;;
  let eadessiftmulti g =
    let order = eadesfas g in
    multifas ~order:(Some order) siftorderfas g
  ;;
  let eadessiftrev g =
    let order = eadesfas g in
    revfas ~order:(Some order) siftorderfas g
  ;;
  let eadesmove g =
    let order = eadesfas g in
    sortingfas ~order:(Some order) moveorderfas g
  ;;
  let eadesmovemulti g =
    let order = eadesfas g in
    multifas ~order:(Some order) moveorderfas g
  ;;
  let eadesmoverev g =
    let order = eadesfas g in
    revfas ~order:(Some order) moveorderfas g
  ;;
end