% === 400G DSP settings ===
dsp_options = struct();
dsp_options.mode = "run_id";
dsp_options.recipe = @dsp_400g_recipe;
dsp_options.append_to_db = false;
% dsp_options.append_mpi_reduction_db = false;
dsp_options.start_occurence = 1;
dsp_options.max_occurences = 1;
dsp_options.debug_plots = false;



dsp_options.database_type = "mysql";
dsp_options.dataBase = "labor_highspeed";

if ismac
    dsp_options.storage_path = "/Volumes/media/labdata/sioe_labor";
else
    dsp_options.storage_path = "W:\labdata\sioe_labor";
end
dsp_options.server = "192.168.178.192";
dsp_options.port = 3306;
dsp_options.user = "silas";
dsp_options.password = "silas";

db = DBHandler("dataBase", [dsp_options.dataBase], ...
    "type", dsp_options.database_type, ...
    "server", dsp_options.server, ...
    "user", dsp_options.user, ...
    "password", dsp_options.password);

%% Select runs

maxRunIds = 1;                   % keep small until the recipe settings are settled

fp = QueryFilter();
fp.where('Runs','fiber_length','EQUALS', 10);
fp.where('Runs','wavelength','EQUALS', 1310);
% fp.where('Runs','bitrate','LESS_THAN', 330e9);
fp.where('Runs','symbolrate','EQUALS', 174e9);
fp.where('Runs','pam_level','EQUALS', 6);
fp.where('Runs','rop_attenuation','EQUALS', 0);
fp.where('Runs','is_mpi','EQUALS', 0);
fp.where('Runs', 'db_mode','EQUALS', 1);

fields = db.getTableFieldNames('Runs');
[dataTable, query] = db.queryDB(fp, fields);
% disp(query);

dataTable = sortrows(dataTable, {'bitrate', 'run_id'});

% dataTable = dataTable(1:maxRunIds, :);

run_ids = dataTable.run_id(:).';

if isempty(run_ids)
    error("investigate_400g_algorithms:MissingRunIds", ...
        "No 400G runs match the current filters.");
end

fprintf("Selected %d run_id(s): %s\n", numel(run_ids), mat2str(run_ids));

%% Parameter sweep

dsp_options.userParameters = struct();

% Examples for parameter loops. DataStorage expands every vector-valued field.
% dsp_options.userParameters.len_tr = 4096*2;
% dsp_options.userParameters.pf_ncoeffs = [1,2,3];

% dsp_options.userParameters.decoding_mode = [db_decoder.sequencedetection];
% dsp_options.userParameters.pf_ncoeffs = [1, 2, 3];
% dsp_options.userParameters.mu_dc = [0, 1e-5, 1e-4];
% dsp_options.userParameters.run_ml_mlse_db = [false, true];
% dsp_options.userParameters.run_mlse_db = [false, true];

wh = DataStorage(dsp_options.userParameters);

n_realizations = (dsp_options.max_occurences - dsp_options.start_occurence + 1);
n_userparams = prod(wh.dim);
n_run_ids = numel(run_ids);
parallel_jobs = n_userparams * n_run_ids;
queried_jobs = n_realizations * n_userparams * n_run_ids;

fprintf("-> [ %d run_id(s) x %d userParam combination(s) = %d job(s) ] x %d realizations = %d total jobs \n", ...
    n_run_ids, n_userparams, parallel_jobs, n_realizations, queried_jobs);

%% Run

[results, wh] = submitJobs(run_ids, dsp_options, processingMode.serial, ...
    "wh", wh, ...
    "waitbar", true);


%% Quick result overview

printBerSummary(wh);
plotBerVsBitrateQuick(wh, dataTable);
plotMlseBerPrecodedVsBitrateByPf(wh, dataTable);

function printBerSummary(wh)
storageNames = fieldnames(wh.sto);
if isempty(storageNames)
    fprintf("No non-empty recipe outputs were stored.\n");
    return
end

fprintf("\nBER summary by stored package:\n");
for storageIdx = 1:numel(storageNames)
    storageName = storageNames{storageIdx};
    values = wh.sto.(storageName)(:).';
    berValues = extractBerValues(values);

    if isempty(berValues)
        fprintf("  %-18s no BER values\n", storageName);
    else
        fprintf("  %-18s min %.3e | median %.3e | n %d\n", ...
            storageName, min(berValues), median(berValues), numel(berValues));
    end
end
end

function berValues = extractBerValues(values)
berValues = [];
for valueIdx = 1:numel(values)
    packageCell = values{valueIdx};
    if isempty(packageCell)
        continue
    end
    if ~iscell(packageCell)
        packageCell = {packageCell};
    end

    for packageIdx = 1:numel(packageCell)
        package = packageCell{packageIdx};
        if isstruct(package) && isfield(package, "metrics")
            metrics = package.metrics;
            if isprop(metrics, "BER")
                berValues(end+1) = metrics.BER; %#ok<AGROW>
            elseif isstruct(metrics) && isfield(metrics, "BER")
                berValues(end+1) = metrics.BER; %#ok<AGROW>
            end
        end
    end
end
berValues = berValues(isfinite(berValues));
end

function plotBerVsBitrateQuick(wh, dataTable)
plotData = buildQuickBerTable(wh, dataTable);
if isempty(plotData)
    fprintf("No BER values available for quick BER-vs-bitrate plot.\n");
    return
end

storageNames = unique(plotData.storage_name, "stable");
decodingModes = [db_decoder.memoryless, db_decoder.sequencedetection];
berMetrics = ["BER", "BER_precoded"];
lineStyles = ["-", ":"];
rawMarkers = [".", "x"];
markers = ["o", "square", "diamond", "^", "v", ">"];

fig = figure(402); clf;
ax = axes(fig); hold(ax, "on");

for storageIdx = 1:numel(storageNames)
    storageName = storageNames(storageIdx);
    for modeIdx = 1:numel(decodingModes)
        decodingMode = decodingModes(modeIdx);
        for metricIdx = 1:numel(berMetrics)
            metricName = berMetrics(metricIdx);
            rowMask = plotData.storage_name == storageName & ...
                plotData.decoding_mode == decodingMode & ...
                plotData.metric_name == metricName;
            if ~any(rowMask)
                continue
            end

            modeData = sortrows(plotData(rowMask, :), "bitrate_Gbps");
            summaryData = groupsummary(modeData, "bitrate_Gbps", "min", "BER");
            summaryData = sortrows(summaryData, "bitrate_Gbps");
            color = quickPlotColor(modeIdx);
            marker = markers(1 + mod(storageIdx - 1, numel(markers)));
            label = sprintf("%s, %s, %s", storageName, ...
                decodingModeLabel(decodingMode), metricName);

            scatter(ax, modeData.bitrate_Gbps, modeData.BER, ...
                12, ...
                "Marker", rawMarkers(metricIdx), ...
                "MarkerEdgeColor", color, ...
                "MarkerEdgeAlpha", 0.25, ...
                "HandleVisibility", "off");

            plot(ax, summaryData.bitrate_Gbps, summaryData.min_BER, ...
                "LineStyle", lineStyles(metricIdx), ...
                "Marker", marker, ...
                "MarkerSize", 5, ...
                "LineWidth", 1.4, ...
                "Color", color, ...
                "DisplayName", label);
        end
    end
end

yline(ax, [2.2e-4, 4.85e-3, 2e-2], ...
    "LineWidth", 1, ...
    "LineStyle", "--", ...
    "Color", [0.25 0.25 0.25], ...
    "HandleVisibility", "off");

xlabel(ax, "Bitrate [Gb/s]");
ylabel(ax, "BER");
title(ax, "Quick BER vs bitrate");
set(ax, "YScale", "log");
grid(ax, "on");
box(ax, "on");
legend(ax, "Location", "best", "Interpreter", "none");

if exist("beautifyBERplot", "file")
    beautifyBERplot("logscale", true, "setcolors", false, ...
        "setmarkers", false, "changemarkers", false);
end
end

function plotData = buildQuickBerTable(wh, dataTable)
storageNames = fieldnames(wh.sto);
if isempty(storageNames)
    plotData = table();
    return
end

runIds = dataTable.run_id(:);
bitrates = dataTable.bitrate(:);

storageCol = strings(0, 1);
runIdCol = zeros(0, 1);
bitrateCol = zeros(0, 1);
decodingCol = db_decoder.empty(0, 1);
metricCol = strings(0, 1);
berCol = zeros(0, 1);

for storageIdx = 1:numel(storageNames)
    storageName = storageNames{storageIdx};
    storageValues = wh.sto.(storageName);
    for linIdx = 1:numel(storageValues)
        [phys, storedValue] = wh.getPhysAndValueByLinIndex(storageName, linIdx);
        metricRows = extractBerMetricRows({storedValue});
        if isempty(metricRows) || ~isfield(phys, "decoding_mode")
            continue
        end

        runId = resolveRunId(phys, runIds);
        bitrate = resolveBitrate(runId, runIds, bitrates);
        if ~isfinite(bitrate)
            continue
        end

        nRows = height(metricRows);
        storageCol(end+1:end+nRows, 1) = string(storageName);
        runIdCol(end+1:end+nRows, 1) = double(runId);
        bitrateCol(end+1:end+nRows, 1) = double(bitrate);
        decodingCol(end+1:end+nRows, 1) = phys.decoding_mode;
        metricCol(end+1:end+nRows, 1) = metricRows.metric_name;
        berCol(end+1:end+nRows, 1) = metricRows.BER;
    end
end

plotData = table(storageCol, runIdCol, bitrateCol, decodingCol, metricCol, berCol, ...
    'VariableNames', ["storage_name", "run_id", "bitrate", ...
    "decoding_mode", "metric_name", "BER"]);
if ~isempty(plotData)
    plotData = plotData(isfinite(plotData.BER) & plotData.BER > 0, :);
    plotData.bitrate_Gbps = plotData.bitrate .* 1e-9;
end
end

function plotMlseBerPrecodedVsBitrateByPf(wh, dataTable)
storageName = "mlse_package";
if ~isfield(wh.sto, storageName)
    fprintf("No %s storage available for BERp plot.\n", storageName);
    return
end

runIds = dataTable.run_id(:);
bitrates = dataTable.bitrate(:);
pfCol = zeros(0, 1);
bitrateCol = zeros(0, 1);
berpCol = zeros(0, 1);

storageValues = wh.sto.(storageName);
for linIdx = 1:numel(storageValues)
    [phys, storedValue] = wh.getPhysAndValueByLinIndex(storageName, linIdx);
    if ~isfield(phys, "pf_ncoeffs")
        continue
    end

    berp = extractMinMetricValue(storedValue, "BER_precoded");
    if ~isfinite(berp) || berp <= 0
        continue
    end

    runId = resolveRunId(phys, runIds);
    bitrate = resolveBitrate(runId, runIds, bitrates);
    if ~isfinite(bitrate)
        continue
    end

    pfCol(end+1, 1) = double(phys.pf_ncoeffs); %#ok<AGROW>
    bitrateCol(end+1, 1) = double(bitrate) * 1e-9; %#ok<AGROW>
    berpCol(end+1, 1) = double(berp); %#ok<AGROW>
end

if isempty(berpCol)
    fprintf("No BERp values available for %s.\n", storageName);
    return
end

plotData = table(pfCol, bitrateCol, berpCol, ...
    'VariableNames', ["pf_ncoeffs", "bitrate_Gbps", "BERp"]);
summaryData = groupsummary(plotData, ["pf_ncoeffs", "bitrate_Gbps"], ...
    "min", "BERp");

fig = figure(403); clf;
ax = axes(fig); hold(ax, "on");
pfValues = [1, 2, 3];
markers = ["o", "square", "diamond"];
colors = lines(numel(pfValues));

for pfIdx = 1:numel(pfValues)
    pfValue = pfValues(pfIdx);
    rowMask = summaryData.pf_ncoeffs == pfValue;
    if ~any(rowMask)
        continue
    end

    curveData = sortrows(summaryData(rowMask, :), "bitrate_Gbps");
    plot(ax, curveData.bitrate_Gbps, curveData.min_BERp, ...
        "LineWidth", 1.3, ...
        "Marker", markers(pfIdx), ...
        "MarkerSize", 5, ...
        "Color", colors(pfIdx, :), ...
        "DisplayName", sprintf("pf\\_ncoeffs = %d", pfValue));
end

xlabel(ax, "Bitrate [Gb/s]");
ylabel(ax, "min BERp");
title(ax, "MLSE BERp vs bitrate");
set(ax, "YScale", "log");
grid(ax, "on");
box(ax, "on");
legend(ax, "Location", "best", "Interpreter", "none");
set(fig, "Position", [1200, 450, 560, 360]);
end

function minValue = extractMinMetricValue(value, metricName)
minValue = NaN;
if isempty(value)
    return
end
if ~iscell(value)
    value = {value};
end

metricValues = NaN(1, numel(value));
for packageIdx = 1:numel(value)
    package = value{packageIdx};
    if ~isstruct(package) || ~isfield(package, "metrics")
        continue
    end
    metricValues(packageIdx) = readMetricValue(package.metrics, metricName);
end

metricValues = metricValues(isfinite(metricValues) & metricValues > 0);
if ~isempty(metricValues)
    minValue = min(metricValues);
end
end

function metricRows = extractBerMetricRows(values)
metricNames = strings(0, 1);
berValues = zeros(0, 1);
requestedMetrics = ["BER", "BER_precoded"];

for valueIdx = 1:numel(values)
    packageCell = values{valueIdx};
    if isempty(packageCell)
        continue
    end
    if ~iscell(packageCell)
        packageCell = {packageCell};
    end

    for packageIdx = 1:numel(packageCell)
        package = packageCell{packageIdx};
        if ~isstruct(package) || ~isfield(package, "metrics")
            continue
        end

        metrics = package.metrics;
        for metricIdx = 1:numel(requestedMetrics)
            metricName = requestedMetrics(metricIdx);
            value = readMetricValue(metrics, metricName);
            if isfinite(value)
                metricNames(end+1, 1) = metricName; %#ok<AGROW>
                berValues(end+1, 1) = value; %#ok<AGROW>
            end
        end
    end
end

metricRows = table(metricNames, berValues, ...
    'VariableNames', ["metric_name", "BER"]);
end

function value = readMetricValue(metrics, metricName)
value = NaN;
fieldName = char(metricName);
if isstruct(metrics) && isfield(metrics, fieldName)
    value = metrics.(fieldName);
elseif isobject(metrics) && isprop(metrics, fieldName)
    value = metrics.(fieldName);
end
end

function runId = resolveRunId(phys, runIds)
if isfield(phys, "run_id")
    runId = phys.run_id;
else
    runId = runIds(1);
end
end

function bitrate = resolveBitrate(runId, runIds, bitrates)
rowIdx = find(double(runIds) == double(runId), 1, "first");
if isempty(rowIdx)
    bitrate = NaN;
else
    bitrate = bitrates(rowIdx);
end
end

function color = quickPlotColor(modeIdx)
colors = [ ...
    0.1059 0.6196 0.4667; ...
    0.8510 0.3725 0.0078];
color = colors(1 + mod(modeIdx - 1, size(colors, 1)), :);
end

function label = decodingModeLabel(decodingMode)
switch decodingMode
    case db_decoder.memoryless
        label = "memoryless";
    case db_decoder.sequencedetection
        label = "sequence detection";
    otherwise
        label = string(decodingMode);
end
end