let transition_rule_4bodies beta prng neighbourhoods =
  let nb_neighbours_max = Array.fold_left
    (fun accu neighbourhood ->
       max accu
         (Array.fold_left (fun sum neighbour -> sum + Array.length neighbourhoods.(neighbour)) 0 neighbourhood))
    0 neighbourhoods in
  let terms = Array.init (2*nb_neighbours_max+1)
    (fun i ->
       let sum_neighbours = i-nb_neighbours_max in
         cosh (beta *. float (sum_neighbours-1)) /.
         cosh (beta *. float (sum_neighbours+1))) in
fun configuration cell ->
  let product accu neighbour =
    let sum_neighbours = ref 0 in
    let nn = neighbourhoods.(neighbour) in
    for i=0 to pred (Array.length nn) do
      let cell' = nn.(i) in
      if cell' <> cell then
        if configuration.(nn.(i)) then incr sum_neighbours else decr sum_neighbours;
    done;
    accu *.
      terms.(!sum_neighbours+nb_neighbours_max)
  in
(*   Array.fold_left product 1. neighbourhoods.(cell) *)
(*   |> barker *)
(*   |> metropolis prng *) (*seems to prevent inlining*)
  metropolis prng (barker (Array.fold_left product 1. neighbourhoods.(cell)))