let quantile_normalization aa =
  assert (is_rectangular aa);
  if Array.length aa = 0 || Array.length aa.(0) = 0
     || Array.length aa.(0) = 1 then
    Array.copy aa
  else
    let num_expts = Float.of_int (Array.length aa.(0)) in
    let num_pts = Array.length aa in

    let comp1 (a,_) (b,_) = Pervasives.compare a b in
    let comp2 (_,a) (_,b) = Pervasives.compare a b in

    let aa = transpose aa in
    let aa = Array.map (Array.mapi ~f:(fun a b -> (a, b))) aa in
    (Array.iter ~f:(Array.sort ~cmp:comp2)) aa;
    let avg i =
      (Array.fold ~f:(fun sum expt -> snd expt.(i) +. sum) ~init:0.0 aa)
      /. num_expts in
    let norms = Array.init num_pts avg in
    let aa = Array.map (Array.mapi ~f:(fun i (idx,_) -> idx, norms.(i))) aa in
    Array.iter ~f:(Array.sort ~cmp:comp1) aa;
    transpose (Array.map ~f:(Array.map ~f:snd) aa)