(** PROSITE-like patterns in sequence, whose elements come from a sub-logic. *)

module type PARAM =
  sig
    val toks_match : Syntax.t_list  (** Syntaxtic representation of the 'match' prefix symbol. *)
    val toks_begin : Syntax.t_list  (** Syntactic representation of the 'begin' symbol. *)
    val toks_end : Syntax.t_list  (** Syntactic representation of the 'end' symbol. *)
    val toks_sep : Syntax.t_list  (** Syntactic representation of the element separator symbol. *)
    val toks_any : Syntax.t_list  (** Syntactic representation of the 'any element' symbol. *)
    val feat_length : bool  (** Whether to generate motif lengths as features. *)
    val feat_dist : bool  (** Whether to generate distribution of elements as features. *)
    val feat_atoms : bool  (** Whether to generate atom featues. *)
    val gen_floating : bool  (** True if 'begin' and 'end' are not compulsory in motifs generated by [gen]. *)
    val gen_len : bool  (** Whether to allow [gen] to generate features about the lengths. *) 
    val gen_link_cond : bool -> int * int -> bool -> bool  (** Condition to be verified on gaps/links in motifs: the 1st/3rd argument respectively tells whether the left/right element is begin/end. *)
    val gen_min_xs : int  (** Threshold number of generated motifs by [gen]. *)
  end

module Make (Param : PARAM) (A : Logic.T) =
  (* A.gen must satisfy: A.gen(f,g,_) = {x}, and A.gen(f,g,A.gen(f,g)) is empty *)
struct
include Logic.Default

module IntervParam = struct let verbose = true end

module Len = Openinterval.DotDot(Intpow.Make)
module Dist = Openinterval.DotDot(Floatpow.Make)

type e = int * int
  (* interval for elastic segments *)
type a = Begin | Atom of A.t | Elastic of e | End
type t = Match of (bool * a list * bool * string)
       | Len of Len.t
       | Dist of (A.t * Dist.t) option

(* private functions *)

let print_debug s =
  print_endline s;
  flush stderr

let a_entails a1 a2 =
  match a1, a2 with
  | Begin, Begin -> true
  | End, End -> true
  | Atom x1, Atom x2 -> A.entails x1 x2
  | _, _ -> false

let e_contains (min1,max1) (min2,max2) = min1 <= min2 & max2 <= max1

let e_append (min1,max1) (min2,max2) = (min1+min2, max1+max2)

let e_union (min1,max1) (min2,max2) = (min min1 min2, max max1 max2)

let rec get_length l =
  let min, max = get_len2 l in
  Len.parse (Syntax.of_list
    (Token.Ident "in" :: Token.Nat min :: Token.Dot :: Token.Dot :: Token.Nat max :: []))
and get_len2 =
  function
  | [] -> 0, 0
  | Begin::l -> get_len2 l
  | End::_ -> 0, 0
  | Atom _::l -> let min, max = get_len2 l in min+1, max+1
  | Elastic (min,max)::l -> let min', max' = get_len2 l in min+min', max+max'

let rec get_dist a l =
  let round x = floor (1000. *. x) /. 1000. in
  let (xmin, ymin), (xmax, ymax) = get_dist2 a l in
  let dmin = round (float xmin /. float ymin) in
  let dmax = round (float xmax /. float ymax) in
  Dist.parse (Syntax.of_list (Token.Ident "in" :: Syntax.toks_of_float (dmin,-3) @ Token.Dot :: Token.Dot :: Syntax.toks_of_float (dmax,-3)))
and get_dist2 a =
  function
  | [] -> (0,0), (0,0)
  | Begin::l -> get_dist2 a l
  | End::_ -> (0,0), (0,0)
  | Atom a'::l ->
      let (xmin,ymin), (xmax,ymax) = get_dist2 a l in
      if A.entails a' a
      then (xmin+1,ymin+1), (xmax+1,ymax+1)
      else (xmin,ymin+1), (xmax,ymax+1)
  | Elastic (min,max)::l ->
      let (xmin,ymin), (xmax,ymax) = get_dist2 a l in
      (xmin,ymin+max), (xmax+max,ymax+max)


(* public functions *)

open Token

let rec motif_cons n x xs =
  if n <= 0
  then xs
  else motif_cons (n-1) x
    ( match x, xs with
    | Begin, Begin::_ -> xs
    | Begin, _ -> Begin::xs
    | End, _ -> [End]
    | Atom a, _ -> Atom a::xs
    | Elastic (min,max), Elastic (min',max')::xs' -> Elastic (min+min',max+max')::xs'
    | Elastic e, _ -> Elastic e::xs)

let motif_append m1 m2 =
  List.fold_right (motif_cons 1) m1 m2

let print_e = function
(*
      | (0,1) -> [Interro]
*)
      | (1,1) -> []
      | (min,max) ->
          if max = min
          then [LeftPar; Nat min; RightPar]
          else [LeftPar; Nat min; Comma; Nat max; RightPar]

let rec print_seq = function
  | l -> print_seq_loop true l
and print_seq_loop first_atom = function
  | [] -> []
  | Begin::l -> Param.toks_begin @ print_seq_loop true l
  | End::l -> Param.toks_end @ print_seq_loop true l
  | x::l -> print_seq_sep first_atom @ print_atom x @ print_seq_loop false l
and print_seq_sep first_atom =
  if first_atom
  then []
  else Param.toks_sep
and print_atom = function
  | Atom a -> A.print a
  | Elastic e ->  Param.toks_any @ print_e e
  | _ -> assert false

let print_comment = function
  | "" -> []
  | c -> [Colon; String c]

let print = function
  | Match (_,l,_, c) -> Param.toks_match @ print_seq l @ print_comment c
  | Len length -> Ident "length" :: PP_tilda :: Len.print length
  | Dist None -> [Ident "dist"]
  | Dist (Some (a,d)) -> Ident "dist" :: PP_tilda :: A.print a @ PP_tilda :: Dist.print d

let parse_repeat = parser
  | [<'LeftPar; 'Nat n; 'RightPar>] -> n
  | [<>] -> 1

let parse_e = parser
      | [<'LeftPar;
          'Nat min ?? "Syntax error: positive integer expected in motif gap, after '('";
          str
         >] ->
          ( match str with parser
          | [<'RightPar>] -> (min,min)
          | [<'Comma;
              'Nat max when max >= min ?? "Syntax error: integer expected in motif gap, after: " ^ Syntax.stringizer [LeftPar; Nat min; Comma];
              'RightPar ?? "Syntax error: missing ')' in motif gap, after: " ^ Syntax.stringizer [LeftPar; Nat min; Comma; Nat max]
             >] -> (min,max)
          | [<>] -> raise (Stream.Error ("Syntax error: ')' or ',' expected after: " ^ Syntax.stringizer [LeftPar; Nat min]))
          )
      | [<>] -> (1,1)
(*
      | [<'Interro>] -> (0,1)
*)

let rec parse_seq = parser
  | [<_ = Syntax.parse_tokens Param.toks_begin;
      l, e = parse_seq_loop1 ?? "Syntax error: motif element, gap or end expected after: " ^ Syntax.stringizer Param.toks_begin
     >] -> true, motif_cons 1 Begin l, e
  | [<l, e = parse_seq_loop1>] -> false, l, e
  | [<>] -> false, [], false
and parse_seq_loop1 = parser
  | [<_ = Syntax.parse_tokens Param.toks_end>] ->
       [End], true
  | [<_ = Syntax.parse_tokens Param.toks_any;
      elast = parse_e;
      l, e = parse_seq_loop2
     >] ->
       motif_cons 1 (Elastic elast) l, e
  | [<a, _ = parse_atom; (*n = parse_repeat;*)
      l, e = parse_seq_loop2
     >] ->
       motif_cons 1 a l, e
  | [<>] -> [], false
and parse_seq_loop2 = parser
  | [<_ = Syntax.parse_tokens Param.toks_end>] -> [End], true
  | [<_ = Syntax.parse_tokens Param.toks_sep;
      l, e = parse_seq_loop1 ?? "Syntax error: motif element, gap or end expected after: " ^ Syntax.stringizer Param.toks_sep
     >] -> l, e
  | [<>] -> [], false
and parse_seq_is = parser
  | [<x, a = parse_atom; (*n = parse_repeat;*)
      l = parse_seq_is ?? "Syntax error: motif element expected after: " ^ Syntax.stringizer (A.print a)
     >] -> motif_cons 1 x l
  | [<>] -> [End]
and parse_atom = parser
  | [<a = A.parse>] ->
       if try A.entails (A.top ()) a with Not_found -> false
       then Elastic (1,1), a
       else Atom a, a

let parse_comment = parser
  | [<'Colon; 'String c>] -> c
  | [<>] -> ""

let parse = parser
  | [<'Ident l when l.[0] = 'l';
      length = Len.parse ?? "Syntax error in motif after '" ^ l ^ "'"
     >] -> Len length
  | [<'Ident d when d.[0] = 'd'; str>] -> (match str with parser
    | [<a = A.parse;
        d = Dist.parse ?? "Syntax error in motif after: " ^ Syntax.stringizer (Ident d :: A.print a)
       >] -> Dist (Some (a,d))
    | [<>] -> Dist None)
  | [<'Ident "is";
      toks = A.parse_compact ?? "Syntax error after 'is'";
      c = parse_comment
     >] -> Match (true, Begin::parse_seq_is (Syntax.of_list toks), true, c)
  | [<b, l, e = Syntax.parse_tokens_and Param.toks_match parse_seq;
      c = parse_comment
     >] -> Match (b,l,e,c)


let rec simpl = function
  | Match (b,l,e, c) ->
      [<simpl_begin (b,e) l; simpl_end (b,e) (List.rev l); simpl_middle (b,e) [] l; simpl_elt (b,e) [] l>]
  | _ -> [<>]
and simpl_begin (b,e) = function
  | _::Elastic (min,max)::l ->
     if min = 0
     then [<'Match (false,l,e,"")>]
     else [<'Match (false,Elastic (min,min)::l,e,"")>]
  | _::l -> [<'Match (false,l,e,"")>]
  | [] -> [<>]
and simpl_end (b,e) = function
  | _::Elastic (min,max)::l ->
     if min = 0
     then [<'Match (b,List.rev l,false,"")>]
     else [<'Match (b,List.rev (Elastic (min,min)::l),false,"")>]
  | _::l -> [<'Match (b,List.rev l,false,"")>]
  | [] -> [<>]
and simpl_middle (b,e) acc = function
  | [] -> [<>]
  | Atom a::l ->
     (*if acc = [] then simpl_middle (b,e) [x] l
     else if l = [] then [<>]
     else*) [<'Match (b,motif_append acc (motif_cons 1 (Elastic (1,1)) l),e,""); simpl_middle (b,e) (acc@[Atom a]) l>]
  | Elastic el::l -> simpl_middle (b,e) (acc@[Elastic el]) l
  | Begin::l -> simpl_middle (b,e) [Begin] l
  | End::_ -> [<>]
and simpl_elt (b,e) acc = function
  | [] -> [<>]
  | Atom a::l -> 
     [<Common.stream_map (fun a' -> Match (b,acc@(Atom a'::l),e,"")) (A.simpl a);
       simpl_elt (b,e) (acc@[Atom a]) l>]
  | Elastic (min,max)::l ->
     if max > min
     then [<'Match (b,acc@(if min=0 then l else Elastic (min,min)::l),e,"");
             simpl_elt (b,e) (acc@[Elastic (min,max)]) l>]
     else simpl_elt (b,e) (acc@[Elastic (min,max)]) l
  | Begin::l -> simpl_elt (b,e) [Begin] l
  | End::_ -> [<>]

(* logical operations *)

(* top undefined to make 'match []', 'len in ..', and 'dist a in ..' incomparable *)

let rec entails f g =
  match f, g with
  | Match (_,lf,_,_), Match (_,lg,_,_) -> match_contains lf lg
  | Len len1, Len len2 -> Len.entails len1 len2
  | Match (true,l,true,_), Len len -> Len.entails (get_length l) len
  | Dist _, Dist None -> true
  | Dist (Some (a1,d1)), Dist (Some (a2,d2)) -> A.entails a1 a2 & A.entails a2 a1 & Dist.entails d1 d2
  | Match (true,l,true,_), Dist None -> true
  | Match (true,l,true,_), Dist (Some (a,d)) -> Dist.entails (get_dist a l) d
  | _, _ -> false
and match_contains lf lg =
  match_begin lf (match lg with [] -> [] | Begin::_ -> lg | _ -> motif_cons 1 (Elastic (0,snd (get_len2 lf))) lg)
and match_begin f g =
  match f, g with
  | _, [] -> true
  | _, [Elastic (0,_)] -> true
  | [], _ -> false
  | Elastic (min1,max1)::f', Elastic (min2,max2)::g' ->
      max1 <= max2 &
      let min3, max3 = max 0 (min2-min1), max2-max1 in
      min3 <= max3 &
      match_begin f' (Elastic (min3, max3)::g')
  | Elastic (_,_)::_, _::_ -> false
  | Begin::f', Begin::g' ->
      match_begin f' g'
  | Begin::f', (Atom _ | End)::_ -> false
  | End::_, End::_ -> true
  | End::_, (Atom _ | Begin)::_ -> false
  | Atom a::f', Atom b::g' ->
      A.entails a b & match_begin f' g'
  | Atom _::_, (Begin | End)::_ -> false
  | _::f', Elastic (0,0)::g' ->
      match_begin f g'
  | Begin::f', Elastic (_,_)::_ ->
      match_begin f' g
  | End::_, Elastic (0,_)::End::_ -> true
  | End::_, Elastic (_,_)::_ -> false
  | Atom a::f', Elastic (0,max)::Atom b::g' ->
      (A.entails a b & match_begin f' g') or (match_begin f' (Elastic (0, max-1)::Atom b::g'))
  | Atom _::f', Elastic (mn,mx)::g' ->
      match_begin f' (Elastic (max 0 (mn-1),mx-1)::g')

let conj f g =
  if entails f g then f
  else if entails g f then g
  else raise Not_found

let add f g = conj f g

let sub f g = if entails f g then raise Not_found else f

let rec features = function
  | Match (true,l,true,_) ->
      (true,Match (false,[],false,"")) ::
      (if Param.feat_length then features (Len (get_length l)) else []) @
      (if Param.feat_dist
       then List.fold_left
        (fun res a -> features (Dist (Some (a,get_dist a l))) @ res)
        []
        (List.fold_left (fun res -> function Atom a -> LSet.add a res | _ -> res) (LSet.empty ()) l)
       else []) @
      (if Param.feat_atoms
      then
	List.fold_left
	  (fun res -> function
	    | Atom a -> List.map (fun (vis,x) -> (vis,Match (false,[Atom x],false,""))) (A.features a) @ res
	    | _ -> res)
	  []
	  l
      else [])
  | Match m -> []
  | Len len -> List.map (fun (vis,x) -> vis,Len x) (Len.features len)
  | Dist None -> []
  | Dist (Some (a,d)) -> (true,Dist None) :: List.map (fun (vis,x) -> vis,Dist (Some (a,x))) (Dist.features d)

(* operation 'gen' *)

type heur = float
type coord = int * int
type elast = int * int
type link = {e : elast; c : coord}
type cell = {
    mutable zs : a LSet.t;
    mutable accs : coord LSet.t;
    mutable succs : link list;
    mutable preds : link list;
    mutable heur_right : heur;
    mutable heur_left : heur;
    mutable heur_end : bool;
  }
type mat_dag = {
    nf : int;
    ng : int;
    arfa : a array;
    arga : a array;
    arfe : elast array;
    arge : elast array;
    mat : cell array array;
    ht_freq : (a, int ref * int ref) Hashtbl.t;
      (* for the nb of z in mat.zs subsumed by some entry key of the hash table *) 
  }

let heur_min =
  (* 0. *)
  1. (* prob *)

let heur_max h1 h2 =
  (* max h1 h2 *)
  min h1 h2 (* prob *)

let heur_better h1 h2 =
  (* h1 >= h2 *)
  h1 <= h2 (* prob *)

let heur_atom md = function
  | z::_ ->
      (* -. log (float !(snd (Hashtbl.find md.ht_freq z)) /. float (md.nf * md.ng)) *)
      (*float (md.nf * md.ng + 1 - !(snd (Hashtbl.find md.ht_freq z))) *)
      sqrt (float !(snd (Hashtbl.find md.ht_freq z)) /. float (md.nf * md.ng)) (* prob *)
  | _ -> assert false

let heur_elastic md (min,max) =
(*  let k = 1. /. float (md.nf * md.ng) in *)
(*  float min *. (k +. k) +. (k /. float (max - min + 1)) *)
(*  float (min+1) /. float (max+1) /. float (max+1) *)
(*  (if max = 0 then 0. else -5.) -. (0.5 *. float (max - min)) *)
  float (max - min + 1) (* prob *)

let heur_join md v1 e v2 =
  let h = heur_elastic md e in
  (* v1 +. h +. v2 *)
  v1 *. h *. v2 (* prob *)

let get_link md (i,j) (i',j') =
  assert (i<i' & j<j');
  let e_i : e ref = let d_i = i' - i - 1 in ref (d_i,d_i) in
  for k = i to i'-1 do e_i := e_append !e_i md.arfe.(k) done;
  let e_j : e ref = let d_j = j' - j - 1 in ref (d_j,d_j) in
  for k = j to j'-1 do e_j := e_append !e_j md.arge.(k) done;
  let e = e_union !e_i !e_j in
  {e = e; c = (i',j')}

let add_link md (i,j) cell cell' l =
  cell.succs <- l::cell.succs;
  cell'.preds <- {l with c=(i,j)}::cell'.preds

let get_add_link md (i,j) cell (i',j') cell' =
  let l = get_link md (i,j) (i',j') in
  add_link md (i,j) cell cell' l;
  l

let remove_link md (i,j) cell (i',j') cell' =
  cell.succs <- List.filter (fun l -> l.c <> (i',j')) cell.succs;
  cell'.preds <- List.filter (fun l -> l.c <> (i,j)) cell'.preds

let gen_match_matdag_succs_cond md (i,j) =
  let i1, j1 = i+1, j+1 in
  let cell = md.mat.(i).(j) in
  let cell_right = md.mat.(i).(j1) in
  let cell_down = md.mat.(i1).(j) in
  let cell_diag = md.mat.(i1).(j1) in

  let my_add_link cell' l =
    if Param.gen_floating or cell'.zs = [End] or cell'.succs <> [] (* cell is reachable from end *)
    then add_link md (i,j) cell cell' l in

  cell.accs <- begin
    let accs1 = LSet.union cell_right.accs cell_down.accs in
    if LSet.is_empty cell_diag.zs then accs1 else LSet.add (i1,j1) accs1 end;
  cell_diag.accs <- []; (* to allow some memory to be freed *)

  if not (LSet.is_empty cell.zs)
  then begin
    let accs_good =
      List.iter
        (fun (i',j') ->
           let cell' = md.mat.(i').(j') in
           if ((i,j) = (0,0) & cell'.preds=[])
           then my_add_link cell' {e = (0,0); c = (i',j')} (* Warning: e is not significant *)
           else
             let d = max (i'-i-1) (j'-j-1) in
             if Param.gen_link_cond (cell.zs=[Begin]) (d,d) (cell'.zs=[End])
             then
               let l = get_link md (i,j) (i',j') in
               if Param.gen_link_cond (cell.zs=[Begin]) l.e (cell'.zs=[End])
               then my_add_link cell' l)
        cell.accs in
    if not Param.gen_floating & cell.zs <> [End] & cell.succs = [] (* cell is not reachable from end *)
    then cell.zs <- LSet.empty () end

let gen_match_matdag lf lg lhs =
  let nf = List.fold_left (fun n -> function Elastic _ -> n | _ -> n+1) 0 lf in
  let ng = List.fold_left (fun n -> function Elastic _ -> n | _ -> n+1) 0 lg in
  let md : mat_dag = {
    nf = nf;
    ng = ng;
    arfa = Array.make (nf+1) Begin;
    arga = Array.make (ng+1) Begin;
    arfe = Array.make nf (0,0);
    arge = Array.make ng (0,0);
    mat = Array.make_matrix (nf+2) (ng+2)
      { zs=LSet.empty ();
        accs=LSet.empty ();
        succs=[]; preds=[]; 
        heur_right=heur_min; heur_left=heur_min; heur_end=true};
    ht_freq = Hashtbl.create nf
  } in
  let cell0 = md.mat.(0).(0) in
(* initializing arfa and arfe *)
  ignore (List.fold_left (fun i ->
    function
    | Elastic e -> md.arfe.(i) <- e; i
    | a -> md.arfa.(i+1) <- a; i+1
  ) 0 lf);
(* initializing arga and arge *)
  ignore (List.fold_left (fun j ->
    function
    | Elastic e -> md.arge.(j) <- e; j
    | a -> md.arga.(j+1) <- a; j+1
  ) 0 lg);
(* computing the fields zs, succs of cells *)
  let zhs =
        List.fold_left
          (fun res lh ->
             List.fold_left (fun res' ->
               function
               | Atom a -> LSet.add a res'
               | _ -> res')
             res lh)
          (LSet.empty ()) lhs in
  for s = md.nf + md.ng downto 2 do
    for i = max 1 (s-md.ng) to min md.nf (s-1) do
      let j = s - i in
    (* computing zs *)
      let zs =
        match md.arfa.(i), md.arga.(j) with
        | Atom x, Atom y ->
              let zhs' =
                List.fold_left
                  (fun res z ->
                     if A.entails x z & A.entails y z & not (List.exists (fun z' -> A.entails z' z) res)
                     then LSet.add z (List.filter (fun z' -> not (A.entails z z')) res)
                     else res)
                  (LSet.empty ()) zhs in
              ( match Common.mapfilter
                (fun z -> try if not (A.entails (A.top ()) z) then Some (Atom z) else None with Not_found -> Some (Atom z))
                (let zs' = A.gen x y zhs' in if zs'=[] then zhs' else zs') with
              | [] -> LSet.empty () 
              | x::_ -> LSet.singleton x)
        | Begin, Begin -> LSet.singleton Begin
        | End, End -> LSet.singleton End
        | _, _ -> LSet.empty () in
      List.iter (fun z ->
        try incr (fst (Hashtbl.find md.ht_freq z))
        with Not_found -> Hashtbl.add md.ht_freq z (ref 1,ref 0)
      ) zs;
    (* computing succs *)
      let cell = {
        zs=zs; 
        accs=LSet.empty ();
        succs=[]; preds=[]; 
        heur_right=heur_min; heur_left=heur_min; heur_end=true} in
      md.mat.(i).(j) <- cell;
      gen_match_matdag_succs_cond md (i,j)
    done
  done;
(* computing frequencies of each element *)
  Hashtbl.iter (fun z (c1,c2) ->
    Hashtbl.iter (fun z' (c1',c2') ->
      if a_entails z' z then c2 := !c2 + !c1'
    ) md.ht_freq
  ) md.ht_freq;
(* computing succs of cell (0,0), taking into account Param.gen_floating *)
  if Param.gen_floating then gen_match_matdag_succs_cond md (0,0)
  else if md.nf>0 & md.ng>0 & md.mat.(1).(1).zs=[Begin] then ignore (get_add_link md (0,0) cell0 (1,1) md.mat.(1).(1))
  else ();
(* in case Param.gen_floating is false, erase unreachable cells from Begin *)
  if not Param.gen_floating then
    for s = 2 to md.nf+md.ng do
      for i = max 1 (s-md.ng) to min md.nf (s-1) do
        let j = s-i in
        let cell = md.mat.(i).(j) in
        if cell.preds = [] & cell.zs <> LSet.empty ()
        then begin
          cell.zs <- LSet.empty ();
          List.iter (fun l ->
            let cell' = md.mat.(fst l.c).(snd l.c) in
            cell'.preds <- List.filter (fun l -> l.c <> (i,j)) cell'.preds
          ) cell.succs;
          cell.succs <- [] end
      done
    done;
(* defining the field 'max_heur' of cells *)
  for s = md.nf+md.ng downto 2 do
    for i = max 1 (s-md.ng) to min md.nf (s-1) do
      let j = s-i in
      let cell = md.mat.(i).(j) in
      if not (LSet.is_empty cell.zs) then
        let h = heur_atom md cell.zs in
        let hr, he =
          List.fold_left
            (fun (hr,he) l ->
               let hr' = heur_join md h l.e md.mat.(fst l.c).(snd l.c).heur_right in
               if heur_better hr' hr or (not Param.gen_floating & he)
               then (hr',false)
               else (hr,he))
            (h,true)
            cell.succs in
        cell.heur_right <- hr;
        cell.heur_end <- he
    done
  done;
  for s = 2 to md.nf+md.ng do
    for i = max 1 (s-md.ng) to min md.nf (s-1) do
      let j = s-i in
      let cell = md.mat.(i).(j) in
      if not (LSet.is_empty cell.zs) then
        let h = heur_atom md cell.zs in
        let hl, hb =
          List.fold_left
            (fun (hl,hb) l ->
               if l.c = (0,0)
               then (hl,false)
               else
                 let hl' = heur_join md h l.e md.mat.(fst l.c).(snd l.c).heur_left in
                 if heur_better hl' hl or (not Param.gen_floating & hb)
                 then (hl',false)
                 else (hl,hb))
            (h,true)
            cell.preds in
        cell.heur_left <- hl;
        if Param.gen_floating & hb
        then let cell0 = md.mat.(0).(0) in cell0.succs <- {e=(0,0); c=(i,j)}::cell0.succs
          (* make this node a possible start for a motif, because every left extension is worse *)
    done
  done;
(* returning the results *)
  md


let gen_insert_motif m (ms, lhs) =
  if List.exists (fun lh -> match_contains lh m) lhs
  then (ms, lhs)
  else (LSet.add m ms, lhs)

let gen_insert_motifs ms' (ms,lhs) =
  List.fold_left
    (fun res m' -> gen_insert_motif m' res)
    (ms,lhs)
    ms'

module Sps = Set.Make
(struct
  type t = heur * heur * a list * cell
  let compare (h1,_,_,_ as x1) (h2,_,_,_ as x2) =
    if x1=x2 then 0
    else if heur_better h1 h2 then 1
    else -1
end)

let rec gen_match_motifs md ps (ms,lhs) =
  if List.length ms >= Param.gen_min_xs
  then ms
  else
    if Sps.is_empty ps
    then ms
    else
      let (h,hm,m,cell as x) = Sps.max_elt ps in
      let ps' = Sps.remove x ps in
         let (ms',lhs') =
           if cell.succs = [] or (Param.gen_floating & cell.heur_end)
           then begin
             gen_insert_motif m (ms,lhs) end
           else (ms,lhs) in
         let ps'' =
           List.fold_left
             (fun res l ->
                let cell' = md.mat.(fst l.c).(snd l.c) in
                let h' = heur_join md hm l.e cell'.heur_right in
                let hm' = heur_join md hm l.e (heur_atom md cell'.zs) in
                let ma = [List.hd cell'.zs] in
                Sps.add (h', hm', m @ (if l.e=(0,0) then [] else [Elastic l.e]) @ ma, cell') res)
             ps'
             cell.succs in
         gen_match_motifs md ps'' (ms',lhs')


let gen_match lf lg lhs =
  let md = gen_match_matdag lf lg lhs in
  let ps =
    List.fold_left
      (fun res l ->
         let cell = md.mat.(fst l.c).(snd l.c) in
         Sps.add (cell.heur_right, heur_atom md cell.zs, [List.hd cell.zs], cell) res)
      Sps.empty
      md.mat.(0).(0).succs in
  gen_match_motifs md ps ([],lhs)


let rec gen f g hs =
  match f, g with
  | Match (_,[],_,_), Match _
  | Match _, Match (_,[],_,_) -> LSet.empty ()
  | Match (bf,lf,ef,_), Match (bg,lg,eg,_) ->
      let lens =
        if bf & ef & bg & eg & Param.gen_len
        then gen (Len (get_length lf)) (Len (get_length lg)) hs
        else LSet.empty () in
      if not (LSet.is_empty lens) then lens else
      let lhs = List.fold_left (fun res -> function Match (_,lh,_,_) -> LSet.add lh res | _ -> res) (LSet.empty ()) hs in
      let motifs = gen_match lf lg lhs in
      List.map
        (fun m ->
           let b = try List.hd m = Begin with _ -> false in
           let e = try List.hd (List.rev m) = End with _ -> false in
           Match (b,m,e,""))
        motifs
  | Match (true,lf,true,_), Len leng ->
      gen (Len (get_length lf)) (Len leng) hs
  | Len lenf, Match (true,lg,true,_) ->
      gen (Len lenf) (Len (get_length lg)) hs
  | Len lenf, Len leng ->
      let lenhs = Common.mapfilter (function Len len -> Some len | _ -> None) hs in
      List.fold_left
        (fun res len -> LSet.add (Len len) res)
        (LSet.empty ())
        (Len.gen lenf leng lenhs) 
  | _, _ -> LSet.empty ()

end
