% Driver for the one-dimensional labs.
%
% The editable script supplies two functions:
%
%     names = method_names()                 a cell array of method labels
%     [x, u] = solve1d(mid, N, f, ua, ub)    method mid, resolution N
%
% solve1d picks its own grid x (uniform, graded, Chebyshev, whatever the method
% wants), returns the solution values on that grid, and the driver measures the
% error against the manufactured solution at exactly those points. f is a
% function handle giving the source term, so a method is free to evaluate it
% anywhere.
%
% The refinement sweep is written as a recursion rather than a for loop. That
% is not a stylistic choice: numbl compiles a call made from inside a loop
% through a generic path whenever the callee itself contains loops that call
% further functions, and the solves being timed here are exactly of that shape.
% The penalty is a factor of five to twenty, it falls unevenly on different
% methods, and it would therefore corrupt every timing in this tutorial.
% Recursion is compiled normally. See CLAUDE.md.
p = lp_params();
kind = p.kind;
names = method_names();
nm = numel(names);
ub = lp_exact1d(kind, [0; 1]);
f = @(xx) drv_rhs(kind, xx);

out = struct();
out.names = names;
out.kind = kind;

if strcmp(p.mode, 'single')
    t0 = tic;
    [x, u] = solve1d(p.mid, p.N, f, ub(1), ub(2));
    out.ms = toc(t0) * 1000;
    uex = lp_exact1d(kind, x(:));
    err = u(:) - uex(:);
    out.x = x(:)';
    out.u = u(:)';
    out.uex = uex(:)';
    out.err = err';
    out.maxerr = max(abs(err));
    out.n = numel(x) - 1;
    xf = linspace(0, 1, p.ngrid)';
    out.xf = xf';
    out.uf = lp_exact1d(kind, xf)';
else
    ns = p.ns(:)';
    nn = numel(ns);
    E = zeros(nm, nn);
    T = zeros(nm, nn);
    Q = zeros(nm, nn);
    nmax = p.nmaxs(:)';
    timed = p.timed;
    ET = drv_methods(1, NaN(2*nm, nn), ns, f, ub(1), ub(2), kind, timed, nmax, nm);
    for k = 1:nm
        E(k, :) = ET(2*k-1, :);
        T(k, :) = ET(2*k, :);
        Q(k, :) = lp_orders(E(k, :), ns);
    end
    out.ns = ns;
    out.E = lp_rows(E);
    out.T = lp_rows(T);
    out.orders = lp_rows(Q);
end

lp_write(out);

function ET = drv_methods(k, ET, ns, f, ua, ub, kind, timed, nmax, nm)
% One method per level of the recursion; row 2k-1 of ET collects its errors and
% row 2k its milliseconds.
if k > nm
    return;
end
% a warm-up call, so that the first timed solve is not also paying for the
% compilation of the method
[xw, uw] = solve1d(k, 8, f, ua, ub);
R = drv_res(1, NaN(2, numel(ns)), k, ns, f, ua, ub, kind, timed, nmax(k));
ET(2*k-1, :) = R(1, :);
ET(2*k, :) = R(2, :);
ET = drv_methods(k+1, ET, ns, f, ua, ub, kind, timed, nmax, nm);
end

function R = drv_res(j, R, mid, ns, f, ua, ub, kind, timed, nmax)
% One resolution per level of the recursion.
if j > numel(ns)
    return;
end
N = ns(j);
if N <= nmax
    t0 = tic;
    [x, u] = solve1d(mid, N, f, ua, ub);
    el = toc(t0);
    if timed > 0 && el < 0.005
        % a sub-millisecond solve says more about the clock than about the
        % method, so take the best of a few repeats
        el = drv_best(4, el, mid, N, f, ua, ub);
    end
    uex = lp_exact1d(kind, x(:));
    R(1, j) = max(abs(u(:) - uex(:)));
    R(2, j) = el * 1000;
end
R = drv_res(j+1, R, mid, ns, f, ua, ub, kind, timed, nmax);
end

function el = drv_best(k, el, mid, N, f, ua, ub)
if k == 0
    return;
end
t0 = tic;
[x, u] = solve1d(mid, N, f, ua, ub);
e2 = toc(t0);
if e2 < el
    el = e2;
end
el = drv_best(k-1, el, mid, N, f, ua, ub);
end

function y = drv_rhs(kind, x)
[~, y] = lp_exact1d(kind, x);
end