let causal_states significance_level discard_cell neighbourhoods past_depth pof_depth =
  let length = Array.length neighbourhoods in
  let plc_indexes = Array.init length (* plc in terms of indexes (instead of states) *)
    (fun i -> if discard_cell i then [||]
(*      else Light_cones.build_lc_index neighbourhoods    pop_depth i) in *)
     else Light_cones.build_lc_index neighbourhoods (succ past_depth) i
          |> (ArrayLabels.sub ~pos:1 ~len:past_depth)) in
  let flc_indexes = Array.init length
    (fun i -> if discard_cell i then [||]
(*      else Light_cones.build_lc_index neighbourhoods (succ future_depth) i *)
(*           |> (ArrayLabels.sub ~pos:1 ~len:future_depth)) in *)
     else Light_cones.build_lc_index neighbourhoods pof_depth i) in
  let plc_identifiers = Light_cones.lc_identifiers plc_indexes in
  let flc_identifiers = Light_cones.lc_identifiers flc_indexes in
(* fun pop_configurations future_configurations -> *)
fun past_configurations pof_configurations ->
(*   let plcs, nb_plcs = plc_identifiers (Light_cones.lcState_of_lcIndex    pop_configurations) in *)
(*   let flcs, nb_flcs = flc_identifiers (Light_cones.lcState_of_lcIndex future_configurations) in *)
  let plcs, nb_plcs = plc_identifiers (Light_cones.lcState_of_lcIndex past_configurations) in
  let flcs, nb_flcs = flc_identifiers (Light_cones.lcState_of_lcIndex  pof_configurations) in
  let unsorted_distributions = Array.make nb_plcs [] in
  Femtolib.Array.iter2 (fun plc flc -> (* for each point (plc,flc), unsorted_distributions.(plc) *)
    if plc <> -1 && flc <> -1 then     (* contains one element, equal to flc. *)
    unsorted_distributions.(plc) <- flc :: unsorted_distributions.(plc)) plcs flcs;
  (*In a bad mix of side effect/return value, returns [fold_left (+) 0 unsorted_distribution]*)
  (* [output.(flc)] will be overwritten with the number of occurences of flc in
     [unsorted_distribution].*)
  let sums unsorted_distribution output =
    List.fold_left
      (fun accu flc -> output.(flc) <- output.(flc) +1; succ accu)
      0 unsorted_distribution in
  let clusters = Clustering.create_clusters () in
  let distribution = Array.make nb_flcs 0 in
  let test = Statistic.Chi_square.close significance_level in
  let add_unsorted_distribution unsorted_distribution =
    for i=0 to pred nb_flcs do distribution.(i) <- 0 done;
    let nb_samples = sums unsorted_distribution distribution in(*!mix of side effect and return*)
    Clustering.add clusters test distribution nb_samples in
(*   let causal_state_of_plcId = Array.map add_unsorted_distribution unsorted_distributions in *)
  let causal_state_of_plcId = Array.make nb_plcs (-2) in (* -2 is a dummy value *)
  let shuffle = Array.init nb_plcs identity in
  Femtolib.Array.shuffle shuffle;
  for i=0 to pred nb_plcs do (* add the distributions in the random order specified by [shuffle]*)
    let i = shuffle.(i) in
    causal_state_of_plcId.(i) <- add_unsorted_distribution unsorted_distributions.(i)
  done;
  (* An array of the same size as the lattice, containing the causal_states. *)
  let causal_states = Array.map (function -1 -> -1 | plc -> causal_state_of_plcId.(plc)) plcs in
  let nbOccurences_of_causalStateId = Clustering.nb_samples clusters in
(*   Printf.printf " (Array.length nbOccurences_of_causalStateId); *)
  causal_states, nbOccurences_of_causalStateId