function plot_FEC_violin(res, varargin)
% plot_fec_violin_from_res  Self-contained violin plot (no dependency on other routines)
%
% Usage:
%   plot_fec_violin_from_res('WDM_20260106_....mat', 'tech', 'FFE', 'fec', 3.8e-3, 'eval_ptr', 2, 'ylim', [-10 -6]);
%   plot_fec_violin_from_res(res,                       'tech', 'FFE', 'fec', 3.8e-3, 'eval_ptr', 2);
%
% Notes:
% - Works with partially-filled res.* cell arrays (missing eval distances / missing realizations).
% - Uses only complete realizations per channel (all ROP points present).
% - Requires violinplot.m on path.

p = inputParser;
p.addParameter('tech', 'FFE', @(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('title_prefix', 'ROP at FEC crossing', @(x)ischar(x)||isstring(x));
p.parse(varargin{:});

if exist('violinplot','file') ~= 2
    error('violinplot.m not found on path. Add it to path first.');
end


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

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

% Determine distances dim robustly from the stored cell arrays (prefer ffe)
dims = size(res.ffe);
if numel(dims) < 4, dims(end+1:4) = 1; end
N_distances = dims(4);
if numel(distances) ~= N_distances
    distances = (1:N_distances).';
end

% Choose eval distance
eval_ptr = p.Results.eval_ptr;
if isempty(eval_ptr)
    eval_ptr = N_distances; % default: last available
end
if eval_ptr < 1 || eval_ptr > N_distances
    error('eval_ptr=%d out of range. Available: 1..%d', eval_ptr, N_distances);
end

% Choose technique field
tech = upper(string(p.Results.tech));
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 one of: FFE, DFE, VNLE, MLSE, DBT.', tech);
end

fec = p.Results.fec;

% ---------- Compute crossings per channel ----------

% ---------- Plot ----------
figure('Name', sprintf('%s violin @ %s', tech, distLabel(distances, eval_ptr)));
hold on;
c = linspecer(N_distances);
for eval_ptr = 1:N_distances
    S_cell = cell(1, N_ch);
    
    for ch = 1:N_ch
        cells = sliceCells4D(C4, ch, eval_ptr, N_rop);    % N_rop x N_realiz cell
        S_cell{ch} = fecCrossings_only(rop, cells, fec); % 1 x K crossings (may be empty)
    end
    
    % Pad to matrix for violinplot: rows=realizations (max K), cols=channels
    Kmax = max(cellfun(@numel, S_cell));
    if isempty(Kmax) || Kmax == 0
        warning('No FEC crossings found for %s at eval_ptr=%d (%s).', tech, eval_ptr, distLabel(distances, eval_ptr));
        return;
    end
    
    S_mat = NaN(Kmax, N_ch);
    for ch = 1:N_ch
        k = numel(S_cell{ch});
        if k > 0
            S_mat(1:k, ch) = S_cell{ch}(:);
        end
    end
    
    catLabels = arrayfun(@(nm) sprintf('%d nm', nm), round(wavelengthplan), 'UniformOutput', false);

    if numel(S_mat)>numel(wavelengthplan)
        S_mat(isnan(S_mat)) = 0;
        v=violinplot(S_mat, catLabels, ...
            'ViolinColor', c(eval_ptr,:), ...
            'ViolinAlpha', 0.10, ...
            'MarkerSize', 20, ...
            'ShowMedian', true, ...
            'EdgeColor', c(eval_ptr,:), ...
            'ShowWhiskers', false, ...
            'ShowData', true, ...
            'ShowBox', false, ...
            'Bandwidth', 0.05);
    else
        channs = 1:numel(wavelengthplan);
        scatter( wavelengthplan ,S_mat,10,c(eval_ptr,:),'Marker','o','LineWidth',2,'DisplayName',distLabel(distances, eval_ptr));
    
        channs = 1:numel(wavelengthplan);
        scatter( round(wavelengthplan(S_mat==0)),S_mat(S_mat==0),10,c(eval_ptr,:),'Marker','o','HandleVisibility','off');
    end

end

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

if ~isempty(p.Results.ylim)
    ylim(p.Results.ylim);
end

end


%% ================= helpers (self-contained) =================

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);

chMax   = dims(1);
ropMax  = dims(2);
rMax    = dims(3);
eMax    = dims(4);

if ch > chMax || eval_ptr > eMax
    return;
end

ropUse = min(N_rop, ropMax);
rUse   = min(N_realiz, rMax);

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

cells2D(1:ropUse, 1:rUse) = tmp;
end

function lbl = distLabel(distances, eval_ptr)
if eval_ptr <= numel(distances)
    d = distances(eval_ptr);
    if isfinite(d)
        lbl = sprintf('%.0f km', d);
        return;
    end
end
lbl = sprintf('eval\\_%d', eval_ptr);
end

function S = fecCrossings_only(rop, cellsNxR, fec)
% Returns 1xK crossing positions for all complete realizations (drops bad ones).
Y = extractCompleteBER(cellsNxR);  % N_rop x K
if isempty(Y), S = []; return; end

% Optional quality gate (same as your previous scripts)
ok = mean(Y,1,'omitnan') <= 0.1;
Y  = Y(:, ok);
if isempty(Y), S = []; return; end

nR = size(Y,2);
S  = nan(1,nR);

rop = rop(:);
for j = 1:nR
    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

S = S(isfinite(S)); % keep only successful crossings
end

function Y = extractCompleteBER(cellSlice)
% Keep only those realization columns where ALL N_rop entries are non-empty
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