plot_pareto.m

function [fig_nnz, fig_time] = plot_pareto(results, labels, varargin)
%PLOT_PARETO  Publication-quality accuracy-vs-cost plots for R/R3-Trust output.
%
%   Produces TWO figures from the same data:
%     (1) Relative Error (%) vs Relative NNZ (%)  -- accuracy vs storage/flops
%     (2) Relative Error (%) vs Time (s)          -- accuracy vs runtime
%   Each iteration of an R_TRUST / R3_TRUST run is one decomposition (one
%   point on the curve). Both figures share the same styling, legend, axis
%   cropping and export options.
%
%   plot_pareto(out) plots a single run.
%   plot_pareto({out1, out2, ...}, {'R-Trust','R3-Trust', ...}) overlays
%   several methods. The out-structs must contain the fields relative_NNZ,
%   relative_error and time (all returned by R_TRUST / R3_TRUST).
%
%   [fig_nnz, fig_time] = plot_pareto(...) returns the two figure handles
%   for any further customisation.
%
%   plot_pareto(..., 'Name', value) options:
%     'YScale'      'linear' (default) | 'log'   -- error axis (both figures).
%     'TimeScale'   'linear' (default) | 'log'   -- time axis (figure 2 only);
%                                                   'log' helps if runtimes
%                                                   span orders of magnitude.
%     'ParetoOnly'  false (default) | true       -- keep only non-dominated
%                                                   points per series.
%     'CommonRange' true (default) | false       -- crop the NNZ figure to the
%                                                   x/y range shared by ALL
%                                                   series (both sweep a
%                                                   comparable storage budget,
%                                                   so the overlap is meaningful).
%     'CommonRangeTime' false (default) | true    -- same cropping for the TIME
%                                                   figure. OFF by default: the
%                                                   methods occupy different
%                                                   runtimes, so cropping would
%                                                   hide most of the slower
%                                                   curve. Leave off to see both.
%     'Pad'         outward axis margin as a fraction (default 0.02) so that
%                   markers on the cropped edges are not clipped.
%     'Width'       figure width  in inches (default 5.0; ~3.4 for 1 column)
%     'Height'      figure height in inches (default 3.8)
%     'FontSize'    base font size (default 16)
%     'LineWidth'   curve line width (default 2)
%     'MarkerSize'  marker size (default 7)
%     'MarkerStep'  show a marker every Nth point (default 1 = all)
%     'Save'        base filename; if given, writes BOTH figures as
%                   <base>_nnz.{pdf,png,fig} and <base>_time.{pdf,png,fig}
%                   (vector PDF, 300 dpi PNG, editable MATLAB .fig).
%
%   Example:
%     plot_pareto({out_RTrust, out_R3Trust, out_ADMM}, ...
%                 {'R-Trust','R3-Trust','ADMM'}, ...
%                 'YScale','log', 'Save','fig_pareto');

% --- normalise inputs ----------------------------------------------------
if ~iscell(results), results = {results}; end
if nargin < 2 || isempty(labels), labels = {}; end
if ischar(labels), labels = {labels}; end

p = inputParser;
addParameter(p, 'YScale',      'linear');
addParameter(p, 'TimeScale',   'linear');
addParameter(p, 'ParetoOnly',  false);
addParameter(p, 'CommonRange', true);    % crop the NNZ figure to the shared range
addParameter(p, 'CommonRangeTime', false); % crop the TIME figure (off: methods have
                                           % different runtimes, show both fully)
addParameter(p, 'Pad',         0.02);
addParameter(p, 'Width',       5.0);
addParameter(p, 'Height',      3.8);
addParameter(p, 'FontSize',    16);
addParameter(p, 'LineWidth',   2);
addParameter(p, 'MarkerSize',  7);
addParameter(p, 'MarkerStep',  1);
addParameter(p, 'Save',        '');
parse(p, varargin{:});
opt = p.Results;

% --- Okabe-Ito colour-blind-safe palette + distinct markers --------------
colors = [0.00 0.45 0.70;    % blue
          0.84 0.37 0.00;    % vermillion
          0.00 0.62 0.45;    % bluish green
          0.90 0.62 0.00;    % orange
          0.80 0.47 0.65;    % reddish purple
          0.34 0.71 0.91;    % sky blue
          0.00 0.00 0.00];   % black
markers = {'o','s','^','d','v','>','<'};

% --- Figure 1: Relative Error vs Relative NNZ ----------------------------
[fig_nnz, ax_nnz] = draw_one(results, labels, opt, colors, markers, ...
    @(o) 100 * [o(:).relative_NNZ], 'Relative NNZ (\%)', 'linear', opt.CommonRange);

% --- Check whether every input series carries a usable time field --------
has_time = true;
for s = 1:numel(results)
    out = results{s};
    if ~isfield(out, 'time') || isempty([out(:).time])
        has_time = false;
        break;
    end
end

% --- Figure 2: Relative Error vs Time (only if all inputs have time) -----
if has_time
    [fig_time, ax_time] = draw_one(results, labels, opt, colors, markers, ...
        @(o) [o(:).time], 'Time (s)', opt.TimeScale, opt.CommonRangeTime);
else
    fig_time = [];
    ax_time  = [];
    warning('plot_pareto:noTime', ...
        'Skipping time vs error figure: the ''time'' field is missing in one or more inputs.');
end

% --- export both ---------------------------------------------------------
if ~isempty(opt.Save)
        % Ensure a "Figures" subfolder exists in the current directory
    fig_dir = fullfile(pwd, 'Figures_Embedding_Radiotherapy');
    if ~exist(fig_dir, 'dir')
        mkdir(fig_dir);
    end

    % Build the full base path inside the Figures folder
    save_base = fullfile(fig_dir, opt.Save);
    save_one(fig_nnz,  ax_nnz,  [save_base '_nnz'],  opt);
    if has_time
        save_one(fig_time, ax_time, [save_base '_time'], opt);
    end
end
end

% ========================================================================
% Draw one figure: Relative Error (y) vs a chosen x quantity
% ========================================================================
function [fig, ax] = draw_one(results, labels, opt, colors, markers, getx, xlab, xscale, docrop)
fig = figure('Units','inches', 'Position',[1 1 opt.Width opt.Height], ...
             'Color','w', 'PaperPositionMode','auto');
ax = axes(fig); hold(ax, 'on');

h = cell(numel(results), 1);
xrange = zeros(numel(results), 2);       % [min x, max x] per series
yrange = zeros(numel(results), 2);       % [min y, max y] per series
for s = 1:numel(results)
    out = results{s};
    x = getx(out);                       % chosen cost axis (NNZ % or time)
    y = 100 * [out(:).relative_error];   % error axis (%)

    if opt.ParetoOnly
        keep = local_pareto(x, y);
        x = x(keep); y = y(keep);
    end
    [x, ord] = sort(x); y = y(ord);      % left-to-right for a clean curve
    xrange(s,:) = [min(x) max(x)];
    yrange(s,:) = [min(y) max(y)];

    c = colors(mod(s-1, size(colors,1)) + 1, :);
    mk = markers{mod(s-1, numel(markers)) + 1};

    % --- decide line vs scatter based on monotonicity of y -----------------
    % After sorting by x, y is monotonic if it is either non-increasing or
    % non-decreasing throughout. A clean Pareto front is monotonic, so we draw
    % a connected line; otherwise the points don't form a curve, so we scatter.
    dy = diff(y);
    is_monotonic = all(dy >= 0) || all(dy <= 0);

    if is_monotonic || numel(y) < 2
        h{s} = plot(ax, x, y, ['-' mk], ...
            'Color', c, 'LineWidth', opt.LineWidth, ...
            'MarkerSize', opt.MarkerSize, 'MarkerFaceColor', c, ...
            'MarkerEdgeColor', 'w');
        if opt.MarkerStep > 1                 % thin markers (MATLAB R2016b+)
            try
                set(h{s}, 'MarkerIndices', 1:opt.MarkerStep:numel(x));
            catch
            end
        end
    else
        h{s} = scatter(ax, x, y, (opt.MarkerSize)^2, ...
            'Marker', mk, 'MarkerFaceColor', c, ...
            'MarkerEdgeColor', 'w', 'LineWidth', 1);
    end
end

% --- styling -------------------------------------------------------------
set(ax, 'YScale', opt.YScale, 'XScale', xscale);
xlabel(ax, xlab, 'Interpreter','latex');
ylabel(ax, 'Relative Error (\%)', 'Interpreter','latex');
set(ax, 'FontSize', opt.FontSize, 'TickLabelInterpreter','latex', ...
        'LineWidth', 1, 'Layer','top', 'Box','on');
grid(ax, 'on'); set(ax, 'GridAlpha', 0.12);

if ~isempty(labels)
    legend(ax, [h{:}], labels, 'Interpreter','latex', ...
           'Location','best', 'Box','off', 'FontSize', opt.FontSize-2);
end

% --- crop to the range shared by all series ------------------------------
% Avoids dead regions where only one curve exists. Limits = intersection of
% the per-series x- and y-spans, with a small outward pad for edge markers.
if docrop && numel(results) > 1
    xlo = max(xrange(:,1));  xhi = min(xrange(:,2));
    ylo = max(yrange(:,1));  yhi = min(yrange(:,2));
    if xhi > xlo
        xlim(ax, pad_range([0.85*xlo 1.15*xhi], xscale, opt.Pad));
    end
    if yhi > ylo
        ylim(ax, pad_range([0.85*ylo 1.15*yhi], opt.YScale, opt.Pad));
    end
    % xlim(ax, [0,2]);
    % ylim(ax, [0,80]);
end
end



% ========================================================================
% Export one figure as vector PDF, raster PNG and editable .fig
% ========================================================================
function save_one(fig, ax, base, opt)
if exist('exportgraphics', 'file')         % R2020a+
    exportgraphics(ax, [base '.pdf'], 'ContentType', 'vector');
    exportgraphics(ax, [base '.png'], 'Resolution', 300);
else                                       % older MATLAB fallback
    set(fig, 'PaperUnits','inches', 'PaperSize',[opt.Width opt.Height], ...
             'PaperPosition',[0 0 opt.Width opt.Height]);
    print(fig, base, '-dpdf', '-painters');
    print(fig, base, '-dpng', '-r300');
end
if exist('savefig', 'file')                % editable MATLAB figure
    savefig(fig, [base '.fig']);
else
    hgsave(fig, [base '.fig']);            % Octave fallback
end
end

% ========================================================================
% Widen [lo hi] by a small fraction so edge markers are not clipped
% ========================================================================
function lim = pad_range(rng, scale, frac)
lo = rng(1); hi = rng(2);
if strcmpi(scale, 'log')
    r = (hi / lo) ^ frac;            % multiplicative pad on a log axis
    lim = [lo / r, hi * r];
else
    d = frac * (hi - lo);            % additive pad on a linear axis
    lim = [lo - d, hi + d];
end
end

% ========================================================================
% Mask of non-dominated points (minimise both x and y)
% ========================================================================
function keep = local_pareto(x, y)
[xs, ord] = sort(x(:));
ys = y(ord);
keep_sorted = false(numel(xs), 1);
best_y = inf;
for i = 1:numel(xs)
    if ys(i) < best_y            % lower error as cost increases => on the front
        keep_sorted(i) = true;
        best_y = ys(i);
    end
end
keep = false(size(x));
keep(ord(keep_sorted)) = true;
end