simplex_sampling.ml

let d = 3
let m = 1000
  
let replicate_with_failures k rng st =
  let rec aux l k =
    if k = 0 then l else
      match rng st with
      | Some x -> aux (x::l) (k-1)
      | None -> aux l (k-1) in
  aux [] k
  
let sample_without_replacement (type elt) ?compare:(compare:elt->elt->int=compare) k (rng:elt CCRandom.t) st=
  let module S = Set.Make(struct type t=elt let compare = compare end) in
  let rec aux s k =
    if k = 0 then
      S.elements s
    else
      let x = rng st in
      if S.mem x s then
        aux s k
      else
        aux (S.add x s) (k-1) in
  aux S.empty k

exception Diff_list_incompatible_size

let diff_list ( - ) l =
  let rec diff_list acc = function
    | [a;b] -> (b - a)::acc
    | a::( b::q as r ) -> diff_list ( (b-a)::acc ) r
    | [_] | [] -> raise Diff_list_incompatible_size
  in
  diff_list [] l  
    
let inner_simplex_uniform m ~len =
  let smap x =
    List.sort compare (m::0::x)
    |> diff_list ( - )
    |> fun x -> Some x in 
  CCRandom.( pure smap <*> sample_without_replacement (len-1) (int_range 1 @@ m-1) )


let sample f k  = replicate_with_failures k @@ f m ~len:d


exception Unexpected_size of int

let map_3 f = function
  | [ a; b; c] -> f a b c
  | l -> raise @@ Unexpected_size (List.length l)

let () =
  let sampler =
    if Array.length Sys.argv >= 2 then
      begin
      match Sys.argv.(1) with
      | "-uniform" -> inner_simplex_uniform
      | "-split_list" | _   -> CCRandom.split_list
      end
    else
      CCRandom.split_list in
  CCRandom.run (sample sampler 50000)
  |> List.iter @@ map_3 @@ Printf.printf "%d %d %d\n"