pdist2 full benchmarking script

function bench_pdist2_full ()
  N_vals = [50, 100, 200, 500, 1000, 2000, 3000, 4000, 5000];
  P = 10;
  n_runs = 7;
  n_warmup = 3;

  n_sizes = numel (N_vals);
  res = zeros (n_sizes, 6);

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

    [t_orig, std_orig] = run_timed (@() pdist2_orig (X, Y), n_runs, n_warmup);
    [t_blocked, std_blocked] = run_timed (@() pdist2_blocked (X, Y), n_runs, n_warmup);

    speedup = t_orig / t_blocked;
    res(idx, :) = [N, t_orig, std_orig, t_blocked, std_blocked, speedup];

    clear X Y;
  endfor

  printf ("\npdist2: Original (meshgrid) vs Blocked (row-wise)\n");
  printf ("Config: P=%d, runs=%d, warmup=%d\n", P, n_runs, n_warmup);
  printf ("%s\n", repmat ("-", 1, 75));
  printf ("%-6s %12s %12s %12s %12s %10s\n", "N", "t_orig", "std_orig", "t_blocked", "std_blocked", "speedup");
  printf ("%s\n", repmat ("-", 1, 75));
  for i = 1:n_sizes
    printf ("%-6d %12.6f %12.6f %12.6f %12.6f %10.2fx\n", ...
            res(i,1), res(i,2), res(i,3), res(i,4), res(i,5), res(i,6));
  endfor
  printf ("%s\n\n", repmat ("-", 1, 75));

  fid = fopen ("bench_pdist2_full.csv", "w");
  fprintf (fid, "N,t_orig,std_orig,t_blocked,std_blocked,speedup\n");
  for i = 1:n_sizes
    fprintf (fid, "%d,%.6f,%.6f,%.6f,%.6f,%.4f\n", ...
             res(i,1), res(i,2), res(i,3), res(i,4), res(i,5), res(i,6));
  endfor
  fclose (fid);
  printf ("Saved: bench_pdist2_full.csv\n");

  h = figure ("visible", "off");
  loglog (N_vals, res(:,2), "ro-", "linewidth", 2, "markersize", 8, "markerfacecolor", "r");
  hold on;
  loglog (N_vals, res(:,4), "bs-", "linewidth", 2, "markersize", 8, "markerfacecolor", "b");
  hold off;
  xlabel ("N"); ylabel ("Time (s)"); title ("pdist2: Original vs Blocked");
  legend ({"Original (meshgrid)", "Blocked (row-wise)"}, "location", "northwest");
  grid on;
  print (h, "bench_pdist2_full.png", "-dpng", "-r150");
  close (h);
  printf ("Saved: bench_pdist2_full.png\n");

  h2 = figure ("visible", "off");
  semilogx (N_vals, res(:,6), "ko-", "linewidth", 2, "markersize", 8, "markerfacecolor", "g");
  hold on;
  plot ([N_vals(1), N_vals(end)], [1, 1], "r--", "linewidth", 1);
  hold off;
  xlabel ("N"); ylabel ("Speedup (t_orig / t_blocked)");
  title ("pdist2: Speedup from Blocked Implementation");
  grid on;
  ylim ([0, max(res(:,6)) * 1.1]);
  print (h2, "bench_pdist2_speedup.png", "-dpng", "-r150");
  close (h2);
  printf ("Saved: bench_pdist2_speedup.png\n");
endfunction

function [t_med, t_std] = 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_med = median (times);
  t_std = std (times);
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
    D(i, :) = sqrt (sum ((X(i, :) - Y) .^ 2, 2))';
  endfor
endfunction