imperativeCCNP.ml

open! DataStruct 
open DissimilarityLib

module IMap = Map.Make(Int) 



module Make = 
   functor (D : sig type t val d : t -> t -> int end) ->
   struct 

     let refine q set =
       let tree = ref IMap.empty in 
       for i = 0 to Array.length set - 1 do 
         let dist = D.d q set.(i) in 
         if IMap.mem dist !tree then 
           tree := IMap.add dist (set.(i) :: IMap.find dist !tree) !tree 
         else 
           tree := IMap.add dist [set.(i)] !tree 
       done;
       IMap.to_seq !tree 
       |> Seq.map snd 
       |> Array.of_seq


     let swap arr i j =
       let old_i = arr.(i) in 
       arr.(i) <- arr.(j);
       arr.(j) <- old_i


     let reverse_short_distance p q subsets = 
       let is_close subset = 
         subset <> [] && D.d q (List.hd subset) <= D.d p q
       in 
       let i = ref 0 in 
       while !i < Array.length subsets && is_close subsets.(!i) do incr i done;
       for j = 0 to !i/2 - 1 do 
         swap subsets j (!i-j-1)
       done


     let choose_pivot inners outers = 
       if inners = [] then 
         (List.hd outers, inners, List.tl outers, true)
       else 
         (List.hd inners, List.tl inners, outers, false)


     let rec recursive_refine p (inners, set, outers) =
       if (inners = [] && outers = []) || Array.length set <= 1 then 
         [set]
       else 
         let (q, ins, outs, q_is_in_outers) = choose_pivot inners outers in 
         let subsets = refine q set in 
         if q_is_in_outers then reverse_short_distance p q subsets;
         let concats = ref [] in 
         for i = Array.length subsets - 1 downto 0 do 
           let ins_i, outs_i = ref ins, ref outs in 
           for j = 0 to i-1 do 
             ins_i := List.rev_append subsets.(j) !ins_i;
           done;
           for j = i+1 to Array.length subsets - 1 do 
             outs_i := List.rev_append subsets.(j) !outs_i;
           done;
           let copoints = recursive_refine p (!ins_i, Array.of_list subsets.(i),!outs_i) in 
           concats := List.append copoints !concats;
         done;
         !concats

     exception Not_Robinson 

     
     let separate_if_separable p copoint =
       let l = Array.length copoint in  
       let x_min = copoint.(0) in 
       let x_max = copoint.(l-1) in 
       let delta = D.d p x_min in 
       let diam = D.d x_min x_max in 
       if diam <= delta then 
         [(x_min, copoint)]
       else 
         let i = ref 0 in 
         while not (
                   D.d x_min copoint.(!i) <= delta  
                   && D.d x_max copoint.(!i+1) <= delta 
                   && D.d copoint.(!i) copoint.(!i+1) >= delta
                 )
         do 
           incr i;
           if !i >= l - 1 then raise Not_Robinson
         done;
         [ (x_min, Array.sub copoint 0 (!i + 1));
           (x_max, Array.sub copoint (!i + 1) (l - !i - 1)) 
         ]
               

     type side = Left | Right 

     let append_all array list_ref =
       for i = 0 to Array.length array - 1 do 
         list_ref := array.(i) :: !list_ref
       done
       

     let sort_by_bipartition p set =
       let l = Array.length set in 
       let sides = Array.make l Right in 
       let lefts = ref [] in 
       let rights = ref [] in 
       let last_undecided = ref (l-1) in 
       let move_right i = 
         append_all (snd set.(i)) rights;
         sides.(i) <- Right;
         decr last_undecided
       in 
       let move_left i = 
         append_all (snd set.(i)) lefts;
         sides.(i) <- Left;
         decr last_undecided
       in
       for i = l-1 downto 0 do 
         let q = fst set.(i) in 
         if i = !last_undecided then move_right i;
         for j = !last_undecided downto 0 do 
           let x = fst set.(j) in 
           if not (D.d q x = D.d p q) then 
             if (D.d q x < D.d p q) = (sides.(i) = Right) then begin
                 for k = !last_undecided downto j+1 do move_left k done;
                 move_right j;
               end
             else begin
                 for k = !last_undecided downto j+1 do move_right k done;
                 move_left j
               end;
         done;
       done;
       List.rev_append !lefts (p :: !rights)


     let rec find_compatible_order set =
       let n = Array.length set in 
       if n <= 1 then set 
       else 
         let p = set.(0) in 
         let set' = Array.sub set 1 (n-1) in  
         let copoints = recursive_refine p ([p], set', []) |> Array.of_list in 
         let represented_copoints = ref [] in 
         for i = 0 to Array.length copoints - 1 do 
           let order_i = find_compatible_order copoints.(i) in 
           represented_copoints := List.rev_append (separate_if_separable p order_i) !represented_copoints;
         done;
         sort_by_bipartition p (Array.of_list (List.rev !represented_copoints))
         |> Array.of_list
   end


let find_compatible_order (type elt) diss = 
  let open Dissimilarity in 
  let module D = struct type t = elt let d = diss.d end in 
  let module M = Make(D) in 
  try 
    let order = M.find_compatible_order (Array.of_list diss.elements) in 
    if Dissimilarity.is_compatible_order diss (Array.to_list order) then Some order
    else None 
  with 
    | M.Not_Robinson -> None