################################################################################
# Preamble to load or create data
if (exist ("hist.var") == 2)
  load hist.var               # Load variables
else
  rand ("state", 1);
  x1 = rand (1e6,1);
  x2 = rand (1e4,1);
  x3 = rand (1e3,1);
  x4 = rand (1e2,1);
  save -binary hist.var x1 x2 x3 x4
endif

[nn, xx] = hist (1:10);       # compile m-file by running once

################################################################################
# Utility function
function retval = discard_minmax (x)
  [~, idx] = min (x);
  x(idx) = [];
  [~, idx] = max (x);
  x(idx) = [];
  retval = x; 
endfunction

################################################################################
# Run benchmarking

N = 50 + 2;   # Number of tests at each size (plus 2 for discarding)

bm.x1e6 = zeros (N,1);
for i = 1:N
  tic;
  [nn, xx] = hist (x1);
  bm.x1e6(i) = toc;
endfor
bm.x1e6 = discard_minmax (bm.x1e6);

bm.x1e4 = zeros (N,1);
for i = 1:N
  tic;
  [nn, xx] = hist (x2);
  bm.x1e4(i) = toc;
endfor
bm.x1e4 = discard_minmax (bm.x1e4);

bm.x1e3 = zeros (N,1);
for i = 1:N
  tic;
  [nn, xx] = hist (x2);
  bm.x1e3(i) = toc;
endfor
bm.x1e3 = discard_minmax (bm.x1e3);

bm.x1e2 = zeros (N,1);
for i = 1:N
  tic;
  [nn, xx] = hist (x3);
  bm.x1e2(i) = toc;
endfor
bm.x1e2 = discard_minmax (bm.x1e2);

################################################################################
# Report results
printf ("1e6: %f\n", mean (bm.x1e6));
printf ("1e4: %f\n", mean (bm.x1e4));
printf ("1e3: %f\n", mean (bm.x1e3));
printf ("1e2: %f\n", mean (bm.x1e2));

