function beautifyBERplot(options)
% BEAUTIFYBERPLOT  Publication-style beautification for BER / NSD plots.
% Style aligned to manual plotting script (lw=0.8, ms=4, white markers).
%
% Usage examples:
%   beautifyBERplot
%   beautifyBERplot("polyfit",true,"polyorder",1)
%   beautifyBERplot("polyfit",true,"fitmethod","pchip")

arguments
    options.logscale   (1,1) logical = 1
    options.setmarkers (1,1) logical = 1
    options.setcolors (1,1) logical = 1
    options.polyfit    (1,1) logical = 0
    options.polyorder  (1,1) double  = 2
    options.fitmethod  (1,1) string  = "polyfit"
end

ax = gca;

% ---------------------------------------------------------
% Extract real data objects from ax.Children
% ---------------------------------------------------------
children = ax.Children;
isLine   = arrayfun(@(h) isa(h,'matlab.graphics.chart.primitive.Line'), children);
isScatter = arrayfun(@(h) isa(h,'matlab.graphics.chart.primitive.Scatter'), children);

% Remove helper lines (e.g. yline, fit overlays)
isConstantLine = arrayfun(@(h) isprop(h,'Tag') && strcmp(h.Tag,'ConstantLine'), children);

dataObjects = children((isLine | isScatter) & ~isConstantLine);
n = numel(dataObjects);

% ---------------------------------------------------------
% Fixed color palette (deterministic across figures)
% ---------------------------------------------------------
if n > 0
    try
        cmap = linspecer(n);
    catch
        cmap = linespecer_fallback(n);
    end
end

markers = {'o','s','o','o','^','v','d','>'};
lw = 0.8;
ms = 2;

% ---------------------------------------------------------
% Apply line/scatter + marker styling
% ---------------------------------------------------------
for i = 1:n
    dataObj = dataObjects(n-i+1);   % preserve plotting order

    if isa(dataObj,'matlab.graphics.chart.primitive.Line')
        dataObj.LineWidth = lw;

        if options.setcolors
            dataObj.Color = cmap(i,:);
        end

        if options.setmarkers
            if strcmp(dataObj.Marker,'none')
                dataObj.Marker = markers{mod(i-1,numel(markers))+1};
            end
            dataObj.MarkerSize = ms;
            if options.setcolors
                dataObj.MarkerEdgeColor = cmap(i,:);
            end
            dataObj.MarkerFaceColor = [1 1 1];
        end
    elseif isa(dataObj,'matlab.graphics.chart.primitive.Scatter')
        if options.setcolors
            dataObj.CData = cmap(i,:);
            dataObj.MarkerEdgeColor = cmap(i,:);
            markerFaceColor = dataObj.MarkerFaceColor;
            hasNoFace = ischar(markerFaceColor) && strcmp(markerFaceColor,'none');
            if ~hasNoFace
                dataObj.MarkerFaceColor = cmap(i,:);
            end
        end

        if options.setmarkers
            if strcmp(dataObj.Marker,'none')
                dataObj.Marker = markers{mod(i-1,numel(markers))+1};
            end
            dataObj.SizeData = ms^2;
        end
    end
end

% ---------------------------------------------------------
% Optional smoothing / fitting overlay
% ---------------------------------------------------------
if options.polyfit && n > 0
    hold on
    for i = 1:n
        dataObj = dataObjects(n-i+1);
        x = dataObj.XData(:);
        y = dataObj.YData(:);
        valid = isfinite(x) & isfinite(y);

        if sum(valid) < options.polyorder+1
            continue
        end

        xf = linspace(min(x(valid)),max(x(valid)),200);

        switch lower(options.fitmethod)
            case "polyfit"
                p  = polyfit(x(valid),y(valid),options.polyorder);
                yf = polyval(p,xf);

            case "smoothingspline"
                try
                    f  = fit(x(valid),y(valid),'smoothingspline');
                    yf = feval(f,xf);
                catch
                    yf = interp1(x(valid),y(valid),xf,'pchip');
                end

            case "loess"
                ys = smooth(x(valid),y(valid),0.2,'loess');
                yf = interp1(x(valid),ys,xf,'linear','extrap');

            case "pchip"
                yf = interp1(x(valid),y(valid),xf,'pchip');

            otherwise
                p  = polyfit(x(valid),y(valid),options.polyorder);
                yf = polyval(p,xf);
        end

        if isa(dataObj,'matlab.graphics.chart.primitive.Line')
            basecol = dataObj.Color;
        else
            basecol = dataObj.CData(1,:);
            if numel(basecol) ~= 3
                basecol = cmap(i,:);
            end
        end
        lightcol = basecol + 0.4*(1-basecol);
        lightcol(lightcol>1)=1;

        plot(xf,yf,'-','Color',lightcol,...
            'LineWidth',0.7,'Marker','none',...
            'HandleVisibility','off');
    end
    hold off
end

% ---------------------------------------------------------
% Axes formatting (matches first script)
% ---------------------------------------------------------
if options.logscale
    set(ax,'YScale','log')
end

grid(ax,'on')
set(ax,'GridLineStyle','--',...
       'GridLineWidth',0.5,...
       'MinorGridLineWidth',0.5)

set(ax,'FontSize',12,...
       'Box','on',...
       'LineWidth',0.8)

set(findall(ax,'-property','Interpreter'),'Interpreter','latex')
set(gcf,'Color','w')

end

% ---------------------------------------------------------
% Fallback palette (only used if linspecer is missing)
% ---------------------------------------------------------
function cmap = linespecer_fallback(n)
cmap = lines(n);
end