pdist & pdist2 benchmarking script

function bench_blocked ()
  N_vals = [500, 1000, 2000, 3000, 4000];
  P = 10;
  n_runs = 5;
  n_warmup = 2;

  res_pdist = zeros (numel (N_vals), 4);
  res_pdist2 = zeros (numel (N_vals), 4);

  for idx = 1:numel (N_vals)
    N = N_vals(idx);
    rand ("seed", 12345);
    X = rand (N, P);
    Y = rand (N, P);

    t_orig = run_timed (@() pdist_orig (X), n_runs, n_warmup);
    t_blocked = run_timed (@() pdist_blocked (X), n_runs, n_warmup);
    res_pdist(idx, :) = [N, t_orig, t_blocked, t_orig / t_blocked];

    t_orig = run_timed (@() pdist2_orig (X, Y), n_runs, n_warmup);
    t_blocked = run_timed (@() pdist2_blocked (X, Y), n_runs, n_warmup);
    res_pdist2(idx, :) = [N, t_orig, t_blocked, t_orig / t_blocked];
  endfor

  print_table ("pdist: Original vs Blocked", {"N", "t_orig", "t_blocked", "speedup"}, res_pdist);
  print_table ("pdist2: Original vs Blocked", {"N", "t_orig", "t_blocked", "speedup"}, res_pdist2);

  save_csv ("bench_pdist_blocked.csv", {"N", "t_orig", "t_blocked", "speedup"}, res_pdist);
  save_csv ("bench_pdist2_blocked.csv", {"N", "t_orig", "t_blocked", "speedup"}, res_pdist2);

  make_plot ("pdist_blocked.png", "pdist: Original vs Blocked", N_vals, ...
             res_pdist(:,2), res_pdist(:,3), "Original", "Blocked");
  make_plot ("pdist2_blocked.png", "pdist2: Original vs Blocked", N_vals, ...
             res_pdist2(:,2), res_pdist2(:,3), "Original", "Blocked");
endfunction

function t = run_timed (fn, n_runs, n_warmup)
  for w = 1:n_warmup, fn (); endfor
  times = zeros (n_runs, 1);
  for r = 1:n_runs, tic; fn (); times(r) = toc; endfor
  t = median (times);
endfunction

function print_table (title, headers, data)
  printf ("\n%s\n%s\n", title, repmat ("-", 1, 60));
  printf ("%-8s %12s %12s %12s\n", headers{:});
  for i = 1:rows (data)
    printf ("%-8d %12.4f %12.4f %12.2fx\n", data(i,1), data(i,2), data(i,3), data(i,4));
  endfor
  printf ("\n");
endfunction

function save_csv (fname, headers, data)
  fid = fopen (fname, "w");
  fprintf (fid, "%s\n", strjoin (headers, ","));
  for i = 1:rows (data)
    fprintf (fid, "%d,%.6f,%.6f,%.4f\n", data(i,1), data(i,2), data(i,3), data(i,4));
  endfor
  fclose (fid);
  printf ("Saved: %s\n", fname);
endfunction

function make_plot (fname, ttl, N_vals, t1, t2, leg1, leg2)
  h = figure ("visible", "off");
  semilogy (N_vals, t1, "ro-", "linewidth", 2, "markersize", 8);
  hold on;
  semilogy (N_vals, t2, "bs-", "linewidth", 2, "markersize", 8);
  hold off;
  xlabel ("N"); ylabel ("Time (s)"); title (ttl);
  legend ({leg1, leg2}, "location", "northwest");
  grid on;
  print (h, fname, "-dpng", "-r150");
  close (h);
  printf ("Saved: %s\n", fname);
endfunction

function Y = pdist_orig (X)
  n = rows (X);
  order = nchoosek (1:n, 2);
  Y = sqrt (sum ((X(order(:,1),:) - X(order(:,2),:)) .^ 2, 2))';
endfunction

function Y = pdist_blocked (X)
  n = rows (X);
  n_pairs = n * (n - 1) / 2;
  Y = zeros (1, n_pairs);
  idx = 1;
  for i = 1:n-1
    j_vals = (i+1):n;
    diff_block = X(i, :) - X(j_vals, :);
    dists = sqrt (sum (diff_block .^ 2, 2))';
    n_j = numel (j_vals);
    Y(idx:idx + n_j - 1) = dists;
    idx = idx + n_j;
  endfor
endfunction

function D = pdist2_orig (X, Y)
  [ix, iy] = meshgrid (1:rows (X), 1:rows (Y));
  D = sqrt (sum ((X(ix(:),:) - Y(iy(:),:)) .^ 2, 2));
  D = reshape (D, rows (Y), rows (X))';
endfunction

function D = pdist2_blocked (X, Y)
  M = rows (X);
  N = rows (Y);
  D = zeros (M, N);
  for i = 1:M
    diff_ij = X(i, :) - Y;
    D(i, :) = sqrt (sum (diff_ij .^ 2, 2))';
  endfor
endfunction