% Offline benchmark of the Pac-Man problem in real MATLAB.
%
% Runs the tutorial's own MATLAB layer, unchanged: the chapter driver and the
% editable script are concatenated into main.m exactly as the page does, and
% executed here instead of in numbl. Only the two best methods from the
% in-browser benchmark are swept (the harmonic basis and the graded-panel
% double layer), with a deeper panel ladder than the page uses, since a native
% LAPACK has no trouble with the larger dense systems.
%
% Because this runs the chapter's own script, whatever that script does badly
% it does badly here too, which is a virtue for comparing MATLAB against numbl
% and a hazard for comparing MATLAB against the Python and Julia ports. The
% shape the script uses, a scalar double loop for the kernel matrix, happens to
% be the fastest of the three shapes available in both MATLAB and numbl, so
% nothing is being given away; if that ever stops being true, the honest fix is
% to change the chapter, not to keep a faster copy here.
%
%   matlab -batch "run('bench/matlab/bench_pacman.m')"
%
% from the repository root writes bench-pacman-matlab.json into the
% repository root (run() changes the working directory while it executes, so
% "the current directory" would be a moving target). Load that file from the
% benchmark chapter's "load offline results" button; because it ran on the
% same machine as your browser, the timings share a clock and can share an
% axis.

here = fileparts(mfilename('fullpath'));
root = fileparts(fileparts(here));
outfile = fullfile(root, 'bench-pacman-matlab.json');

% the deep basis ladders are legitimately rank deficient at the 1e-12 level;
% Chapter 9 discusses why, and the warning would repeat once per solve
warning('off', 'MATLAB:rankDeficientMatrix');

kind = 'both';
% ladders: zeros skip a method (rows: SW+SOR, basis, panels, Chebyshev)
nsm = [0, 0, 0, 0, 0, 0;
       4, 8, 16, 32, 48, 0;
       4, 8, 12, 16, 20, 24;
       0, 0, 0, 0, 0, 0];

work = tempname;
mkdir(work);
libdir = fullfile(root, 'src', 'matlab');
libfiles = dir(fullfile(libdir, '*.m'));
for k = 1:numel(libfiles)
    copyfile(fullfile(libdir, libfiles(k).name), work);
end
driver = fileread(fullfile(root, 'src', 'chapters', 'ch11', 'driver.m'));
script = fileread(fullfile(root, 'src', 'chapters', 'ch11', 'pacman.m'));
fid = fopen(fullfile(work, 'main.m'), 'w');
fwrite(fid, sprintf('%s\n%s', driver, script));
fclose(fid);

params = struct('kind', kind, 'mode', 'sweep', 'nsm', nsm);
fid = fopen(fullfile(work, 'params.json'), 'w');
fwrite(fid, jsonencode(params));
fclose(fid);

old = cd(work);
cleanup = onCleanup(@() cd(old));
run('main.m');
res = jsondecode(fileread('out.json'));
cd(old);

% out.json holds one row per method; keep the two that actually ran
slots = {'fd2', 'mps', 'bie', 'spectral'};
keep = [2, 3];
methods = cell(1, numel(keep));
for q = 1:numel(keep)
    k = keep(q);
    ns = pick_row(res.NS, k);
    ok = isfinite(ns) & ns > 0;
    m = struct();
    m.name = res.names{k};
    m.slot = slots{k};
    m.ns = ns(ok);
    m.E = pick_row(res.E, k);
    m.E = m.E(ok);
    m.T = pick_row(res.T, k);
    m.T = m.T(ok);
    methods{q} = m;
end

out = struct();
out.schema = 'laplace-dirichlet-bench-1';
out.label = sprintf('MATLAB %s, %d threads', version('-release'), maxNumCompThreads);
out.kind = kind;
out.date = char(datetime('now', 'Format', 'yyyy-MM-dd'));
out.methods = methods;

fid = fopen(outfile, 'w');
fwrite(fid, jsonencode(out));
fclose(fid);
fprintf('wrote %s\n', outfile);

function r = pick_row(M, k)
% jsondecode gives a matrix when the rows are rectangular and a cell of
% columns when they are not; accept both.
if iscell(M)
    r = M{k}(:)';
else
    r = M(k, :);
end
end