(* ---------- *)
(* Exercice 1 *)
(* Alignements *)

(* 1.1 Segment optimal *)

let somme_tableau (x : float array) : float =
  Array.fold_left (+.) 0. x
;;

let carre xi = xi *. xi;;

let somme_tableau_carre (x : float array) : float =
  Array.fold_left (+.) 0. (Array.map carre x)
;;

let somme_produit (x : float array) (y : float array) : float =
  somme_tableau (Array.map2 ( *. ) x y)
;;

let segment_opt (x : float array) (y : float array) (p : int) (q : int) : float * float =
  (* extraire x et y, et calculer q - p + 1 *)
  let x = Array.sub x p (q - p + 1)
  and y = Array.sub y p (q - p + 1) in
  let q_mp_p1 = float_of_int (q - p + 1) in
  (* les sommes et produits etc *)
  let somme_xicarre = somme_tableau_carre x
  and somme_xi_carre = carre (somme_tableau x)
  and somme_xiyi = somme_produit x y
  and somme_xi = somme_tableau x
  and somme_yi = somme_tableau y in
  (* d'abord a *)
  let nom_a = (q_mp_p1 *. somme_xiyi) -. (somme_xi *. somme_yi)
  and denom_a = (q_mp_p1 *. somme_xicarre) -. somme_xi_carre in
  let a = nom_a /. denom_a in
  (* ensuite b *)
  let nom_b = somme_yi -. (a *. somme_xi)
  and denom_b = q_mp_p1 in
  let b = nom_b /. denom_b in
  a, b
;;


let cout_un_segment (x : float array) (y : float array) (p : int) (q : int) : float =
  let a, b = segment_opt x y p q in
  let x = Array.sub x p (q - p + 1)
  and y = Array.sub y p (q - p + 1) in
  let yi_maxi_mb = Array.map2 (fun xi yi -> yi -. (a *. xi) -. b) x y in
  somme_tableau_carre yi_maxi_mb
;;

(* Exemples *)

let x = [|1.0; 2.0; 3.0; 4.0; 5.0; 6.0; 6.4; 6.5; 7.0; 8.0; 9.0; 10.|];;

let y = [|1.1; 2.1; 2.5; 2.7; 3.3; 4.2; 2.5; 2.0; 1.0; 1.5; 1.9; 2.4|];;


#load "graphics.cma";;
open Graphics;;
let rayon = 5;;

let h = 50.0;;

let dilate xi yi =
  (int_of_float (h *. xi), int_of_float (h *. yi))
;;

let trace_point xi yi : unit =
  let x, y = dilate xi yi in
  fill_circle x y rayon
;;

(* Trace le nuage de points *)
let montre_points x y : unit =
  open_graph "";
  clear_graph ();
  let n = Array.length x in
  for i = 0 to (n - 1) do
    trace_point x.(i) y.(i);
  done;
;;

(* montre_points x y;; *)


let penalite = 1.0;;

let cout_min (x : float array) (y : float array) (pen : float) : float =
  let n = Array.length x in
  assert (n = Array.length y);
  (* on stocke E_q dans erreurs.(q+1) pour q = 0, ..., n-1 *)
  let erreurs = Array.make n infinity in
  erreurs.(0) <- pen; (* 1 ou 2 points *)
  erreurs.(1) <- pen; (* 1 unique segment, sans erreur *)
  for q = 2 to n-1 do
    erreurs.(q) <- erreurs.(q-1) +. pen;
    for p = q-1 downto 0 do
      let eps = cout_un_segment x y p q in
      let err_pm1 = if p = 0 then 0.0 else erreurs.(p - 1) in
      let nouvelle_erreur = pen +. eps +. err_pm1 in
      if nouvelle_erreur < erreurs.(q) then begin
        erreurs.(q) <- nouvelle_erreur;
      end;
    done;
  done;
  erreurs.(n-1)
;;

cout_min x y penalite;;

(* 1.3 *)

let cout_min_2 (x : float array) (y : float array) (pen : float) : float array * int array =
  let n = Array.length x in
  assert (n = Array.length y);
  (* on stocke E_q dans erreurs.(q+1) pour q = 0, ..., n-1 *)
  let erreurs = Array.make n infinity in
  erreurs.(0) <- pen; (* 1 ou 2 points *)
  erreurs.(1) <- pen; (* 1 unique segment, sans erreur *)
  let antecedent = Array.make n 0 in
  for q = 2 to n-1 do
    erreurs.(q) <- erreurs.(q-1) +. pen;
    for p = q-1 downto 0 do
      let eps = cout_un_segment x y p q in
      let err_pm1 = if p = 0 then 0.0 else erreurs.(p - 1) in
      let nouvelle_erreur = pen +. eps +. err_pm1 in
      if nouvelle_erreur < erreurs.(q) then begin
        erreurs.(q) <- nouvelle_erreur;
        antecedent.(q) <- p;
      end;
    done;
  done;
  erreurs, antecedent
;;

let parcours a =
  let n = Array.length a in
  let i = ref (n - 1) in
  while !i >= 0 do
  let j = a.(!i) - 1 in
  print_int (!i - j);
  print_newline ();
  i := j
  done;
;;

let err, an = cout_min_2 x y penalite;;
parcours an;;

let cout_min_memo x y penalite =
  let n = Array.length x in
  let tab_cout = Array.make n infinity in
  tab_cout.(0) <- penalite;
  tab_cout.(1) <- penalite;
  let rec cout_min_rec q =
    if q < 0 then
      0.
    else if tab_cout.(q) <> infinity then
      tab_cout.(q)
    else begin
      let err_min = ref (penalite +. (cout_min_rec (q - 1))) in
      (* Initialisation de la recherche du minimum :
      cas où le point est seul sur son segment (p = q) *)
      for p = q - 1 downto 0 do
        let eps = cout_un_segment x y p q in
        let err = penalite +. eps +. cout_min_rec (p - 1) in
        if err < !err_min then
          err_min := err;
      done;
      tab_cout.(q) <- !err_min;
      !err_min;
    end;
  in
  cout_min_rec (n - 1)
;;

let err = cout_min_memo x y penalite;;