pdist & pdist2 benchmarking script (switch statement placement)

function bench_switch ()
  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_inside = run_timed (@() pdist_switch_inside (X), n_runs, n_warmup);
    t_outside = run_timed (@() pdist_switch_outside (X), n_runs, n_warmup);
    res_pdist(idx, :) = [N, t_inside, t_outside, t_inside / t_outside];

    t_inside = run_timed (@() pdist2_switch_inside (X, Y), n_runs, n_warmup);
    t_outside = run_timed (@() pdist2_switch_outside (X, Y), n_runs, n_warmup);
    res_pdist2(idx, :) = [N, t_inside, t_outside, t_inside / t_outside];
  endfor

  print_table ("pdist: Switch Inside vs Outside", {"N", "t_inside", "t_outside", "ratio"}, res_pdist);
  print_table ("pdist2: Switch Inside vs Outside", {"N", "t_inside", "t_outside", "ratio"}, res_pdist2);

  save_csv ("bench_pdist_switch.csv", {"N", "t_inside", "t_outside", "ratio"}, res_pdist);
  save_csv ("bench_pdist2_switch.csv", {"N", "t_inside", "t_outside", "ratio"}, res_pdist2);

  make_plot ("pdist_switch.png", "pdist: Switch Inside vs Outside", N_vals, ...
             res_pdist(:,2), res_pdist(:,3), "Inside", "Outside");
  make_plot ("pdist2_switch.png", "pdist2: Switch Inside vs Outside", N_vals, ...
             res_pdist2(:,2), res_pdist2(:,3), "Inside", "Outside");
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.3fx\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_switch_inside (X)
  n = rows (X);
  n_pairs = n * (n - 1) / 2;
  Y = zeros (1, n_pairs);
  block_size = 1024;
  idx = 1;
  for i = 1:n-1
    j_vals = (i+1):n;
    n_j = numel (j_vals);
    for j_start = 1:block_size:n_j
      j_end = min (j_start + block_size - 1, n_j);
      j_block = j_vals(j_start:j_end);
      switch "euclidean"
        case "euclidean"
          diff_block = X(i, :) - X(j_block, :);
          dists = sqrt (sum (diff_block .^ 2, 2))';
      endswitch
      Y(idx:idx + numel(j_block) - 1) = dists;
      idx = idx + numel (j_block);
    endfor
  endfor
endfunction

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

function D = pdist2_switch_inside (X, Y)
  M = rows (X);
  N = rows (Y);
  D = zeros (M, N);
  block_size = 1024;
  for i_start = 1:block_size:M
    i_end = min (i_start + block_size - 1, M);
    block_rows = i_start:i_end;
    n_block = numel (block_rows);
    for bi = 1:n_block
      i = block_rows(bi);
      switch "euclidean"
        case "euclidean"
          diff_ij = X(i, :) - Y;
          D(i, :) = sqrt (sum (diff_ij .^ 2, 2))';
      endswitch
    endfor
  endfor
endfunction

function D = pdist2_switch_outside (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