% Driver for the benchmark chapter: one hard problem, every method.
%
% The editable script supplies
%
%     names = method_names()
%     v = solve_pac(mid, n, g, px, py)
%     d = dof_pac(mid, n)
%
% and must return the computed solution at the points (px, py), whatever grid,
% basis, or boundary discretisation it used internally. n is the method's own
% resolution parameter: intervals per unit length for the difference scheme,
% basis terms for the expansion, refinement depth for the panels, points per
% direction for the collocation. dof_pac reports the number of unknowns behind
% that parameter, which the page shows so that the parameters stay comparable.
%
% g(x, y) is the Dirichlet data: the exact solution restricted to the
% boundary. The domain is the Pac-Man
%
%     0 < r < 1,   0 < theta < 3*pi/2,
%
% and the exact solutions are those of the L-shaped chapter (lp_lexact), whose
% corner sits at the same origin with the same interior angle.
p = lp_params();
kind = p.kind;
names = method_names();
nm = numel(names);
g = @(xx, yy) lp_lexact(kind, xx, yy);

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

% the boundary, drawn over the field pictures: the edge theta = 0, the arc,
% and the edge theta = 3*pi/2, walked as one closed polyline
ta = (0:120)' * (3*pi/2) / 120;
e1 = linspace(0, 1, 15)';
e2 = linspace(-1, 0, 15)';
out.cx = [e1; cos(ta); zeros(15, 1)]';
out.cy = [zeros(15, 1); sin(ta); e2]';

if strcmp(p.mode, 'field')
    % a picture on a fixed grid of spacing 1/32, masked to the open domain;
    % pixels on the two straight edges are left out so that the boundary
    % methods are never evaluated exactly on the curve they integrate over
    m = 34;
    xs = (-m:m)' / 32;
    [X, Y] = meshgrid(xs, xs);
    R2 = X.^2 + Y.^2;
    onedge = (abs(Y) < 1e-12 & X > -1e-12) | (abs(X) < 1e-12 & Y < 1e-12);
    inside = R2 < 1 - 1e-9 & ~(X > 1e-12 & Y < -1e-12) & ~onedge;
    px = X(inside);
    py = Y(inside);
    t0 = tic;
    v = solve_pac(p.mid, p.n, g, px, py);
    out.ms = toc(t0) * 1000;
    ve = g(px, py);
    out.maxerr = max(abs(v - ve));
    out.unknowns = dof_pac(p.mid, p.n);

    U = NaN(2*m+1, 2*m+1);
    E = NaN(2*m+1, 2*m+1);
    U(inside) = v;
    E(inside) = abs(v - ve);
    out.x = xs';
    out.y = xs';
    out.U = lp_rows(U);
    out.Err = lp_rows(E);
else
    % errors at the nodes of a fixed lattice of spacing 1/8, kept at least
    % 0.15 from the arc and 0.25 from the straight edges, and reaching within
    % 0.125 of the re-entrant corner (but excluding the corner itself, which
    % is a boundary point)
    [X, Y] = meshgrid((-8:8)' / 8);
    R2 = X.^2 + Y.^2;
    keep = R2 <= 0.7225 + 1e-9 & R2 > 1e-12 & ~(X > 1e-12 & Y < -1e-12) ...
        & ~(X > 1e-12 & Y < 0.25 - 1e-9) & ~(Y < -1e-12 & X > -0.25 + 1e-9);
    px = X(keep);
    py = Y(keep);
    ve = g(px, py);

    % Each method gets its own list of resolutions. numbl's jsondecode turns
    % a JSON array of arrays into a cell of column vectors where MATLAB gives
    % a matrix; the offline harness runs this same file in real MATLAB, so
    % both shapes are accepted.
    NSc = p.nsm;
    if iscell(NSc)
        nn = numel(NSc{1});
        NS = zeros(nm, nn);
        for k = 1:nm
            row = NSc{k};
            NS(k, :) = row(:)';
        end
    else
        NS = NSc;
        nn = size(NS, 2);
    end
    E = NaN(nm, nn);
    T = NaN(nm, nn);
    Q = NaN(nm, nn);
    DOF = NaN(nm, nn);
    for k = 1:nm
        % a ladder of zeros skips the method entirely, which is how the
        % offline harness runs a subset of the methods through this driver
        nw = 0;
        for j = 1:nn
            if isfinite(NS(k, j)) && NS(k, j) > 0
                nw = NS(k, j);
                break;
            end
        end
        if nw == 0
            continue;
        end
        vw = solve_pac(k, nw, g, px(1), py(1));   % warm-up
        for j = 1:nn
            n = NS(k, j);
            if ~isfinite(n) || n <= 0
                continue;
            end
            % Keep the fastest of several runs. Three runs at least whenever a solve
            % is cheap enough to allow it, many more when it is very cheap, and one
            % when it is expensive: the budget below caps the wall clock per point.
            % The minimum over runs is what makes these numbers comparable at all,
            % here and in the offline scripts, which apply the same rule. A single
            % measurement of a fifty-millisecond solve reports whatever else the
            % machine happened to be doing at the time.
            el = Inf;
            tot = 0;
            rep = 0;
            while rep == 0 || (rep < 50 && tot < 0.6 && (rep < 3 || tot < 0.1))
                t0 = tic;
                v = solve_pac(k, n, g, px, py);
                dt = toc(t0);
                el = min(el, dt);
                tot = tot + dt;
                rep = rep + 1;
            end
            T(k, j) = el * 1000;
            E(k, j) = max(abs(v - ve));
            DOF(k, j) = dof_pac(k, n);
        end
        Q(k, :) = lp_orders(E(k, :), NS(k, :));
    end
    out.ns = NS(1, :);
    out.NS = lp_rows(NS);
    out.DOF = lp_rows(DOF);
    out.E = lp_rows(E);
    out.T = lp_rows(T);
    out.orders = lp_rows(Q);
end

lp_write(out);