function WDM_model_10km_queue(options)
%WDM_model_10km_queue_memopt
% Queue/pipeline version with:
% (a) reduced variable lifetime / fewer unnecessary copies (clear large temporaries early)
% (b) worker-side memory monitoring (RSS logging to per-worker files)

%%% Run parameters
arguments
    options.num_channels = 16;
    options.channel_spacing = 400e9;
    options.fiber_length_km = 10;
    options.rand_key = 1;
    options.num_realiz = 1;
    options.fwm_mitigation_technique = "co";
end

%% Add the imdd_simulation framework to the path
if ispc
    addpath(genpath('C:\Users\Silas\Documents\MATLAB\imdd_simulation'));
else
    addpath(genpath('/work_beegfs/sutef391/imdd_simulation'));
end

warning('off','MATLAB:datetime:AmbiguousTimeZone');

%% How many workers?
cpus = str2double(getenv('SLURM_CPUS_PER_TASK'));
if ~isfinite(cpus) || cpus < 1, cpus = max(1, feature('numcores')); end

% Use a per-job, node-local JobStorageLocation to avoid stale locks on $HOME
tmpbase = getenv('TMPDIR');
if isempty(tmpbase), tmpbase = tempdir; end
jsl = fullfile(tmpbase, sprintf('matlab_jobstorage_%s_%s', ...
    getenv('USER'), getenv('SLURM_JOB_ID')));
if ~exist(jsl,'dir'); mkdir(jsl); end

% Start pool once
c = parcluster('local');
c.NumWorkers = cpus;
c.JobStorageLocation = jsl;

p = gcp('nocreate');
if isempty(p) || p.NumWorkers ~= cpus
    if ~isempty(p), delete(p); end
    p = parpool(c, cpus);
end
fprintf('parpool up with %d workers; JobStorage=%s\n', p.NumWorkers, c.JobStorageLocation);

%% result filename (timestamp + optional job id)
t      = datetime('now','TimeZone','local','Format','yyyyMMdd_HHmmss');
jobid  = getenv('SLURM_JOB_ID'); if isempty(jobid), jobid = 'nojid'; end
host   = getenv('HOSTNAME'); if isempty(host), host = 'localhost'; end

foldname = sprintf('%dkm_%dch_%dghz_%s', options.fiber_length_km(end), options.num_channels, options.channel_spacing.*1e-9, options.fwm_mitigation_technique);
if ispc
    output_root = fullfile('C:\Users\Silas\Documents\MATLAB\Datensätze\FWM_2025\',foldname,'\');
else
    output_root = fullfile('/work_beegfs/sutef391/results_WDM',foldname,'\');
end
if ~exist(output_root,'dir'), mkdir(output_root); end

fname = sprintf('WDM_%s_%s_%s_%dkm_%dch_%dghz_%s.mat', char(t), host, jobid, options.fiber_length_km(end), options.num_channels, options.channel_spacing.*1e-9, options.fwm_mitigation_technique);

% Worker memory logs directory
memlog_dir = fullfile(output_root, 'memlogs');
if ~exist(memlog_dir,'dir'), mkdir(memlog_dir); end

%% Settings
s.num_realiz        = options.num_realiz;
s.wavelengthplan    = calcWavelengthPlan(options.num_channels, options.channel_spacing, 1310);
link_length         = options.fiber_length_km;
s.pmd               = 0;%0.1;
s.gamma             = 0;%0.0023;

s.M                 = 4;
fsym                = 112e9;
fdac                = 2*fsym;
fadc                = 120000000000;
s.random_key        = options.rand_key;

% Laser / s.Modulator
vbias_rel           = 0.5;
u_pi                = 4.6;
vbias               = -vbias_rel*u_pi;
laser_linewidth     = 0e6;

% EQ SETTINGS
dfe_order   = [0 0 0];
len_tr      = 4096*2;
mu_ffe1     = 0.0001;
mu_ffe2     = 0.0008;
mu_ffe3     = 0.001;
mu_dc       = 0.005;
mu_ffe      = [mu_ffe1 mu_ffe3 mu_ffe3];
mu_dfe      = 0.0004;

% DB Stuff
db_precode      = 0;
db_encode       = 0;
duob_mode       = db_mode.no_db;
apply_pulsef    = 0;

rcalpha     = 0.05;
Pform       = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha);

N = numel(s.wavelengthplan);
f_plan = physconst('lightspeed')./(s.wavelengthplan.*1e-9);
margin = 2e12; % some THz left and right
f_span = (max(f_plan)+margin)-(min(f_plan)-margin);
f_nyq  = f_span/2;

kover = 4;
upsample_required = f_nyq./(fdac*kover/2);
upsample_pow      = 2^nextpow2(upsample_required);

s.f_opt     = fdac*kover*upsample_pow;
s.f_opt_nyq = s.f_opt/2;

s.rop = -12:0.75:-0.75;
% s.rop = -12:1:0;
%% ---------- Intermediate evaluation points (distance dimension) ----------
segment_length = 1; % km (must match fiber loop below)

eval_dist_km = [0];
eval_dist_km = eval_dist_km(eval_dist_km <= link_length);

if link_length > 0
    if isempty(eval_dist_km) || eval_dist_km(end) ~= link_length
        eval_dist_km = unique([eval_dist_km link_length], 'stable');
    end
end

eval_seg = eval_dist_km ./ segment_length;
if any(abs(eval_seg - round(eval_seg)) > 1e-12)
    error('eval_dist_km must be integer multiples of segment_length=%g km.', segment_length);
end
eval_seg = round(eval_seg);
nEval    = numel(eval_seg);
% -------------------------------------------------------------------------

%% Preallocate outputs (eval distance as 4th dimension)
output_ffe  = cell(N, length(s.rop), s.num_realiz, nEval);
output_dfe  = cell(N, length(s.rop), s.num_realiz, nEval);
output_vnle = cell(N, length(s.rop), s.num_realiz, nEval);
output_mlse = cell(N, length(s.rop), s.num_realiz, nEval);
output_dbt  = cell(N, length(s.rop), s.num_realiz, nEval);

s.p_launch = 3;
s.p = options.fwm_mitigation_technique;

switch s.p
    case "co"
        pol_rot = 100.*ones(1,N);
        d_local = 0;
    case "pair"
        pol_rot = repmat([100,100,0,0],1,N/4);
        d_local = 0;
    case "alt"
        pol_rot = repmat([100,0,100,0],1,N/4);
        d_local = 0;
    case "seg"
        pol_rot = 100.*ones(1,N);
        d_local = 3;
    otherwise
        error('Unknown fwm_mitigation_technique: %s', string(s.p));
end

for realiz = 1:s.num_realiz

    % Per-realization TX storage (needed later for DSP: Symbols/Tx_bits)
    signal_cell = cell(1,N);
    Symbols     = cell(1,N);
    Tx_bits     = cell(1,N);

    % -------- Job queue containers --------
    F = parallel.FevalFuture.empty(0,1);
    meta = struct('l',{},'ri',{},'realiz',{},'eval_ptr',{});
    % -------------------------------------

    %% ---------- TX per channel ----------
    for l = 1:N

        [Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ...
            "fsym",fsym,"M",s.M,"order",17,"useprbs",0, ...
            "fs_out",fdac, ...
            "applyclipping",0,"clipfactor",1.5, ...
            "applypulseform",apply_pulsef,"pulseformer",Pform, ...
            "randkey",s.random_key+l+realiz, ...
            "db_precode",db_precode,"db_encode",db_encode, ...
            "mrds_code",0,"mrds_blocklength",512,"duobinary_mode",duob_mode ...
            ).process();

        Lp_awg = Filter('filtdegree',3,"f_cutoff",56e9,"fs",fdac*kover, ...
            "filterType",filtertypes.gaussian,"active",true);

        El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",1,"kover",kover, ...
            "bit_resolution",6,"upsampling_method","samplehold","precomp_sinc_rolloff",0, ...
            "H_lpf",Lp_awg,"dac_max",0.6,"dac_min",-0.6).process(Digi_sig);

        % Digi_sig not needed after AWG
        clear Digi_sig

        % Electrical Driver Amplifier
        El_sig = El_sig.normalize("mode","oneone");
        scaling = 0.6*(u_pi/2-abs(vbias-u_pi/2));
        El_sig = El_sig .* scaling;

        % E/O Conversion
        Eml_out = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs, ...
            "lambda",s.wavelengthplan(l),"bias",vbias,"u_pi",u_pi, ...
            "linewidth",laser_linewidth,"randomkey",s.random_key+l+realiz,"alpha",0.8).process(El_sig);

        % El_sig not needed after EML
        clear El_sig

        signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",pol_rot(l)).process(Eml_out);

        % Eml_out not needed after pol controller
        clear Eml_out Lp_awg
    end

    disp('Signal generated for all channels.');

    %% ---------- WDM mux + launch ----------
    Opt_sig_wdm = Optical_Multiplex("fs_in",fdac*kover,"fs_out",upsample_pow*fdac*kover, ...
        "lambda_center",1310,"random_key",0,"filtype",1,"B",120e9).process(signal_cell);

    % IMPORTANT: signal_cell is not needed anymore after multiplex (Symbols/Tx_bits remain)
    clear signal_cell

    Opt_sig_wdm = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ...
        "amplification_db",s.p_launch+10*log10(N)).process(Opt_sig_wdm);

    %% ---------- Fiber propagation ----------
    Opt_sig_wdm_fib = Opt_sig_wdm;

    % Opt_sig_wdm no longer needed as separate handle/object
    clear Opt_sig_wdm

    nSegments = link_length/segment_length;
    if abs(nSegments - round(nSegments)) > 1e-12
        error('fiber_length_km=%g must be an integer multiple of segment_length=%g km.', link_length, segment_length);
    end
    nSegments = round(nSegments);

    zdw = 1310;
    randomize_D = true;

    % Guard for 0 km: avoid calling getDispersionVector(0,...) if it doesn't support it
    if nSegments > 0
        Dvec = getDispersionVector(nSegments, d_local, zdw, randomize_D, s.random_key+realiz);
    else
        Dvec = [];
    end

    eval_ptr = 1;
    % =================== Queue throttle (prevents OOM) ===================
    % Limit how many futures can be outstanding (running + queued + finished-not-yet-fetched).
    % Start conservatively; you can raise to 4*NumWorkers if stable.
    maxInFlight = max(2, p.NumWorkers);
    % =====================================================================

    % -------- Evaluate at 0 km (BTB) if requested --------
    if eval_ptr <= nEval && eval_seg(eval_ptr) == 0
        disp('0 km before demux.');
        Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1, ...
            "fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_wdm_fib);

        cnt = 0;
        total = length(s.rop)*N;
        fprintf('total of %d jobs to enqueue at 0 km\n', total);

        for ri = 1:length(s.rop)
            for l = 1:N
                cnt = cnt + 1;

                % ---- Throttle before enqueueing more futures ----
                [F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt] = ...
                    throttle_inflight(F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, maxInFlight);

                % Enqueue one job per (l,ri)
                F(end+1,1) = parfeval(p, @rx_job, 5, ...
                    Opt_sig_wdm_demux{l}, ...
                    s.rop(ri), ...
                    Symbols{l}, Tx_bits{l}, ...
                    s, l, ri, realiz, eval_ptr, ...
                    fdac, kover, fadc, fsym, ...
                    len_tr, mu_dc, mu_ffe, mu_dfe, dfe_order, duob_mode, ...
                    memlog_dir);

                fprintf('Enqueued job %d/%d for realiz %d/%d, l=%d/%d, ri=%d/%d at 0 km (inflight=%d/%d)\n', ...
                    cnt, total, realiz, s.num_realiz, l, N, ri, length(s.rop), numel(F), maxInFlight);

                meta(end+1) = struct('l',l,'ri',ri,'realiz',realiz,'eval_ptr',eval_ptr);
            end
        end

        clear Opt_sig_wdm_demux
        eval_ptr = eval_ptr + 1;
    end
    % -----------------------------------------------------


    for seg = 1:nSegments

        fprintf('Realiz %d/%d: Segment %d/%d \n', realiz, s.num_realiz, seg, nSegments);

        Opt_sig_wdm_fib = DP_Fiber("L",segment_length,"D",Dvec(seg),"Dpmd",s.pmd,"Ds",0.07, ...
            "beat_len",10,"corr_len",100,"dz",1,"manakov",0, ...
            "gamma",s.gamma,"lambda",zdw,"n_waveplates",10,"SS_dphimax",0.01, ...
            "SS_dzmax",50,"SS_dzmin",10,"X_alpha",0.3,"X_beta",0,"rng",1).process(Opt_sig_wdm_fib);

        % -------- Evaluate at intermediate distance (enqueue jobs) --------
        if eval_ptr <= nEval && seg == eval_seg(eval_ptr)

            Opt_sig_wdm_demux = Optical_Demultiplex("attenuation",0,"B",200e9,"filtype",1, ...
                "fs_out",fdac*kover,"fs_in",fdac*kover*upsample_pow,"lambda_center",1310).process(Opt_sig_wdm_fib);

            cnt = 0;
            total = length(s.rop)*N;
            fprintf('total of %d jobs to enqueue at 0 km\n', total);

            for ri = 1:length(s.rop)
                for l = 1:N

                    cnt = cnt + 1;

                    % ---- Throttle before enqueueing more futures ----
                    [F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt] = ...
                        throttle_inflight(F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, maxInFlight);

                    F(end+1,1) = parfeval(p, @rx_job, 5, ...
                        Opt_sig_wdm_demux{l}, ...
                        s.rop(ri), ...
                        Symbols{l}, Tx_bits{l}, ...
                        s, l, ri, realiz, eval_ptr, ...
                        fdac, kover, fadc, fsym, ...
                        len_tr, mu_dc, mu_ffe, mu_dfe, dfe_order, duob_mode, ...
                        memlog_dir);

                    fprintf('Enqueued job %d/%d for realiz %d/%d, l=%d/%d, ri=%d/%d at 0 km (inflight=%d/%d)\n', ...
                        cnt, total, realiz, s.num_realiz, l, N, ri, length(s.rop), numel(F), maxInFlight);

                    meta(end+1) = struct('l',l,'ri',ri,'realiz',realiz,'eval_ptr',eval_ptr);
                end
            end

            clear Opt_sig_wdm_demux
            eval_ptr = eval_ptr + 1;
        end
        % ----------------------------------------------------------------

        % Non-blocking harvest (your existing line)
        [F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt] = ...
            collect_done(F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, false);

    end

    while ~isempty(F)
        [F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt] = ...
            collect_done(F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, true);
    end

    %% Save results (per realization)
    res = struct();
    res.settings      = s;
    res.eval_dist_km  = eval_dist_km;
    res.ffe           = output_ffe;
    res.dfe           = output_dfe;
    res.vnle          = output_vnle;
    res.mlse          = output_mlse;
    res.dbt           = output_dbt;

    save(fullfile(output_root, fname), 'res', '-v7.3');
    fprintf('Saved results to: %s\n', fullfile(output_root, fname));
    disp(datetime('now','TimeZone','local','Format','yyyyMs.Mdd_HHmmss'));

    % Per-realization large arrays that are no longer needed
    clear Opt_sig_wdm_fib Symbols Tx_bits Dvec

end

end % end main


%% ========================= Local functions =========================

function [ffe_results, dfe_results, vnle_results, mlse_results, dbt_results] = rx_job( ...
    Opt_sig_chan, rop_db, Symbols_l, Tx_bits_l, s, l, ri, realiz, eval_ptr, ...
    fdac, kover, fadc, fsym, len_tr, mu_dc, mu_ffe, mu_dfe, dfe_order, duob_mode, memlog_dir)

% NOTE: keep plotting OFF in workers
debug_plots = 0;

% Create per-worker logfile (avoid contention)
logfile = make_worker_logfile(memlog_dir);

log_mem(logfile, 'job_start', l, ri, realiz, eval_ptr, rop_db, Opt_sig_chan);

%%%%%% ROP %%%%%%
Opt_sig_wdm_rx = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power", ...
    "amplification_db", rop_db).process(Opt_sig_chan);

% Opt_sig_chan no longer needed after amp
clear Opt_sig_chan
log_mem(logfile, 'after_amp', l, ri, realiz, eval_ptr, rop_db, Opt_sig_wdm_rx);

%%%%%% PD Square Law %%%%%%
assert(fdac*kover==Opt_sig_wdm_rx.fs,'Sampling Frequencies do not match! Check previous steps');
PD_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20, ...
    "nep",1.8e-11,"randomkey",s.random_key + l + realiz).process(Opt_sig_wdm_rx);

% Opt_sig_wdm_rx not needed after PD
clear Opt_sig_wdm_rx
log_mem(logfile, 'after_pd', l, ri, realiz, eval_ptr, rop_db, PD_sig);

%%%%%% Low-pass RX (PD, El. Connectors and Scope) %%%%%%
rx_bwl = 100e9;
PD_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover, ...
    "filterType",filtertypes.butterworth,"active",true).process(PD_sig);
log_mem(logfile, 'after_rx_lpf', l, ri, realiz, eval_ptr, rop_db, PD_sig);

%%%%%% Low-pass Scope %%%%%%
Lp_scpe = Filter('filtdegree',4,"f_cutoff",80e9,"fs",fadc, ...
    "filterType",filtertypes.butterworth,"active",true);

%%%%%% Scope %%%%%%
Scpe_sig = Scope("fsimu",fdac*kover,"fadc",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',0,'H_lpf',Lp_scpe).process(PD_sig);

% PD_sig no longer needed after scope
clear PD_sig Lp_scpe
log_mem(logfile, 'after_scope', l, ri, realiz, eval_ptr, rop_db, Scpe_sig);

Scpe_sig_2sps = Scpe_sig.resample("fs_out",2*fsym);

% Scpe_sig no longer needed after resample
clear Scpe_sig
log_mem(logfile, 'after_resample', l, ri, realiz, eval_ptr, rop_db, Scpe_sig_2sps);

[~, Scpe_cell, ~, ~] = Scpe_sig_2sps.tsynch("reference", Symbols_l, "fs_ref", fsym, "debug_plots", debug_plots);

% Scpe_sig_2sps no longer needed after sync
clear Scpe_sig_2sps
log_mem(logfile, 'after_sync', l, ri, realiz, eval_ptr, rop_db);

Rx_sig = Scpe_cell{1}.normalize("mode","rms");

% Scpe_cell no longer needed
clear Scpe_cell
log_mem(logfile, 'after_rxsig', l, ri, realiz, eval_ptr, rop_db, Rx_sig);

% -------------------- FFE --------------------
ffe_order   = [50, 0, 0];
eq_ffe      = 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);

ffe_results = ffe(eq_ffe,s.M,Rx_sig,Symbols_l,Tx_bits_l, ...
    "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ...
    "eth_style_symbol_mapping",0);
clear eq_ffe
log_mem(logfile, 'after_ffe', l, ri, realiz, eval_ptr, rop_db);

% -------------------- DFE --------------------
eq_dfe      = EQ("Ne",ffe_order,"Nb",[2,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);

dfe_results = ffe(eq_dfe,s.M,Rx_sig,Symbols_l,Tx_bits_l, ...
    "precode_mode",duob_mode,'showAnalysis',0,"postFFE",[], ...
    "eth_style_symbol_mapping",0);
clear eq_dfe
log_mem(logfile, 'after_dfe', l, ri, realiz, eval_ptr, rop_db);

% -------------------- VNLE + MLSE --------------------
pf_ncoeffs  = 1;
ffe_order3  = [50, 5, 5];
eq_v        = EQ("Ne",ffe_order3,"Nb",dfe_order, ...
    "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",1);
pf_         = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1);

mlse_ = MLSE("duobinary_output",0,'M',s.M,'trellis_states',PAMmapper(s.M,0).levels);

[vnle_results, mlse_results] = vnle_postfilter_mlse(eq_v, pf_, mlse_, s.M, Rx_sig, Symbols_l, Tx_bits_l, ...
    "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", [], "eth_style_symbol_mapping", 0);

clear eq_v pf_ mlse_
log_mem(logfile, 'after_vnle_mlse', l, ri, realiz, eval_ptr, rop_db);

% -------------------- DB target --------------------
mlse_db_ = MLSE("DIR",[1,1],"duobinary_output",0,"M",s.M,'trellis_states',PAMmapper(s.M,0).levels);
ffe_order = [50, 5, 5];
eq_       = EQ("Ne",ffe_order,"Nb",dfe_order,"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",1);
dbt_results = duobinary_target(eq_,mlse_db_, s.M, Rx_sig, Symbols_l, Tx_bits_l, ...
    "precode_mode", duob_mode, 'showAnalysis', 0, "postFFE", []);

clear mlse_db_
log_mem(logfile, 'job_end', l, ri, realiz, eval_ptr, rop_db);

% Rx_sig no longer needed
clear Rx_sig Symbols_l Tx_bits_l
end


function [F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt] = ...
    collect_done(F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, block)

if nargin < 9, block = false; end
if isempty(F), return; end

if block
    timeout_first = Inf;
else
    timeout_first = 0;
end

first = true;
while ~isempty(F)
    if first
        timeout = timeout_first;
        first = false;
    else
        timeout = 0;
    end

    try
        [k, ffe_r, dfe_r, vnle_r, mlse_r, dbt_r] = fetchNext(F, timeout);
    catch
        break;
    end

    if isempty(k)
        break;
    end

    m = meta(k);

    output_ffe{m.l, m.ri, m.realiz, m.eval_ptr}  = ffe_r;
    output_dfe{m.l, m.ri, m.realiz, m.eval_ptr}  = dfe_r;
    output_vnle{m.l, m.ri, m.realiz, m.eval_ptr} = vnle_r;
    output_mlse{m.l, m.ri, m.realiz, m.eval_ptr} = mlse_r;
    output_dbt{m.l, m.ri, m.realiz, m.eval_ptr}  = dbt_r;

    F(k) = [];
    meta(k) = [];
end
end


%% ========================= Memory logging helpers =========================

function logfile = make_worker_logfile(memlog_dir)
% One file per worker process to avoid write contention.
pid = get_pid_safe();
t = getCurrentTask();
if isempty(t)
    tid = -1;
else
    tid = t.ID;
end
logfile = fullfile(memlog_dir, sprintf('memlog_pid%d_task%d.txt', pid, tid));
end

function log_mem(logfile, tag, l, ri, realiz, eval_ptr, rop_db, varargin)
% Append one line with RSS/VmSize plus optional "largest variable" info.
% Uses /proc on Linux where available.
ts = datetime('now','TimeZone','local','Format','yyyy-MM-dd HH:mm:ss.SSS');

[rssMB, vmsMB] = proc_mem_mb();

% Optional: include size of a specific variable/object if provided
extra = "";
if ~isempty(varargin)
    try
        x = varargin{1}; %#ok<NASGU>
        w = whos('x');
        extra = sprintf(' | x_bytes=%d', w.bytes);
    catch
        extra = " | x_bytes=NA";
    end
end

line = sprintf('%s | %s | l=%d ri=%d realiz=%d eval=%d rop=%.3f | RSS=%.1fMB Vm=%.1fMB%s\n', ...
    char(ts), tag, l, ri, realiz, eval_ptr, rop_db, rssMB, vmsMB, extra);

fid = fopen(logfile, 'a');
if fid ~= -1
    fwrite(fid, line);
    fclose(fid);
end
end

function [rssMB, vmsMB] = proc_mem_mb()
% RSS/VmSize from /proc (Linux). Falls back to NaN if unavailable.
rssMB = NaN; vmsMB = NaN;

if isunix
    try
        txt = fileread('/proc/self/status');
        rssMB = parse_kb_field(txt, 'VmRSS:') / 1024;
        vmsMB = parse_kb_field(txt, 'VmSize:') / 1024;
        return;
    catch
        % fall through
    end
end

% Fallback (Windows): memory() sometimes works
if ispc
    try
        m = memory;
        % MemUsedMATLAB is bytes
        rssMB = double(m.MemUsedMATLAB) / 1024^2;
        vmsMB = NaN;
    catch
    end
end
end

function kb = parse_kb_field(txt, key)
kb = NaN;
idx = strfind(txt, key);
if isempty(idx), return; end
i = idx(1) + length(key);
% read until end of line
j = find(txt(i:end)==newline, 1, 'first') + i - 2;
val = strtrim(txt(i:j));
% format: "123456 kB"
parts = split(val);
kb = str2double(parts{1});
end

function pid = get_pid_safe()
pid = -1;
try
    pid = feature('getpid');
catch
    % no-op
end
end

function [F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt] = ...
    throttle_inflight(F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt, maxInFlight)
%THROTTLE_INFLIGHT
% Enforces a HARD cap on outstanding futures (running + queued + finished-not-yet-fetched).
% Blocks until numel(F) < maxInFlight, then returns.

while numel(F) >= maxInFlight

    % Wait until at least one finishes OR fails
    % (wait('finished',1) only wakes for finished, so we also handle failures via periodic checks)
    try
        wait(F, 'finished', 1);
    catch
        % If wait fails for some reason, fall back to a short pause
        pause(0.1);
    end

    % Harvest finished and failed futures (robust snapshot)
    n = numel(F);

    % Safety: meta must match F length
    if numel(meta) ~= n
        % Keep them aligned (rare; indicates earlier mismatch)
        m = min(numel(meta), n);
        F = F(1:m);
        meta = meta(1:m);
        n = m;
    end

    states = cell(1,n);
    for k = 1:n
        states{k} = F(k).State;
    end

    doneFinished = strcmp(states, 'finished');
    doneFailed   = strcmp(states, 'failed');
    done         = doneFinished | doneFailed;

    idxDone = find(done);

    if isempty(idxDone)
        pause(0.05);
        continue;
    end

    % Fetch outputs for finished; mark failed as empty results
    for jj = 1:numel(idxDone)
        ii = idxDone(jj);
        m  = meta(ii);

        if doneFinished(ii)
            [ffe_r, dfe_r, vnle_r, mlse_r, dbt_r] = fetchOutputs(F(ii));
        else
            % Failed
            try
                err = F(ii).Error;
                warning('rx_job failed: realiz=%d eval=%d l=%d ri=%d | %s', ...
                    m.realiz, m.eval_ptr, m.l, m.ri, err.message);
            catch
            end
            ffe_r = []; dfe_r = []; vnle_r = []; mlse_r = []; dbt_r = [];
        end

        output_ffe{m.l, m.ri, m.realiz, m.eval_ptr}  = ffe_r;
        output_dfe{m.l, m.ri, m.realiz, m.eval_ptr}  = dfe_r;
        output_vnle{m.l, m.ri, m.realiz, m.eval_ptr} = vnle_r;
        output_mlse{m.l, m.ri, m.realiz, m.eval_ptr} = mlse_r;
        output_dbt{m.l, m.ri, m.realiz, m.eval_ptr}  = dbt_r;
    end

    % Remove harvested futures/meta entries using numeric indices (safe)
    F(idxDone) = [];
    meta(idxDone) = [];
end

end