let transition_rule_4bodies_optimized beta prng neighbourhoods =
  let nb_neighbours = 8 in
  let terms = Array.init (2*nb_neighbours+1)
    (fun i ->
       let sum_neighbours = i-nb_neighbours in
         cosh (beta *. float (sum_neighbours-1)) /.
         cosh (beta *. float (sum_neighbours+1))) in
fun configuration cell ->
  let n = neighbourhoods.(cell) in
  assert (Array.length n = 8);
  let cell0 = configuration.(n.(0)) in
  let cell1 = configuration.(n.(1)) in
  let cell2 = configuration.(n.(2)) in
  let cell3 = configuration.(n.(3)) in
  let cell4 = configuration.(n.(4)) in
  let cell5 = configuration.(n.(5)) in
  let cell6 = configuration.(n.(6)) in
  let cell7 = configuration.(n.(7)) in
(*   let cell_centre = if configuration.(cell) then 1 else -1 in ### *)
  let term c1 c2 c3 =
(*     let sum_neighbours = ref cell_centre in *)
(*     if c1 then incr sum_neighbours else decr sum_neighbours; *)
    let sum_neighbours = ref (if c1 then 1 else -1) in
    if c2 then incr sum_neighbours else decr sum_neighbours;
    if c3 then incr sum_neighbours else decr sum_neighbours;
    terms.(!sum_neighbours+nb_neighbours)
 in
(*   term cell0 cell1 cell2 *. term cell1 cell3 cell4 *. term cell4 cell7 cell6 *. term cell2 cell5 cell6 *)
(*   |> barker *)
(*   |> metropolis prng *) (*seems to prevent inlining*)
  metropolis prng
    (barker
       (term cell0 cell1 cell2 *. term cell1 cell3 cell4 *. term cell4 cell7 cell6 *. term cell2 cell5 cell6))