function plot_FEC_violin(res, varargin)
% plot_FEC_violin  Violin plot of ROP-at-FEC crossings per wavelength/channel
% Uses the Bechtold community violinplot (MATLAB Central 45134).
%
% Requirements:
% - community violinplot must be on path, either as:
%     (a) violinplot.m   (shadows MATLAB official), OR
%     (b) violinplot_community.m (renamed + function name adjusted)
%
% Usage:
%   plot_FEC_violin(res, 'tech','VNLE', 'fec',3.8e-3, 'eval_ptr',[], 'ylim',[-10 -6], 'bandwidth',0.15);

% -------- args --------
p = inputParser;
p.addParameter('tech',      'VNLE', @(x)ischar(x)||isstring(x));
p.addParameter('fec',       3.8e-3, @(x)isnumeric(x)&&isscalar(x));
p.addParameter('eval_ptr',  [],     @(x)isnumeric(x)&&isscalar(x));
p.addParameter('ylim',      [],     @(x)isnumeric(x)&&numel(x)==2);
p.addParameter('bandwidth', [],     @(x)isnumeric(x)&&isscalar(x));
p.parse(varargin{:});
opt = p.Results;

tech = upper(string(opt.tech));
fec  = opt.fec;

% -------- pick result field --------
switch tech
    case "FFE",  C4 = res.ffe;
    case "DFE",  C4 = res.dfe;
    case "VNLE", C4 = res.vnle;
    case "MLSE", C4 = res.mlse;
    case "DBT",  C4 = res.dbt;
    otherwise, error('Unknown tech "%s". Use FFE/DFE/VNLE/MLSE/DBT.', tech);
end

rop            = res.settings.rop(:);
wavelengthplan = res.settings.wavelengthplan(:);
distances      = res.eval_dist_km(:);

N_ch  = numel(wavelengthplan);
N_rop = numel(rop);

dims = size(C4);
if numel(dims) < 4, dims(end+1:4) = 1; end
N_dist = dims(4);

eval_ptr = opt.eval_ptr;
if isempty(eval_ptr), eval_ptr = N_dist; end
eval_ptr = max(1, min(N_dist, eval_ptr));

% -------- compute crossings per channel + "completeness" --------
S_cell   = cell(1, N_ch);  % crossings values (finite only)
K_total  = zeros(1, N_ch); % number of complete realizations for that channel
K_cross  = zeros(1, N_ch); % number of realizations that crossed

for ch = 1:N_ch
    cells = sliceCells4D(C4, ch, eval_ptr, N_rop); % N_rop x N_realiz (may contain [])
    [S_cell{ch}, K_total(ch), K_cross(ch)] = crossings_and_counts(rop, cells, fec);
end

% -------- build vector+category representation (NO NaNs in xAll) --------
catLabels = arrayfun(@(nm)sprintf('%d nm', round(nm)), wavelengthplan, 'UniformOutput', false);
catsAll = categorical(strings(0,1), catLabels, 'Ordinal', false);
xAll    = zeros(0,1);

for ch = 1:N_ch
    v = S_cell{ch}(:);
    if ~isempty(v)
        xAll    = [xAll; v];
        catsAll = [catsAll; repmat( ...
            categorical(string(catLabels{ch}), catLabels, 'Ordinal', false), ...
            numel(v), 1)];
    end
end


% -------- plot --------
figure('Name', sprintf('%s FEC violin @ %s', tech, distLabel(distances, eval_ptr)));
hold on;

% Violin only if we have any data at all
if ~isempty(xAll)
    violinplot_community(xAll, catsAll,'Bandwidth',    0.15, ...                 % KDE bandwidth (critical!)
    'ViolinColor',  [0.2 0.4 0.8], ...         % or Nx3 for per-channel colors
    'ViolinAlpha',  0.15, ...
    'EdgeColor',    [0.3 0.3 0.3], ...
    'MarkerSize',   14, ...
    'ShowData',     true, ...
    'ShowMean',     false, ...
    'ShowMedian',   true, ...
    'ShowBox',      false, ...
    'ShowWhiskers', false);
else
    warning('No FEC crossings found at eval_ptr=%d (%s). Plotting only markers.', eval_ptr, distLabel(distances, eval_ptr));
    set(gca,'XTick',1:N_ch,'XTickLabel',catLabels);
end

% Marker X positions (violin groups are at 1..N_ch)
xpos = 1:N_ch;

% -------- overlay mean markers (black/red/blue) --------
yBlue = -9;

for ch = 1:N_ch
    if K_total(ch) == 0 || K_cross(ch) == 0
        % no complete realizations OR none crossed -> blue X at -9
        scatter(xpos(ch), yBlue, 80, 'x', 'LineWidth', 2, 'MarkerEdgeColor', [0 0 1], 'HandleVisibility','off');
        continue;
    end

    mu = mean(S_cell{ch}, 'omitnan');

    if K_cross(ch) < K_total(ch)
        % some complete realizations did not cross -> red X
        scatter(xpos(ch), mu, 80, 'x', 'LineWidth', 2, 'MarkerEdgeColor', [1 0 0], 'HandleVisibility','off');
    else
        % all complete realizations crossed -> black X
        scatter(xpos(ch), mu, 80, 'x', 'LineWidth', 2, 'MarkerEdgeColor', [0 0 0], 'HandleVisibility','off');
    end
end

title(sprintf('%s | BER %.2e | %s', tech, fec, distLabel(distances, eval_ptr)));
ylabel('ROP at FEC crossing');
grid on; box on;

if ~isempty(opt.ylim)
    ylim(opt.ylim);
end

end


% ================= helpers =================

function cells2D = sliceCells4D(C4, ch, eval_ptr, N_rop)
dims = size(C4);
if numel(dims) < 4, dims(end+1:4) = 1; end
N_realiz = dims(3);

cells2D = cell(N_rop, N_realiz);

if ch > dims(1) || eval_ptr > dims(4)
    return;
end

ropUse = min(N_rop, dims(2));
rUse   = min(N_realiz, dims(3));

tmp = squeeze(C4(ch, 1:ropUse, 1:rUse, eval_ptr));
tmp = reshape(tmp, ropUse, rUse);
cells2D(1:ropUse, 1:rUse) = tmp;
end

function [S, K_total, K_cross] = crossings_and_counts(rop, cellsNxR, fec)
% counts "complete" realizations and how many of those crossed
% S returns only finite crossings (no NaN).

Y = extractCompleteBER(cellsNxR);   % N_rop x K_total
K_total = size(Y,2);

if K_total == 0
    S = [];
    K_cross = 0;
    return;
end

% optional quality gate (keep your style)
ok = mean(Y,1,'omitnan') <= 0.1;
Y  = Y(:,ok);
K_total = size(Y,2);

if K_total == 0
    S = [];
    K_cross = 0;
    return;
end

rop = rop(:);
S = nan(1, K_total);

for j = 1:K_total
    y = Y(:,j);
    above = (y > fec);
    idx = find(above(1:end-1) & ~above(2:end), 1, 'first');
    if ~isempty(idx)
        x1 = rop(idx);   y1 = y(idx);
        x2 = rop(idx+1); y2 = y(idx+1);
        if isfinite(y1) && isfinite(y2) && (y2 ~= y1)
            t    = (fec - y1) / (y2 - y1);
            S(j) = x1 + t*(x2 - x1);
        end
    end
end

K_cross = sum(isfinite(S));
S = S(isfinite(S));
end

function Y = extractCompleteBER(cellSlice)
if isempty(cellSlice), Y = []; return; end
nR = size(cellSlice,2);
keep = false(1,nR);

for r = 1:nR
    col = cellSlice(:,r);
    keep(r) = all(cellfun(@(c) ~isempty(c), col));
end

if ~any(keep), Y = []; return; end
Y = cellfun(@(c) c.metrics.BER, cellSlice(:,keep), 'UniformOutput', true);
end

function lbl = distLabel(distances, eval_ptr)
if eval_ptr <= numel(distances) && isfinite(distances(eval_ptr))
    lbl = sprintf('%.0f km', distances(eval_ptr));
else
    lbl = sprintf('eval_%d', eval_ptr);
end
end