PCA Benchmarking

## Benchmark script for pca.m vectorization optimization
## Compares pca_original.m vs optimized pca.m

addpath("inst");
addpath("src");

## Configuration
num_runs = 10;

## Test configurations: [n, d]
test_configs = [
  1000, 20;
  2000, 30;
  5000, 50;
  8000, 50;
];

printf("=== PCA Benchmark ===\n\n");

## Initialize results
results = {};
result_idx = 1;

for cfg_idx = 1:rows(test_configs)
  n = test_configs(cfg_idx, 1);
  d = test_configs(cfg_idx, 2);

  printf("Dataset: n=%d, d=%d\n", n, d);

  ## Generate test data with 10% missing values
  rand("seed", 42);
  randn("seed", 42);
  X = randn(n, d);
  missing_count = round(0.1 * numel(X));
  missing_idx = randperm(numel(X), missing_count);
  X(missing_idx) = NaN;

  ## Also generate weights for weighted test
  weights = rand(n, 1) + 0.5;  # positive weights

  ## Test 1: With missing values (rows='pairwise') - coeff only
  printf("  Testing with missing values (pairwise)...\n");

  original_times = zeros(1, num_runs);
  for run = 1:num_runs
    tic;
    coeff1 = pca_original(X, "Algorithm", "eig", "Rows", "pairwise");
    original_times(run) = toc;
  endfor
  original_mean = mean(original_times);

  optimized_times = zeros(1, num_runs);
  for run = 1:num_runs
    tic;
    coeff2 = pca(X, "Algorithm", "eig", "Rows", "pairwise");
    optimized_times(run) = toc;
  endfor
  optimized_mean = mean(optimized_times);

  speedup = original_mean / optimized_mean;

  ## Validate correctness
  coeff_err = norm(coeff1 - coeff2);
  if (coeff_err > 1e-10)
    printf("    WARNING: Results differ! coeff_err=%e\n", coeff_err);
  endif

  printf("    original=%.4fs, optimized=%.4fs, speedup=%.2fx\n", ...
         original_mean, optimized_mean, speedup);

  results{result_idx, 1} = n;
  results{result_idx, 2} = d;
  results{result_idx, 3} = "missing_pairwise";
  results{result_idx, 4} = original_mean;
  results{result_idx, 5} = optimized_mean;
  results{result_idx, 6} = speedup;
  result_idx += 1;

  ## Test 2: With missing values and weights - coeff only
  printf("  Testing with missing values + weights...\n");

  original_times = zeros(1, num_runs);
  for run = 1:num_runs
    tic;
    coeff1 = pca_original(X, "Algorithm", "eig", "Rows", "pairwise", "Weights", weights);
    original_times(run) = toc;
  endfor
  original_mean = mean(original_times);

  optimized_times = zeros(1, num_runs);
  for run = 1:num_runs
    tic;
    coeff2 = pca(X, "Algorithm", "eig", "Rows", "pairwise", "Weights", weights);
    optimized_times(run) = toc;
  endfor
  optimized_mean = mean(optimized_times);

  speedup = original_mean / optimized_mean;

  printf("    original=%.4fs, optimized=%.4fs, speedup=%.2fx\n", ...
         original_mean, optimized_mean, speedup);

  results{result_idx, 1} = n;
  results{result_idx, 2} = d;
  results{result_idx, 3} = "missing_weighted";
  results{result_idx, 4} = original_mean;
  results{result_idx, 5} = optimized_mean;
  results{result_idx, 6} = speedup;
  result_idx += 1;

  ## Test 3: Without missing values (baseline)
  printf("  Testing without missing values...\n");

  X_clean = randn(n, d);

  original_times = zeros(1, num_runs);
  for run = 1:num_runs
    tic;
    [coeff1, score1, latent1] = pca_original(X_clean);
    original_times(run) = toc;
  endfor
  original_mean = mean(original_times);

  optimized_times = zeros(1, num_runs);
  for run = 1:num_runs
    tic;
    [coeff2, score2, latent2] = pca(X_clean);
    optimized_times(run) = toc;
  endfor
  optimized_mean = mean(optimized_times);

  speedup = original_mean / optimized_mean;

  printf("    original=%.4fs, optimized=%.4fs, speedup=%.2fx\n", ...
         original_mean, optimized_mean, speedup);

  results{result_idx, 1} = n;
  results{result_idx, 2} = d;
  results{result_idx, 3} = "no_missing";
  results{result_idx, 4} = original_mean;
  results{result_idx, 5} = optimized_mean;
  results{result_idx, 6} = speedup;
  result_idx += 1;

  printf("\n");
endfor

## Save results to CSV
fid = fopen("pca_benchmark_results.csv", "w");
fprintf(fid, "n,d,test_type,original_time,optimized_time,speedup\n");
for i = 1:size(results, 1)
  fprintf(fid, "%d,%d,%s,%.6f,%.6f,%.2f\n", ...
          results{i, 1}, results{i, 2}, results{i, 3}, ...
          results{i, 4}, results{i, 5}, results{i, 6});
endfor
fclose(fid);

printf("Results saved to pca_benchmark_results.csv\n\n");

## Summary
speedups = cell2mat(results(:, 6));
printf("=== Summary ===\n");
printf("Average speedup: %.2fx\n", mean(speedups));
printf("Max speedup: %.2fx\n", max(speedups));
printf("Min speedup: %.2fx\n", min(speedups));