diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index 229bce9..f88eadb 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -526,13 +526,10 @@ classdef Signal end y_margin = 0.05 * y_range; - ylim([y_min - y_margin, y_max + y_margin]); - - % Set ticks automatically, avoid overpopulation - try - yticks(round(linspace(y_min, y_max, min(10, max(4, ceil(y_range/10)))))); - end - grid on; + ylim([y_min - y_margin, y_max + y_margin]); + + yticks(-200:10:200); + grid on; % Add legend if not already present if isempty(get(gca, 'Legend')) diff --git a/Functions/EQ_recipes/mpi_recipe_dev.m b/Functions/EQ_recipes/mpi_recipe_dev.m index 962ec62..1911f6e 100644 --- a/Functions/EQ_recipes/mpi_recipe_dev.m +++ b/Functions/EQ_recipes/mpi_recipe_dev.m @@ -15,12 +15,13 @@ end %% Execution toggles -run_conventional_ffe = 0; -run_hpf = 1; -run_a2_tracked_levels = 0; -run_a2_residual = 0; -run_a1 = 0; -run_tracking_adaptive = 0; +run_conventional_ffe = 1; +run_ma = 0; +run_hpf = 0; +run_a2_tracked_levels = 1; +run_a2_residual = 1; +run_a1 = 1; +run_tracking_adaptive = 1; plot_output_signals = 0; %% Shared fixed EQ settings @@ -75,12 +76,68 @@ if run_conventional_ffe "plain_ffe", "baseline", block_update); output.(char(storageName)) = ffe_results; + eq_noise_ffe = equalized_signal - Symbols; + + if plot_output_signals plotEqSignals(equalized_signal,Symbols,options,400,-1); end end -%% Conventional FFE + +%% Moving Average Filter -> Conventional FFE +if run_ma + eq_ffe = FFE_plain( ... + "sps", eq_sps, ... + "order", eq_order, ... + "decide", false, ... + "adaption_technique", eq_adaption, ... + "len_tr", eq_len_tr, ... + "epochs_tr", eq_epochs_tr, ... + "mu_tr", eq_mu_tr, ... + "dd_mode", true, ... + "epochs_dd", eq_epochs_dd, ... + "mu_dd", 0.0012, ... + "optmize_mus", false, ... + "plot_mu_optimization", options.debug_plots, ... + "save_debug", eq_save_debug); + + %%% APPLY MA filter + movavg_window_length = options.userParameters.bufferlen; + movavg_window_symmetry = "causal"; % "causal" or "noncausal" + + switch movavg_window_symmetry + case "causal" + dc_estimate = movmean(Scpe_sig.signal, ... + [movavg_window_length-1,0],"Endpoints","shrink"); + case "noncausal" + dc_estimate = movmean(Scpe_sig.signal, ... + movavg_window_length,"Endpoints","shrink"); + otherwise + error("mpi_recipe_dev:InvalidMovAvgSymmetry", ... + "movavg_window_symmetry must be 'causal' or 'noncausal'."); + end + + Scpe_sig_filt = Scpe_sig; + Scpe_sig_filt.signal = Scpe_sig.signal - dc_estimate; + clear dc_estimate; + + storageName = "MovAvg_plus_conventional_FFE"; + [ffe_results,equalized_signal] = runFfe(eq_ffe, "Moving-average + conventional FFE", ... + Scpe_sig_filt, Symbols, Tx_bits, options); + clear Scpe_sig_filt; + ffe_results = attachMpiReductionConfig(ffe_results, eq_ffe, storageName, ... + "MovAvg_plus_conventional_FFE", char(movavg_window_symmetry), block_update); + output.(char(storageName)) = ffe_results; + + if plot_output_signals + plotEqSignals(equalized_signal,Symbols,options,400,-1); + end +end + + + +%% High Pass Filter -> Conventional FFE if run_hpf eq_ffe = FFE_plain( ... "sps", eq_sps, ... @@ -99,11 +156,12 @@ if run_hpf %%% APPLY HPF in MHZ order f = Filter("f_cutoff",options.userParameters.hpf,"filterType","butterworth","lowpass",0,"filtdegree",4,"fs",options.fsym*2); - Scpe_sig = f.process(Scpe_sig); + Scpe_sig_filt = f.process(Scpe_sig); storageName = "HPF_plus_conventional_FFE"; [ffe_results,equalized_signal] = runFfe(eq_ffe, "Conventional FFE", ... - Scpe_sig, Symbols, Tx_bits, options); + Scpe_sig_filt, Symbols, Tx_bits, options); + clear Scpe_sig_filt; ffe_results = attachMpiReductionConfig(ffe_results, eq_ffe, storageName, ... "HPF_plus_conventional_FFE", "hpf", block_update); output.(char(storageName)) = ffe_results; @@ -191,7 +249,7 @@ if run_a1 "epochs_dd", eq_epochs_dd, ... "mu_dd", 0.012, ... "dc_smoothing_a1", 1, ... - "dc_avg_bufferlength_a1", 2048, ... + "dc_avg_bufferlength_a1", 2048, ... %was 2048 "dc_avg_update_blocklength_a1", block_update, ... "save_debug", eq_save_debug); @@ -238,7 +296,19 @@ if run_tracking_adaptive "dc_tracking", "adaptive", block_update); output.(char(storageName)) = ffe_results; + % fignum = 106; + % + % eq_noise = equalized_signal - Symbols; + % dn = sprintf("FFE DCT; SIR: %d dB",options.dataTable.sir); + % showEQNoisePSD(eq_noise, "fignum", fignum, "displayname", dn,"colormode","diverging"); + % ylim([-70 -30]); + % + % dn = sprintf("FFE only; SIR: %d dB",options.dataTable.sir); + % showEQNoisePSD(eq_noise_ffe, "fignum", fignum+1, "displayname", dn,"colormode","diverging"); + % ylim([-70 -30]); + if plot_output_signals + plotEqSignals(equalized_signal,Symbols,options,450,-1); end end diff --git a/Functions/EQ_visuals/showEQNoisePSD.m b/Functions/EQ_visuals/showEQNoisePSD.m index 1a16831..ba613da 100644 --- a/Functions/EQ_visuals/showEQNoisePSD.m +++ b/Functions/EQ_visuals/showEQNoisePSD.m @@ -5,8 +5,14 @@ arguments options.fignum (1,1) double = NaN % Default to NaN if not provided options.displayname (1,:) char = '' % Default to an empty string if not provided options.color = []; + options.colormode (1,1) string {mustBeMember(options.colormode, ["diverging", "qualitative", "paired"])} = "diverging"; end + persistent nextGroupIdx + if isempty(nextGroupIdx) + nextGroupIdx = 0; + end + % Determine the figure number to use or create a new figure if isnan(options.fignum) fig = figure; % Create a new figure and get its handle @@ -15,16 +21,16 @@ end end hold on ax = gca; - % N = numel(ax.Children); - if isempty(options.color) - N = sum(arrayfun(@(x) strcmp(x.LineStyle, '-'), ax.Children)); - cmap = linspecer(8); - options.color = cmap(mod(N, size(cmap, 1)) + 1, :); + linesBefore = findall(ax, 'Type', 'Line'); + + useAutomaticColor = isempty(options.color); + if useAutomaticColor + options.color = selectEQNoisePSDColor(ax, options.colormode); end % Ensure the figure is ready before calling spectrum eq_noise = eq_noise - mean(eq_noise.signal); - eq_noise.spectrum("displayname", options.displayname, "fignum", fig.Number, "normalizeTo0dB", 0,"color",options.color,"fft_length",4096*2); + eq_noise.spectrum("displayname", options.displayname, "fignum", fig.Number, "normalizeTo0dB", 0,"color",options.color,"fft_length",4096); if ~isnan(options.postfilter_taps) % Hold on to the figure for further plotting @@ -44,7 +50,88 @@ end legend('show'); end + if useAutomaticColor + nextGroupIdx = nextGroupIdx + 1; + newLines = getNewLines(ax, linesBefore); + tagEQNoisePSDLines(newLines, nextGroupIdx); + recolorEQNoisePSDLines(ax, options.colormode); + end + xlim([-eq_noise.fs/2* 1e-9 eq_noise.fs/2* 1e-9]); % ylim([-15, 0]); end + +function color = selectEQNoisePSDColor(ax, colormode) +if colormode == "qualitative" + N = sum(arrayfun(@(x) strcmp(x.LineStyle, '-'), ax.Children)); + cmap = linspecer(8); + color = cmap(mod(N, size(cmap, 1)) + 1, :); +else + color = [0 0 0]; +end +end + +function newLines = getNewLines(ax, linesBefore) +linesAfter = findall(ax, 'Type', 'Line'); +newLines = linesAfter(~ismember(linesAfter, linesBefore)); +end + +function tagEQNoisePSDLines(lines, groupIdx) +for lineIdx = 1:numel(lines) + setappdata(lines(lineIdx), 'showEQNoisePSDGroupIdx', groupIdx); +end +end + +function recolorEQNoisePSDLines(ax, colormode) +if colormode == "qualitative" + return +end + +lines = findall(ax, 'Type', 'Line'); +groupIdx = NaN(size(lines)); + +for lineIdx = 1:numel(lines) + if isappdata(lines(lineIdx), 'showEQNoisePSDGroupIdx') + groupIdx(lineIdx) = getappdata(lines(lineIdx), 'showEQNoisePSDGroupIdx'); + end +end + +lines = lines(~isnan(groupIdx)); +groupIdx = groupIdx(~isnan(groupIdx)); + +if isempty(lines) + return +end + +groups = unique(groupIdx, 'stable'); +groups = sort(groups); +cols = getEQNoisePSDColormap(colormode, numel(groups)); + +for groupColorIdx = 1:numel(groups) + sameGroup = groupIdx == groups(groupColorIdx); + set(lines(sameGroup), 'Color', cols(groupColorIdx, :)); +end +end + +function cols = getEQNoisePSDColormap(colormode, numColors) +switch colormode + case "paired" + cols = cbrewer2('paired', numColors); + otherwise + cols = getDivergingWithoutBrightCenter(numColors); +end +end + +function cols = getDivergingWithoutBrightCenter(numColors) +numBaseColors = max(6, numColors); +centerExcludeFraction = 0.3; + +baseCols = cbrewer2('RdYlBu', numBaseColors); +centerIdx = (numBaseColors + 1) / 2; +excludeHalfWidth = centerExcludeFraction * numBaseColors / 2; +keepIdx = find(abs((1:numBaseColors) - centerIdx) > excludeHalfWidth); +sampleIdx = round(linspace(1, numel(keepIdx), numColors)); + +cols = baseCols(keepIdx(sampleIdx), :); +end diff --git a/projects/Diss/MPI_revisit/algorithms/PLOT_mpi_reduction_db_delay_regimes_pam_comparison.m b/projects/Diss/MPI_revisit/algorithms/PLOT_mpi_reduction_db_delay_regimes_pam_comparison.m new file mode 100644 index 0000000..e85a22c --- /dev/null +++ b/projects/Diss/MPI_revisit/algorithms/PLOT_mpi_reduction_db_delay_regimes_pam_comparison.m @@ -0,0 +1,514 @@ +%% BER over SIR by MPI delay/coherence regime and PAM level +% 1) gather block_update = 1 data from the database +% 2) clean data and assign path-length regimes +% 3) plot PAM4/PAM6/PAM8 BER-over-SIR curves in the same figure + +clear; clc; + +%% 1) Gather data + +studyName = "block_update_sweep"; +selectedBlockUpdate = 1; +selectedPamLevels = [4 6 8]; + +algorithmSelection = table( ... + ["dc_tracking"; ... + "plain_ffe"; ... + "a2_tracked_levels"; ... + "a2_residual"; ... + "a1_moving_average"], ... + [true; true; false; false; false], ... + 'VariableNames', ["algorithm", "enabled"]); +selectedAlgorithms = algorithmSelection.algorithm(algorithmSelection.enabled); + +useBoundedLines = true; +usePolyfit = true; +polyfitOrderMax = 4; +boundaryPolyfitOrderMax = 4; +maxBerForPlot = 0.1; +scatterKeepFraction = 0.2; +boundMode = "fitStd"; % "fitStd", "directStd", "fitCI", or "directCI" +confidenceLevel = 0.95; +regimeNames = ["0-1 m", "10-100 m", "300 m", "1000 m"]; + +db = DBHandler( ... + "dataBase", "labor", ... + "type", "mysql", ... + "server", "192.168.178.192", ... + "user", "silas", ... + "password", "silas"); +db.refresh(); + +fp = QueryFilter(); +fp.where('Runs', 'fiber_length', 'EQUALS', 0); +fp.where('Runs', 'db_mode', 'EQUALS', '"no_db"'); +fp.where('Runs', 'v_bias', 'EQUALS', 2.65); +fp.where('MpiReductionResults', 'study_name', 'EQUALS', char(studyName)); +fp.where('MpiReductionResults', 'block_update', 'EQUALS', selectedBlockUpdate); + +selectedFields = { ... + 'Runs.run_id', ... + 'Runs.sir', ... + 'Runs.pam_level', ... + 'Runs.symbolrate', ... + 'Runs.interference_path_length', ... + 'Runs.power_mpi_interference', ... + 'Runs.power_pd_in', ... + 'Runs.v_bias', ... + 'Runs.loop_id', ... + 'MpiReductionResults.occurrence_idx', ... + 'MpiReductionResults.storage_name', ... + 'MpiReductionResults.algorithm', ... + 'MpiReductionResults.algorithm_variant', ... + 'MpiReductionResults.eq_class', ... + 'MpiReductionResults.block_update', ... + 'MpiReductionResults.BER', ... + 'MpiReductionResults.BER_precoded', ... + 'MpiReductionResults.SNR', ... + 'MpiReductionResults.GMI', ... + 'MpiReductionResults.AIR'}; +selectedFields = selectedFields(:); + +[rawData, query] = db.queryDB(fp, selectedFields); +disp(query); +fprintf("Fetched %d MPI reduction rows.\n", height(rawData)); + +%% 2) Clean data and assign regimes + +data = rawData; +numericFields = ["run_id", "sir", "pam_level", "symbolrate", ... + "interference_path_length", "occurrence_idx", "block_update", ... + "power_mpi_interference", "BER", "BER_precoded", "SNR", "GMI", "AIR"]; +for fieldIdx = 1:numel(numericFields) + fieldName = numericFields(fieldIdx); + if ismember(fieldName, string(data.Properties.VariableNames)) + data.(char(fieldName)) = numericColumn(data.(char(fieldName))); + end +end + +stringFields = ["storage_name", "algorithm", "algorithm_variant", "eq_class"]; +for fieldIdx = 1:numel(stringFields) + fieldName = stringFields(fieldIdx); + if ismember(fieldName, string(data.Properties.VariableNames)) + data.(char(fieldName)) = stringColumn(data.(char(fieldName))); + end +end + +data(data.run_id == 3866, :) = []; +data(data.run_id == 3865, :) = []; +data(data.run_id == 3796, :) = []; +data(data.run_id == 3797, :) = []; +data(data.run_id == 3798, :) = []; +data(data.run_id == 4002, :) = []; +data(data.run_id == 4200, :) = []; +data(data.run_id == 4199, :) = []; +data(data.loop_id == 91, :) = []; +data(data.loop_id == 59, :) = []; + +data = data(isfinite(data.BER) & data.BER > 0 & data.BER < maxBerForPlot, :); +data.sir_exact = -7 - data.power_mpi_interference; + +data.path_regime = strings(height(data), 1); +data.path_regime(data.interference_path_length == 0 | ... + data.interference_path_length == 1) = "0-1 m"; +data.path_regime(data.interference_path_length >= 10 & ... + data.interference_path_length <= 100) = "10-100 m"; +data.path_regime(data.interference_path_length == 300) = "300 m"; +data.path_regime(data.interference_path_length == 1000) = "1000 m"; +data = data(data.path_regime ~= "", :); + +data = data(ismember(data.pam_level, selectedPamLevels), :); +data = data(ismember(data.algorithm, selectedAlgorithms), :); +selectedAlgorithms = selectedAlgorithms(ismember(selectedAlgorithms, unique(data.algorithm, "stable"))); + +if isempty(data) + warning("PLOT_mpi_reduction_db_delay_regimes_pam_comparison:NoRows", ... + "No rows remain after BER/study/block/PAM/algorithm/regime filtering."); + return +end + +data.clean_keep = true(height(data), 1); +groupId = findgroups(data.path_regime, data.pam_level, data.algorithm, data.sir); +for curGroup = unique(groupId(isfinite(groupId))).' + rowMask = groupId == curGroup; + berValues = data.BER(rowMask); + if nnz(rowMask) > 3 + data.clean_keep(rowMask) = ~isoutlier(berValues); + end +end +cleanData = data(data.clean_keep, :); + +groupVars = ["path_regime", "pam_level", "algorithm", "sir"]; +summaryTable = groupsummary(cleanData, groupVars, ... + {"mean", "min", "max"}, "BER"); +sirExactTable = groupsummary(cleanData, groupVars, "median", "sir_exact"); +summaryTable = sortrows(summaryTable, groupVars); +sirExactTable = sortrows(sirExactTable, groupVars); +summaryTable.sir_exact = sirExactTable.median_sir_exact; +summaryTable = addBerIntervalBounds(summaryTable, cleanData, groupVars, confidenceLevel); + +fprintf("Cleaned to %d rows across %d regime/PAM/algorithm/SIR groups.\n", ... + height(cleanData), height(summaryTable)); +disp(groupcounts(cleanData, ["path_regime", "pam_level", "algorithm"])); + +%% 3) Plot PAM curves in one delay-regime comparison figure + + + +for regimeIdx = 1:numel(regimeNames) + regimeName = regimeNames(regimeIdx); + figure(); hold on; + + for algIdx = 1:numel(selectedAlgorithms) + algorithmName = selectedAlgorithms(algIdx); + algColor = algorithmColor(algorithmName); + displayName = algorithmDisplayName(algorithmName); + + for pamIdx = 1:numel(selectedPamLevels) + pamLevel = selectedPamLevels(pamIdx); + [marker, lineStyle] = pamStyle(pamLevel); + + rawMask = cleanData.path_regime == regimeName & ... + cleanData.pam_level == pamLevel & ... + cleanData.algorithm == algorithmName; + curveMask = summaryTable.path_regime == regimeName & ... + summaryTable.pam_level == pamLevel & ... + summaryTable.algorithm == algorithmName; + + if ~any(curveMask) + continue + end + + scatterRows = downsampleRows(cleanData(rawMask, :), scatterKeepFraction); + % scatter(scatterRows.sir_exact, scatterRows.BER, ... + % 14, ... + % "Marker", marker, ... + % "MarkerEdgeColor", algColor, ... + % "MarkerFaceColor", algColor, ... + % "MarkerEdgeAlpha", 0.35, ... + % "MarkerFaceAlpha", 0.35, ... + % "HandleVisibility", "off"); + + sirValues = summaryTable.sir_exact(curveMask).'; + meanBer = summaryTable.mean_BER(curveMask).'; + [boundCenterBer, boundLowerBer, boundUpperBer] = ... + selectBerBounds(summaryTable(curveMask, :), boundMode); + valid = isfinite(sirValues) & isfinite(meanBer) & meanBer > 0; + + if useBoundedLines && exist("boundedline", "file") && any(valid) + [xBand, centerBand, yBounds] = berIntervalBounds( ... + sirValues, boundCenterBer, boundLowerBer, boundUpperBer, ... + boundMode, boundaryPolyfitOrderMax); + + [hl, hp] = boundedline(xBand, centerBand, yBounds, ... + 'alpha', 'transparency', 0.08, ... + 'cmap', algColor, ... + 'nan', 'fill', ... + 'orientation', 'vert'); + set(hl, "LineStyle", "none", "Marker", "none", "HandleVisibility", "off"); + set(hp, "LineStyle", "-", "HandleVisibility", "off", "Marker", "none"); + end + + plot(sirValues(valid), meanBer(valid), ... + "LineStyle", lineStyle, ... + "Marker", marker, ... + "MarkerSize", 3.5, ... + "LineWidth", 1, ... + "Color", algColor, ... + "MarkerFaceColor", "w", ... + "MarkerEdgeColor", algColor, ... + "DisplayName", sprintf("%s, PAM %.0f", displayName, pamLevel)); + + if usePolyfit + fitMask = valid & meanBer > 0; + if nnz(fitMask) >= 2 + fitOrder = min(polyfitOrderMax, nnz(fitMask) - 1); + fitCoeff = polyfit(sirValues(fitMask), log10(meanBer(fitMask)), fitOrder); + xFit = linspace(min(sirValues(fitMask)), max(sirValues(fitMask)), 300); + yFit = 10 .^ polyval(fitCoeff, xFit); + + % plot(xFit, yFit, ... + % "LineStyle", "--", ... + % "LineWidth", 1.1, ... + % "Color", algColor, ... + % "HandleVisibility", "off"); + end + end + end + end + + yline(2.2e-4, "LineWidth", 1, "LineStyle", "--", "HandleVisibility", "off"); + yline(3.8e-3, "LineWidth", 1, "LineStyle", "--", "HandleVisibility", "off"); + yline(2e-2, "LineWidth", 1, "LineStyle", "--", "HandleVisibility", "off"); + + title(sprintf("Interference path regime %s", regimeName)); + xlabel("SIR (dB)"); + ylabel("BER"); + set(gca, "YScale", "log"); + ylim([9e-5, maxBerForPlot]); + xlim([15, 45]); + grid on; + box on; + legend("Location", "northeast", "Interpreter", "none"); + + if exist("beautifyBERplot", "file") + beautifyBERplot("logscale", true, "setcolors", false, ... + "setmarkers", false, "changemarkers", false); + end +end + +% mat2tikz_improved("C:/Users/Silas/Documents/6971e0b65b380ca6d71c837f/04_Experimental_Evaluation/tikz/mpi/ber_vs_sir_delay_regimes_pam_comparison.tikz","cleanfigure",1); + +%% Local helpers + +function values = numericColumn(values) + if iscell(values) + values = string(values); + end + if isstring(values) || ischar(values) + values = str2double(values); + end + values = double(values); +end + +function values = stringColumn(values) + if iscell(values) + values = string(values); + elseif ischar(values) + values = string(values); + end + values = strip(string(values)); +end + +function rows = downsampleRows(rows, keepFraction) + if keepFraction >= 1 || height(rows) <= 1 + return + end + + keepEvery = max(1, round(1 / keepFraction)); + keepIdx = 1:keepEvery:height(rows); + rows = rows(keepIdx, :); +end + +function summaryTable = addBerIntervalBounds(summaryTable, cleanData, groupVars, confidenceLevel) + nGroups = height(summaryTable); + ciCenter = NaN(nGroups, 1); + ciLower = NaN(nGroups, 1); + ciUpper = NaN(nGroups, 1); + stdCenter = NaN(nGroups, 1); + stdLower = NaN(nGroups, 1); + stdUpper = NaN(nGroups, 1); + + for groupIdx = 1:nGroups + rowMask = true(height(cleanData), 1); + for varIdx = 1:numel(groupVars) + varName = groupVars(varIdx); + rowMask = rowMask & cleanData.(char(varName)) == summaryTable.(char(varName))(groupIdx); + end + + [ciCenter(groupIdx), ciLower(groupIdx), ciUpper(groupIdx)] = ... + logBerMeanConfidenceInterval(cleanData.BER(rowMask), confidenceLevel); + [stdCenter(groupIdx), stdLower(groupIdx), stdUpper(groupIdx)] = ... + logBerMeanStdInterval(cleanData.BER(rowMask)); + end + + summaryTable.ci_center_BER = ciCenter; + summaryTable.ci_lower_BER = ciLower; + summaryTable.ci_upper_BER = ciUpper; + summaryTable.std_center_BER = stdCenter; + summaryTable.std_lower_BER = stdLower; + summaryTable.std_upper_BER = stdUpper; +end + +function [centerBer, lowerBer, upperBer] = logBerMeanConfidenceInterval(berValues, confidenceLevel) + berValues = berValues(isfinite(berValues) & berValues > 0); + if isempty(berValues) + centerBer = NaN; + lowerBer = NaN; + upperBer = NaN; + return + end + + logBer = log10(berValues(:)); + centerLog = mean(logBer, "omitnan"); + + if numel(logBer) < 2 + centerBer = 10 .^ centerLog; + lowerBer = centerBer; + upperBer = centerBer; + return + end + + alpha = 1 - confidenceLevel; + if exist("tinv", "file") + tCritical = tinv(1 - alpha/2, numel(logBer) - 1); + else + tCritical = 1.96; + end + + halfWidth = tCritical * std(logBer, 0, "omitnan") / sqrt(numel(logBer)); + centerBer = 10 .^ centerLog; + lowerBer = 10 .^ (centerLog - halfWidth); + upperBer = 10 .^ (centerLog + halfWidth); +end + +function [centerBer, lowerBer, upperBer] = logBerMeanStdInterval(berValues) + berValues = berValues(isfinite(berValues) & berValues > 0); + if isempty(berValues) + centerBer = NaN; + lowerBer = NaN; + upperBer = NaN; + return + end + + logBer = log10(berValues(:)); + centerLog = mean(logBer, "omitnan"); + stdLog = std(logBer, 0, "omitnan"); + + centerBer = 10 .^ centerLog; + lowerBer = 10 .^ (centerLog - stdLog); + upperBer = 10 .^ (centerLog + stdLog); +end + +function [centerBer, lowerBer, upperBer] = selectBerBounds(curveSummary, boundMode) + if boundMode == "fitCI" || boundMode == "directCI" + centerBer = curveSummary.ci_center_BER.'; + lowerBer = curveSummary.ci_lower_BER.'; + upperBer = curveSummary.ci_upper_BER.'; + else + centerBer = curveSummary.std_center_BER.'; + lowerBer = curveSummary.std_lower_BER.'; + upperBer = curveSummary.std_upper_BER.'; + end +end + +function [xBand, centerBand, yBounds] = berIntervalBounds(sirValues, centerBer, lowerBer, upperBer, boundMode, maxOrder) + if boundMode == "directCI" || boundMode == "directStd" + [xBand, centerBand, yBounds] = directBerBounds(sirValues, centerBer, lowerBer, upperBer); + return + end + + [xBand, centerBand, yBounds] = fittedBerBounds(sirValues, centerBer, lowerBer, upperBer, maxOrder); +end + +function [xBand, centerBand, yBounds] = directBerBounds(sirValues, centerBer, lowerBer, upperBer) + valid = isfinite(sirValues) & isfinite(centerBer) & isfinite(lowerBer) & ... + isfinite(upperBer) & centerBer > 0 & lowerBer > 0 & upperBer > 0; + + xBand = sirValues(valid).'; + centerBand = centerBer(valid).'; + lowerBand = lowerBer(valid).'; + upperBand = upperBer(valid).'; + + [xBand, orderIdx] = sort(xBand(:)); + centerBand = centerBand(orderIdx); + lowerBand = lowerBand(orderIdx); + upperBand = upperBand(orderIdx); + + lowerTmp = min(lowerBand, upperBand); + upperBand = max(lowerBand, upperBand); + lowerBand = lowerTmp; + + centerBand = min(max(centerBand, lowerBand), upperBand); + yBounds = [max(centerBand - lowerBand, 0), max(upperBand - centerBand, 0)]; +end + +function [xBand, centerBand, yBounds] = fittedBerBounds(sirValues, meanBer, minBer, maxBer, maxOrder) + valid = isfinite(sirValues) & isfinite(meanBer) & isfinite(minBer) & ... + isfinite(maxBer) & meanBer > 0 & minBer > 0 & maxBer > 0; + + x = sirValues(valid); + yMean = meanBer(valid); + yMin = minBer(valid); + yMax = maxBer(valid); + + if numel(x) < 2 + xBand = x(:); + centerBand = yMean(:); + yBounds = [max(yMean(:) - yMin(:), 0), max(yMax(:) - yMean(:), 0)]; + return + end + + [x, orderIdx] = sort(x(:)); + yMean = yMean(orderIdx); + yMin = yMin(orderIdx); + yMax = yMax(orderIdx); + + xBand = linspace(min(x), max(x), 300).'; + fitOrder = min(maxOrder, numel(unique(x)) - 1); + + if fitOrder < 1 + centerBand = interp1(x, yMean, xBand, "linear", "extrap"); + lowerBand = interp1(x, yMin, xBand, "linear", "extrap"); + upperBand = interp1(x, yMax, xBand, "linear", "extrap"); + else + centerBand = fitLogBer(x, yMean, xBand, fitOrder); + lowerBand = fitLogBer(x, yMin, xBand, fitOrder); + upperBand = fitLogBer(x, yMax, xBand, fitOrder); + end + + lowerTmp = min(lowerBand, upperBand); + upperBand = max(lowerBand, upperBand); + lowerBand = lowerTmp; + + centerBand = min(max(centerBand, lowerBand), upperBand); + yBounds = [max(centerBand - lowerBand, 0), max(upperBand - centerBand, 0)]; +end + +function yFit = fitLogBer(x, y, xFit, fitOrder) + coeff = polyfit(x, log10(y), fitOrder); + yFit = 10 .^ polyval(coeff, xFit); +end + +function label = algorithmDisplayName(algorithmName) + algorithmName = string(algorithmName); + switch algorithmName + case "plain_ffe" + label = "FFE only"; + case "a2_tracked_levels" + label = "ACT"; + case "a2_residual" + label = "L-DCA"; + case "a1_moving_average" + label = "DCA"; + case "dc_tracking" + label = "DCT"; + otherwise + label = algorithmName; + end +end + +function color = algorithmColor(algorithmName) + algorithmName = string(algorithmName); + switch algorithmName + case "plain_ffe" + color = [0.3467 0.5360 0.6907]; + case "a2_tracked_levels" + color = [0.9153 0.2816 0.2878]; + case "a2_residual" + color = [0.4416 0.7490 0.4322]; + case "a1_moving_average" + color = [1.0000 0.5984 0.2000]; + case "dc_tracking" + color = [0.6769 0.4447 0.7114]; + otherwise + color = [0 0 0]; + end +end + +function [marker, lineStyle] = pamStyle(pamLevel) + switch pamLevel + case 4 + lineStyle = "-"; + marker = "o"; + case 6 + lineStyle = "--"; + marker = "square"; + case 8 + lineStyle = ":"; + marker = "diamond"; + otherwise + lineStyle = "-."; + marker = "^"; + end +end diff --git a/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m b/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m index b7db98d..78ea07f 100644 --- a/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m +++ b/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m @@ -3,7 +3,7 @@ dsp_options = struct(); dsp_options.mode = "run_id"; dsp_options.recipe = @mpi_recipe_dev; dsp_options.append_to_db = false; -dsp_options.start_occurence = 5; +dsp_options.start_occurence = 1; dsp_options.max_occurences = 15; dsp_options.debug_plots = false; @@ -13,7 +13,7 @@ mpi_reduction_study_name = "pam6_8_sweep"; mpi_reduction_writer_path = fullfile(fileparts(mfilename('fullpath')),"db"); addpath(mpi_reduction_writer_path); -dsp_options.database_type = "mysql"; +dsp_options.database_type = "mysql"; dsp_options.dataBase = "labor"; @@ -51,7 +51,7 @@ for i = 1 fp.where('Runs', 'symbolrate', 'EQUALS', B); % 72 96 112 fp.where('Runs', 'fiber_length', 'EQUALS', 0); fp.where('Runs', 'interference_path_length', 'EQUALS', 1000); - fp.where('Runs', 'sir', 'LESS_EQUAL', 20); + fp.where('Runs', 'sir', 'LESS_EQUAL', 50); fp.where('Runs', 'db_mode', 'EQUALS', '"no_db"'); % fp.where('Runs', 'is_mpi', 'EQUALS', 1); fp.where('Runs', 'pam_level', 'EQUALS', M); @@ -62,14 +62,20 @@ for i = 1 [~, sortIdx] = sort(dataTable.sir, 'descend'); dataTable = dataTable(sortIdx, :); - dataTable = dataTable(1,:); + % desired_sir = [15,22,26,30,38]; + % desired_sir = [22,30,38]; + % % Filter dataTable to only keep rows with sir in desired_sir + % isDesired = ismember(dataTable.sir, desired_sir); + % dataTable = dataTable(isDesired, :); + % dataTable = dataTable(1,:); run_ids = dataTable.run_id; %% === Warehouse setup === dsp_options.userParameters = struct(); dsp_options.userParameters.block_update = 1;%[1,2,4,8,16,32,64,112,224,448,448*2,1024,2048,4096,8192,16384];%%logspace(-3.8,-1,22);%[linspace(2,4096,22)]; - dsp_options.userParameters.hpf = [2e6:2e6:10e6,20e6:10e6:100e6]; + % dsp_options.userParameters.hpf = [2e6:0.5e6:10e6]; + % dsp_options.userParameters.bufferlen = [112,224,448,448*2,1024,2048,4096,8192,16384]; wh = DataStorage(dsp_options.userParameters); %% @@ -85,7 +91,7 @@ for i = 1 %% === Run === % submitJobs returns the raw per-job results and the filled Warehouse. For a % single run_id and remove_dc = [0, 1], results is a 2-by-1 cell array. - [results, wh] = submitJobs(run_ids, dsp_options, processingMode.parallel, ... + [results, wh] = submitJobs(run_ids, dsp_options, processingMode.serial, ... "wh", wh, ... "waitbar", true); @@ -109,9 +115,9 @@ storageNames = fieldnames(wh.sto); x_base = dataTable.sir(:).'; x_base = dsp_options.userParameters.block_update; -x_base = dsp_options.userParameters.hpf; - -figure(2026); clf; hold on +x_base = dsp_options.userParameters.bufferlen; +% x_base = dsp_options.userParameters.hpf; +figure(2029); hold on for storage_idx = 1:numel(storageNames) storageName = storageNames{storage_idx}; result = wh.sto.(storageName); diff --git a/projects/Diss/MPI_revisit/parallelization_analysis/PLOT_mpi_reduction_db_parallelization_vs_sir.m b/projects/Diss/MPI_revisit/parallelization_analysis/PLOT_mpi_reduction_db_parallelization_vs_sir.m index ff04768..a7ca35f 100644 --- a/projects/Diss/MPI_revisit/parallelization_analysis/PLOT_mpi_reduction_db_parallelization_vs_sir.m +++ b/projects/Diss/MPI_revisit/parallelization_analysis/PLOT_mpi_reduction_db_parallelization_vs_sir.m @@ -8,7 +8,7 @@ clear; clc; studyName = "block_update_sweep"; pathLengthToPlot = 1000; -selectedPamLevels = 6; % set [] to use all PAM levels in the query result +selectedPamLevels = 4; % set [] to use all PAM levels in the query result selectedAlgorithms = []; % set [] to use all algorithms in the query result useBoundedLines = true; diff --git a/workerError.mat b/workerError.mat index a369927..34c4a9f 100644 Binary files a/workerError.mat and b/workerError.mat differ