start with 400G analysis work

Classes/04_DSP/Equalizer/ML_MLSE_DUOBINARY.m

@@ -0,0 +1,409 @@
+classdef ML_MLSE_DUOBINARY < handle
+    % ---------------------------------------------------------------------
+    % W. Lanneer and Y. Lefevre,
+    % “Machine Learning-Based Pre-Equalizers for Maximum Likelihood
+    %  Sequence Estimation in High-Speed PONs,” EUSIPCO 2023
+    % ---------------------------------------------------------------------
+    % This implementation reproduces the closed-loop ML-based
+    % pre-equalizer training for MLSE, supporting both training and
+    % detection (decision-directed) modes.
+    % ---------------------------------------------------------------------
+
+    properties
+        sps
+        order
+        e
+        e_tr
+        error
+
+        len_tr
+        mu_tr
+        epochs_tr
+
+        dd_mode
+        mu_dd
+        epochs_dd
+        adaptive_mu
+
+        constellation
+        L
+        alpha
+        DIR
+        DIR_flip
+        trellis_states
+        traceback_depth
+        delta
+
+        % Internal variables
+        S
+        Nf
+        nStates
+        nFeasible
+        combs
+        first_sym
+        last_sym
+        valid
+        valid_to_idx
+        valid_from_idx
+        w
+
+        % Fast lookup
+        nSym
+        key_table
+        trans_index
+        true_to_state_idx
+
+        % Debug metrics
+        ber = []
+        ber_dd = []
+        ce  = ones(1,1)
+    end
+
+    methods
+        function obj = ML_MLSE_DUOBINARY(options)
+            arguments(Input)
+                options.sps = 2;
+                options.order = 15;
+                options.len_tr = 4096;
+                options.mu_tr = 0.001;
+                options.epochs_tr = 5;
+                options.dd_mode = 1;
+                options.mu_dd = 1e-5;
+                options.epochs_dd = 5;
+                options.adaptive_mu = 1;
+                options.delta = 0;
+                options.traceback_depth = 1024;
+                options.L = 1;
+            end
+
+            fn = fieldnames(options);
+            for n = 1:numel(fn)
+                obj.(fn{n}) = options.(fn{n});
+            end
+
+            obj.e = zeros(obj.order,1);
+            obj.error = 0;
+        end
+
+        % ==============================================================
+        % PROCESS
+        % ==============================================================
+        function [X,X_viterbi] = process(obj, X, D)
+            % Normalize input RMS
+            X = X.normalize("mode","rms");
+            obj.constellation = sort(unique(D.signal),'ascend');
+            obj.nSym = numel(obj.constellation);
+
+            if length(X)/length(D) ~= obj.sps
+                warning('Signal length does not fit to reference!');
+            end
+
+            % --- Parameters
+            obj.S = obj.nSym;
+            obj.Nf = obj.order * obj.sps;
+            obj.nStates = obj.S^obj.L;
+            obj.nFeasible = obj.nStates * obj.S;
+
+            % --- Trellis mapping
+            obj.trellis_states = reshape(obj.constellation,1,[]);
+            pre_comb_mat  = repmat(obj.trellis_states, obj.L, 1);
+            pre_comb_cell = mat2cell(pre_comb_mat, ones(1,obj.L), size(pre_comb_mat,2));
+            obj.combs     = fliplr(combvec(pre_comb_cell{:}).');
+            obj.first_sym = obj.combs(:,1);
+            obj.last_sym  = obj.combs(:,end);
+            obj.nStates   = size(obj.combs,1);
+
+            % --- Valid transitions
+            obj.valid = false(obj.nStates);
+            for from = 1:obj.nStates
+                for to = 1:obj.nStates
+                    if all(obj.combs(to,2:end) == obj.combs(from,1:end-1))
+                        obj.valid(to,from) = true;
+                    end
+                end
+            end
+            [obj.valid_to_idx,obj.valid_from_idx] = find(obj.valid);
+
+            % --- Initialize weights
+            if isempty(obj.w) || any(size(obj.w) ~= [obj.Nf+1,obj.nFeasible])
+                % obj.w = randn(obj.Nf+1,obj.nFeasible);
+                obj.w = zeros(obj.Nf+1,obj.nFeasible);
+            end
+
+            % --- Fast lookup tables
+            [~, sym_idx_mat] = ismember(obj.combs, obj.constellation);
+            key_vals = 1 + sum((sym_idx_mat - 1) .* (obj.nSym .^ (0:obj.L-1)), 2);
+            max_key = obj.nSym^obj.L;
+            obj.key_table = zeros(max_key,1,'uint32');
+            obj.key_table(key_vals) = 1:obj.nStates;
+
+            obj.trans_index = sparse(obj.nStates,obj.nStates);
+            for i = 1:length(obj.valid_from_idx)
+                f = obj.valid_from_idx(i);
+                t = obj.valid_to_idx(i);
+                obj.trans_index(t,f) = i;
+            end
+
+            % ==============================================================
+            % TRAINING
+            % ==============================================================
+            fprintf('\n--- Training mode ---\n');
+            obj.equalize(X.signal, D.signal, obj.mu_tr, obj.epochs_tr, obj.len_tr, true);
+            obj.e_tr = obj.e;
+
+            % ==============================================================
+            % DECISION-DIRECTED / TESTING
+            % ==============================================================
+            fprintf('--- Decision-directed / detection mode ---\n');
+            [y, y_vit] = obj.equalize(X.signal, D.signal, obj.mu_dd, obj.epochs_dd, X.length, false);
+
+            X_viterbi = X;
+            X.signal = y;
+            X_viterbi.signal = y_vit;
+        end
+
+        % ==============================================================
+        % EQUALIZE
+        % ==============================================================
+        function [y,y_ref] = equalize(obj,x,d,mu,epochs,N,training)
+            debug = 0;
+            showPlots = 0;
+            y = zeros(N,1);
+            nSymbols = ceil(N/obj.sps);
+
+            for epoch = 1:epochs
+                pm = zeros(obj.nStates,1);
+                pred = zeros(nSymbols,obj.nStates,'uint32');
+                pm_sto = nan(obj.nStates,nSymbols,'like',pm);
+                CE_accum = 0;
+
+                start_sample = 1;
+                end_sample   = N;
+                start_symbol = 1 + floor((start_sample - 1)/obj.sps);
+
+                % --- initialize true state
+                if numel(d) >= obj.L && start_symbol >= obj.L
+                    init_seq = d(start_symbol-obj.L+1:start_symbol);
+                    key_init = obj.seq2key(init_seq);
+                    true_to_state_idx = obj.key_table(key_init);
+                    if true_to_state_idx==0, true_to_state_idx=1; end
+                else
+                    true_to_state_idx = uint32(1);
+                end
+
+                for sample = start_sample:obj.sps:end_sample
+                    symbol = (sample - start_sample)/obj.sps + 1;
+                    sym_idx = start_symbol + (symbol - 1);
+
+                    % --- Observation window (with delta)
+                    i1 = sample - obj.Nf + 1 + obj.delta;
+                    i2 = sample + obj.delta;
+                    buf = x(max(1,i1):min(length(x),i2));
+                    padL = max(0,1 - i1);
+                    padR = max(0,i2 - length(x));
+                    yk = [zeros(padL,1); buf(:); zeros(padR,1)];
+                    yk = [yk;1];
+
+                    % --- Branch metrics
+                    c_hat = (yk.' * obj.w).';
+                    pm = pm - min(pm);
+                    v_tilde = pm(obj.valid_from_idx) + c_hat;
+
+                    % --- allocate once
+                    if epoch==1 && symbol==1
+                        obj.true_to_state_idx = ones(ceil(N/obj.sps),1,'uint32');
+                    end
+
+                    % --- previous "to" becomes "from"
+                    if symbol>1
+                        true_from_state_idx = obj.true_to_state_idx(symbol-1);
+                    else
+                        true_from_state_idx = 1;
+                    end
+
+                    % --- compute or reuse "to" state
+                    if epoch==1
+                        if sym_idx>=obj.L
+                            key_to = obj.seq2key(d(sym_idx-obj.L+1:sym_idx));
+                            state_idx = obj.key_table(key_to);
+                            if state_idx==0
+                                state_idx = true_from_state_idx;
+                            end
+                            obj.true_to_state_idx(symbol) = state_idx;
+                        else
+                            obj.true_to_state_idx(symbol) = true_from_state_idx;
+                        end
+                    end
+                    true_to_state_idx = obj.true_to_state_idx(symbol);
+
+                    % --- fast Dirac creation
+                    dirac = zeros(obj.nFeasible,1);
+                    trans_idx = obj.trans_index(true_to_state_idx,true_from_state_idx);
+                    if trans_idx~=0
+                        dirac(trans_idx)=1;
+                    end
+
+                    % --- ensure valid (from,to)
+                    if ~any(dirac)
+                        mask  = obj.valid_from_idx==true_from_state_idx & ...
+                            obj.valid_to_idx  ==true_to_state_idx;
+                        if any(mask)
+                            dirac(mask) = 1;
+                        else
+                            idx = find(obj.valid_from_idx==true_from_state_idx,1,'first');
+                            dirac(idx) = 1;
+                            obj.true_to_state_idx(symbol) = obj.valid_to_idx(idx);
+                        end
+                    end
+
+                    % ===================================================================
+                    % TRAINING MODE (weight update)
+                    % ===================================================================
+                    if training
+                        % --- Softmax and CE
+                        v_shift = -(v_tilde - min(v_tilde));
+                        v_shift = min(v_shift,100);
+                        expv = exp(v_shift);
+                        p = expv./(sum(expv)+eps);
+                        CE_symbol(symbol) = -log(p(dirac==1)+eps);
+
+                        % --- CE smoothing and adaptive μ
+                        if sym_idx>obj.L
+                            CE_smooth(symbol)=0.01*CE_symbol(symbol)+0.99*CE_symbol(symbol-1);
+                        else
+                            CE_smooth(symbol)=CE_symbol(symbol);
+                        end
+                        CE_accum=CE_accum+CE_symbol(symbol);
+
+                        % --- Gradient update
+                        dmp=(dirac-p)';
+                        dL_Dw=(yk).*dmp;
+                        if sym_idx>=obj.L
+                            if obj.adaptive_mu
+                                mu_eff=CE_smooth(symbol);
+                                mu_eff=max(min(mu_eff,0.2),1e-4);
+                            else
+                                mu_eff=mu;
+                            end
+                            obj.w=obj.w - mu_eff.*dL_Dw;
+                        end
+                    end
+
+                    % ===================================================================
+                    % DECODING MODE (Viterbi only)
+                    % ===================================================================
+                    % Compare-Select (always executed)
+                    vmat=inf(obj.nStates,obj.nStates);
+                    vmat(obj.valid)=v_tilde;
+                    [pm_next,pred(symbol,:)]=min(vmat,[],2);
+                    pm_next=pm_next-min(pm_next);
+                    pm=pm_next;
+                    pm_sto(:,symbol)=pm;
+                end
+
+                % --- Traceback
+                [~,s_end]=min(pm);
+                vpath=zeros(symbol,1,'uint32');
+                vpath(symbol)=s_end;
+                for n=symbol:-1:2
+                    vpath(n-1)=pred(n,vpath(n));
+                end
+
+                y_ref=d(start_symbol:end);
+                y=obj.first_sym(vpath);
+
+                % --- BER/CE reporting and plots
+                if training
+                    err=sum(y~=y_ref(1:length(y)));
+                    ser=err/length(y);
+                    [ber, ~] = obj.calculateDuobinaryBer(y, y_ref);
+                    if isfinite(ber)
+                        fprintf('Epoch %d - BER: %.2e\n',epoch,ber);
+                        obj.ber(epoch)=ber;
+                    else
+                        fprintf('Epoch %d - SER: %.2e\n',epoch,ser);
+                        obj.ber(epoch)=ser;
+                    end
+                    obj.ce(epoch)=CE_accum/symbol;
+
+                    if debug && mod(epoch,10)==1 && showPlots
+                        figure(10);clf
+                        subplot(3,2,1:2);
+                        imagesc(obj.w);axis xy;colorbar;title('Filter W');
+                        subplot(3,2,3);
+                        vtilde_mat=NaN(obj.nStates,obj.nStates);
+                        vtilde_mat(obj.valid)=v_tilde;
+                        imagesc(vtilde_mat);axis xy;colorbar;title('Path Metrics (v\_tilde)');
+                        subplot(3,2,4);
+                        plot(1:symbol,pm_sto);title('Path Metric Evolution');
+                        subplot(3,2,5);hold on;
+                        scatter(1:symbol,CE_symbol,1,'.');
+                        scatter(1:symbol,CE_smooth,1,'.');
+                        title('Cross Entropy');
+                        subplot(3,2,6);hold on;
+                        yyaxis left
+                        scatter(1:length(obj.ce),obj.ce,10,'s','filled');
+                        ylabel('Cross Entropy');
+                        yyaxis right
+                        scatter(1:length(obj.ber),obj.ber,10,'d','filled');
+                        set(gca,'YScale','log');
+                        ylabel('BER (log)');
+                        xlabel('Epoch');grid on;
+                        title('Convergence');
+                        drawnow;
+                    end
+                else
+                    [ber, ser] = obj.calculateDuobinaryBer(y, y_ref);
+                    if isfinite(ber)
+                        fprintf('DD epoch %d - BER: %.2e\n',epoch,ber);
+                        obj.ber_dd(epoch)=ber;
+                    else
+                        fprintf('DD epoch %d - SER: %.2e\n',epoch,ser);
+                        obj.ber_dd(epoch)=ser;
+                    end
+                end
+            end
+        end
+
+    end
+
+    methods (Access=private)
+        % ==============================================================
+        % Helper: convert a detected precoded sequence back to PAM data
+        % ==============================================================
+        function data = invertDuobinaryPrecoder(obj, precoded)
+            encoded = Duobinary().encode(precoded,"M",obj.S);
+            encoded_signal = Signal(encoded);
+            decoded_signal = Duobinary().decode(encoded_signal,"M",obj.S);
+            data = decoded_signal.signal;
+        end
+
+        % ==============================================================
+        % Helper: calculate BER after inverting the duobinary precoder
+        % ==============================================================
+        function [ber, ser] = calculateDuobinaryBer(obj, detected, reference)
+            n = min(numel(detected), numel(reference));
+            detected = detected(1:n);
+            reference = reference(1:n);
+
+            ser = sum(detected ~= reference) / n;
+
+            detected_data = obj.invertDuobinaryPrecoder(detected);
+            reference_data = obj.invertDuobinaryPrecoder(reference);
+            mapper = PAMmapper(obj.S, 0);
+            detected_bits = mapper.demap(detected_data);
+            reference_bits = mapper.demap(reference_data);
+            [~,~,ber,~] = calc_ber(reference_bits, detected_bits, ...
+                "skip_front",10,"skip_end",10,"returnErrorLocation",1);
+        end
+
+        % ==============================================================
+        % Helper: Sequence → key (always scalar)
+        % ==============================================================
+        function key = seq2key(obj, seq)
+            [~, idx] = ismember(flip(seq), obj.constellation);
+            pow = (obj.nSym .^ (0:obj.L-1)).';
+            key = 1 + sum((idx(:) - 1) .* pow);
+        end
+    end
+end

Classes/04_DSP/Equalizer/VNLE.m

@@ -30,6 +30,9 @@ classdef VNLE < handle
         optmize_mus = 0;
         mu_optimization
         mu_optimization_iter = 0;
+        mu_optimization_len
+        plot_mu_optimization = 0
+        mu_optimization_fignum = 3020;
 
         x_norm
         ce
@@ -55,6 +58,9 @@ classdef VNLE < handle
                 options.decide = false;
                 options.save_debug = 0;
                 options.optmize_mus = 0;
+                options.mu_optimization_len = 2^15;
+                options.plot_mu_optimization = 0;
+                options.mu_optimization_fignum = 3020;
 
             end
 
@@ -110,7 +116,6 @@ classdef VNLE < handle
             lbdesc = [num2str(obj.order),' tap FFE'];
             X  = X.logbookentry(lbdesc); % append to logbook
 
-            N = X;
             N = X - D;
 
 
@@ -138,13 +143,18 @@ classdef VNLE < handle
                            ones(1,obj.ce(3))*mu(3) ]);
             end
 
+            x = x(:);
+            d = d(:);
             x = [zeros(floor(obj.order(1)/2),1); x; zeros(obj.order(1),1)];
+            n_symbols = floor(N / obj.sps);
+            y = zeros(n_symbols,1);
+            d_hat = zeros(n_symbols,1);
 
             if showviz
-                f = figure(111);
+                figure(111);
                 subplot(2,2,1:2);
                 hold on
-                a = scatter(1:numel(x),x,1,'.');
+                scatter(1:numel(x),x,1,'.');
                 a2 = scatter(1,1,1,'.');
                 a3 = scatter(1,1,2,'.');
                 a4 = xline(1);
@@ -215,40 +225,50 @@ classdef VNLE < handle
             mu_range = [1e-5, 1e-2];
             mu_dc_range = [1e-5, 1e-1];
 
-            vars = [optimizableVariable("mu_tr",mu_range,"Transform","log"), ...
+            [x_opt,d_opt,N_opt] = obj.optimizationSignals(x,d);
-                optimizableVariable("mu_dd",mu_range,"Transform","log")];
+
+            vars = obj.muOptimizableVariables("mu_tr",mu_range);
+            vars = [vars, obj.muOptimizableVariables("mu_dd",mu_range)];
             optimize_mu_dc = obj.mu_dc ~= 0;
             if optimize_mu_dc
                 vars = [vars, optimizableVariable("mu_dc",mu_dc_range,"Transform","log")];
             end
 
             obj.mu_optimization_iter = 0;
-            obj.mu_optimization = bayesopt(@(p)obj.muObjective(p,x,d),vars, ...
+            fprintf("VNLE mu opt uses %d samples / %d symbols\n",N_opt,numel(d_opt));
-                "MaxObjectiveEvaluations",10, ...
+            obj.mu_optimization = bayesopt(@(p)obj.muObjective(p,x_opt,d_opt),vars, ...
+                "MaxObjectiveEvaluations",20, ...
                 "AcquisitionFunctionName","expected-improvement-plus", ...
                 "IsObjectiveDeterministic",false, ...
                 "Verbose",0, ...
                 "PlotFcn",[]);
-            obj.mu_tr = obj.mu_optimization.XAtMinObjective.mu_tr;
+
-            obj.mu_dd = obj.mu_optimization.XAtMinObjective.mu_dd;
+            obj.mu_tr = obj.muVectorFromParams(obj.mu_optimization.XAtMinObjective,"mu_tr");
+            obj.mu_dd = obj.muVectorFromParams(obj.mu_optimization.XAtMinObjective,"mu_dd");
             if optimize_mu_dc
                 obj.mu_dc = obj.mu_optimization.XAtMinObjective.mu_dc;
             end
 
-            objective_db = 10*log10(obj.mu_optimization.MinObjective);
             if optimize_mu_dc
-                fprintf("\nVNLE mu opt done: mu_tr=%9.3e, mu_dd=%9.3e, mu_dc=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB\n", ...
+                fprintf("\nVNLE mu opt done: mu_tr=[%s], mu_dd=[%s], mu_dc=%9.3e, BER=%9.3e\n", ...
-                    obj.mu_tr,obj.mu_dd,obj.mu_dc,obj.mu_optimization.MinObjective,objective_db);
+                    obj.formatMuVector(obj.mu_tr),obj.formatMuVector(obj.mu_dd), ...
+                    obj.mu_dc,obj.mu_optimization.MinObjective);
             else
-                fprintf("\nVNLE mu opt done: mu_tr=%9.3e, mu_dd=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB\n", ...
+                fprintf("\nVNLE mu opt done: mu_tr=[%s], mu_dd=[%s], BER=%9.3e\n", ...
-                    obj.mu_tr,obj.mu_dd,obj.mu_optimization.MinObjective,objective_db);
+                    obj.formatMuVector(obj.mu_tr),obj.formatMuVector(obj.mu_dd), ...
+                    obj.mu_optimization.MinObjective);
+            end
+
+            if obj.plot_mu_optimization
+                obj.plotMuOptimization();
             end
         end
 
         function objective = muObjective(obj,params,x,d)
-            old_debug = obj.save_debug;
+            old_state = obj.captureObjectiveState();
-            old_mu_dc = obj.mu_dc;
+            cleanup = onCleanup(@()obj.restoreObjectiveState(old_state));
-            obj.save_debug = 1;
+
+            obj.save_debug = 0;
             optimize_mu_dc = ismember("mu_dc",string(params.Properties.VariableNames));
             if optimize_mu_dc
                 obj.mu_dc = params.mu_dc;
@@ -257,34 +277,217 @@ classdef VNLE < handle
             obj.e = zeros(sum(obj.ce),1);
             obj.e_dc = 0;
             obj.debug_struct = struct();
-            obj.equalize(x,d,params.mu_tr,obj.epochs_tr,obj.len_tr,1,0);
+            muTrCandidate = obj.muVectorFromParams(params,"mu_tr");
-            obj.equalize(x,d,params.mu_dd,obj.epochs_dd,numel(x),0,0);
+            muDdCandidate = obj.muVectorFromParams(params,"mu_dd");
+            N_tr = min(obj.len_tr,numel(x));
+            obj.equalize(x,d,muTrCandidate,obj.epochs_tr,N_tr,1,0);
+            [signal,~] = obj.equalize(x,d,muDdCandidate,obj.epochs_dd,numel(x),0,0);
 
-            objective = mean(obj.debug_struct.error(end,:),"omitnan");
+            [ber,errors] = obj.berObjective(signal,d);
+            objective = ber;
             if ~isfinite(objective)
                 objective = inf;
             end
 
-            objective_db = 10*log10(objective);
             obj.mu_optimization_iter = obj.mu_optimization_iter + 1;
             if optimize_mu_dc
-                fprintf("\rVNLE mu opt %02d: mu_tr=%9.3e, mu_dd=%9.3e, mu_dc=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB", ...
+                fprintf("\rVNLE mu opt %02d: mu_tr=[%s], mu_dd=[%s], mu_dc=%9.3e, BER=%9.3e, errors=%d", ...
-                    obj.mu_optimization_iter,params.mu_tr,params.mu_dd,params.mu_dc,objective,objective_db);
+                    obj.mu_optimization_iter,obj.formatMuVector(muTrCandidate), ...
+                    obj.formatMuVector(muDdCandidate),params.mu_dc,ber,errors);
             else
-                fprintf("\rVNLE mu opt %02d: mu_tr=%9.3e, mu_dd=%9.3e, MSE=%9.3e, MSE_dB=%7.2f dB", ...
+                fprintf("\rVNLE mu opt %02d: mu_tr=[%s], mu_dd=[%s], BER=%9.3e, errors=%d", ...
-                    obj.mu_optimization_iter,params.mu_tr,params.mu_dd,objective,objective_db);
+                    obj.mu_optimization_iter,obj.formatMuVector(muTrCandidate), ...
+                    obj.formatMuVector(muDdCandidate),ber,errors);
             end
-            obj.save_debug = old_debug;
+
-            obj.mu_dc = old_mu_dc;
+            clear cleanup
+        end
+
+        function [x_opt,d_opt,N_opt] = optimizationSignals(obj,x,d,opt_len)
+            if nargin < 4
+                opt_len = obj.mu_optimization_len;
+            end
+
+            N_available = min(numel(x),numel(d) * obj.sps);
+
+            if isempty(opt_len) || opt_len <= 0 || isinf(opt_len)
+                N_opt = N_available;
+            else
+                N_opt = min(N_available,max(obj.len_tr,opt_len));
+            end
+
+            N_opt = obj.sps * floor(N_opt / obj.sps);
+            N_opt = max(obj.sps,N_opt);
+
+            n_symbols = N_opt / obj.sps;
+            x_opt = x(1:N_opt);
+            d_opt = d(1:n_symbols);
+        end
+
+        function [ber,errors] = berObjective(~,signal,d)
+            M = numel(unique(d));
+            mapper = PAMmapper(M,0);
+            eq_signal_sd = Signal(signal);
+            eq_signal_hd = mapper.quantize(eq_signal_sd);
+            tx_symbols = Signal(d);
+            rx_bits = mapper.demap(eq_signal_hd);
+            tx_bits = mapper.demap(tx_symbols);
+            skip_front = min(1000,max(0,floor(numel(rx_bits.signal) / 4)));
+            [~,errors,ber,~] = calc_ber(rx_bits.signal,tx_bits.signal, ...
+                "skip_front",skip_front, ...
+                "skip_end",0, ...
+                "returnErrorLocation",1);
+        end
+
+        function vars = muOptimizableVariables(obj,prefix,mu_range)
+            vars = optimizableVariable.empty;
+            activeOrders = find(obj.order > 0);
+
+            for idx = 1:numel(activeOrders)
+                orderIdx = activeOrders(idx);
+                varName = sprintf("%s_%d",prefix,orderIdx);
+                vars = [vars, optimizableVariable(varName,mu_range,"Transform","log")]; %#ok<AGROW>
+            end
+        end
+
+        function mu = muVectorFromParams(obj,params,prefix)
+            mu = zeros(1,3);
+            for orderIdx = 1:numel(mu)
+                varName = sprintf("%s_%d",prefix,orderIdx);
+                if ismember(varName,string(params.Properties.VariableNames))
+                    mu(orderIdx) = params.(varName);
+                elseif numel(obj.(char(prefix))) >= orderIdx
+                    mu(orderIdx) = obj.(char(prefix))(orderIdx);
+                else
+                    mu(orderIdx) = obj.(char(prefix))(1);
+                end
+            end
+        end
+
+        function plotMuOptimization(obj)
+            if isempty(obj.mu_optimization)
+                return
+            end
+
+            X = obj.mu_optimization.XTrace;
+            objective = obj.mu_optimization.ObjectiveTrace;
+            objective = objective(:);
+            valid = isfinite(objective);
+
+            if isempty(X) || ~any(valid)
+                return
+            end
+
+            var_names = X.Properties.VariableNames;
+            n_vars = numel(var_names);
+            eval_idx = (1:numel(objective)).';
+            objective_plot = obj.positiveObjectiveForLogPlot(objective);
+            best_plot = obj.positiveObjectiveForLogPlot(cummin(objective));
+
+            figure(obj.mu_optimization_fignum);
+            clf;
+            t = tiledlayout(2,2,"TileSpacing","compact","Padding","compact");
+            title(t,"VNLE Bayesian mu optimization");
+
+            nexttile;
+            h_candidate = semilogy(eval_idx,objective_plot,"o-","DisplayName","candidate");
+            obj.addOptimizationDataTips(h_candidate,X,objective,objective_plot,eval_idx,var_names);
+            hold on;
+            h_best_trace = semilogy(eval_idx,best_plot,"k-","LineWidth",1.2,"DisplayName","best so far");
+            h_best_trace.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow("Evaluation",eval_idx);
+            h_best_trace.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow("Best BER",best_plot);
+            grid on;
+            xlabel("Evaluation");
+            ylabel("BER");
+            legend("Location","best");
+
+            if n_vars < 2
+                return
+            end
+
+            pairs = nchoosek(1:n_vars,2);
+            n_pair_plots = min(size(pairs,1),3);
+            [~,best_idx] = min(objective);
+
+            for pair_idx = 1:n_pair_plots
+                nexttile;
+                x_name = var_names{pairs(pair_idx,1)};
+                y_name = var_names{pairs(pair_idx,2)};
+                x_data = X.(x_name);
+                y_data = X.(y_name);
+                c_data = log10(objective_plot);
+
+                h_scatter = scatter(log10(x_data),log10(y_data),35,c_data,"filled");
+                obj.addOptimizationDataTips(h_scatter,X,objective,objective_plot,eval_idx,var_names);
+                hold on;
+                h_best = plot(log10(x_data(best_idx)),log10(y_data(best_idx)),"kp", ...
+                    "MarkerSize",12, ...
+                    "MarkerFaceColor","y", ...
+                    "DisplayName","best");
+                obj.addOptimizationDataTips(h_best,X(best_idx,:),objective(best_idx),objective_plot(best_idx),eval_idx(best_idx),var_names);
+                grid on;
+                xlabel("log10(" + string(x_name) + ")");
+                ylabel("log10(" + string(y_name) + ")");
+                cb = colorbar;
+                cb.Label.String = "log10(BER)";
+                title(string(x_name) + " vs " + string(y_name));
+            end
+        end
+
+        function addOptimizationDataTips(~,plot_handle,X,objective,objective_plot,eval_idx,var_names)
+            plot_handle.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow("Evaluation",eval_idx);
+            plot_handle.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow("BER",objective);
+            plot_handle.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow("BER shown",objective_plot);
+
+            for var_idx = 1:numel(var_names)
+                var_name = var_names{var_idx};
+                plot_handle.DataTipTemplate.DataTipRows(end+1) = dataTipTextRow(var_name,X.(var_name));
+            end
+        end
+
+        function objective_plot = positiveObjectiveForLogPlot(~,objective)
+            objective_plot = objective;
+            positive_values = objective(isfinite(objective) & objective > 0);
+
+            if isempty(positive_values)
+                floor_value = 1e-12;
+            else
+                floor_value = min(positive_values) / 10;
+            end
+
+            objective_plot(~isfinite(objective_plot) | objective_plot <= 0) = floor_value;
+        end
+
+        function state = captureObjectiveState(obj)
+            state.e = obj.e;
+            state.e_dc = obj.e_dc;
+            state.error = obj.error;
+            state.mu_dc = obj.mu_dc;
+            state.save_debug = obj.save_debug;
+            state.debug_struct = obj.debug_struct;
+        end
+
+        function restoreObjectiveState(obj,state)
+            obj.e = state.e;
+            obj.e_dc = state.e_dc;
+            obj.error = state.error;
+            obj.mu_dc = state.mu_dc;
+            obj.save_debug = state.save_debug;
+            obj.debug_struct = state.debug_struct;
+        end
+
+        function s = formatMuVector(~,mu)
+            s = strtrim(sprintf("%9.3e ",mu));
         end
 
         %% Functions needed During Adaption
          function x_in_vnle_format = calcVNLENonlinVecs(~,x_in_block,I_2,I_3,N_,norm_)
             % These are the second and third order input signal products of the VNLE EQ
             % ∑ h1 x_in(k-n1) + ∑∑ h2 x_in(k-n1)*x_in(k-n2) + ∑∑∑ h3 x_in(k-n1)*x_in(k-n2)*x_in(k-n3)
+            x_in_block = x_in_block(:);
             l1=length(x_in_block);
-            l2=length(I_2);
+            l2=size(I_2,1);
-            l3=length(I_3);
+            l3=size(I_3,1);
             final_length = l1+l2+l3;
 
             x_in_vnle_format = zeros(final_length,1);
@@ -382,6 +585,7 @@ classdef VNLE < handle
         end 
 
         function powerNorm  = calcPowerNormalization(~,v)
+            v = v(:);
 
             powerNorm(1) = sqrt(mean(abs(v   ).^2));
             powerNorm(2) = sqrt(mean(abs(v.^2).^2));

Classes/DataBaseHandler/Metricstruct.m

@@ -13,6 +13,9 @@ classdef Metricstruct
         BER_precoded (1,1) double {mustBeNumeric, mustBeNonnegative} = 0
         numBitErr_precoded (1,1) double {mustBeInteger, mustBeNonnegative} = 0
 
+        % BER_DB_memoryless (1,1) double {mustBeNumeric} = NaN
+        % BER_DB_sequencedetection (1,1) double {mustBeNumeric} = NaN
+
         SNR (1,1) double {mustBeNumeric} = NaN
         SNR_level (:,1) double {mustBeNumeric} = []
         STD (1,1) double {mustBeNumeric} = NaN

Functions/EQ_blocks/duobinary_signaling.m

@@ -21,10 +21,13 @@ arguments
     tx_symbols
     tx_bits
     options.precode_mode db_mode
+    options.decoding_mode db_decoder = db_decoder.sequencedetection;
+
     options.showAnalysis = 0
     options.eth_style_symbol_mapping = 0
     options.postFFE = []
     options.database = []
+    
 end
 
 %% Process signals through equalizer
@@ -35,6 +38,12 @@ if ~isempty(options.postFFE)
     [eq_signal, eq_noise] = options.postFFE.process(eq_signal, tx_symbols);
 end
 
+
+
+run_both_detection_schemes = 1;
+if options.decoding_mode == db_decoder.sequencedetection || run_both_detection_schemes %MLSE
+
+
         if isa(mlse_,'MLSE_viterbi')
             [mlse_signal] = mlse_.process(eq_signal);
         else
@@ -49,29 +58,52 @@ else
             % ref_sym_dec = Duobinary().decode(ref_sym_dbenc); %ref_sym wieder zurück!
 
             mlse_.trellis_states = PAMmapper(M,0).levels;
-    mlse_.trellis_state_mode = 1;
+            % mlse_.scale_mode = 2;
-    [mlse_signal,LLR,GMI_MLSE] = mlse_.process(eq_signal,ref_sym_dpc);
+            if M == 6
+                mlse_.trellis_exclusion = 0;
+            end
+            mlse_.debug = 0;
+            [mlse_signal,LLR,GMI_MLSE] = mlse_.process(eq_signal.normalize("mode","rms"),ref_sym_dpc);
         end
 
-
-
-% tx_symbols_ = Duobinary().decode(tx_symbols);
-% [mlse_signal,~,GMI_MLSE] = mlse_.process(eq_signal,tx_symbols);
-
         % Apply duobinary encoding and decoding
-mlse_signal = Duobinary().encode(mlse_signal);
+        mlse_signal = Duobinary().encode(mlse_signal,"M",M);
-mlse_signal = Duobinary().decode(mlse_signal);
+
+        mlse_signal = Duobinary().decode(mlse_signal,"M",M);
 
         % Demap symbols to bits
         rx_bits = PAMmapper(M, 0, "eth_style", options.eth_style_symbol_mapping).demap(mlse_signal);
 
-%% Calculate BER and metrics
+        [bits_db, ~, ber_sequencedetection, ~] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", 100, "skip_end", 150, "returnErrorLocation", 1);
-[bits_db, errors_db, ber_db, error_pos] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", 100, "skip_end", 150, "returnErrorLocation", 1);
+
+end
+
+if options.decoding_mode == db_decoder.memoryless || run_both_detection_schemes
+
+        
+        db_ref_constellation = unique(tx_symbols.signal);
+        eq_signal_hd = PAMmapper(M, 0).quantize(eq_signal,'custom_const',db_ref_constellation.');
+        eq_signal_hd = Duobinary().decode(eq_signal_hd,"M",M);
+
+        % Demap
+        rx_bits = PAMmapper(M, 0, "eth_style", options.eth_style_symbol_mapping).demap(eq_signal_hd);
+
+        [bits_db, ~, ber_memoryless, ~] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", 100, "skip_end", 150, "returnErrorLocation", 1);
+
+end
+
+
+
+% Calculate BER and metrics
 
 % Calculate performance metrics after duobinary FFE!
 [snr, snr_lvl] = calc_snr(tx_symbols.signal, eq_noise.signal); %SNR of duobinary sequence - not directly comparable to 
-[gmi] = calc_air(eq_signal, tx_symbols, "skip_front", 10000, "skip_end", 10000);
+% [gmi] = calc_air(eq_signal, tx_symbols, "skip_front", 10000, "skip_end", 10000); % Not working for Duobinary == Channel with memory
+% air = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi ./ log2(double(M));
+if options.decoding_mode == db_decoder.sequencedetection || run_both_detection_schemes
+    [gmi] = GMI_MLSE;
     air = tx_symbols.fs .* floor(log2(double(M))*10)/10 .* gmi ./ log2(double(M));
+end
 [evm_total, evm_lvl] = calc_evm(eq_signal, tx_symbols);
 [std_total, std_lvl] = calc_std(eq_signal, tx_symbols);
 [std_rxraw_total, std_rxraw_lvl] = calc_std(rx_signal.resample("fs_out", tx_symbols.fs), tx_symbols);
@@ -91,9 +123,12 @@ db_results.metrics.result_id = NaN;
 db_results.metrics.run_id = NaN;
 db_results.metrics.eqParam_id = NaN;
 db_results.metrics.date_of_processing = datetime('now');
-db_results.metrics.BER = ber_db;
+db_results.metrics.BER = ber_sequencedetection; % THIS IS A CONVENTION
+db_results.metrics.BER_precoded = ber_memoryless; % THIS IS A CONVENTION
+% db_results.metrics.BER_DB_memoryless = ber_memoryless;
+% db_results.metrics.BER_DB_sequencedetection = ber_sequencedetection;
 db_results.metrics.numBits = bits_db;
-db_results.metrics.numBitErr = errors_db;
+% db_results.metrics.numBitErr = NaN;
 db_results.metrics.SNR = snr;
 db_results.metrics.SNR_level = snr_lvl;
 db_results.metrics.STD = std_total;

Functions/EQ_recipes/dsp_400g_recipe.m

@@ -0,0 +1,449 @@
+function output = dsp_400g_recipe(Scpe_sig_raw, Symbols, Tx_bits, options)
+%dsp_400g_recipe Run 400G equalizer schemes from a synchronized scope signal.
+
+arguments
+    Scpe_sig_raw
+    Symbols
+    Tx_bits
+    options.fsym
+    options.M
+    options.duob_mode
+    options.dataTable table
+    options.userParameters struct = struct()
+    options.preprocess_mode string = "auto"
+    options.tx_pulseformer = []
+    options.debug_plots (1,1) logical = false
+end
+
+p = defaultRecipeParameters();
+p.preprocess_mode = options.preprocess_mode;
+p = applyUserParameters(p, options.userParameters);
+
+output = struct();
+
+Scpe_sig = preprocessSignal(Scpe_sig_raw, Symbols, options.fsym, ...
+    "mode", p.preprocess_mode, ...
+    "tx_pulseformer", options.tx_pulseformer, ...
+    "debug_plots", options.debug_plots);
+[Scpe_sig, Symbols, Tx_bits] = alignDspInputs(Scpe_sig, Symbols, Tx_bits, p.eq_sps);
+if isempty(p.ml_mlse_len_tr)
+    p.ml_mlse_len_tr = floor(length(Scpe_sig) / 4);
+end
+
+if p.plot_input_signal
+    showLevelScatter(Scpe_sig.normalize("mode","rms"), Symbols, ...
+        "fsym", options.fsym, ...
+        "fignum", p.input_plot_fignum, ...
+        "normalize", true);
+end
+
+if options.duob_mode ~= db_mode.db_encoded
+    if p.run_ffe
+        eq_ffe = FFE("epochs_tr", p.epochs_tr, ...
+            "epochs_dd", p.epochs_dd, ...
+            "len_tr", p.len_tr, ...
+            "mu_dd", p.ffe_mu_dd, ...
+            "mu_tr", p.ffe_mu_tr, ...
+            "order", p.ffe_order(1), ...
+            "sps", p.eq_sps, ...
+            "decide", false, ...
+            "optmize_mus", p.optimize_mus, ...
+            "dd_mode", p.dd_mode, ...
+            "adaption_technique", p.ffe_adaption, ...
+            "dc_tracking_mu", p.mu_dc);
+
+        [ffe_results, equalized_signal] = runFfe(eq_ffe, "FFE", ...
+            Scpe_sig, Symbols, Tx_bits, options);
+        ffe_results.config.equalizer_structure = equalizer_structure.ffe;
+        ffe_results.recipe_config = collectRecipeConfig("ffe", eq_ffe, p, options);
+        output.ffe_package = ffe_results;
+
+        if p.plot_output_signals
+            plotEqSignals(equalized_signal, Symbols, options, p.output_plot_fignum, -1);
+        end
+    end
+
+
+    if p.run_vnle
+        %     eq_vnle = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2, ...
+        % "mu_dd",[0.0004 0.0005 0.0006],"mu_tr",[0.0001 0.0008 0.001], ...
+        % "order",[150,5,5],"sps",2,"decide",1, ...
+        % "optmize_mus",1,"mu_optimization_len",2^15);
+
+        eq_vnle = EQ("Ne", p.vnle_ffe_order, ...
+            "Nb", p.vnle_dfe_order, ...
+            "training_length", p.len_tr, ...
+            "training_loops", p.epochs_tr, ...
+            "dd_loops", p.epochs_dd, ...
+            "K", p.eq_K, ...
+            "DCmu", p.mu_dc, ...
+            "DDmu", [p.eq_mu_ffe p.mu_dfe], ...
+            "DFEmu", p.dfe_mu_feedback, ...
+            "FFEmu", 0, ...
+            "plotfinal", 0, ...
+            "ideal_dfe", false);
+
+        [vnle_results, equalized_signal] = runFfe(eq_vnle, "VNLE", ...
+            Scpe_sig, Symbols, Tx_bits, options);
+        vnle_results.config.equalizer_structure = equalizer_structure.ffe;
+        ffe_results.recipe_config = collectRecipeConfig("vnle", eq_ffe, p, options);
+        output.vnle_package = vnle_results;
+
+        if p.plot_output_signals
+            plotEqSignals(equalized_signal, Symbols, options, p.output_plot_fignum, -1);
+        end
+    end
+
+
+    if p.run_dfe
+        eq_dfe = EQ("Ne", p.dfe_ffe_order, ...
+            "Nb", p.dfe_feedback_order, ...
+            "training_length", p.len_tr, ...
+            "training_loops", p.epochs_tr, ...
+            "dd_loops", p.epochs_dd, ...
+            "K", p.eq_K, ...
+            "DCmu", p.mu_dc, ...
+            "DDmu", [p.eq_mu_ffe p.mu_dfe], ...
+            "DFEmu", p.dfe_mu_feedback, ...
+            "FFEmu", 0, ...
+            "plotfinal", 0, ...
+            "ideal_dfe", false);
+
+        dfe_results = ffe(eq_dfe, options.M, Scpe_sig, Symbols, Tx_bits, ...
+            "precode_mode", options.duob_mode, ...
+            "showAnalysis", options.debug_plots, ...
+            "postFFE", [], ...
+            "eth_style_symbol_mapping", 0);
+        dfe_results.config.equalizer_structure = equalizer_structure.dfe;
+        dfe_results.metrics.print("description", resultDescription("DFE", options));
+        dfe_results.recipe_config = collectRecipeConfig("dfe", eq_dfe, p, options);
+        output.dfe_package = dfe_results;
+    end
+
+    if p.run_vnle_mlse
+        eq_vnle = EQ("Ne", p.vnle_ffe_order, ...
+            "Nb", p.vnle_dfe_order, ...
+            "training_length", p.len_tr, ...
+            "training_loops", p.epochs_tr, ...
+            "dd_loops", p.epochs_dd, ...
+            "K", p.eq_K, ...
+            "DCmu", p.mu_dc, ...
+            "DDmu", [p.eq_mu_ffe p.mu_dfe], ...
+            "DFEmu", p.dfe_mu_feedback, ...
+            "FFEmu", 0, ...
+            "plotfinal", 0, ...
+            "ideal_dfe", false);
+        pf = Postfilter("ncoeff", p.pf_ncoeffs, "useBurg", true);
+        mlse = buildMlse(options.M, options.duob_mode, p, "pf_mlse");
+
+        [vnle_results, mlse_results] = vnle_postfilter_mlse(eq_vnle, pf, mlse, ...
+            options.M, Scpe_sig, Symbols, Tx_bits, ...
+            "precode_mode", options.duob_mode, ...
+            "showAnalysis", options.debug_plots, ...
+            "postFFE", [], ...
+            "eth_style_symbol_mapping", 0);
+
+        vnle_results.config.equalizer_structure = equalizer_structure.vnle;
+        vnle_results.recipe_config = collectRecipeConfig("vnle", eq_vnle, p, options);
+        mlse_results.config.equalizer_structure = equalizer_structure.vnle_pf_mlse;
+        mlse_results.recipe_config = collectRecipeConfig("vnle_pf_mlse", eq_vnle, p, options);
+        vnle_results.metrics.print("description", resultDescription("VNLE", options));
+        mlse_results.metrics.print("description", resultDescription("VNLE + PF + MLSE", options));
+
+        output.vnle_package = vnle_results;
+        output.mlse_package = mlse_results;
+    end
+
+    if p.run_dbtgt
+        eq_dbtgt = EQ("Ne", p.dbtgt_ffe_order, ...
+            "Nb", p.dbtgt_dfe_order, ...
+            "training_length", p.len_tr, ...
+            "training_loops", p.epochs_tr, ...
+            "dd_loops", p.epochs_dd, ...
+            "K", p.eq_K, ...
+            "DCmu", p.mu_dc, ...
+            "DDmu", [p.eq_mu_ffe p.mu_dfe], ...
+            "DFEmu", p.dfe_mu_feedback, ...
+            "FFEmu", 0, ...
+            "plotfinal", 0, ...
+            "ideal_dfe", true);
+        mlse_db = buildMlse(options.M, options.duob_mode, p, "db_target");
+
+        dbtgt_results = runDuobinaryTarget(eq_dbtgt, mlse_db, ...
+            Scpe_sig, Symbols, Tx_bits, options, p);
+        dbtgt_results.config.equalizer_structure = equalizer_structure.vnle_db_mlse;
+        dbtgt_results.recipe_config = collectRecipeConfig("vnle_db_mlse", eq_dbtgt, p, options);
+        dbtgt_results.metrics.print("description", resultDescription("VNLE DB target + MLSE", options));
+        output.dbtgt_package = dbtgt_results;
+    end
+
+    if p.run_ml_mlse
+        ml_mlse_equalizer = ML_MLSE("epochs_tr", p.ml_mlse_epochs_tr, ...
+            "epochs_dd", p.ml_mlse_epochs_dd, ...
+            "len_tr", p.ml_mlse_len_tr, ...
+            "mu_dd", p.ml_mlse_mu_dd, ...
+            "mu_tr", p.ml_mlse_mu_tr, ...
+            "order", p.ml_mlse_order, ...
+            "sps", p.eq_sps, ...
+            "traceback_depth", p.ml_mlse_traceback_depth, ...
+            "L", p.ml_mlse_L, ...
+            "delta", p.ml_mlse_delta, ...
+            "adaptive_mu", p.ml_mlse_adaptive_mu);
+
+        ml_mlse_results = ml_mlse(ml_mlse_equalizer, options.M, ...
+            Scpe_sig, Symbols, Tx_bits, ...
+            "precode_mode", options.duob_mode);
+        ml_mlse_results.config.equalizer_structure = equalizer_structure.ml_mlse;
+        ml_mlse_results.recipe_config = collectRecipeConfig("ml_mlse", ml_mlse_equalizer, p, options);
+        output.mlmlse_package = ml_mlse_results;
+    end
+else
+    if p.run_ml_mlse_db
+
+        ml_mlse_db_equalizer = ML_MLSE_DUOBINARY("epochs_tr", p.ml_mlse_epochs_tr, ...
+            "epochs_dd", p.ml_mlse_epochs_dd, ...
+            "len_tr", p.ml_mlse_len_tr, ...
+            "mu_dd", p.ml_mlse_mu_dd, ...
+            "mu_tr", p.ml_mlse_mu_tr, ...
+            "order", p.ml_mlse_order, ...
+            "sps", p.eq_sps, ...
+            "traceback_depth", p.ml_mlse_traceback_depth, ...
+            "L", p.ml_mlse_L, ...
+            "delta", p.ml_mlse_delta, ...
+            "adaptive_mu", p.ml_mlse_adaptive_mu);
+
+        % Use a precoded reference for the detector; the received waveform remains encoded.
+        ref_sym = PAMmapper(options.M,0).map(Tx_bits);
+        Symbols_precoded = Duobinary().precode(ref_sym); % precoded
+
+        
+        ml_mlse_db_results = ml_mlse(ml_mlse_db_equalizer, options.M, ...
+            Scpe_sig, Symbols_precoded, Tx_bits, ...
+            "precode_mode", db_mode.db_precoded);
+
+        ml_mlse_db_results.config.equalizer_structure = equalizer_structure.ml_mlse;
+        ml_mlse_db_results.config.comment = 'function: ML-based MLSE; duobinary encoded';
+        ml_mlse_db_results.recipe_config = collectRecipeConfig("ml_mlse_db", ml_mlse_db_equalizer, p, options);
+        output.mlmlse_db_package = ml_mlse_db_results;
+    end
+
+    if p.run_mlse_db
+        eq_db_enc = EQ("Ne", p.dbtgt_ffe_order, ...
+            "Nb", p.dbtgt_dfe_order, ...
+            "training_length", p.len_tr, ...
+            "training_loops", p.epochs_tr, ...
+            "dd_loops", p.epochs_dd, ...
+            "K", p.eq_K, ...
+            "DCmu", p.mu_dc, ...
+            "DDmu", [p.eq_mu_ffe p.mu_dfe], ...
+            "DFEmu", p.dfe_mu_feedback, ...
+            "FFEmu", 0, ...
+            "plotfinal", 0, ...
+            "ideal_dfe", true);
+        mlse_db_enc = MLSE("DIR", [1,1], ...
+            "duobinary_output", 0, ...
+            "M", options.M, ...
+            "trellis_states", PAMmapper(options.M,0).levels);
+
+        if isempty(p.decoding_mode)
+            mlse_db_results = duobinary_signaling(eq_db_enc, mlse_db_enc, ...
+                options.M, Scpe_sig, Symbols, Tx_bits, ...
+                "precode_mode", options.duob_mode, ...
+                "showAnalysis", options.debug_plots, ...
+                "postFFE", []);
+            mlse_db_results.metrics.print("description", resultDescription(["DB Encoded; MLSE"], options));
+        else
+            mlse_db_results = duobinary_signaling(eq_db_enc, mlse_db_enc, ...
+                options.M, Scpe_sig, Symbols, Tx_bits, ...
+                "precode_mode", options.duob_mode, ...
+                "showAnalysis", options.debug_plots, ...
+                "postFFE", [],"decoding_mode",p.decoding_mode);
+            mlse_db_results.metrics.print("description", resultDescription(["DB Encoded "+string(p.decoding_mode)], options));
+        end
+        mlse_db_results.config.equalizer_structure = equalizer_structure.db_encoded;
+        mlse_db_results.recipe_config = collectRecipeConfig("mlse_db", eq_db_enc, p, options);
+        output.mlse_db_package = mlse_db_results;
+        
+    end
+end
+
+end
+
+function p = defaultRecipeParameters()
+p = struct();
+
+p.run_ffe = false;
+p.run_vnle = false;
+p.run_dfe = false;
+p.run_vnle_mlse = false;
+p.run_dbtgt = false;
+p.run_ml_mlse = true;     % non-encoded and precoded branches
+p.run_ml_mlse_db = true;  % db_encoded: ML-based MLSE
+p.run_mlse_db = true;     % db_encoded: conventional MLSE
+
+
+p.preprocess_mode = "auto";
+p.plot_input_signal = false;
+p.plot_output_signals = false;
+p.input_plot_fignum = 400;
+p.output_plot_fignum = 410;
+
+p.eq_sps = 2;
+p.eq_K = 2;
+p.len_tr = 4096*2;
+p.epochs_tr = 5;
+p.epochs_dd = 5;
+p.dd_mode = true;
+p.optimize_mus = true;
+
+p.ffe_order = [50, 0, 0];
+p.dfe_ffe_order = [50, 5, 5];
+p.vnle_ffe_order = [50, 5, 5];
+p.dbtgt_ffe_order = [50, 5, 5];
+p.vnle_dfe_order = [0, 0, 0];
+p.dbtgt_dfe_order = [0, 0, 0];
+p.dfe_feedback_order = [2, 0, 0];
+
+p.eq_mu_ffe = [0.0001, 0.0008, 0.001];
+p.ffe_mu_tr = 0.4;
+p.ffe_mu_dd = 0.1;
+p.ffe_adaption = "nlms";
+p.mu_dfe = 0.0004;
+p.mu_dc = 1.021e-05;
+p.dfe_mu_feedback = 0.005;
+
+p.pf_ncoeffs = 1;
+p.use_viterbi = false;
+p.mlse_scale_mode = 2;
+p.mlse_trellis_state_mode = 2;
+p.dbtgt_trellis_state_mode = 3;
+p.decoding_mode = [];
+
+p.ml_mlse_mu_tr = 0.03;
+p.ml_mlse_mu_dd = 0.03;
+p.ml_mlse_epochs_tr = 100;
+p.ml_mlse_epochs_dd = 1;
+p.ml_mlse_len_tr = [];
+p.ml_mlse_order = 11;
+p.ml_mlse_traceback_depth = 256;
+p.ml_mlse_L = 1;
+p.ml_mlse_delta = 4;
+p.ml_mlse_adaptive_mu = false;
+end
+
+function p = applyUserParameters(p, userParameters)
+if isempty(userParameters)
+    return
+end
+
+paramNames = fieldnames(userParameters);
+for paramIdx = 1:numel(paramNames)
+    paramName = paramNames{paramIdx};
+    if ~isfield(p, paramName)
+        warning("dsp_400g_recipe:UnknownUserParameter", ...
+            "Ignoring unknown user parameter '%s'.", paramName);
+        continue
+    end
+    p.(paramName) = userParameters.(paramName);
+end
+end
+
+function [Scpe_sig, Symbols, Tx_bits] = alignDspInputs(Scpe_sig, Symbols, Tx_bits, sps)
+nSymbols = min(length(Symbols), floor(length(Scpe_sig) / sps));
+if nSymbols <= 0
+    error("dsp_400g_recipe:EmptyAlignedSignal", ...
+        "No overlapping samples remain after preprocessing and synchronization.");
+end
+
+Scpe_sig.signal = real(Scpe_sig.signal(1:sps*nSymbols));
+Symbols.signal = Symbols.signal(1:nSymbols,:);
+
+% if isprop(Tx_bits, "signal")
+%     Tx_bits.signal = Tx_bits.signal(1:nSymbols,:); <- THIS IS WRONG!! 
+% end
+end
+
+function [ffe_results, equalized_signal] = runFfe(eq_ffe, description, ...
+        Scpe_sig, Symbols, Tx_bits, options)
+[ffe_results, equalized_signal] = ffe(eq_ffe, options.M, Scpe_sig, Symbols, Tx_bits, ...
+    "precode_mode", options.duob_mode, ...
+    "showAnalysis", options.debug_plots, ...
+    "postFFE", [], ...
+    "eth_style_symbol_mapping", 0);
+
+ffe_results.metrics.print("description", resultDescription(description, options));
+end
+
+function dbtgt_results = runDuobinaryTarget(eq_dbtgt, mlse_db, ...
+        Scpe_sig, Symbols, Tx_bits, options, p)
+if isempty(p.decoding_mode)
+    dbtgt_results = duobinary_target(eq_dbtgt, mlse_db, options.M, ...
+        Scpe_sig, Symbols, Tx_bits, ...
+        "precode_mode", options.duob_mode, ...
+        "showAnalysis", options.debug_plots, ...
+        "postFFE", []);
+else
+    dbtgt_results = duobinary_target(eq_dbtgt, mlse_db, options.M, ...
+        Scpe_sig, Symbols, Tx_bits, ...
+        "precode_mode", options.duob_mode, ...
+        "showAnalysis", options.debug_plots, ...
+        "postFFE", [], ...
+        "decoding_mode", p.decoding_mode);
+end
+end
+
+function mlse = buildMlse(M, duobMode, p, mode)
+if p.use_viterbi
+    mlse = MLSE_viterbi("duobinary_output", 0, ...
+        "M", M, ...
+        "trellis_states", PAMmapper(M,0).levels);
+    return
+end
+
+if duobMode == db_mode.no_db && M == 6
+    trellisExclusion = true;
+else
+    trellisExclusion = false;
+end
+
+switch string(mode)
+    case "db_target"
+        mlse = MLSE("DIR", [1,1], ...
+            "duobinary_output", 0, ...
+            "M", M, ...
+            "trellis_states", PAMmapper(M,0).levels, ...
+            "scale_mode", p.mlse_scale_mode, ...
+            "trellis_exclusion", trellisExclusion, ...
+            "trellis_state_mode", p.dbtgt_trellis_state_mode);
+    otherwise
+        mlse = MLSE("duobinary_output", 0, ...
+            "M", M, ...
+            "trellis_states", PAMmapper(M,0).levels, ...
+            "scale_mode", p.mlse_scale_mode, ...
+            "trellis_exclusion", trellisExclusion, ...
+            "trellis_state_mode", p.mlse_trellis_state_mode);
+end
+end
+
+function config = collectRecipeConfig(storageName, eqObject, p, options)
+config = struct();
+config.storage_name = char(storageName);
+config.recipe = "dsp_400g_recipe";
+config.eq_class = class(eqObject);
+config.run_id = options.dataTable.run_id;
+config.user_parameters = options.userParameters;
+config.parameters = p;
+end
+
+function description = resultDescription(prefix, options)
+dt = options.dataTable;
+description = sprintf('%s; run %d; PAM-%d; %.0f GBd; %.0f km', ...
+    prefix, dt.run_id, dt.pam_level, dt.symbolrate * 1e-9, dt.fiber_length);
+end
+
+function plotEqSignals(equalized_signal, Symbols, options, fignum, output_scale)
+showLevelScatter(equalized_signal .* output_scale, Symbols, ...
+    "fsym", options.fsym, ...
+    "fignum", fignum + 1, ...
+    "normalize", true);
+end

Functions/EQ_recipes/mpi_recipe_dev.m

@@ -22,7 +22,7 @@ run_a2_tracked_levels       = 0;
 run_a2_residual             = 0;
 run_a1                      = 0;
 run_tracking_adaptive       = 1;
-plot_output_signals         = 1;
+plot_output_signals         = options.debug_plots;
 
 %% Shared fixed EQ settings
 eq_sps = 2;
@@ -42,11 +42,12 @@ end
 Scpe_sig = preprocessSignal(Scpe_sig_raw, Symbols, options.fsym, ...
     "mode", "auto", ...
     "debug_plots", options.debug_plots);
+[Scpe_sig, Symbols, Tx_bits] = alignDspInputs(Scpe_sig, Symbols, Tx_bits, eq_sps);
 output = struct();
 
 %%
 if plot_output_signals
-    showLevelScatter(Scpe_sig, Symbols, ...
+    showLevelScatter(Scpe_sig.normalize("mode","rms"), Symbols, ...
         "fsym", options.fsym, ...
         "fignum", 400, ...
         "normalize", true);
@@ -283,9 +284,9 @@ if run_tracking_adaptive
         "dc_tracking_persistence_gain", 0, ...
         "dc_tracking_buffer_len", block_update, ...
         "optmize_mus", false, ...
-        "optimize_dc_tracking_params", false, ...
+        "optimize_dc_tracking_params", true, ...
         "dc_tracking_optimization_len", 2^15, ...
-        "dc_tracking_optimization_max_evals", 30, ...
+        "dc_tracking_optimization_max_evals", 20, ...
         "plot_mu_optimization", options.debug_plots, ...
         "save_debug", eq_save_debug);
 
@@ -296,6 +297,7 @@ if run_tracking_adaptive
         "dc_tracking", "adaptive", block_update);
     output.(char(storageName)) = ffe_results;
 
+    if options.debug_plots
         fignum = 106;
 
         eq_noise = equalized_signal - Symbols;
@@ -303,9 +305,10 @@ if run_tracking_adaptive
         showEQNoisePSD(eq_noise, "fignum", fignum, "displayname", dn,"colormode","diverging");
         ylim([-70 -30]);
         %
-    % dn = sprintf("FFE only; SIR: %d dB",options.dataTable.sir);
+        dn = sprintf("FFE only; SIR: %d dB",options.dataTable.sir);
-    % showEQNoisePSD(eq_noise_ffe, "fignum", fignum+1, "displayname", dn,"colormode","diverging");
+        showEQNoisePSD(eq_noise_ffe, "fignum", fignum+1, "displayname", dn,"colormode","diverging");
-    % ylim([-70 -30]);
+        ylim([-70 -30]);
+    end
 
     if plot_output_signals
 
@@ -387,6 +390,22 @@ function description = resultDescription(prefix,options)
     description = sprintf('%s; SIR %g dB',prefix,sir);
 end
 
+function [Scpe_sig, Symbols, Tx_bits] = alignDspInputs(Scpe_sig, Symbols, Tx_bits, sps)
+    nSymbols = min(length(Symbols), floor(length(Scpe_sig) / sps));
+    if nSymbols <= 0
+        error("mpi_recipe_dev:EmptyAlignedSignal", ...
+            "No overlapping samples remain after preprocessing and synchronization.");
+    end
+
+    Scpe_sig.signal = Scpe_sig.signal(1:sps*nSymbols);
+    Scpe_sig.signal = real(Scpe_sig.signal);
+    Symbols.signal = Symbols.signal(1:nSymbols,:);
+
+    if isprop(Tx_bits, "signal")
+        Tx_bits.signal = Tx_bits.signal(1:nSymbols,:);
+    end
+end
+
 function plotEqSignals(equalized_signal,Symbols,options,fignum,output_scale)
 
     showLevelScatter(equalized_signal .* output_scale, Symbols, ...

Functions/EQ_visuals/showMarkovDiagram.m

@@ -8,7 +8,8 @@ sequence = sequence(2+mod(1,length(sequence)):end); %filtered sequences often ha
 
 x = sequence;
 
-levels = sort(unique(x)).';              % or provide known 1x6 level values
+levels = sort(unique(x)).';
+
 
 [~,ix] = min(abs(x - levels),[],2);
 x = levels(ix);

Functions/Metrics/calc_air.m

@@ -29,8 +29,18 @@ end
 % CALC AIR
 %%% new implementation of AIR
 constellation = unique(reference_signal);
+% map reference symbols to constellation indices
 reference_idx = arrayfun(@(x) find(constellation == x, 1), reference_signal);
-ach_inf_rate = air_garcia_implementation(constellation',test_signal',reference_idx');
+
+% compute probability mass function of symbol indices
+[unique_idx,~,ic] = unique(reference_idx);
+counts = accumarray(ic,1);
+pmf = zeros(size(constellation));
+pmf(unique_idx) = counts / sum(counts);
+
+% ensure pmf corresponds to full set of constellation indices 1:N
+% (pmf already aligned because unique_idx are indices into constellation)
+ach_inf_rate = air_garcia_implementation(constellation',test_signal',reference_idx',pmf);
 
 
 function [data_,reference_]=trimseq(data,reference,skipstart,skip_end)

Tests/04_DSP/Equalizer/ML_MLSE_DUOBINARY_test.m

@@ -0,0 +1,84 @@
+classdef ML_MLSE_DUOBINARY_test < IMDDTestCase
+    methods (Test, TestTags = {'unit', 'fast', 'dsp', 'ml_mlse', 'duobinary'})
+        function constructorStoresDuobinaryConfiguration(testCase)
+            eq = ML_MLSE_DUOBINARY();
+
+            testCase.verifyEqual(class(eq), 'ML_MLSE_DUOBINARY');
+            testCase.verifyEqual(eq.sps, 2);
+            testCase.verifyEqual(eq.order, 15);
+            testCase.verifyEqual(eq.ber, []);
+            testCase.verifyEqual(eq.ber_dd, []);
+        end
+
+        function processReturnsPrecodeDomainAndFiniteDiagnostics(testCase)
+            [~, ~, txPrecoded, txEncoded] = makePam4Fixture(512);
+
+            eq = ML_MLSE_DUOBINARY( ...
+                "sps", 1, ...
+                "order", 1, ...
+                "len_tr", length(txEncoded), ...
+                "epochs_tr", 1, ...
+                "epochs_dd", 1, ...
+                "mu_tr", 0.01, ...
+                "mu_dd", 0.01, ...
+                "adaptive_mu", false, ...
+                "L", 1);
+
+            [detected, detectedViterbi] = eq.process(txEncoded, txPrecoded);
+
+            testCase.verifyEqual(length(detected), length(txEncoded));
+            testCase.verifyEqual(length(detectedViterbi), length(txEncoded));
+            testCase.verifyTrue(all(ismembertol(detected.signal, unique(txPrecoded.signal), 1e-12)));
+            testCase.verifySize(eq.ber, [1 1]);
+            testCase.verifySize(eq.ber_dd, [1 1]);
+            testCase.verifyTrue(isfinite(eq.ber(1)));
+            testCase.verifyTrue(isfinite(eq.ber_dd(1)));
+
+        end
+
+        function encodeDecodeInvertsPrecodingAwayFromInitialState(testCase)
+            [~, txSymbols, txPrecoded, ~] = makePam4Fixture(512);
+
+            encoded = Duobinary().encode(txPrecoded, "M", 4);
+            decoded = Duobinary().decode(encoded, "M", 4);
+
+            testCase.verifyEqual(decoded.signal(11:end-10), ...
+                txSymbols.signal(11:end-10), "AbsTol", 1e-12);
+        end
+
+        function resultWrapperReportsBothPrecodedAndOriginalBer(testCase)
+            [txBits, ~, txPrecoded, txEncoded] = makePam4Fixture(32000);
+
+            eq = ML_MLSE_DUOBINARY( ...
+                "sps", 1, ...
+                "order", 1, ...
+                "len_tr", length(txEncoded), ...
+                "epochs_tr", 1, ...
+                "epochs_dd", 1, ...
+                "mu_tr", 0.01, ...
+                "mu_dd", 0.01, ...
+                "adaptive_mu", false, ...
+                "L", 1);
+
+            results = ml_mlse(eq, 4, txEncoded, txPrecoded, txBits, ...
+                "precode_mode", db_mode.db_precoded);
+
+            testCase.verifyTrue(isfinite(results.metrics.BER));
+            testCase.verifyTrue(isfinite(results.metrics.BER_precoded));
+            testCase.verifyGreaterThanOrEqual(results.metrics.numBits, 0);
+            testCase.verifyGreaterThanOrEqual(results.metrics.numBitErr, 0);
+            testCase.verifyGreaterThanOrEqual(results.metrics.numBitErr_precoded, 0);
+        end
+    end
+end
+
+function [txBits, txSymbols, txPrecoded, txEncoded] = makePam4Fixture(nSymbols)
+    pattern = [0 0; 0 1; 1 1; 1 0];
+    bits = repmat(pattern, ceil(nSymbols / size(pattern, 1)), 1);
+    bits = bits(1:nSymbols, :);
+
+    txBits = Informationsignal(bits);
+    txSymbols = PAMmapper(4, 0).map(txBits);
+    txPrecoded = Duobinary().precode(txSymbols);
+    txEncoded = Duobinary().encode(txPrecoded, "M", 4);
+end

projects/Advanced_DSP_for_400G_IMDD_experiments/Auswertung_JLT/power_vs_wavelength.m

@@ -5,12 +5,12 @@ db = DBHandler("dataBase", [dataBase], "type", database_type);
 
 fp = QueryFilter();
 fp.where('power_state_info', 'pam_level','EQUALS', 4);
-fp.where('power_state_info', 'db_mode','EQUALS', 1);
+fp.where('power_state_info', 'db_mode','EQUALS', 0);
 % fp.where('power_state_info', 'fiber_length','EQUALS', 1);
 fp.where('power_state_info', 'is_mpi','EQUALS', 0);
 
 fields =  db.getTableFieldNames('power_state_info');
-% [dataTable,~] = db.queryDB(fp, fields);
+[dataTable,~] = db.queryDB(fp, fields);
 
 fiber_len       = unique(dataTable.fiber_length);
 cnt = 0;

projects/Diss/400G_revisit/PLOT_BER_BEST_ALGOS_WITH_DUOBINARY_SIGNALING.m

@@ -0,0 +1,343 @@
+%% 400G BER over bitrate: best normal algorithms plus duobinary signaling
+% Normal algorithms are reduced to the best BER per gross rate across
+% db_mode 0/1, pre-emphasis on/off, and BER/BER_precoded result variants.
+% Duobinary signaling uses db_mode = 2 and only the sequence-detection BER
+% stored in the BER field.
+
+clear; clc;
+
+%% 1) Query data
+
+selectedPamLevel = 8;
+selectedFiberLengthKm = 10;
+selectedWavelengthNm = 1310;
+selectedRopAttenuation = 0;     % set [] to use all ROP attenuation values
+selectedIsMpi = 0;              % set [] to use all entries
+
+normalDbModes = [double(db_mode.no_db), double(db_mode.db_precoded)];
+duobinaryDbMode = double(db_mode.db_encoded);
+
+maxBerForPlot = 0.5;
+showRawEntries = false;
+showBestLine = true;
+
+algoStyles = defaultAlgorithmStyles();
+
+db = DBHandler( ...
+    "dataBase", "labor_highspeed", ...
+    "type", "mysql", ...
+    "server", "192.168.178.192", ...
+    "user", "silas", ...
+    "password", "silas");
+db.refresh();
+
+fp = QueryFilter();
+fp.where('Runs', 'fiber_length', 'EQUALS', selectedFiberLengthKm);
+fp.where('Runs', 'pam_level', 'EQUALS', selectedPamLevel);
+fp.where('Runs', 'wavelength', 'EQUALS', selectedWavelengthNm);
+if ~isempty(selectedRopAttenuation)
+    fp.where('Runs', 'rop_attenuation', 'EQUALS', selectedRopAttenuation);
+end
+if ~isempty(selectedIsMpi)
+    fp.where('Runs', 'is_mpi', 'EQUALS', selectedIsMpi);
+end
+
+selectedFields = db.getTableFieldNames('dashboard_ungrouped_alltime');
+selectedFields = appendMissingFields(selectedFields, ...
+    {'Runs.precomp_amp'; 'Runs.is_mpi'});
+selectedFields = selectedFields(:);
+
+[rawData, query] = db.queryDB(fp, selectedFields);
+disp(query);
+fprintf("Fetched %d 400G result rows.\n", height(rawData));
+
+%% 2) Clean data and build the five plotted curves
+
+data = rawData;
+numericFields = ["result_id", "run_id", "eq_id", "bitrate", "grossrate", ...
+    "symbolrate", "pam_level", "wavelength", "fiber_length", "db_mode", ...
+    "rop_attenuation", "precomp_amp", "is_mpi", "numBits", "numBitErr", ...
+    "BER", "numBitErr_precoded", "BER_precoded", "STD", "STDrx", ...
+    "GMI", "AIR", "NGMI", "EVM", "Alpha"];
+for fieldIdx = 1:numel(numericFields)
+    fieldName = numericFields(fieldIdx);
+    if ismember(fieldName, string(data.Properties.VariableNames))
+        data.(char(fieldName)) = numericColumn(data.(char(fieldName)));
+    end
+end
+
+if ~ismember("precomp_amp", string(data.Properties.VariableNames))
+    warning("plot_best_algos:NoPrecompAmp", ...
+        "Runs.precomp_amp was not returned. Falling back to pre_emphasis = (db_mode == 0).");
+    data.pre_emphasis = data.db_mode == double(db_mode.no_db);
+else
+    data.pre_emphasis = derivePreEmphasis(data.precomp_amp, data.db_mode);
+end
+
+normalRows = data(ismember(data.db_mode, normalDbModes), :);
+normalPlotData = buildNormalMetricRows(normalRows);
+normalPlotData = normalPlotData(isfinite(normalPlotData.BER_plot) & ...
+    normalPlotData.BER_plot > 0 & normalPlotData.BER_plot < maxBerForPlot, :);
+
+duobinaryRows = data(data.db_mode == duobinaryDbMode, :);
+if ismember("equalizer_structure", string(duobinaryRows.Properties.VariableNames))
+    duobinaryRows = duobinaryRows( ...
+        equalizerMask(duobinaryRows.equalizer_structure, ...
+        equalizer_structure.db_encoded), :);
+end
+duobinaryPlotData = buildDuobinarySignalingRows(duobinaryRows);
+duobinaryPlotData = duobinaryPlotData(isfinite(duobinaryPlotData.BER_plot) & ...
+    duobinaryPlotData.BER_plot > 0 & ...
+    duobinaryPlotData.BER_plot < maxBerForPlot, :);
+
+plotData = [normalPlotData; duobinaryPlotData];
+if isempty(plotData)
+    warning("plot_best_algos:NoRows", ...
+        "No rows remain after length/PAM/wavelength/BER filtering.");
+    return
+end
+
+plotData.bitrate_Gbps = plotData.bitrate .* 1e-9;
+plotData.grossrate_Gbps = plotData.grossrate .* 1e-9;
+
+fprintf("Remaining candidate BER rows: %d\n", height(plotData));
+disp(groupcounts(plotData, ["algorithm_key", "db_mode", "pre_emphasis", "precode"]));
+
+bestPlotData = bestBerByAlgorithmAndGrossRate(plotData);
+fprintf("Keeping %d best-BER rows across algorithm/gross-rate groups.\n", ...
+    height(bestPlotData));
+disp(groupcounts(bestPlotData, "algorithm_key"));
+
+%% 3) Plot one figure with five lines
+
+availableStyles = algoStyles(hasAlgorithmRows(bestPlotData, algoStyles), :);
+if isempty(availableStyles)
+    warning("plot_best_algos:NoSelectedAlgorithms", ...
+        "None of the configured algorithm styles match the queried rows.");
+    return
+end
+
+fig = figure(); clf;
+ax = axes(fig); hold(ax, "on");
+
+for styleIdx = 1:height(availableStyles)
+    style = availableStyles(styleIdx, :);
+    rowMask = bestPlotData.algorithm_key == style.algorithm_key;
+    if ~any(rowMask)
+        continue
+    end
+
+    algoData = sortrows(bestPlotData(rowMask, :), "grossrate_Gbps");
+
+    if showRawEntries
+        scatter(ax, algoData.grossrate_Gbps, algoData.BER_plot, ...
+            9, ...
+            "Marker", ".", ...
+            "MarkerEdgeColor", style.color, ...
+            "MarkerFaceColor", style.color, ...
+            "MarkerEdgeAlpha", 0.25, ...
+            "MarkerFaceAlpha", 0.25, ...
+            "HandleVisibility", "off");
+    end
+
+    if showBestLine
+        plot(ax, algoData.grossrate_Gbps, algoData.BER_plot, ...
+            "LineStyle", style.lineStyle, ...
+            "Marker", style.marker, ...
+            "MarkerSize", 5, ...
+            "LineWidth", 1.5, ...
+            "Color", style.color, ...
+            "MarkerFaceColor", style.markerFaceColor, ...
+            "MarkerEdgeColor", style.color, ...
+            "DisplayName", style.name);
+    end
+end
+
+yline(ax, [2.2e-4, 4.85e-3, 2e-2], ...
+    "LineWidth", 1, ...
+    "LineStyle", "--", ...
+    "Color", [0.25 0.25 0.25], ...
+    "HandleVisibility", "off");
+
+title(ax, sprintf("PAM-%d, %.0f km, %.0f nm", ...
+    selectedPamLevel, selectedFiberLengthKm, selectedWavelengthNm));
+xlabel(ax, "Gross rate [Gb/s]");
+ylabel(ax, "BER");
+set(ax, "YScale", "log");
+ylim(ax, [1e-5, maxBerForPlot]);
+grid(ax, "on");
+box(ax, "on");
+
+xTicks = unique(bestPlotData.grossrate_Gbps(isfinite(bestPlotData.grossrate_Gbps)));
+if ~isempty(xTicks)
+    xticks(ax, xTicks);
+    xlim(ax, [min(xTicks), max(xTicks)]);
+end
+
+legend(ax, "Location", "best", "Interpreter", "none");
+if exist("beautifyBERplot", "file")
+    beautifyBERplot("logscale", true, "setcolors", false, ...
+        "setmarkers", false, "changemarkers", false);
+end
+
+set(fig, "Position", 1e3 .* [0.1000 0.5500 0.7200 0.4200]);
+
+%% Local helpers
+
+function fields = appendMissingFields(fields, extraFields)
+fields = cellstr(fields);
+extraFields = cellstr(extraFields);
+for idx = 1:numel(extraFields)
+    if ~any(strcmp(fields, extraFields{idx}))
+        fields{end+1, 1} = extraFields{idx}; %#ok<AGROW>
+    end
+end
+end
+
+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 styles = defaultAlgorithmStyles()
+styles = table( ...
+    ["vnle"; ...
+     "vnle_pf_mlse"; ...
+     "vnle_db_mlse"; ...
+     "ml_mlse"; ...
+     "db_encoded"], ...
+    [equalizer_structure.vnle; ...
+     equalizer_structure.vnle_pf_mlse; ...
+     equalizer_structure.vnle_db_mlse; ...
+     equalizer_structure.ml_mlse; ...
+     equalizer_structure.db_encoded], ...
+    ["VNLE"; ...
+     "VNLE + PF + MLSE"; ...
+     "VNLE DBt. + MLSE"; ...
+     "ML pre-EQ + Viterbi"; ...
+     "Duobinary signaling"], ...
+    ["o"; "square"; "diamond"; "^"; "v"], ...
+    ["-"; "-"; "-"; "-"; "-"], ...
+    ["w"; "w"; "w"; "w"; "w"], ...
+    [clr.Paired.red; ...
+     clr.Paired.green; ...
+     clr.Paired.blue; ...
+     clr.Paired.purple; ...
+     clr.Paired.orange], ...
+    'VariableNames', ["algorithm_key", "eq", "name", "marker", ...
+    "lineStyle", "markerFaceColor", "color"]);
+end
+
+function preEmphasis = derivePreEmphasis(precompAmp, dbMode)
+preEmphasis = false(size(dbMode));
+
+validPrecomp = isfinite(precompAmp);
+preEmphasis(validPrecomp) = precompAmp(validPrecomp) > -45;
+
+missingPrecomp = ~validPrecomp;
+preEmphasis(missingPrecomp) = dbMode(missingPrecomp) == double(db_mode.no_db);
+end
+
+function plotData = buildNormalMetricRows(data)
+baseRows = data(isfinite(data.BER), :);
+baseRows.precode = false(height(baseRows), 1);
+baseRows.BER_plot = baseRows.BER;
+baseRows.algorithm_key = algorithmKeyFromEqualizer(baseRows.equalizer_structure);
+baseRows = baseRows(baseRows.algorithm_key ~= "", :);
+
+if ismember("BER_precoded", string(data.Properties.VariableNames))
+    precodedRows = data(isfinite(data.BER_precoded), :);
+    precodedRows.precode = true(height(precodedRows), 1);
+    precodedRows.BER_plot = precodedRows.BER_precoded;
+    precodedRows.algorithm_key = algorithmKeyFromEqualizer( ...
+        precodedRows.equalizer_structure);
+    precodedRows = precodedRows(precodedRows.algorithm_key ~= "", :);
+    plotData = [baseRows; precodedRows];
+else
+    warning("plot_best_algos:NoPrecodedBer", ...
+        "BER_precoded was not returned. Plotting only BER rows for normal algorithms.");
+    plotData = baseRows;
+end
+end
+
+function plotData = buildDuobinarySignalingRows(data)
+plotData = data(isfinite(data.BER), :);
+plotData.precode = false(height(plotData), 1);
+plotData.BER_plot = plotData.BER;
+plotData.algorithm_key = repmat("db_encoded", height(plotData), 1);
+end
+
+function algorithmKey = algorithmKeyFromEqualizer(equalizerColumn)
+eqNumeric = equalizerNumeric(equalizerColumn);
+algorithmKey = strings(size(eqNumeric));
+
+algorithmKey(eqNumeric == enumValue(equalizer_structure.vnle)) = "vnle";
+algorithmKey(eqNumeric == enumValue(equalizer_structure.vnle_pf_mlse)) = ...
+    "vnle_pf_mlse";
+algorithmKey(eqNumeric == enumValue(equalizer_structure.vnle_db_mlse)) = ...
+    "vnle_db_mlse";
+algorithmKey(eqNumeric == enumValue(equalizer_structure.ml_mlse)) = "ml_mlse";
+end
+
+function mask = equalizerMask(equalizerColumn, eqValue)
+eqNumeric = equalizerNumeric(equalizerColumn);
+mask = eqNumeric == enumValue(eqValue);
+end
+
+function eqNumeric = equalizerNumeric(equalizerColumn)
+if isa(equalizerColumn, "equalizer_structure")
+    eqNumeric = double(equalizerColumn);
+elseif isnumeric(equalizerColumn)
+    eqNumeric = double(equalizerColumn);
+else
+    equalizerString = string(equalizerColumn);
+    eqNumeric = str2double(equalizerString);
+
+    enumNames = ["vnle", "ffe", "dfe", "vnle_pf_mlse", ...
+        "vnle_db_mlse", "db_encoded", "ml_mlse"];
+    enumValues = [ ...
+        enumValue(equalizer_structure.vnle), ...
+        enumValue(equalizer_structure.ffe), ...
+        enumValue(equalizer_structure.dfe), ...
+        enumValue(equalizer_structure.vnle_pf_mlse), ...
+        enumValue(equalizer_structure.vnle_db_mlse), ...
+        enumValue(equalizer_structure.db_encoded), ...
+        enumValue(equalizer_structure.ml_mlse)];
+
+    for idx = 1:numel(enumNames)
+        missingNumeric = isnan(eqNumeric);
+        eqNumeric(missingNumeric & equalizerString == enumNames(idx)) = ...
+            enumValues(idx);
+    end
+end
+end
+
+function value = enumValue(enumEntry)
+value = double(enumEntry);
+end
+
+function bestData = bestBerByAlgorithmAndGrossRate(data)
+groupVars = ["algorithm_key", "grossrate_Gbps"];
+groupId = findgroups(data(:, groupVars));
+keepIdx = NaN(max(groupId), 1);
+
+for curGroup = 1:max(groupId)
+    rowIdx = find(groupId == curGroup);
+    [~, localBestIdx] = min(data.BER_plot(rowIdx));
+    keepIdx(curGroup) = rowIdx(localBestIdx);
+end
+
+bestData = sortrows(data(keepIdx, :), groupVars);
+end
+
+function keep = hasAlgorithmRows(data, algoStyles)
+keep = false(height(algoStyles), 1);
+for idx = 1:height(algoStyles)
+    keep(idx) = any(data.algorithm_key == algoStyles.algorithm_key(idx));
+end
+end

projects/Diss/400G_revisit/PLOT_BER_VS_ALGO.m

@@ -0,0 +1,294 @@
+%% 400G BER over bitrate from labor_highspeed.dashboard_ungrouped_alltime
+% 1) gather all BER entries for one PAM format and fiber length
+% 2) derive pre-emphasis and precoding groups
+% 3) plot one bitrate-vs-BER tile per equalizer structure
+
+clear; clc;
+
+%% 1) Gather data
+
+selectedPamLevel = 4;
+selectedFiberLengthKm = 10;
+selectedWavelength =1310;         % set [] to use all wavelengths
+selectedRopAttenuation = [];     % set [] to use all ROP attenuation values
+selectedIsMpi = [];              % set [] to use all entries
+
+maxBerForPlot = 0.5;
+showRawEntries = true;
+showMedianLine = true;
+
+eqStyles = defaultEqualizerStyles();
+comboStyles = defaultCombinationStyles();
+
+db = DBHandler( ...
+    "dataBase", "labor_highspeed", ...
+    "type", "mysql");
+db.refresh();
+
+fp = QueryFilter();
+fp.where('Runs', 'fiber_length', 'EQUALS', selectedFiberLengthKm);
+fp.where('Runs', 'pam_level', 'EQUALS', selectedPamLevel);
+if ~isempty(selectedWavelength)
+    fp.where('Runs', 'wavelength', 'EQUALS', selectedWavelength);
+end
+if ~isempty(selectedRopAttenuation)
+    fp.where('Runs', 'rop_attenuation', 'EQUALS', selectedRopAttenuation);
+end
+if ~isempty(selectedIsMpi)
+    fp.where('Runs', 'is_mpi', 'EQUALS', selectedIsMpi);
+end
+
+selectedFields = db.getTableFieldNames('dashboard_ungrouped_alltime');
+selectedFields = [selectedFields; {'Runs.precomp_amp'; 'Runs.is_mpi'}];
+selectedFields = selectedFields(:);
+
+[rawData, query] = db.queryDB(fp, selectedFields);
+disp(query);
+fprintf("Fetched %d 400G result rows.\n", height(rawData));
+
+%% 2) Clean data and derive analysis groups
+
+data = rawData;
+numericFields = ["result_id", "run_id", "eq_id", "bitrate", "grossrate", ...
+    "symbolrate", "pam_level", "wavelength", "fiber_length", "db_mode", ...
+    "rop_attenuation", "precomp_amp", "is_mpi", "numBits", "numBitErr", ...
+    "BER", "numBitErr_precoded", "BER_precoded", "STD", "STDrx", ...
+    "GMI", "AIR", "NGMI", "EVM", "Alpha"];
+for fieldIdx = 1:numel(numericFields)
+    fieldName = numericFields(fieldIdx);
+    if ismember(fieldName, string(data.Properties.VariableNames))
+        data.(char(fieldName)) = numericColumn(data.(char(fieldName)));
+    end
+end
+
+if ~ismember("precomp_amp", string(data.Properties.VariableNames))
+    warning("analyze_db:NoPrecompAmp", ...
+        "Runs.precomp_amp was not returned. Falling back to pre_emphasis = (db_mode == 0).");
+    data.pre_emphasis = data.db_mode == 0;
+else
+    data.pre_emphasis = derivePreEmphasis(data.precomp_amp, data.db_mode);
+end
+
+plotData = buildBerMetricRows(data);
+plotData = plotData(isfinite(plotData.BER_plot) & ...
+    plotData.BER_plot > 0 & plotData.BER_plot < maxBerForPlot, :);
+
+if isempty(plotData)
+    warning("analyze_db:NoRows", ...
+        "No rows remain after fiber/PAM/BER filtering.");
+    return
+end
+
+plotData.bitrate_Gbps = plotData.bitrate .* 1e-9;
+plotData.grossrate_Gbps = plotData.grossrate .* 1e-9;
+
+fprintf("Remaining plotted BER rows: %d\n", height(plotData));
+disp(groupcounts(plotData, ["equalizer_structure", "pre_emphasis", "precode"]));
+
+bestPlotData = bestBerByBitrateAndGroup(plotData);
+fprintf("Keeping %d best-BER rows across bitrate/EQ/pre-emphasis/precode groups.\n", ...
+    height(bestPlotData));
+
+%% 3) Plot bitrate versus BER
+
+availableEqStyles = eqStyles(hasEqualizerRows(bestPlotData, eqStyles), :);
+if isempty(availableEqStyles)
+    warning("analyze_db:NoSelectedEqualizers", ...
+        "None of the configured equalizer styles match the queried rows.");
+    return
+end
+
+fig = figure(401); clf;
+tiledlayout(1, height(availableEqStyles), ...
+    "TileSpacing", "compact", ...
+    "Padding", "compact");
+
+for eqIdx = 1:height(availableEqStyles)
+    eqStyle = availableEqStyles(eqIdx, :);
+    ax = nexttile; hold(ax, "on");
+    eqMask = equalizerMask(bestPlotData.equalizer_structure, eqStyle.eq);
+
+    for comboIdx = 1:height(comboStyles)
+        comboStyle = comboStyles(comboIdx, :);
+        rowMask = eqMask & ...
+            bestPlotData.pre_emphasis == comboStyle.pre_emphasis & ...
+            bestPlotData.precode == comboStyle.precode;
+
+        if ~any(rowMask)
+            continue
+        end
+
+        comboData = sortrows(bestPlotData(rowMask, :), "bitrate_Gbps");
+        comboColor = emphasisColor(eqStyle.color, comboStyle.pre_emphasis);
+        label = sprintf("%s, %s", eqStyle.name, comboStyle.name);
+
+        if showRawEntries
+            scatter(ax, comboData.bitrate_Gbps, comboData.BER_plot, ...
+                5, ...
+                "Marker", '.', ...
+                "MarkerEdgeColor", comboColor, ...
+                "MarkerFaceColor", comboColor, ...
+                "MarkerEdgeAlpha", 0.25, ...
+                "MarkerFaceAlpha", 0.25, ...
+                "HandleVisibility", "off");
+        end
+
+        if showMedianLine
+            summaryTable = summarizeBerByBitrate(comboData);
+            plot(ax, summaryTable.bitrate_Gbps, summaryTable.median_BER_plot, ...
+                "LineStyle", comboStyle.lineStyle, ...
+                "Marker", eqStyle.marker, ...
+                "MarkerSize", 4, ...
+                "LineWidth", 1.4, ...
+                "Color", comboColor, ...
+                "MarkerFaceColor", comboStyle.markerFaceColor, ...
+                "MarkerEdgeColor", comboColor, ...
+                "DisplayName", label);
+        end
+    end
+
+    yline(ax, [2.2e-4, 4.85e-3, 2e-2], ...
+        "LineWidth", 1, ...
+        "LineStyle", "--", ...
+        "Color", [0.25 0.25 0.25], ...
+        "HandleVisibility", "off");
+
+    title(ax, sprintf("PAM-%d, %s", selectedPamLevel, eqStyle.name));
+    xlabel(ax, "Gross rate [Gb/s]");
+    ylabel(ax, "BER");
+    set(ax, "YScale", "log");
+    ylim(ax, [1e-5, maxBerForPlot]);
+    grid(ax, "on");
+    box(ax, "on");
+
+    xTicks = unique(bestPlotData.bitrate_Gbps(isfinite(bestPlotData.bitrate_Gbps)));
+    if ~isempty(xTicks)
+        xticks(ax, xTicks);
+        xlim(ax, [min(xTicks), max(xTicks)]);
+    end
+
+    legend(ax, "Location", "best", "Interpreter", "none");
+    if exist("beautifyBERplot", "file")
+        beautifyBERplot("logscale", true, "setcolors", false, ...
+            "setmarkers", false, "changemarkers", false);
+    end
+end
+
+set(fig, "Position", 1e3 .* [0.1000 0.5500 1.4113 0.3200]);
+
+%% 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 styles = defaultEqualizerStyles()
+    styles = table( ...
+        [equalizer_structure.vnle; ...
+         equalizer_structure.vnle_pf_mlse; ...
+         equalizer_structure.vnle_db_mlse; ...
+         equalizer_structure.ml_mlse], ...
+        ["VNLE"; ...
+         "VNLE + PF + MLSE"; ...
+         "VNLE DBt. + MLSE"; ...
+         "ML pre-EQ + Viterbi"], ...
+        ["o"; "square"; "diamond"; "^"], ...
+        [clr.Paired.red; ...
+         clr.Paired.green; ...
+         clr.Paired.blue; ...
+         clr.Paired.purple], ...
+        'VariableNames', ["eq", "name", "marker", "color"]);
+end
+
+function styles = defaultCombinationStyles()
+    styles = table( ...
+        [true; true; false; false], ...
+        [true; false; true; false], ...
+        ["w/ pre-emph., w/ precode"; ...
+         "w/ pre-emph., w/o precode"; ...
+         "w/o pre-emph., w/ precode"; ...
+         "w/o pre-emph., w/o precode"], ...
+        ["--"; "-"; "--"; "-"], ...
+        ["w"; "w"; "none"; "none"], ...
+        'VariableNames', ["pre_emphasis", "precode", "name", ...
+        "lineStyle", "markerFaceColor"]);
+end
+
+function preEmphasis = derivePreEmphasis(precompAmp, dbMode)
+    preEmphasis = false(size(dbMode));
+
+    validPrecomp = isfinite(precompAmp);
+    % In the 400G measurement scripts, -50 dB is the low/no-pre-emphasis
+    % setting, while -38/-37/-34 dB are the active pre-emphasis settings.
+    preEmphasis(validPrecomp) = precompAmp(validPrecomp) > -45;
+
+    missingPrecomp = ~validPrecomp;
+    preEmphasis(missingPrecomp) = dbMode(missingPrecomp) == 0;
+end
+
+function plotData = buildBerMetricRows(data)
+    baseRows = data(isfinite(data.BER), :);
+    baseRows.precode = false(height(baseRows), 1);
+    baseRows.BER_plot = baseRows.BER;
+
+    if ismember("BER_precoded", string(data.Properties.VariableNames))
+        precodedRows = data(isfinite(data.BER_precoded), :);
+        precodedRows.precode = true(height(precodedRows), 1);
+        precodedRows.BER_plot = precodedRows.BER_precoded;
+        plotData = [baseRows; precodedRows];
+    else
+        warning("analyze_db:NoPrecodedBer", ...
+            "BER_precoded was not returned. Plotting only precode = 0 rows.");
+        plotData = baseRows;
+    end
+end
+
+function mask = equalizerMask(equalizerColumn, eqValue)
+    if isa(equalizerColumn, "equalizer_structure")
+        mask = equalizerColumn == eqValue;
+    elseif isnumeric(equalizerColumn)
+        mask = double(equalizerColumn) == double(int32(eqValue));
+    else
+        mask = string(equalizerColumn) == string(eqValue);
+    end
+end
+
+function keep = hasEqualizerRows(data, eqStyles)
+    keep = false(height(eqStyles), 1);
+    for idx = 1:height(eqStyles)
+        keep(idx) = any(equalizerMask(data.equalizer_structure, eqStyles.eq(idx)));
+    end
+end
+
+function summaryTable = summarizeBerByBitrate(data)
+    summaryTable = groupsummary(data, "bitrate_Gbps", "median", "BER_plot");
+    summaryTable = sortrows(summaryTable, "bitrate_Gbps");
+end
+
+function bestData = bestBerByBitrateAndGroup(data)
+    groupVars = ["equalizer_structure", "pre_emphasis", "precode", "bitrate_Gbps"];
+    groupId = findgroups(data(:, groupVars));
+    keepIdx = NaN(max(groupId), 1);
+
+    for curGroup = 1:max(groupId)
+        rowIdx = find(groupId == curGroup);
+        [~, localBestIdx] = min(data.BER_plot(rowIdx));
+        keepIdx(curGroup) = rowIdx(localBestIdx);
+    end
+
+    bestData = sortrows(data(keepIdx, :), groupVars);
+end
+
+function color = emphasisColor(baseColor, preEmphasis)
+    if preEmphasis
+        color = 0.65 .* baseColor + 0.35;
+    else
+        color = baseColor;
+    end
+end

projects/Diss/400G_revisit/PLOT_DUOBINARY_DETECTION_BER_VS_BITRATE.m

@@ -0,0 +1,242 @@
+%% Duobinary transmission: BER over bitrate for detection algorithms
+% DB transmission means the sequence was precoded and encoded at the Tx
+% (Runs.db_mode = db_mode.db_encoded). In this stored-result convention:
+%   BER          -> VNLE + MLSE
+%   BER_precoded -> VNLE + memoryless detection
+
+clear; clc;
+
+%% 1) Query data
+
+selectedPamLevels = [4, 6, 8];
+selectedFiberLengthKm = 10;
+selectedWavelengthNm = 1310;
+selectedRopAttenuation = 0;
+selectedIsMpi = 0;
+selectedDbMode = db_mode.db_encoded;
+selectedEqualizerStructure = equalizer_structure.db_encoded;
+
+maxBerForPlot = 0.5;
+showRawEntries = false;
+showBestLine = true;
+
+db = DBHandler( ...
+    "dataBase", "labor_highspeed", ...
+    "type", "mysql", ...
+    "server", "192.168.178.192", ...
+    "user", "silas", ...
+    "password", "silas");
+db.refresh();
+
+fp = QueryFilter();
+fp.where('Runs', 'fiber_length', 'EQUALS', selectedFiberLengthKm);
+fp.where('Runs', 'wavelength', 'EQUALS', selectedWavelengthNm);
+fp.where('Runs', 'rop_attenuation', 'EQUALS', selectedRopAttenuation);
+fp.where('Runs', 'is_mpi', 'EQUALS', selectedIsMpi);
+fp.where('Runs', 'db_mode', 'EQUALS', double(selectedDbMode));
+
+selectedFields = db.getTableFieldNames('dashboard_ungrouped_alltime');
+selectedFields = appendMissingFields(selectedFields, {'Runs.is_mpi'});
+selectedFields = selectedFields(:);
+
+[rawData, query] = db.queryDB(fp, selectedFields);
+disp(query);
+fprintf("Fetched %d duobinary result rows.\n", height(rawData));
+
+%% 2) Clean and reshape BER metrics
+
+data = rawData;
+numericFields = ["result_id", "run_id", "eq_id", "bitrate", "grossrate", ...
+    "symbolrate", "pam_level", "wavelength", "fiber_length", "db_mode", ...
+    "rop_attenuation", "is_mpi", "numBits", "numBitErr", "BER", ...
+    "numBitErr_precoded", "BER_precoded", "STD", "STDrx", ...
+    "GMI", "AIR", "NGMI", "EVM", "Alpha"];
+for fieldIdx = 1:numel(numericFields)
+    fieldName = numericFields(fieldIdx);
+    if ismember(fieldName, string(data.Properties.VariableNames))
+        data.(char(fieldName)) = numericColumn(data.(char(fieldName)));
+    end
+end
+
+data = data(ismember(data.pam_level, selectedPamLevels), :);
+if ismember("equalizer_structure", string(data.Properties.VariableNames)) && ...
+        ~isempty(selectedEqualizerStructure)
+    data = data(equalizerMask(data.equalizer_structure, selectedEqualizerStructure), :);
+end
+
+plotData = buildDetectionMetricRows(data);
+plotData = plotData(isfinite(plotData.BER_plot) & ...
+    plotData.BER_plot > 0 & plotData.BER_plot < maxBerForPlot, :);
+
+if isempty(plotData)
+    warning("plot_duobinary_detection:NoRows", ...
+        "No rows remain after duobinary/PAM/wavelength/BER filtering.");
+    return
+end
+
+plotData.bitrate_Gbps = plotData.bitrate .* 1e-9;
+plotData = sortrows(plotData, ...
+    ["pam_level", "wavelength", "detection_type", "bitrate_Gbps", "run_id"]);
+
+fprintf("Remaining plotted BER rows: %d\n", height(plotData));
+disp(groupcounts(plotData, ["pam_level", "wavelength", "detection_type"]));
+
+%% 3) Plot BER versus bitrate
+
+detectionStyles = defaultDetectionStyles();
+availablePamLevels = selectedPamLevels(ismember(selectedPamLevels, unique(plotData.pam_level).'));
+
+fig = figure(430); clf;
+tiledlayout(1, numel(availablePamLevels), ...
+    "TileSpacing", "compact", ...
+    "Padding", "compact");
+
+for pamIdx = 1:numel(availablePamLevels)
+    pamLevel = availablePamLevels(pamIdx);
+    ax = nexttile; hold(ax, "on");
+    pamMask = plotData.pam_level == pamLevel;
+
+    for styleIdx = 1:height(detectionStyles)
+        style = detectionStyles(styleIdx, :);
+        rowMask = pamMask & plotData.detection_type == style.detection_type;
+        if ~any(rowMask)
+            continue
+        end
+
+        detectionData = sortrows(plotData(rowMask, :), "bitrate_Gbps");
+        summaryTable = summarizeBerByBitrate(detectionData);
+
+        if showRawEntries
+            scatter(ax, detectionData.bitrate_Gbps, detectionData.BER_plot, ...
+                9, ...
+                "Marker", ".", ...
+                "MarkerEdgeColor", style.color, ...
+                "MarkerFaceColor", style.color, ...
+                "MarkerEdgeAlpha", 0.25, ...
+                "MarkerFaceAlpha", 0.25, ...
+                "HandleVisibility", "off");
+        end
+
+        if showBestLine
+            plot(ax, summaryTable.bitrate_Gbps, summaryTable.BER_plot, ...
+                "LineStyle", style.lineStyle, ...
+                "Marker", style.marker, ...
+                "MarkerSize", 5, ...
+                "LineWidth", 1.4, ...
+                "Color", style.color, ...
+                "MarkerFaceColor", style.markerFaceColor, ...
+                "MarkerEdgeColor", style.color, ...
+                "DisplayName", style.name);
+        end
+    end
+
+    yline(ax, [2.2e-4, 4.85e-3, 2e-2], ...
+        "LineWidth", 1, ...
+        "LineStyle", "--", ...
+        "Color", [0.25 0.25 0.25], ...
+        "HandleVisibility", "off");
+
+    title(ax, sprintf("PAM-%d, %.0f nm", pamLevel, selectedWavelengthNm));
+    xlabel(ax, "Bitrate [Gb/s]");
+    ylabel(ax, "BER");
+    set(ax, "YScale", "log");
+    ylim(ax, [1e-5, maxBerForPlot]);
+    grid(ax, "on");
+    box(ax, "on");
+
+    xTicks = unique(plotData.bitrate_Gbps(pamMask & isfinite(plotData.bitrate_Gbps)));
+    if ~isempty(xTicks)
+        xticks(ax, xTicks);
+        xlim(ax, [min(xTicks), max(xTicks)]);
+    end
+
+    legend(ax, "Location", "best", "Interpreter", "none");
+    if exist("beautifyBERplot", "file")
+        beautifyBERplot("logscale", true, "setcolors", false, ...
+            "setmarkers", false, "changemarkers", false);
+    end
+end
+
+set(fig, "Position", 1e3 .* [0.1000 0.5500 1.4113 0.3200]);
+
+%% Local helpers
+
+function fields = appendMissingFields(fields, extraFields)
+fields = cellstr(fields);
+extraFields = cellstr(extraFields);
+for idx = 1:numel(extraFields)
+    if ~any(strcmp(fields, extraFields{idx}))
+        fields{end+1, 1} = extraFields{idx}; %#ok<AGROW>
+    end
+end
+end
+
+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 plotData = buildDetectionMetricRows(data)
+baseRows = data(isfinite(data.BER), :);
+baseRows.detection_type = repmat("VNLE + MLSE", height(baseRows), 1);
+baseRows.BER_plot = baseRows.BER;
+
+if ismember("BER_precoded", string(data.Properties.VariableNames))
+    memorylessRows = data(isfinite(data.BER_precoded), :);
+    memorylessRows.detection_type = repmat("VNLE + memoryless", ...
+        height(memorylessRows), 1);
+    memorylessRows.BER_plot = memorylessRows.BER_precoded;
+    plotData = [baseRows; memorylessRows];
+else
+    warning("plot_duobinary_detection:NoPrecodedBer", ...
+        "BER_precoded was not returned. Plotting only VNLE + MLSE rows.");
+    plotData = baseRows;
+end
+end
+
+function styles = defaultDetectionStyles()
+styles = table( ...
+    ["VNLE + MLSE"; "VNLE + memoryless"], ...
+    ["VNLE + MLSE"; "VNLE + memoryless"], ...
+    ["o"; "square"], ...
+    ["-"; "--"], ...
+    ["w"; "none"], ...
+    [clr.Paired.blue; clr.Paired.orange], ...
+    'VariableNames', ["detection_type", "name", "marker", ...
+    "lineStyle", "markerFaceColor", "color"]);
+end
+
+function mask = equalizerMask(equalizerColumn, eqValue)
+if isa(equalizerColumn, "equalizer_structure")
+    mask = equalizerColumn == eqValue;
+elseif isnumeric(equalizerColumn)
+    mask = double(equalizerColumn) == enumValue(eqValue);
+else
+    equalizerString = string(equalizerColumn);
+    numericEqualizer = str2double(equalizerString);
+    mask = equalizerString == string(eqValue) | numericEqualizer == enumValue(eqValue);
+end
+end
+
+function value = enumValue(enumEntry)
+value = double(enumEntry);
+end
+
+function summaryTable = summarizeBerByBitrate(data)
+groupId = findgroups(data.bitrate_Gbps);
+keepIdx = NaN(max(groupId), 1);
+
+for curGroup = 1:max(groupId)
+    rowIdx = find(groupId == curGroup);
+    [~, localBestIdx] = min(data.BER_plot(rowIdx));
+    keepIdx(curGroup) = rowIdx(localBestIdx);
+end
+
+summaryTable = data(keepIdx, :);
+summaryTable = sortrows(summaryTable, "bitrate_Gbps");
+end

projects/Diss/400G_revisit/investigate_400g_algorithms.m

@@ -0,0 +1,361 @@
+% === 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','EQUALS', 330e9);
+fp.where('Runs','pam_level','EQUALS', 4);
+fp.where('Runs','rop_attenuation','EQUALS', 0);
+fp.where('Runs','is_mpi','EQUALS', 0);
+fp.where('Runs', 'db_mode','EQUALS', 2);
+
+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();
+dsp_options.userParameters.run_ml_mlse_db = true;
+dsp_options.userParameters.run_mlse_db = true;
+
+% Enable/disable equalizer branches.
+% dsp_options.userParameters.run_ffe = false;
+% dsp_options.userParameters.run_vnle = false;
+% dsp_options.userParameters.run_dfe = false;
+% dsp_options.userParameters.run_vnle_mlse = true;
+% dsp_options.userParameters.run_dbtgt = true;
+
+% Examples for parameter loops. DataStorage expands every vector-valued field.
+% dsp_options.userParameters.len_tr = 4096*2;
+% dsp_options.userParameters.pf_ncoeffs = 1;
+
+% dsp_options.userParameters.decoding_mode = [db_decoder.memoryless,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);
+
+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", "median", "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.median_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 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

projects/Diss/MPI_revisit/algorithms/PLOT_mpi_reduction_db_delay_regimes.m

@@ -9,7 +9,7 @@ clear; clc;
 
 studyName = "block_update_sweep";
 selectedBlockUpdate = 1;    % set [] to pool all block_update values
-selectedPamLevels = 6;      % set [] to use all PAM levels in the query result
+selectedPamLevels = 4;      % set [] to use all PAM levels in the query result
 
 algorithmSelection = table( ...
     ["plain_ffe"; ...
@@ -17,7 +17,7 @@ algorithmSelection = table( ...
      "a2_residual"; ...
      "a1_moving_average"; ...
      "dc_tracking"], ...
-    [true; true; true; true; true], ...
+    [true; false; false; false; true], ...
     'VariableNames', ["algorithm", "enabled"]);
 selectedAlgorithms = algorithmSelection.algorithm(algorithmSelection.enabled);
 
@@ -286,6 +286,81 @@ for regimeIdx = 1:numel(regimeNames)
     % mat2tikz_improved("C:/Users/Silas/Documents/6971e0b65b380ca6d71c837f/04_Experimental_Evaluation/tikz/mpi/ber_vs_sir_" + regimeName + ".tikz","cleanfigure",1);
 end
 
+%% 4) Plot BER spread over SIR by delay/coherence regime and algorithm
+
+spreadTable = berSpreadSummary(cleanData, groupVars);
+fprintf("Calculated BER spread for %d regime/PAM/algorithm/SIR groups.\n", ...
+    height(spreadTable));
+
+figure(); clf; hold on;
+
+for regimeIdx = 1:numel(regimeNames)
+    regimeName = regimeNames(regimeIdx);
+    lineStyle = pathRegimeLineStyle(regimeName);
+
+    for algIdx = 1:numel(selectedAlgorithms)
+        algorithmName = selectedAlgorithms(algIdx);
+        algColor = algorithmColor(algorithmName);
+        marker = algorithmMarker(algorithmName, algorithmMarkers);
+        displayName = algorithmDisplayName(algorithmName);
+
+        curveMask = spreadTable.path_regime == regimeName & ...
+            spreadTable.algorithm == algorithmName;
+
+        if ~any(curveMask)
+            continue
+        end
+
+        curveTable = sortrows(spreadTable(curveMask, :), ...
+            ["pam_level", "sir_exact"]);
+
+        for pamIdx = 1:numel(selectedPamLevels)
+            pamLevel = selectedPamLevels(pamIdx);
+            pamMask = curveTable.pam_level == pamLevel;
+            if ~any(pamMask)
+                continue
+            end
+
+            x = curveTable.sir_exact(pamMask).';
+            y = curveTable.std_log10_BER(pamMask).';
+            valid = isfinite(x) & isfinite(y);
+
+            if ~any(valid)
+                continue
+            end
+
+            if isscalar(selectedPamLevels)
+                legendText = sprintf("%s, %s", displayName, regimeName);
+            else
+                legendText = sprintf("%s, %s, PAM %.0f", ...
+                    displayName, regimeName, pamLevel);
+            end
+
+            plot(x(valid), y(valid), ...
+                "LineStyle", lineStyle, ...
+                "Marker", marker, ...
+                "MarkerSize", 3.5, ...
+                "LineWidth", 1, ...
+                "Color", algColor, ...
+                "MarkerFaceColor", "w", ...
+                "MarkerEdgeColor", algColor, ...
+                "DisplayName", legendText);
+        end
+    end
+end
+
+xlabel("SIR (dB)");
+ylabel("Std. dev. of log_{10}(BER)");
+title("BER spread over SIR by delay regime and algorithm");
+xlim([15, 45]);
+grid on;
+box on;
+legend("Location", "northeast", "Interpreter", "none");
+
+if exist("beautifyBERplot", "file")
+    beautifyBERplot("logscale", false, "setcolors", false, "setmarkers", false);
+end
+
 %% Local helpers
 
 function values = numericColumn(values)
@@ -487,6 +562,38 @@ function yFit = fitLogBer(x, y, xFit, fitOrder)
     yFit = 10 .^ polyval(coeff, xFit);
 end
 
+function spreadTable = berSpreadSummary(cleanData, groupVars)
+    summaryGroups = groupsummary(cleanData, groupVars);
+    sirExactTable = groupsummary(cleanData, groupVars, "median", "sir_exact");
+    summaryGroups = sortrows(summaryGroups, groupVars);
+    sirExactTable = sortrows(sirExactTable, groupVars);
+    summaryGroups.sir_exact = sirExactTable.median_sir_exact;
+
+    stdLogBer = NaN(height(summaryGroups), 1);
+    stdBer = NaN(height(summaryGroups), 1);
+    for groupIdx = 1:height(summaryGroups)
+        rowMask = true(height(cleanData), 1);
+        for varIdx = 1:numel(groupVars)
+            varName = groupVars(varIdx);
+            rowMask = rowMask & cleanData.(char(varName)) == ...
+                summaryGroups.(char(varName))(groupIdx);
+        end
+
+        berValues = cleanData.BER(rowMask);
+        berValues = berValues(isfinite(berValues) & berValues > 0);
+        if isempty(berValues)
+            continue
+        end
+
+        stdLogBer(groupIdx) = std(log10(berValues), 0, "omitnan");
+        stdBer(groupIdx) = std(berValues, 0, "omitnan");
+    end
+
+    spreadTable = summaryGroups;
+    spreadTable.std_log10_BER = stdLogBer;
+    spreadTable.std_BER = stdBer;
+end
+
 function label = algorithmDisplayName(algorithmName)
     algorithmName = string(algorithmName);
     switch algorithmName
@@ -532,3 +639,18 @@ function marker = algorithmMarker(algorithmName, algorithmMarkers)
     end
     marker = algorithmMarkers{mod(markerIdx - 1, numel(algorithmMarkers)) + 1};
 end
+
+function lineStyle = pathRegimeLineStyle(regimeName)
+    switch string(regimeName)
+        case "0-1 m"
+            lineStyle = "-";
+        case "10-100 m"
+            lineStyle = "--";
+        case "300 m"
+            lineStyle = ":";
+        case "1000 m"
+            lineStyle = "-.";
+        otherwise
+            lineStyle = "-";
+    end
+end

projects/Diss/MPI_revisit/parallelization_analysis/PLOT_mpi_reduction_db_parallelization_vs_sir.m

@@ -132,7 +132,7 @@ summaryTable = addBerStdBounds(summaryTable, cleanData, groupVars);
 
 blockUpdates = unique(cleanData.block_update(isfinite(cleanData.block_update))).';
 blockUpdates = sort(blockUpdates);
-helperBlockUpdates = unique(blockUpdates([1, end]), "stable");
+helperBlockUpdates = unique(blockUpdates([1, end-3]), "stable");
 
 requiredSirRows = table();
 for pamIdx = 1:numel(selectedPamLevels)
@@ -214,14 +214,14 @@ for pamIdx = 1:numel(selectedPamLevels)
                     sirValues, boundCenterBer, boundLowerBer, boundUpperBer, ...
                     boundMode, boundaryPolyfitOrderMax);
 
-                [hl, hp] = boundedline(xBand, centerBand, yBounds, ...
+                % [hl, hp] = boundedline(xBand, centerBand, yBounds, ...
-                    'alpha', 'transparency', 0.1, ...
+                %     'alpha', 'transparency', 0.1, ...
-                    'cmap', algColor, ...
+                %     'cmap', algColor, ...
-                    'nan', 'fill', ...
+                %     'nan', 'fill', ...
-                    'orientation', 'vert');
+                %     'orientation', 'vert');
-                set(hl, "LineStyle", "none", 'LineWidth', 1, "Marker", "none", ...
+                % set(hl, "LineStyle", "none", 'LineWidth', 1, "Marker", "none", ...
-                    "HandleVisibility", "off", "DisplayName", char(displayName));
+                %     "HandleVisibility", "off", "DisplayName", char(displayName));
-                set(hp, "LineStyle", "-", "HandleVisibility", "off", "Marker", "none");
+                % set(hp, "LineStyle", "-", "HandleVisibility", "off", "Marker", "none");
             end
 
             plot(sirValues(valid), meanBer(valid), ...
@@ -242,11 +242,11 @@ for pamIdx = 1:numel(selectedPamLevels)
                     xFit = linspace(min(sirValues(fitMask)), max(sirValues(fitMask)), 300);
                     yFit = 10 .^ polyval(fitCoeff, xFit);
 
-                    % plot(xFit, yFit, ...
+                    plot(xFit, yFit, ...
-                    %     "LineStyle", "--", ...
+                        "LineStyle", "--", ...
-                    %     "LineWidth", 1.1, ...
+                        "LineWidth", 1.1, ...
-                    %     "Color", algColor, ...
+                        "Color", algColor, ...
-                    %     "HandleVisibility", "off");
+                        "HandleVisibility", "off");
                 end
             end
         end

projects/Diss/MPI_revisit/simulation/PLOT_mpi_simulation_linewidth_vs_sir.m

@@ -0,0 +1,608 @@
+%% BER over SIR from saved MPI simulation warehouses grouped by linewidth
+% Loads a saved simulation warehouse and plots one BER-vs-SIR curve per
+% linewidth for each algorithm.
+
+% clear;
+% clc;
+
+%% Load data
+resultFile = "";
+if resultFile == ""
+    resultDir = fullfile(fileparts(mfilename("fullpath")), "results");
+    files = dir(fullfile(resultDir, "mpi_simulation_*.mat"));
+    if isempty(files)
+        error("PLOT_mpi_simulation_linewidth_vs_sir:NoResultFiles", ...
+            "No mpi_simulation_*.mat files found in %s.", resultDir);
+    end
+    [~, newestIdx] = max([files.datenum]);
+    resultFile = fullfile(files(newestIdx).folder, files(newestIdx).name);
+end
+
+loaded = load(resultFile, "wh", "simulation_config");
+wh = loaded.wh;
+fprintf("Loaded MPI simulation warehouse:\n%s\n", resultFile);
+wh.showInfo;
+
+%% Plot settings
+useBoundedLines = true;
+usePolyfit = true;
+polyfitOrderMax = 4;
+boundaryPolyfitOrderMax = 4;
+fecBerThreshold = 3.8e-3;
+maxBerForPlot = 0.1;
+crossingSirWindow = [15 40];
+boundMode = "fitStd";
+algorithmMarkers = {'o','square','diamond','^','v','>','<','pentagram'};
+
+%% Collect and clean
+cleanData = collectMpiSimulationBerRows(wh);
+cleanData = normalizeBerColumnName(cleanData);
+cleanData = cleanData(isfinite(cleanData.BER) & cleanData.BER > 0 & ...
+    cleanData.BER < maxBerForPlot, :);
+
+if isempty(cleanData)
+    warning("PLOT_mpi_simulation_linewidth_vs_sir:NoRows", ...
+        "No valid BER rows remain for plotting.");
+    return
+end
+
+if all(~isfinite(cleanData.laser_linewidth))
+    if isfield(loaded, "simulation_config") && isfield(loaded.simulation_config, "laser_linewidth")
+        cleanData.laser_linewidth(:) = loaded.simulation_config.laser_linewidth;
+    else
+        cleanData.laser_linewidth(:) = 0;
+    end
+end
+
+cleanData.clean_keep = true(height(cleanData), 1);
+groupId = findgroups(cleanData.storage_name, cleanData.block_update, ...
+    cleanData.laser_linewidth, cleanData.sir);
+for curGroup = unique(groupId(isfinite(groupId))).'
+    rowMask = groupId == curGroup;
+    berValues = cleanData.BER(rowMask);
+    if nnz(rowMask) > 3
+        cleanData.clean_keep(rowMask) = ~isoutlier(berValues);
+    end
+end
+cleanData = cleanData(cleanData.clean_keep, :);
+
+groupVars = ["storage_name", "algorithm", "block_update", "laser_linewidth", "sir"];
+summaryTable = groupsummary(cleanData, groupVars, {"mean", "min", "max"}, "BER");
+summaryTable = sortrows(summaryTable, groupVars);
+summaryTable.sir_exact = summaryTable.sir;
+summaryTable = addBerStdBounds(summaryTable, cleanData, groupVars);
+
+selectedBlockUpdates = unique(cleanData.block_update(isfinite(cleanData.block_update))).';
+selectedAlgorithms = unique(cleanData.storage_name, "stable").';
+selectedLinewidths = unique(cleanData.laser_linewidth(isfinite(cleanData.laser_linewidth))).';
+selectedLinewidths = sort(selectedLinewidths);
+
+requiredSirRows = table();
+for blockIdx = 1:numel(selectedBlockUpdates)
+    blockUpdate = selectedBlockUpdates(blockIdx);
+    for algIdx = 1:numel(selectedAlgorithms)
+        storageName = selectedAlgorithms(algIdx);
+        algorithmName = string(summaryTable.algorithm(find(summaryTable.storage_name == storageName, 1, "first")));
+        for linewidthIdx = 1:numel(selectedLinewidths)
+            laserLinewidth = selectedLinewidths(linewidthIdx);
+            curveMask = summaryTable.storage_name == storageName & ...
+                summaryTable.block_update == blockUpdate & ...
+                summaryTable.laser_linewidth == laserLinewidth;
+
+            sirValues = summaryTable.sir_exact(curveMask).';
+            meanBer = summaryTable.mean_BER(curveMask).';
+            [~, ~, requiredSir, fitOrder] = fitBerAtFec( ...
+                sirValues, meanBer, polyfitOrderMax, fecBerThreshold, crossingSirWindow);
+
+            newRow = table(storageName, algorithmName, blockUpdate, ...
+                laserLinewidth, requiredSir, fitOrder, nnz(isfinite(sirValues) & isfinite(meanBer)), ...
+                'VariableNames', {'storage_name', 'algorithm', 'block_update', ...
+                'laser_linewidth', 'required_sir', 'fit_order', 'n_points'});
+            requiredSirRows = [requiredSirRows; newRow]; %#ok<AGROW>
+        end
+    end
+end
+
+fprintf("Cleaned to %d simulation BER rows across %d storage/block/linewidth/SIR groups.\n", ...
+    height(cleanData), height(summaryTable));
+disp(groupcounts(cleanData, ["storage_name", "block_update", "laser_linewidth"]));
+disp(requiredSirRows);
+
+%% Plot BER vs SIR, linewidth as curve family
+for blockIdx = 1:numel(selectedBlockUpdates)
+    blockUpdate = selectedBlockUpdates(blockIdx);
+
+    figure();
+    clf;
+    tiledlayout(numel(selectedAlgorithms), 1, "TileSpacing", "compact");
+
+    for algIdx = 1:numel(selectedAlgorithms)
+        storageName = selectedAlgorithms(algIdx);
+        nexttile;
+        hold on;
+
+        for linewidthIdx = 1:numel(selectedLinewidths)
+            laserLinewidth = selectedLinewidths(linewidthIdx);
+            lineColor = linewidthColor(linewidthIdx, numel(selectedLinewidths));
+            marker = algorithmMarker(storageName, algorithmMarkers);
+            displayName = sprintf("%s, %s", ...
+                algorithmDisplayName(storageName), linewidthLabel(laserLinewidth));
+
+            rawMask = cleanData.storage_name == storageName & ...
+                cleanData.block_update == blockUpdate & ...
+                cleanData.laser_linewidth == laserLinewidth;
+            curveMask = summaryTable.storage_name == storageName & ...
+                summaryTable.block_update == blockUpdate & ...
+                summaryTable.laser_linewidth == laserLinewidth;
+
+            if ~any(curveMask)
+                continue
+            end
+
+            scatterRows = cleanData(rawMask, :);
+            scatter(scatterRows.sir, scatterRows.BER, ...
+                22, ...
+                "Marker", ".", ...
+                "MarkerEdgeColor", lineColor, ...
+                "MarkerFaceColor", lineColor, ...
+                "HandleVisibility", "off");
+
+            sirValues = summaryTable.sir_exact(curveMask).';
+            meanBer = summaryTable.mean_BER(curveMask).';
+            boundCenterBer = summaryTable.std_center_BER(curveMask).';
+            boundLowerBer = summaryTable.std_lower_BER(curveMask).';
+            boundUpperBer = summaryTable.std_upper_BER(curveMask).';
+            valid = isfinite(sirValues) & isfinite(meanBer) & meanBer > 0;
+
+            if useBoundedLines && exist("boundedline", "file") && any(valid)
+                [xBand, centerBand, yBounds] = berStdBounds( ...
+                    sirValues, boundCenterBer, boundLowerBer, boundUpperBer, ...
+                    boundMode, boundaryPolyfitOrderMax);
+
+                [hl, hp] = boundedline(xBand, centerBand, yBounds, ...
+                    'alpha', 'transparency', 0.08, ...
+                    'cmap', lineColor, ...
+                    'nan', 'fill', ...
+                    'orientation', 'vert');
+                set(hl, "LineStyle", "none", "LineWidth", 1, "Marker", "none", ...
+                    "HandleVisibility", "off", "DisplayName", displayName);
+                set(hp, "LineStyle", "-", "HandleVisibility", "off", "Marker", "none");
+            end
+
+            plot(sirValues(valid), meanBer(valid), ...
+                "LineStyle", "-", ...
+                "Marker", marker, ...
+                "MarkerSize", 3, ...
+                "LineWidth", 1, ...
+                "Color", lineColor, ...
+                "MarkerFaceColor", "w", ...
+                "MarkerEdgeColor", lineColor, ...
+                "DisplayName", displayName, ...
+                "HandleVisibility", "on");
+
+            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", lineColor, ...
+                        "HandleVisibility", "off");
+                end
+            end
+        end
+
+        yline(2.2e-4, "LineWidth", 1, "LineStyle", "--", "HandleVisibility", "off");
+        yline(fecBerThreshold, "LineWidth", 1, "LineStyle", "--", "HandleVisibility", "off");
+        yline(2e-2, "LineWidth", 1, "LineStyle", "--", "HandleVisibility", "off");
+
+        title(sprintf("%s, block update %.0f", algorithmDisplayName(storageName), blockUpdate));
+        xlabel("SIR (dB)");
+        ylabel("BER");
+        set(gca, "YScale", "log");
+        ylim([9e-5, maxBerForPlot]);
+        xlim(crossingSirWindow);
+        grid on;
+        box on;
+        legend("Location", "northeast", "Interpreter", "none");
+
+        if exist("beautifyBERplot", "file")
+            beautifyBERplot("logscale", true, "setcolors", false, ...
+                "setmarkers", false, "changemarkers", false);
+        end
+    end
+end
+
+%% Plot required SIR at FEC over linewidth
+for blockIdx = 1:numel(selectedBlockUpdates)
+    blockUpdate = selectedBlockUpdates(blockIdx);
+
+    figure();
+    clf;
+    hold on;
+
+    for algIdx = 1:numel(selectedAlgorithms)
+        storageName = selectedAlgorithms(algIdx);
+        algColor = algorithmColor(storageName);
+        marker = algorithmMarker(storageName, algorithmMarkers);
+        displayName = algorithmDisplayName(storageName);
+        rowMask = requiredSirRows.storage_name == storageName & ...
+            requiredSirRows.block_update == blockUpdate;
+
+        x = requiredSirRows.laser_linewidth(rowMask).';
+        y = requiredSirRows.required_sir(rowMask).';
+        [x, sortIdx] = sort(x);
+        y = y(sortIdx);
+        valid = isfinite(x) & isfinite(y);
+
+        if ~any(valid)
+            continue
+        end
+
+        plot(x(valid), y(valid), ...
+            "LineWidth", 1.4, ...
+            "LineStyle", "-", ...
+            "Marker", marker, ...
+            "MarkerSize", 5, ...
+            "Color", algColor, ...
+            "MarkerFaceColor", "w", ...
+            "MarkerEdgeColor", algColor, ...
+            "DisplayName", char(displayName));
+    end
+
+    set(gca, "XScale", "log");
+    xticks(selectedLinewidths);
+    xticklabels(arrayfun(@linewidthLabel, selectedLinewidths, "UniformOutput", false));
+    ylim(crossingSirWindow);
+    grid on;
+    box on;
+    xlabel("Laser linewidth");
+    ylabel(sprintf("Required SIR at BER = %.1e (dB)", fecBerThreshold));
+    title(sprintf("Required SIR over linewidth, block update %.0f", blockUpdate));
+    legend("Location", "best", "Interpreter", "none");
+
+    if exist("beautifyBERplot", "file")
+        beautifyBERplot("logscale", false, "setcolors", false, ...
+            "setmarkers", false, "changemarkers", false);
+    end
+end
+
+%% Local helpers
+function data = collectMpiSimulationBerRows(wh)
+storageNames = string(fieldnames(wh.sto));
+data = table();
+
+for storageIdx = 1:numel(storageNames)
+    storageName = storageNames(storageIdx);
+    storage = wh.sto.(char(storageName));
+
+    for linIdx = 1:numel(storage)
+        package = storage{linIdx};
+        ber = extractPackageBer(package);
+        if ~isfinite(ber)
+            continue
+        end
+
+        [physValues, physNames] = wh.getPhysIndicesByLinIndex(linIdx);
+        phys = struct();
+        for physIdx = 1:numel(physNames)
+            phys.(char(physNames{physIdx})) = physValues{physIdx};
+        end
+
+        algorithm = extractPackageAlgorithm(package, storageName);
+        newRow = table( ...
+            storageName, ...
+            algorithm, ...
+            readPhysValue(phys, "sir", NaN), ...
+            readPhysValue(phys, "block_update", NaN), ...
+            readPhysValue(phys, "random_key", NaN), ...
+            readPhysValue(phys, "laser_linewidth", NaN), ...
+            ber, ...
+            'VariableNames', {'storage_name', 'algorithm', 'sir', ...
+            'block_update', 'random_key', 'laser_linewidth', 'BER'});
+        data = [data; newRow]; %#ok<AGROW>
+    end
+end
+end
+
+function data = normalizeBerColumnName(data)
+variableNames = string(data.Properties.VariableNames);
+if ismember("BER", variableNames)
+    return
+end
+
+if ismember("ber", variableNames)
+    data.Properties.VariableNames(variableNames == "ber") = {'BER'};
+    return
+end
+
+error("PLOT_mpi_simulation_linewidth_vs_sir:MissingBerColumn", ...
+    "Could not find a BER or ber column in the simulation BER table.");
+end
+
+function value = readPhysValue(phys, name, defaultValue)
+if isfield(phys, name)
+    value = phys.(name);
+else
+    value = defaultValue;
+end
+end
+
+function ber = extractPackageBer(package)
+ber = NaN;
+if isempty(package)
+    return
+end
+
+if iscell(package)
+    package = package{1};
+end
+
+if isstruct(package) && isfield(package, "metrics")
+    metrics = package.metrics;
+    if isstruct(metrics) && isfield(metrics, "BER")
+        ber = metrics.BER;
+    elseif isobject(metrics) && isprop(metrics, "BER")
+        ber = metrics.BER;
+    end
+end
+end
+
+function algorithm = extractPackageAlgorithm(package, fallbackName)
+algorithm = fallbackName;
+if isempty(package)
+    return
+end
+
+if iscell(package)
+    package = package{1};
+end
+
+if isstruct(package) && isfield(package, "mpi_reduction_config") && ...
+        isfield(package.mpi_reduction_config, "algorithm")
+    algorithm = string(package.mpi_reduction_config.algorithm);
+end
+end
+
+function summaryTable = addBerStdBounds(summaryTable, cleanData, groupVars)
+nGroups = height(summaryTable);
+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
+
+    [stdCenter(groupIdx), stdLower(groupIdx), stdUpper(groupIdx)] = ...
+        logBerMeanStdInterval(cleanData.BER(rowMask));
+end
+
+summaryTable.std_center_BER = stdCenter;
+summaryTable.std_lower_BER = stdLower;
+summaryTable.std_upper_BER = stdUpper;
+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 [xBand, centerBand, yBounds] = berStdBounds(sirValues, centerBer, lowerBer, upperBer, boundMode, maxOrder)
+if 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 [xFit, yFit, requiredSir, fitOrder] = fitBerAtFec( ...
+    sirValues, meanBer, polyfitOrderMax, fecBerThreshold, crossingSirWindow)
+
+xFit = NaN;
+yFit = NaN;
+requiredSir = NaN;
+fitOrder = NaN;
+
+valid = isfinite(sirValues) & isfinite(meanBer) & meanBer > 0 & ...
+    sirValues >= crossingSirWindow(1) & sirValues <= crossingSirWindow(2);
+if nnz(valid) < 2
+    return
+end
+
+sirValues = sirValues(valid);
+meanBer = meanBer(valid);
+[sirValues, sortIdx] = sort(sirValues);
+meanBer = meanBer(sortIdx);
+
+fitOrder = min(polyfitOrderMax, nnz(valid) - 1);
+fitCoeff = polyfit(sirValues, log10(meanBer), fitOrder);
+xFit = linspace(max(min(sirValues), crossingSirWindow(1)), ...
+    min(max(sirValues), crossingSirWindow(2)), 300);
+yFit = 10 .^ polyval(fitCoeff, xFit);
+
+thresholdMask = isfinite(yFit) & yFit <= fecBerThreshold;
+if ~any(thresholdMask)
+    return
+end
+
+firstThresholdIdx = find(thresholdMask, 1, "first");
+if firstThresholdIdx == 1
+    requiredSir = xFit(firstThresholdIdx);
+    return
+end
+
+xPair = xFit(firstThresholdIdx - 1:firstThresholdIdx);
+yPair = log10(yFit(firstThresholdIdx - 1:firstThresholdIdx));
+if all(isfinite(yPair)) && diff(yPair) ~= 0
+    requiredSir = interp1(yPair, xPair, log10(fecBerThreshold), ...
+        "linear", "extrap");
+else
+    requiredSir = xFit(firstThresholdIdx);
+end
+
+if requiredSir < crossingSirWindow(1) || requiredSir > crossingSirWindow(2)
+    requiredSir = NaN;
+end
+end
+
+function label = algorithmDisplayName(algorithmName)
+algorithmName = string(algorithmName);
+switch algorithmName
+    case {"plain_ffe", "conventional_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", "conventional_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 color = linewidthColor(linewidthIdx, nLinewidths)
+if nLinewidths <= 1
+    color = [0.3467 0.5360 0.6907];
+    return
+end
+
+if exist("cbrewer2", "file")
+    colors = cbrewer2("Set1", max(nLinewidths, 3));
+else
+    colors = lines(nLinewidths);
+end
+
+color = colors(linewidthIdx, :);
+end
+
+function label = linewidthLabel(laserLinewidth)
+if abs(laserLinewidth) >= 1e6
+    label = sprintf("%.3g MHz", laserLinewidth * 1e-6);
+elseif abs(laserLinewidth) >= 1e3
+    label = sprintf("%.3g kHz", laserLinewidth * 1e-3);
+else
+    label = sprintf("%.3g Hz", laserLinewidth);
+end
+end
+
+function marker = algorithmMarker(algorithmName, algorithmMarkers)
+algorithmOrder = ["conventional_ffe", "dc_tracking", "a2_tracked_levels", ...
+    "a2_residual", "a1_moving_average"];
+markerIdx = find(algorithmOrder == string(algorithmName), 1);
+if isempty(markerIdx)
+    markerIdx = 1;
+end
+marker = algorithmMarkers{mod(markerIdx - 1, numel(algorithmMarkers)) + 1};
+end

projects/Diss/MPI_revisit/simulation/mpi_simulation_worker.m

@@ -0,0 +1,214 @@
+function output = mpi_simulation_worker(userParameters, simulation_config)
+%MPI_SIMULATION_WORKER Build one decorrelated MPI point and run the DSP recipe.
+
+arguments
+    userParameters struct
+    simulation_config struct
+end
+
+config = applyDefaults(simulation_config);
+
+[Scpe_sig_raw, Symbols, Tx_bits, dataTable] = buildMpiScopeSignal(userParameters, config);
+
+output = config.recipe(Scpe_sig_raw, Symbols, Tx_bits, ...
+    "fsym", config.fsym, ...
+    "M", config.M, ...
+    "duob_mode", config.duob_mode, ...
+    "dataTable", dataTable, ...
+    "userParameters", userParameters, ...
+    "debug_plots", config.debug_plots);
+end
+
+function config = applyDefaults(config)
+defaults = struct( ...
+    "M", 4, ...
+    "fsym", 112e9, ...
+    "mpi_path_meter", 1000, ...
+    "laser_linewidth", 150e3, ...
+    "recipe", @mpi_recipe_dev, ...
+    "debug_plots", false, ...
+    "fdac", 256e9, ...
+    "fadc", 256e9, ...
+    "kover", 16, ...
+    "random_key", 1, ...
+    "rcalpha", 0.05, ...
+    "duob_mode", db_mode.no_db, ...
+    "vbias_rel", 0.5, ...
+    "u_pi", 3, ...
+    "laser_wavelength", 1293, ...
+    "link_length_km", 1, ...
+    "rop", -9, ...
+    "rx_bwl", 80e9, ...
+    "scope_bwl", 110e9, ...
+    "alpha", 0);
+
+names = fieldnames(defaults);
+for nameIdx = 1:numel(names)
+    name = names{nameIdx};
+    if ~isfield(config, name) || isempty(config.(name))
+        config.(name) = defaults.(name);
+    end
+end
+end
+
+function [Scpe_sig, Symbols, Tx_bits, dataTable] = buildMpiScopeSignal(userParameters, config)
+sir = readUserParameter(userParameters, "sir", 30);
+randomKey = readUserParameter(userParameters, "random_key", config.random_key);
+laserLinewidth = readUserParameter(userParameters, "laser_linewidth", config.laser_linewidth);
+config.laser_linewidth = laserLinewidth;
+
+Pform = Pulseformer( ...
+    "fsym", config.fsym, ...
+    "fdac", 4 * config.fsym, ...
+    "pulse", "rrc", ...
+    "pulselength", 16, ...
+    "alpha", config.rcalpha);
+
+vbias = -config.vbias_rel * config.u_pi;
+mainSource = buildOpticalSource(Pform, config, randomKey, randomKey + 1, vbias);
+interferenceSource = buildOpticalSource(Pform, config, randomKey + 100000, ...
+    randomKey + 100001, vbias);
+
+Symbols = mainSource.Symbols;
+Tx_bits = mainSource.Tx_bits;
+Opt_sig = applyDecorrelatedMpi(mainSource.Opt_sig, interferenceSource.Opt_sig, ...
+    sir, config.link_length_km);
+
+Rx_sig = Amplifier( ...
+    "amp_mode", "ideal_no_noise", ...
+    "gain_mode", "output_power", ...
+    "amplification_db", config.rop).process(Opt_sig);
+
+Rx_sig = Photodiode( ...
+    "fsimu", config.fdac * config.kover, ...
+    "dark_current", 2e-08, ...
+    "responsivity", 1, ...
+    "temperature", 20, ...
+    "nep", 1.8e-11).process(Rx_sig);
+Rx_sig.signal = real(Rx_sig.signal);
+
+Rx_sig = Filter( ...
+    "filtdegree", 4, ...
+    "f_cutoff", config.rx_bwl, ...
+    "fs", config.fdac * config.kover, ...
+    "filterType", filtertypes.butterworth, ...
+    "active", true).process(Rx_sig);
+Rx_sig.signal = real(Rx_sig.signal);
+
+scopeFilter = Filter( ...
+    "filtdegree", 4, ...
+    "f_cutoff", config.scope_bwl, ...
+    "fs", config.fadc, ...
+    "filterType", filtertypes.butterworth, ...
+    "active", true);
+
+Scpe_sig = Scope( ...
+    "fsimu", config.fdac * config.kover, ...
+    "fadc", config.fadc, ...
+    "delay", 0, ...
+    "fixed_delay", 0, ...
+    "filtertype", filtertypes.butterworth, ...
+    "samplingdelay", 0, ...
+    "rand_samplingdelay", 0, ...
+    "freq_offset", 0, ...
+    "samp_jitter", 0, ...
+    "adcresolution", 8, ...
+    "quantbuffer", 0.1, ...
+    "block_dc", 1, ...
+    "lpf_active", 1, ...
+    "H_lpf", scopeFilter).process(Rx_sig);
+Scpe_sig.signal = real(Scpe_sig.signal);
+
+HighpassFilter = Filter( ...
+    "filtdegree", 6, ...
+    "f_cutoff", 1e6, ...
+    "fs", config.fadc, ...
+    "filterType", filtertypes.butterworth, ...
+    "active", true, "lowpass",0);
+
+Scpe_sig = HighpassFilter.process(Scpe_sig);
+
+dataTable = table( ...
+    sir, ...
+    config.fsym, ...
+    config.M, ...
+    config.mpi_path_meter, ...
+    laserLinewidth, ...
+    randomKey, ...
+    'VariableNames', {'sir', 'symbolrate', 'pam_level', ...
+    'interference_path_length', 'laser_linewidth', 'random_key'});
+end
+
+function source = buildOpticalSource(Pform, config, sourceRandomKey, laserRandomKey, vbias)
+[Digi_sig, Symbols, Tx_bits] = PAMsource( ...
+    "fsym", config.fsym, ...
+    "M", config.M, ...
+    "order", 18, ...
+    "useprbs", 0, ...
+    "fs_out", config.fdac, ...
+    "applyclipping", 0, ...
+    "clipfactor", 1.5, ...
+    "applypulseform", 1, ...
+    "pulseformer", Pform, ...
+    "randkey", sourceRandomKey, ...
+    "duobinary_mode", config.duob_mode, ...
+    "mrds_code", 0, ...
+    "mrds_blocklength", 512).process();
+
+El_sig = M8199A("kover", config.kover).process(Digi_sig);
+
+
+El_sig = Filter("f_cutoff",65e9,"filterType","butterworth","filtdegree",4,"fs",El_sig.fs).process(El_sig);
+
+El_sig = El_sig.normalize("mode", "oneone");
+
+scaling = 0.6 * (config.u_pi / 2 - abs(vbias - config.u_pi / 2));
+El_sig = El_sig .* scaling;
+
+Opt_sig = EML( ...
+    "mode", eml_mode.im_cosinus, ...
+    "power", 3, ...
+    "fsimu", El_sig.fs, ...
+    "lambda", config.laser_wavelength, ...
+    "bias", vbias, ...
+    "u_pi", config.u_pi, ...
+    "linewidth", config.laser_linewidth, ...
+    "randomkey", laserRandomKey, ...
+    "alpha", config.alpha).process(El_sig);
+
+source = struct( ...
+    "Opt_sig", Opt_sig, ...
+    "Symbols", Symbols, ...
+    "Tx_bits", Tx_bits);
+end
+
+function Opt_sig = applyDecorrelatedMpi(main_sig, interference_sig, sir, linkLengthKm)
+interference_sig = Amplifier( ...
+    "amp_mode", "ideal_no_noise", ...
+    "gain_mode", "output_power", ...
+    "amplification_db", main_sig.power - sir).process(interference_sig);
+
+if numel(main_sig.signal) ~= numel(interference_sig.signal)
+    minLength = min(numel(main_sig.signal), numel(interference_sig.signal));
+    main_sig.signal = main_sig.signal(1:minLength);
+    interference_sig.signal = interference_sig.signal(1:minLength);
+end
+combined_sig = main_sig + interference_sig;
+
+Opt_sig = Fiber( ...
+    "fsimu", combined_sig.fs, ...
+    "fiber_length", linkLengthKm, ...
+    "alpha", 0.3, ...
+    "D", 0, ...
+    "lambda0", 1310, ...
+    "gamma", 0, ...
+    "Dslope", 0.07).process(combined_sig);
+end
+
+function value = readUserParameter(userParameters, name, defaultValue)
+if isfield(userParameters, name)
+    value = userParameters.(name);
+else
+    value = defaultValue;
+end
+end

projects/Diss/MPI_revisit/simulation/run_mpi_simulation_recipe.m

@@ -0,0 +1,476 @@
+%RUN_MPI_SIMULATION_RECIPE Queue MPI simulation points and store DSP packages.
+
+clear;
+% clc;
+
+%% Fixed experiment-like simulation setup
+simulation_config = struct();
+simulation_config.M = 4;
+simulation_config.fsym = 112e9;
+simulation_config.mpi_path_meter = 1000;
+simulation_config.laser_linewidth = 150e3;%[150e3 150e3 250e3 500e3 750e3 1e6 5e6 10e6 20e6 50e6];
+simulation_config.laser_linewidths = simulation_config.laser_linewidth;
+simulation_config.recipe = @mpi_recipe_dev;
+simulation_config.debug_plots = true;
+simulation_config.waitbar = true;
+simulation_config.processing_mode = processingMode.parallel;
+simulation_config.num_workers = 0;
+simulation_config.random_keys = 1:10;
+
+simulation_config.fdac = 256e9;
+simulation_config.fadc = 256e9;
+simulation_config.kover = 16;
+simulation_config.random_key = 1;
+simulation_config.rcalpha = 0.05;
+simulation_config.duob_mode = db_mode.no_db;
+simulation_config.vbias_rel = 0.5;
+simulation_config.u_pi = 3;
+simulation_config.laser_wavelength = 1310;
+simulation_config.link_length_km = 0.5;
+simulation_config.rop = -7.5;
+simulation_config.rx_bwl = 70e9;
+simulation_config.scope_bwl = 70e9;
+simulation_config.alpha = 0;
+
+%% Warehouse sweep setup
+sweep_params = struct();
+sweep_params.sir = 15:3:45;
+sweep_params.block_update = 1;
+sweep_params.random_key = simulation_config.random_keys;
+sweep_params.laser_linewidth = simulation_config.laser_linewidths;
+
+wh = DataStorage(sweep_params);
+
+fprintf("Requested %d MPI simulation job(s).\n", wh.getLastLinIndice());
+
+%% Run queued simulations
+[results, wh] = submitMpiSimulationJobs(wh, simulation_config, ...
+    "mode", simulation_config.processing_mode, ...
+    "waitbar", simulation_config.waitbar, ...
+    "numWorkers", simulation_config.num_workers);
+
+%% Save result artifact
+result_dir = fullfile(fileparts(mfilename("fullpath")), "results");
+if ~exist(result_dir, "dir")
+    mkdir(result_dir);
+end
+
+timestamp = string(datetime("now", "Format", "yyyyMMdd_HHmmss"));
+result_file = fullfile(result_dir, "mpi_simulation_" + timestamp + ".mat");
+save(result_file, "wh", "simulation_config", "results");
+fprintf("Saved MPI simulation warehouse to:\n%s\n", result_file);
+
+%% BER inspection
+plotMpiSimulationBer(wh);
+
+function plotMpiSimulationBer(wh)
+storageNames = string(fieldnames(wh.sto));
+if isempty(storageNames)
+    warning("run_mpi_simulation_recipe:NoStorage", ...
+        "The warehouse does not contain any stored DSP packages.");
+    return
+end
+
+useBoundedLines = true;
+usePolyfit = true;
+polyfitOrderMax = 4;
+boundaryPolyfitOrderMax = 4;
+fecBerThreshold = 3.8e-3;
+maxBerForPlot = 0.1;
+boundMode = "fitStd";
+algorithmMarkers = {'o','square','diamond','^','v','>','<','pentagram'};
+
+cleanData = collectMpiSimulationBerRows(wh);
+cleanData = normalizeBerColumnName(cleanData);
+cleanData = cleanData(isfinite(cleanData.BER) & cleanData.BER > 0 & ...
+    cleanData.BER < maxBerForPlot, :);
+
+if isempty(cleanData)
+    warning("run_mpi_simulation_recipe:NoBerRows", ...
+        "No valid BER rows remain for plotting.");
+    return
+end
+
+cleanData.clean_keep = true(height(cleanData), 1);
+groupId = findgroups(cleanData.storage_name, cleanData.block_update, cleanData.sir);
+for curGroup = unique(groupId(isfinite(groupId))).'
+    rowMask = groupId == curGroup;
+    berValues = cleanData.BER(rowMask);
+    if nnz(rowMask) > 3
+        cleanData.clean_keep(rowMask) = ~isoutlier(berValues);
+    end
+end
+cleanData = cleanData(cleanData.clean_keep, :);
+
+groupVars = ["storage_name", "algorithm", "block_update", "sir"];
+summaryTable = groupsummary(cleanData, groupVars, {"mean", "min", "max"}, "BER");
+summaryTable = sortrows(summaryTable, groupVars);
+summaryTable.sir_exact = summaryTable.sir;
+summaryTable = addBerStdBounds(summaryTable, cleanData, groupVars);
+
+selectedBlockUpdates = unique(cleanData.block_update(isfinite(cleanData.block_update))).';
+selectedAlgorithms = unique(cleanData.storage_name, "stable").';
+
+fprintf("Cleaned to %d simulation BER rows across %d storage/block/SIR groups.\n", ...
+    height(cleanData), height(summaryTable));
+disp(groupcounts(cleanData, ["storage_name", "block_update"]));
+
+figure();
+clf;
+tiledlayout(numel(selectedBlockUpdates), 1, "TileSpacing", "compact");
+
+for blockIdx = 1:numel(selectedBlockUpdates)
+    blockUpdate = selectedBlockUpdates(blockIdx);
+    nexttile; hold on;
+
+    for algIdx = 1:numel(selectedAlgorithms)
+        storageName = selectedAlgorithms(algIdx);
+        algColor = algorithmColor(storageName);
+        marker = algorithmMarker(storageName, algorithmMarkers);
+        displayName = algorithmDisplayName(storageName);
+
+        rawMask = cleanData.storage_name == storageName & ...
+            cleanData.block_update == blockUpdate;
+        curveMask = summaryTable.storage_name == storageName & ...
+            summaryTable.block_update == blockUpdate;
+
+        if ~any(curveMask)
+            continue
+        end
+
+        scatterRows = cleanData(rawMask, :);
+        scatter(scatterRows.sir, scatterRows.BER, ...
+            26, ...
+            "Marker", ".", ...
+            "MarkerEdgeColor", algColor, ...
+            "MarkerFaceColor", algColor, ...
+            "HandleVisibility", "off");
+
+        sirValues = summaryTable.sir_exact(curveMask).';
+        meanBer = summaryTable.mean_BER(curveMask).';
+        boundCenterBer = summaryTable.std_center_BER(curveMask).';
+        boundLowerBer = summaryTable.std_lower_BER(curveMask).';
+        boundUpperBer = summaryTable.std_upper_BER(curveMask).';
+        valid = isfinite(sirValues) & isfinite(meanBer) & meanBer > 0;
+
+        if useBoundedLines && exist("boundedline", "file") && any(valid)
+            [xBand, centerBand, yBounds] = berStdBounds( ...
+                sirValues, boundCenterBer, boundLowerBer, boundUpperBer, ...
+                boundMode, boundaryPolyfitOrderMax);
+
+            [hl, hp] = boundedline(xBand, centerBand, yBounds, ...
+                'alpha', 'transparency', 0.1, ...
+                'cmap', algColor, ...
+                'nan', 'fill', ...
+                'orientation', 'vert');
+            set(hl, "LineStyle", "none", "LineWidth", 1, "Marker", "none", ...
+                "HandleVisibility", "off", "DisplayName", char(displayName));
+            set(hp, "LineStyle", "-", "HandleVisibility", "off", "Marker", "none");
+        end
+
+        plot(sirValues(valid), meanBer(valid), ...
+            "LineStyle", "-", ...
+            "Marker", marker, ...
+            "MarkerSize", 3, ...
+            "LineWidth", 1, ...
+            "Color", algColor, ...
+            "MarkerFaceColor", "w", ...
+            "MarkerEdgeColor", algColor, ...
+            "DisplayName", char(displayName), ...
+            "HandleVisibility", "on");
+
+        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
+
+    yline(2.2e-4, "LineWidth", 1, "LineStyle", "--", "HandleVisibility", "off");
+    yline(fecBerThreshold, "LineWidth", 1, "LineStyle", "--", "HandleVisibility", "off");
+    yline(2e-2, "LineWidth", 1, "LineStyle", "--", "HandleVisibility", "off");
+
+    title(sprintf("Simulated MPI, block update %.0f", blockUpdate));
+    xlabel("SIR (dB)");
+    ylabel("BER");
+    set(gca, "YScale", "log");
+    ylim([9e-5, maxBerForPlot]);
+    xlim([min(cleanData.sir) - 1, max(cleanData.sir) + 1]);
+    grid on;
+    box on;
+    legend("Location", "northeast", "Interpreter", "none");
+
+    if exist("beautifyBERplot", "file")
+        beautifyBERplot("logscale", true, "setcolors", false, "setmarkers", false);
+    end
+end
+end
+
+function data = collectMpiSimulationBerRows(wh)
+storageNames = string(fieldnames(wh.sto));
+data = table();
+
+for storageIdx = 1:numel(storageNames)
+    storageName = storageNames(storageIdx);
+    storage = wh.sto.(char(storageName));
+
+    for linIdx = 1:numel(storage)
+        package = storage{linIdx};
+        ber = extractPackageBer(package);
+        if ~isfinite(ber)
+            continue
+        end
+
+        [physValues, physNames] = wh.getPhysIndicesByLinIndex(linIdx);
+        phys = struct();
+        for physIdx = 1:numel(physNames)
+            phys.(char(physNames{physIdx})) = physValues{physIdx};
+        end
+
+        algorithm = extractPackageAlgorithm(package, storageName);
+        newRow = table( ...
+            storageName, ...
+            algorithm, ...
+            readPhysValue(phys, "sir", NaN), ...
+            readPhysValue(phys, "block_update", NaN), ...
+            readPhysValue(phys, "random_key", NaN), ...
+            readPhysValue(phys, "laser_linewidth", NaN), ...
+            ber, ...
+            'VariableNames', {'storage_name', 'algorithm', 'sir', ...
+            'block_update', 'random_key', 'laser_linewidth', 'BER'});
+        data = [data; newRow]; %#ok<AGROW>
+    end
+end
+end
+
+function data = normalizeBerColumnName(data)
+variableNames = string(data.Properties.VariableNames);
+if ismember("BER", variableNames)
+    return
+end
+
+if ismember("ber", variableNames)
+    data.Properties.VariableNames(variableNames == "ber") = {'BER'};
+    return
+end
+
+error("run_mpi_simulation_recipe:MissingBerColumn", ...
+    "Could not find a BER or ber column in the simulation BER table.");
+end
+
+function value = readPhysValue(phys, name, defaultValue)
+if isfield(phys, name)
+    value = phys.(name);
+else
+    value = defaultValue;
+end
+end
+
+function ber = extractPackageBer(package)
+ber = NaN;
+if isempty(package)
+    return
+end
+
+if iscell(package)
+    package = package{1};
+end
+
+if isstruct(package) && isfield(package, "metrics")
+    metrics = package.metrics;
+    if isstruct(metrics) && isfield(metrics, "BER")
+        ber = metrics.BER;
+    elseif isobject(metrics) && isprop(metrics, "BER")
+        ber = metrics.BER;
+    end
+end
+end
+
+function algorithm = extractPackageAlgorithm(package, fallbackName)
+algorithm = fallbackName;
+if isempty(package)
+    return
+end
+
+if iscell(package)
+    package = package{1};
+end
+
+if isstruct(package) && isfield(package, "mpi_reduction_config") && ...
+        isfield(package.mpi_reduction_config, "algorithm")
+    algorithm = string(package.mpi_reduction_config.algorithm);
+end
+end
+
+function summaryTable = addBerStdBounds(summaryTable, cleanData, groupVars)
+nGroups = height(summaryTable);
+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
+
+    [stdCenter(groupIdx), stdLower(groupIdx), stdUpper(groupIdx)] = ...
+        logBerMeanStdInterval(cleanData.BER(rowMask));
+end
+
+summaryTable.std_center_BER = stdCenter;
+summaryTable.std_lower_BER = stdLower;
+summaryTable.std_upper_BER = stdUpper;
+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 [xBand, centerBand, yBounds] = berStdBounds(sirValues, centerBer, lowerBer, upperBer, boundMode, maxOrder)
+if 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", "conventional_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", "conventional_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 = algorithmMarker(algorithmName, algorithmMarkers)
+algorithmOrder = ["conventional_ffe", "dc_tracking", "a2_tracked_levels", ...
+    "a2_residual", "a1_moving_average"];
+markerIdx = find(algorithmOrder == string(algorithmName), 1);
+if isempty(markerIdx)
+    markerIdx = 1;
+end
+marker = algorithmMarkers{mod(markerIdx - 1, numel(algorithmMarkers)) + 1};
+end

projects/Diss/MPI_revisit/simulation/submitMpiSimulationJobs.m

@@ -0,0 +1,87 @@
+function [results, wh] = submitMpiSimulationJobs(wh, simulation_config, options)
+%SUBMITMPISIMULATIONJOBS Execute MPI simulation warehouse points via runBatch.
+
+arguments
+    wh DataStorage
+    simulation_config struct
+    options.mode = processingMode.serial
+    options.waitbar (1,1) logical = true
+    options.numWorkers (1,1) double {mustBeNonnegative, mustBeInteger} = 0
+    options.idleTimeout (1,1) double {mustBePositive} = 300
+    options.cancelExistingQueue (1,1) logical = true
+end
+
+nJobs = wh.getLastLinIndice();
+jobs = repmat(struct("args", {{}}, "label", "", "meta", struct()), 1, nJobs);
+
+for linIdx = 1:nJobs
+    userParameters = buildUserParameters(wh, linIdx);
+    jobs(linIdx).args = {userParameters, simulation_config};
+    jobs(linIdx).label = buildJobLabel(userParameters, linIdx);
+    jobs(linIdx).meta.lin_idx = linIdx;
+    jobs(linIdx).meta.userParameters = userParameters;
+end
+
+results = runBatch(@mpi_simulation_worker, jobs, ...
+    "mode", options.mode, ...
+    "waitbar", options.waitbar, ...
+    "waitbarMessage", "Processing MPI simulations...", ...
+    "numWorkers", options.numWorkers, ...
+    "idleTimeout", options.idleTimeout, ...
+    "cancelExistingQueue", options.cancelExistingQueue, ...
+    "resultHandler", @storeResult, ...
+    "errorHandler", @handleError);
+
+    function userParameters = buildUserParameters(storageWh, linIdx)
+        userParameters = struct();
+        if isempty(storageWh.getDimension())
+            return
+        end
+
+        [values, names] = storageWh.getPhysIndicesByLinIndex(linIdx);
+        for paramIdx = 1:numel(names)
+            userParameters.(char(names{paramIdx})) = values{paramIdx};
+        end
+    end
+
+    function label = buildJobLabel(userParameters, linIdx)
+        label = sprintf("MPI sim job %d", linIdx);
+        if isfield(userParameters, "sir")
+            label = sprintf("%s, SIR %g dB", label, userParameters.sir);
+        end
+        if isfield(userParameters, "block_update")
+            label = sprintf("%s, block %g", label, userParameters.block_update);
+        end
+    end
+
+    function storeResult(val, job, ~)
+        if isempty(val) || ~isstruct(val)
+            return
+        end
+
+        storageNames = fieldnames(val);
+        for storageIdx = 1:numel(storageNames)
+            storageName = storageNames{storageIdx};
+            if isempty(val.(storageName))
+                continue
+            end
+
+            ensureStorage(storageName);
+            wh.addValueToStorageByLinIdx(val.(storageName), storageName, job.meta.lin_idx);
+        end
+    end
+
+    function ensureStorage(storageName)
+        if ~isfield(wh.sto, storageName)
+            wh.addStorage(storageName);
+        end
+    end
+
+    function handleError(ME, job, ~)
+        fprintf("[%s] ERROR [%s]: %s\n", job.label, ME.identifier, ME.message);
+        for st = ME.stack'
+            fprintf("    %s:%d   (%s)\n", st.file, st.line, st.name);
+        end
+        fprintf("Full report:\n%s\n", getReport(ME, "extended"));
+    end
+end

projects/IMDD_base_system/minimal_example.m

@@ -133,14 +133,19 @@ Scpe_sig.signal = Scpe_sig.signal(1:2*length(Symbols));
 if 1
     % -------------------- FFE --------------------
     ffe_order   = [150, 0, 0];
-    eq_      = EQ("Ne",ffe_order,"Nb",[2,0,0], ...
+    eq_      = EQ("Ne",ffe_order,"Nb",[0,0,0], ...
         "training_length",len_tr,"training_loops",5,"dd_loops",5, ...
         "K",2,"DCmu",mu_dc,"DDmu",[mu_ffe mu_dfe],"DFEmu",0.005, ...
         "FFEmu",0,"plotfinal",0,"ideal_dfe",0);
     
     % mu_tr_rls = 9.805e-01; mu_dd_rls = 0.999989348903919;
     eq_ = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",1.933e-04,"order",ffe_order(1),"sps",2,"decide",0,"optmize_mus",1,"dd_mode",1,"adaption_technique","nlms");
+    eq_ = FFE_plain("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",1.933e-04,"order",ffe_order(1),"sps",2,"decide",0,"optmize_mus",1,"dd_mode",1,"adaption_technique","nlms");
     
+    eq_ = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2, ...
+        "mu_dd",[0.0004 0.0005 0.0006],"mu_tr",[0.0001 0.0008 0.001], ...
+        "order",[150,5,5],"sps",2,"decide",1, ...
+        "optmize_mus",1,"mu_optimization_len",2^15);
     
     % eq_ = FFE("epochs_tr",4,"epochs_dd",5,"len_tr",4096,"mu_dd",0.01,"mu_tr",0.01,"order",50,"sps",2,"decide",0, "adaption",adaption_method.nlms,"dd_mode",1);
     % eq_ = FFE_DFE("epochs_tr",5,"epochs_dd",5,"len_tr",512,"ffe_mu_dd",1e-4,"dfe_mu_dd",5e-4,"ffe_mu_tr",0,"dfe_mu_tr",0,"ffe_order",99,"dfe_order",99,"sps",2,"decide",0);
@@ -149,10 +154,10 @@ if 1
         "precode_mode",duob_mode,'showAnalysis',1,"postFFE",[], ...
         "eth_style_symbol_mapping",0);
     
-    output.ffe_results.metrics.print("description",'DFE');
+    output.ffe_results.metrics.print("description",'VNLE');
 end
 %%
-
+if 0
     % -------------------- VNLE + MLSE --------------------
     pf_ncoeffs  = 1;
     ffe_order3  = [50, 5, 5];
@@ -168,9 +173,9 @@ mlse_ = MLSE("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels);
         "precode_mode", duob_mode, 'showAnalysis', 1, "postFFE", [], "eth_style_symbol_mapping", 0);
     
     output.mlse_results.metrics.print("description",'MLSE');
-
+end
 %%
-
+if 0
     % -------------------- DB target --------------------
     mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",M,'trellis_states',PAMmapper(M,0).levels);
     ffe_order = [50, 5, 5];
@@ -180,6 +185,6 @@ output.dbt_results = duobinary_target(eq_,mlse_db_, M, Scpe_sig, Symbols, Tx_bit
         "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", []);
     
     output.dbt_results.metrics.print("description",'Duobinary');
-
+end