diff --git a/Classes/00_signals/Signal.m b/Classes/00_signals/Signal.m index b8113c8..f6dd186 100644 --- a/Classes/00_signals/Signal.m +++ b/Classes/00_signals/Signal.m @@ -765,7 +765,7 @@ classdef Signal pkpos = sort(pkpos); - if mean(w) > 10 || mean(p) > 10 + if mean(w) > 10 || mean(p) > 15 return else sequenceFound = 1; diff --git a/Classes/02_optical/EML.m b/Classes/02_optical/EML.m index 32a902c..4a97686 100644 --- a/Classes/02_optical/EML.m +++ b/Classes/02_optical/EML.m @@ -15,6 +15,7 @@ classdef EML u_pi randomkey randomstream + alpha %on instance creation field @@ -38,6 +39,7 @@ classdef EML options.lambda; options.power; options.linewidth = 0; + options.alpha = 0; options.ampl_imbal = 0; options.pha_imbal = 0; options.bias; @@ -101,9 +103,21 @@ classdef EML end %modulate the laserfield with the electrical signal - opt_out = obj.externalmodulation(laserfield,elec_in); + laserfield = obj.externalmodulation(laserfield,elec_in); + % add chirp + opt_out = obj.chirp(laserfield); + end + + function chirped_field = chirp(obj,laserfield) + + % Chirp + p = abs(laserfield.^2); + derv_p = [0; diff(p)]; + delta_phi = derv_p./(4*pi*p).*obj.alpha; + delta_phi = cumsum(delta_phi); + chirped_field = laserfield.*exp(1i*2*pi*delta_phi); end diff --git a/Classes/DataBaseHandler/Metricstruct.m b/Classes/DataBaseHandler/Metricstruct.m index ef00675..fe9a818 100644 --- a/Classes/DataBaseHandler/Metricstruct.m +++ b/Classes/DataBaseHandler/Metricstruct.m @@ -8,9 +8,9 @@ classdef Metricstruct date_of_processing (1,1) datetime = datetime('now') numBits (1,1) double {mustBeInteger, mustBeNonnegative} = 0 - BER (1,1) double {mustBeNumeric, mustBeNonnegative, mustBeLessThanOrEqual(BER,1)} = 0 + BER (1,1) double {mustBeNumeric, mustBeNonnegative} = 0 numBitErr (1,1) double {mustBeInteger, mustBeNonnegative} = 0 - BER_precoded (1,1) double {mustBeNumeric, mustBeNonnegative, mustBeLessThanOrEqual(BER_precoded,1)} = 0 + BER_precoded (1,1) double {mustBeNumeric, mustBeNonnegative} = 0 numBitErr_precoded (1,1) double {mustBeInteger, mustBeNonnegative} = 0 SNR (1,1) double {mustBeNumeric} = NaN diff --git a/projects/IMDD_base_system/imdd_it.m b/projects/IMDD_base_system/imdd_it.m index bd2461d..748708f 100644 --- a/projects/IMDD_base_system/imdd_it.m +++ b/projects/IMDD_base_system/imdd_it.m @@ -10,12 +10,12 @@ if 1 uloops.laser_wavelength = [1310]; uloops.M = [4]; uloops.link_length = [1]; % 1,2,3,5,6,8,10 - uloops.interference_attenuation = [0,3,6,9,12,15,18,21,24,27,30,45]; + % uloops.interference_attenuation = [0,3,6,9,12,15,18,21,24,27,30,45]; wh = DataStorage(uloops); wh.addStorage("ber"); % wh = submit_simulations(wh,"parallel",0,"simulation_mode",0); - wh = submit_handle(@imdd_model,wh,"parallel",0); + wh = submit_handle(@imdd_model,wh,"serial",0); end diff --git a/projects/IMDD_base_system/imdd_model.m b/projects/IMDD_base_system/imdd_model.m index d400350..ac6bd0c 100644 --- a/projects/IMDD_base_system/imdd_model.m +++ b/projects/IMDD_base_system/imdd_model.m @@ -92,34 +92,34 @@ end if doub_mode ~= db_mode.db_encoded if precomp == 0 && db_precode == 1 - doub_mode = db_mode.db_precoded; - + doub_mode = db_mode.db_precoded; + db_precode = 1; % preceded data (in my measurement set, this corresponds to low precomp too!) - discard_precode = 0; % + discard_precode = 0; % emulate_precode = 0; legendentry = 'low precomp; precoded'; disp('low precomp; precoded') elseif precomp == 1 && db_precode == 1 - doub_mode = db_mode.db_emulate; - + doub_mode = db_mode.db_emulate; + db_precode = 0; % preceded data (in my measurement set, this corresponds to low precomp too!) - discard_precode = 0; % + discard_precode = 0; % emulate_precode = 1; legendentry = 'high precomp; precoded'; disp('high precomp; precoded') elseif precomp == 0 && db_precode == 0 doub_mode = db_mode.db_discard; - + db_precode = 1; % preceded data (in my measurement set, this corresponds to low precomp too!) - discard_precode = 1; % + discard_precode = 1; % emulate_precode = 0; legendentry = 'no precomp; not precoded'; disp('no precomp; not precoded') elseif precomp == 1 && db_precode == 0 - doub_mode = db_mode.no_db; - + doub_mode = db_mode.no_db; + db_precode = 0; % preceded data (in my measurement set, this corresponds to low precomp too!) - discard_precode = 0; % + discard_precode = 0; % emulate_precode = 0; legendentry = 'high precomp; not precoded'; disp('high precomp; not precoded') @@ -137,135 +137,73 @@ end f_nyquist = fsym/2; -%%% run the simulation or measurement or ... -if simulation_mode - Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); - rcalpha = 1; - Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha); - db_precode = 0; - db_encode = 0; - apply_pulsef = 1; - [Digi_sig,Symbols,Tx_bits] = PAMsource(... - "fsym",fsym,"M",M,"order",18,"useprbs",0,... - "fs_out",fdac,... - "applyclipping",0,"clipfactor",1.5,... - "applypulseform",apply_pulsef,"pulseformer",Pform,... - "randkey",random_key,... - "db_precode",db_precode,"db_encode",db_encode,... - "mrds_code",0,"mrds_blocklength",512).process(); +Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",16,"alpha",rcalpha); +rcalpha = 1; +Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha); - Digi_sig.spectrum("displayname",'Digi Spectrum','fignum',10,'normalizeTo0dB',1); +db_precode = 0; +db_encode = 0; +apply_pulsef = 1; +[Digi_sig,Symbols,Tx_bits] = PAMsource(... + "fsym",fsym,"M",M,"order",18,"useprbs",0,... + "fs_out",fdac,... + "applyclipping",0,"clipfactor",1.5,... + "applypulseform",apply_pulsef,"pulseformer",Pform,... + "randkey",random_key,... + "db_precode",db_precode,"db_encode",db_encode,... + "mrds_code",0,"mrds_blocklength",512).process(); + +Digi_sig.spectrum("displayname",'Digi Spectrum','fignum',10,'normalizeTo0dB',1); - %%%%% AWG - % El_sig = M8199A("kover",kover).process(Digi_sig); - El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",1).process(Digi_sig); - % El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',0); - % El_sig = El_sig.setPower(0,"dBm"); +%%%%% AWG +% El_sig = M8199A("kover",kover).process(Digi_sig); +El_sig = AWG("fdac",fdac,"f_cutoff",fsym,"lpf_active",0,"kover",kover,"bit_resolution",12,"upsampling_method","samplehold","precomp_sinc_rolloff",1).process(Digi_sig); +% El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',0); +% El_sig = El_sig.setPower(0,"dBm"); - %%%%% Low-pass el. components %%%%%% - tx_bwl = tx_bw_nyquist.*f_nyquist; - % tx_bwl = 80e9; - El_sig = Filter('filtdegree',4,"f_cutoff",tx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); - % El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',1); +%%%%% Low-pass el. components %%%%%% +tx_bwl = tx_bw_nyquist.*f_nyquist; +% tx_bwl = 80e9; +El_sig = Filter('filtdegree',4,"f_cutoff",tx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(El_sig); +% El_sig.spectrum("displayname",'Digi Spectrum','fignum',100,'normalizeTo0dB',1); - %%%%% Electrical Driver Amplifier %%%%%% - El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",3).process(El_sig); - El_sig = El_sig.normalize("mode","oneone"); +%%%%% Electrical Driver Amplifier %%%%%% +El_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","gain","amplification_db",3).process(El_sig); +El_sig = El_sig.normalize("mode","oneone"); - %%%%% MODULATE E/O CONVERSION %%%%%% - [Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",laser_wavelength,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+1).process(El_sig); +%%%%% MODULATE E/O CONVERSION %%%%%% +[Opt_sig] = EML("mode",eml_mode.im_cosinus,"power",3,"fsimu",El_sig.fs,"lambda",laser_wavelength,"bias",vbias,"u_pi",u_pi,"linewidth",laser_linewidth,"randomkey",random_key+1).process(El_sig); - Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); +Opt_sig = Fiber("fsimu",Opt_sig.fs,"fiber_length",link_length/1000,"alpha",0.3,"D",0,"lambda0",1310,"gamma",0,"Dslope",0.07).process(Opt_sig); - %%%%%% ROP %%%%%% - Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig); +%%%%%% ROP %%%%%% +Rx_sig = Amplifier("amp_mode","ideal_no_noise","gain_mode","output_power","amplification_db",rop).process(Opt_sig); - %%%%%% PD Square Law %%%%%% - Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11).process(Rx_sig); +%%%%%% PD Square Law %%%%%% +Rx_sig = Photodiode("fsimu",fdac*kover,"dark_current",2e-08,"responsivity",1,"temperature",20,"nep",1.8e-11).process(Rx_sig); - %%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%% - rx_bwl = rx_bw_nyquist.*f_nyquist; - % rx_bwl = 80e9; - Rx_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(Rx_sig); +%%%%%% Low-pass RX (PD, El. Connectors and Scope %%%%%% +rx_bwl = rx_bw_nyquist.*f_nyquist; +% rx_bwl = 80e9; +Rx_sig = Filter('filtdegree',4,"f_cutoff",rx_bwl,"fs",fdac*kover,"filterType",filtertypes.butterworth,"active",true).process(Rx_sig); - % %%%%%% Low-pass Scope %%%%%% - Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); +% %%%%%% Low-pass Scope %%%%%% +Lp_scpe = Filter('filtdegree',4,"f_cutoff",110e9,"fs",fadc,"filterType",filtertypes.butterworth,"active",true); - % Rx_sig.spectrum("displayname",'Analog Rx Spectrum','fignum',100,'normalizeTo0dB',1); +% Rx_sig.spectrum("displayname",'Analog Rx Spectrum','fignum',100,'normalizeTo0dB',1); - %%%%%% 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',1,'H_lpf',Lp_scpe).process(Rx_sig); +%%%%%% 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',1,'H_lpf',Lp_scpe).process(Rx_sig); - Scpe_cell{1} = Scpe_sig; +Scpe_cell{1} = Scpe_sig; -else - profile on - basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; - database = DBHandler("pathToDB",[basePath,'silas_labor.db']); - profile off - basePath = 'C:\Users\Silas\Documents\MATLAB\Datensätze\sioe_labor\'; - useGui = 0; - % db = DBHandler("pathToDB",[basePath,'silas_labor.db']); - filterParams = database.tables; - % filterParams.Runs.run_id = 2958; % no db - % filterParams.Runs.run_id = 2937; % no db - filterParams.Configurations = struct( ... - 'bitrate', bitrate, ... - 'db_mode', db_precode+db_encode, ... - 'fiber_length', link_length, ... - 'interference_attenuation', [], ... - 'interference_path_length', [], ... - 'is_mpi', is_mpi, ... - 'pam_level', M, ... - 'precomp_amp', [], ... - 'rop_attenuation', 0, ... - 'symbolrate', [], ... - 'v_awg', [], ... - 'v_bias', [], ... - 'wavelength', laser_wavelength ... - ); - - selectedFields = {'Runs.run_id','Runs.tx_bits_path', 'Runs.tx_symbols_path', 'Runs.rx_sync_path','Runs.rx_raw_path',... - 'Configurations.db_mode','Configurations.pam_level','Configurations.bitrate','Configurations.symbolrate','Configurations.fiber_length','Configurations.wavelength','Configurations.precomp_amp','Measurements.power_rop','Configurations.v_bias',... - 'Configurations.interference_attenuation'}; - - [dataTable,sql_query] = database.queryDB(filterParams, selectedFields); - [~, uniqueIdx] = unique(dataTable.run_id); % Get unique run_id indices - dataTable = dataTable(uniqueIdx,:); % Extract unique configurations for each run_id - fprintf('Found %d entries for requested Configuration. IDs are: %s \n \n',size(dataTable,1),jsonencode(dataTable.run_id(1:min(size(dataTable,1),100)))); - - Tx_bits = load([basePath, char(dataTable.tx_bits_path(end))]); - Tx_bits = Tx_bits.Bits; - - Symbols = load([basePath, char(dataTable.tx_symbols_path(end))]); - Symbols = Symbols.Symbols; - - Scpe_load = load([basePath, char(dataTable.rx_sync_path(end))]); - Scpe_cell = Scpe_load.S; - - - % Raw_signal = load([basePath, char(dataTable.rx_raw_path(1))]); - % Raw_signal.Scpe_sig_raw.plot("displayname",'0db atten','fignum',10101) - % Raw_signal = Raw_signal.Scpe_sig_raw; - % - % Raw_signal = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.55,"fs",Raw_signal.fs,"filterType",filtertypes.gaussian,"active",true).process(Raw_signal); - % - % Scpe_cell{1}.eye(fsym,M,"displayname",'eye','fignum',227); - % - % Raw_signal.spectrum("normalizeTo0dB",0,"fignum",11,"fft_length",2^12); - % Raw_signal.move_it_spectrum("fignum",334); - % Raw_signal.move_it_spectrum("fignum",334); - - fsym = Symbols.fs; - -end if db_precode Symbols_precoded = Symbols; @@ -291,74 +229,74 @@ for occ = 1%:proc_occ Scpe_sig = Filter('filtdegree',4,"f_cutoff",Symbols.fs.*0.5,"fs",Scpe_sig.fs,"filterType",filtertypes.gaussian,"active",true).process(Scpe_sig); Scpe_sig = Scpe_sig - mean(Scpe_sig.signal); - % + % % Pform = Pulseformer("fsym",Scpe_sig.fs,"fdac",2*fsym,"pulse","rrc","pulselength",16,"alpha",rcalpha,"matched",0); - % + % % Scpe_sig_matched = Pform.process(Scpe_sig); - % + % % Scpe_sig.spectrum("normalizeTo0dB",0,"fignum",336,"displayname","scope "); % Scpe_sig_matched.spectrum("normalizeTo0dB",0,"fignum",336,"displayname","matched"); %%% EQUALIZING - - - % eq_mlse = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0,"dc_buffer_len",1,"mu_dc",0.05); - % eq_mlse = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0); - % eq_mlse = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0,"dc_buffer_len",512,"mu_dc",0.05); - - mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3]; - vnle_order=[vnle_order1,vnle_order2,vnle_order3]; - - % %%%%% VNLE + DFE %%%% - if 0 - - eq_vnle_dfe = EQ("Ne",vnle_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); - eq_2 = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",2001,"sps",1,"decide",0); - - [result] = vnle(eq_vnle_dfe,M,Scpe_sig,Symbols,Tx_bits,"precode_mode",doub_mode,"showAnalysis",1,"postFFE",[]); - vnle_dfe_package{occ} = result; - - end - %%%%% VNLE + PF + MLSE %%%% - if 1 - - % len_tr = length(Symbols)-1000; - eq_vnle_ = EQ("Ne",vnle_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); - % eq_vnle_ = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",vnle_order,"sps",2,"decide",0); - pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); - mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); - - [result] = vnle_postfilter_mlse(eq_vnle_,pf_,mlse_,M,Scpe_sig,Symbols,Tx_bits,"precode_mode",doub_mode,'showAnalysis',1); - vnle_pf_package{occ} = result; - - end - %%%%% Duobinary Targeting %%%% - if 1 + % eq_mlse = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0,"dc_buffer_len",1,"mu_dc",0.05); + % eq_mlse = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0); + % eq_mlse = FFE_DCremoval("epochs_tr",5,"epochs_dd",5,"len_tr",len_tr,"mu_dd",mu_ffe(1),"mu_tr",0,"order",ffe_order(1),"sps",2,"decide",0,"dc_buffer_len",512,"mu_dc",0.05); - mlse_db = MLSE_viterbi("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); - eq_db = EQ("Ne",vnle_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); + mu_ffe = [mu_ffe1 mu_ffe2 mu_ffe3]; + vnle_order=[vnle_order1,vnle_order2,vnle_order3]; - [result] = duobinary_target(eq_db, mlse_db, M, Scpe_sig, Symbols, Tx_bits, "precode_mode", doub_mode,'showAnalysis',0); - dbtgt_package{occ} = result; - + % %%%%% VNLE + DFE %%%% + if 0 - end - - %%%%%% %db signaling => db encoded %%%%% - if 0 - mlse_db_enc = MLSE_viterbi("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); - eq_db_enc = EQ("Ne",vnle_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); - [result] = duobinary_signaling(eq_db_enc, mlse_db_enc,M, Scpe_sig ,Symbols, Tx_bits); - dbenc_package{occ} = result; - end + eq_vnle_dfe = EQ("Ne",vnle_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); + eq_2 = FFE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",1e-4,"mu_tr",0,"order",2001,"sps",1,"decide",0); + + [result] = vnle(eq_vnle_dfe,M,Scpe_sig,Symbols,Tx_bits,"precode_mode",doub_mode,"showAnalysis",1,"postFFE",[]); + vnle_dfe_package{occ} = result; + + end + %%%%% VNLE + PF + MLSE %%%% + if 1 + + % len_tr = length(Symbols)-1000; + eq_vnle_ = EQ("Ne",vnle_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); + % eq_vnle_ = VNLE("epochs_tr",5,"epochs_dd",5,"len_tr",4096*2,"mu_dd",[0.0004 0.0005 0.0006],"mu_tr",0,"order",vnle_order,"sps",2,"decide",0); + pf_ = Postfilter("ncoeff",pf_ncoeffs,"useBurg",1); + mlse_ = MLSE_viterbi("duobinary_output",0,'M',M,'trellis_states',PAMmapper(M,0).levels); + + [result] = vnle_postfilter_mlse(eq_vnle_,pf_,mlse_,M,Scpe_sig,Symbols,Tx_bits,"precode_mode",doub_mode,'showAnalysis',1); + vnle_pf_package{occ} = result; + + end + + + %%%%% Duobinary Targeting %%%% + if 1 + + mlse_db = MLSE_viterbi("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); + eq_db = EQ("Ne",vnle_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); + + [result] = duobinary_target(eq_db, mlse_db, M, Scpe_sig, Symbols, Tx_bits, "precode_mode", doub_mode,'showAnalysis',0); + dbtgt_package{occ} = result; + + + end + + %%%%%% %db signaling => db encoded %%%%% + if 0 + mlse_db_enc = MLSE_viterbi("DIR",[1,1],"duobinary_output",0,"M",M,"trellis_states",PAMmapper(M,0).levels); + eq_db_enc = EQ("Ne",vnle_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); + [result] = duobinary_signaling(eq_db_enc, mlse_db_enc,M, Scpe_sig ,Symbols, Tx_bits); + dbenc_package{occ} = result; + end % autoArrangeFigures; disp('- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ') fprintf('\n') - + end diff --git a/projects/IMDD_base_system/simulation_bwl_2.m b/projects/IMDD_base_system/simulation_bwl_2.m index dec8036..e336ad2 100644 --- a/projects/IMDD_base_system/simulation_bwl_2.m +++ b/projects/IMDD_base_system/simulation_bwl_2.m @@ -59,7 +59,7 @@ Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",1 db_precode = 0; db_encode = 0; -duob_mode = db_mode.db_precoded; +duob_mode = db_mode.no_db; apply_pulsef = 1; [Digi_sig,Symbols,Tx_bits] = PAMsource(... "fsym",fsym,"M",M,"order",18,"useprbs",0,... @@ -72,13 +72,12 @@ apply_pulsef = 1; Digi_sig.spectrum("displayname",'Digi Spectrum','fignum',10,'normalizeTo0dB',1); -%% proof of concept -Symbols_db = Duobinary().encode(Symbols); -mim_decoded = Duobinary().decode(Symbols_db,"M",M); -rx_bits_mim_decoded = PAMmapper(M,0,"eth_style",0).demap(mim_decoded); -rx_bits_mim_decoded_.signal = circshift(rx_bits_mim_decoded.signal,0); -[~,~,ber_mim_decode,~] = calc_ber(rx_bits_mim_decoded_.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); -fprintf('BER mim: %.2e \n',ber_mim_decode); +% %% proof of concept memoryless inverse mapping (direct db targeting and decoding) +% Symbols_db = Duobinary().encode(Symbols); +% mim_decoded = Duobinary().decode(Symbols_db,"M",M); +% rx_bits_mim_decoded = PAMmapper(M,0,"eth_style",0).demap(mim_decoded); +% [~,~,ber_mim_decode,~] = calc_ber(rx_bits_mim_decoded_.signal,Tx_bits.signal,"skip_front",100,"skip_end",150,"returnErrorLocation",1); +% fprintf('BER mim: %.2e \n',ber_mim_decode); %% diff --git a/projects/WDM/WDM_auswertung.m b/projects/WDM/WDM_auswertung.m index 60651c1..a6a666c 100644 --- a/projects/WDM/WDM_auswertung.m +++ b/projects/WDM/WDM_auswertung.m @@ -1,207 +1,290 @@ +%% Robust eval/plot script for res.* (handles missing/partial dims) +% res.ffe stored like: output_ffe{ch, rop, realiz, eval_distance} +% BUT: some dimensions may be singleton or not fully filled (e.g. only 0 km, only 1 realization) -%saved like this: output_ffe{ch,rop,realiz,eval_distance} +% -------------------- Basic metadata -------------------- +rop = res.settings.rop(:); +wavelengthplan = res.settings.wavelengthplan(:); +distances = res.eval_dist_km(:); % can be [0], [2 10], etc. +N_ch = numel(wavelengthplan); +N_rop = numel(rop); + +% -------------------- Determine dims robustly -------------------- +% Ensure we have 4D addressing even if MATLAB collapses trailing singletons dims = size(res.ffe); -rop = res.settings.rop; -wavelengthplan = res.settings.wavelengthplan; -N_distances = dims(4); %usually 2km and 10km are evaualted -distances = res.eval_dist_km; %[2,10] -N_ch = length(wavelengthplan); %16 channels -figure(); hold on; -cols = cbrewer2('set2',N_ch); % one color per wavelength (Ch) +if numel(dims) < 4 + dims(end+1:4) = 1; +end +N_realiz = dims(3); +N_distances = dims(4); -fec = 2.2e-4; -fec = 3.8e-3; -Sffe = cell(1,N_ch); -Svnle = cell(1,N_ch); -Smlse = cell(1,N_ch); -Sdbt = cell(1,N_ch); - -% Choose your quantile band. For your old style, use 0.04/0.99: -qLow = 0.0; % lower quantile (e.g., 0.04 for old script) -qHigh = 1; % upper quantile (e.g., 0.99 for old script) -cols = linspecer(N_ch); % one color per wavelength (Ch) -cols = cbrewer2('RdBu',N_ch); - -for eval_ptr = 2 - for ch = 1:N_ch - - % Slice 12x50 cell arrays - ffe_cells = reshape(squeeze(res.ffe(ch,:,:,eval_ptr)),length(rop),[]); - vnle_cells = reshape(squeeze(res.vnle(ch,:,:,eval_ptr)),length(rop),[]); - mlse_cells = reshape(squeeze(res.mlse(ch,:,:,eval_ptr)),length(rop),[]); - dbt_cells = reshape(squeeze(res.dbt(ch,:,:,eval_ptr)),length(rop),[]); - - [Sffe{ch}, noX_ffe] = fecCrossings(rop, ffe_cells, fec); - - [Svnle{ch}, noX_ffe] = fecCrossings(rop, vnle_cells, fec); - - [Smlse{ch}, noX_ffe] = fecCrossings(rop, mlse_cells, fec); - - [Sdbt{ch}, noX_ffe] = fecCrossings(rop, dbt_cells, fec); - - % Extract BER matrices using only complete realizations (12/12 ROP filled) - ffe_mat = extractCompleteBER(ffe_cells); % 12 x K_ffe - vnle_mat = extractCompleteBER(vnle_cells); % 12 x K_vnle - mlse_mat = extractCompleteBER(mlse_cells); % 12 x K_mlse - mlse_alpha_mat = extractCompleteAlphas(mlse_cells); % 12 x K_mlse - dbt_mat = extractCompleteBER(dbt_cells); % 12 x K_dbt - - showLegend = 1; % one legend entry per technique - - % Plot shaded band + mean line with boundedline - % plotBandMeanBL(rop, ffe_mat, cols(ch,:), sprintf('FFE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '--s', showLegend); - % scatter(Sffe,fec.*ones(size(Sffe)),20,'v','MarkerFaceColor','black'); - - % plotBandMeanBL(rop, vnle_mat, cols(ch,:), sprintf('VNLE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '--x', showLegend); - - % plotBandMeanBL(rop, mlse_mat, cols(ch,:), sprintf('VNLE+PF+MLSE @ %.1d nm',round(wavelengthplan(ch))), qLow, qHigh, '-o', showLegend); - - % plotBandMeanBL(rop, dbt_mat, cols(ch,:), sprintf('DBt.+MLSE @ %.1d nm',round(wavelengthplan(l))), qLow, qHigh, '--v', showLegend); - - set(gca,'XScale','linear','YScale','log','TickLabelInterpreter','latex','FontSize',11); - yline([3.8e-3, 2.2e-4], 'HandleVisibility','off','LineWidth',1.5); - - end +% If res.eval_dist_km length differs from dims(4), trust dims(4) +if numel(distances) ~= N_distances + % best-effort fallback + distances = (1:N_distances).'; end -ylabel('BER'); -xlabel('ROP'); -title('BER vs. ROP'); -xlim([min(rop) max(rop)]); -ylim([1e-5 0.3]); -grid on; -legend show; +% -------------------- Plot settings -------------------- +fec = 3.8e-3; % choose one; you can switch to 2.2e-4 if needed -%% VIOLIN +qLow = 0.00; +qHigh = 1.00; + +% Color per wavelength/channel +try + cols = ccbrewer2('Set1', N_ch); +catch + cols = linspecer(N_ch); +end + +% Containers for FEC crossings per channel (per eval distance) +Sffe = cell(N_distances, N_ch); +Sdfe = cell(N_distances, N_ch); +Svnle = cell(N_distances, N_ch); +Smlse = cell(N_distances, N_ch); +Sdbt = cell(N_distances, N_ch); + +% -------------------- Choose which eval distances to plot -------------------- +% If you want all: eval_list = 1:N_distances; +% If you want a specific distance (e.g., last): eval_list = N_distances; +eval_list = 1;%1:N_distances; + +% -------------------- BER vs ROP figures (one figure per eval distance) -------------------- +for eval_ptr = eval_list + + figure('Name', sprintf('BER vs ROP @ %s', distLabel(distances, eval_ptr))); + hold on; -S_cell = Sdbt; -S_cell = Smlse; -S_cell = {Svnle,Smlse,Sdbt}; -S_cell = {Svnle}; -figure(5); hold on; -for i = 1:length(S_cell) - % Pad to rectangular matrix: rows = realizations, cols = wavelengths - Kmax = max(cellfun(@numel, S_cell{i})); - S_mat = NaN(Kmax, N_ch); for ch = 1:N_ch - k = numel(S_cell{i}{ch}); - if k > 0 - S_mat(1:k, ch) = S_cell{i}{ch}; + + % ---- Extract cell slices robustly: [N_rop x N_realiz] ---- + ffe_cells = sliceCells4D(res.ffe , ch, eval_ptr, N_rop); + dfe_cells = sliceCells4D(res.dfe , ch, eval_ptr, N_rop); + vnle_cells = sliceCells4D(res.vnle, ch, eval_ptr, N_rop); + mlse_cells = sliceCells4D(res.mlse, ch, eval_ptr, N_rop); + dbt_cells = sliceCells4D(res.dbt , ch, eval_ptr, N_rop); + + % ---- FEC crossings (works with incomplete columns) ---- + [Sffe{eval_ptr,ch}, ~] = fecCrossings(rop, ffe_cells, fec); + [Sdfe{eval_ptr,ch}, ~] = fecCrossings(rop, dfe_cells, fec); + [Svnle{eval_ptr,ch}, ~] = fecCrossings(rop, vnle_cells, fec); + [Smlse{eval_ptr,ch}, ~] = fecCrossings(rop, mlse_cells, fec); + [Sdbt{eval_ptr,ch}, ~] = fecCrossings(rop, dbt_cells, fec); + + % ---- Extract BER matrices using only complete realizations ---- + ffe_mat = extractCompleteBER(ffe_cells); % N_rop x K + dfe_mat = extractCompleteBER(dfe_cells); + vnle_mat = extractCompleteBER(vnle_cells); + mlse_mat = extractCompleteBER(mlse_cells); + dbt_mat = extractCompleteBER(dbt_cells); + + showLegend = (ch == 1); % one legend entry per technique + + % Plot bands + mean (only if non-empty) + if ~isempty(ffe_mat) + plotBandMeanBL(rop, ffe_mat, cols(ch,:), sprintf('FFE @ %dnm', round(wavelengthplan(ch))), qLow, qHigh, '--s', showLegend); end + % if ~isempty(dfe_mat) + % plotBandMeanBL(rop, dfe_mat, cols(ch,:), sprintf('DFE @ %dnm', round(wavelengthplan(ch))), qLow, qHigh, '-o', showLegend); + % end + % if ~isempty(vnle_mat) + % plotBandMeanBL(rop, vnle_mat, cols(ch,:), sprintf('VNLE @ %dnm', round(wavelengthplan(ch))), qLow, qHigh, '--x', showLegend); + % end + % if ~isempty(mlse_mat) + % plotBandMeanBL(rop, mlse_mat, cols(ch,:), sprintf('VNLE+PF+MLSE @ %dnm', round(wavelengthplan(ch))), qLow, qHigh, '-o', showLegend); + % end + % if ~isempty(dbt_mat) + % plotBandMeanBL(rop, dbt_mat, cols(ch,:), sprintf('DBt+MLSE @ %dnm', round(wavelengthplan(ch))), qLow, qHigh, '--v', showLegend); + % end + end - % --- Violin plot over wavelengths (columns) --- - - cols=linspecer(3); - catLabels = arrayfun(@(nm) sprintf('%d nm', nm), wavelengthplan, 'UniformOutput', false); - vs = violinplot(S_mat, catLabels, ... - 'ViolinColor', cols(i,:), ... - 'ViolinAlpha', 0.10, ... - 'MarkerSize', 20, ... - 'ShowMedian', true, ... - 'EdgeColor', cols(i,:), ... - 'ShowWhiskers', false, ... - 'ShowData', true, ... - 'ShowBox', false, ... - 'Bandwidth', 0.05); - - ylim([floor(min(S_mat,[],'all')), ceil(max(S_mat,[],'all'))]) - % ylim([-8 0]); - ylabel('ROP at FEC crossing'); - title(sprintf('RROP to cross BER %.2e', fec)); - grid on; box on; - + set(gca,'XScale','linear','YScale','log','TickLabelInterpreter','latex','FontSize',11); + yline([3.8e-3, 2.2e-4], 'HandleVisibility','off','LineWidth',1.5); + ylabel('BER'); + xlabel('ROP [dB]'); + title(sprintf('BER vs. ROP @ %s', distLabel(distances, eval_ptr))); + xlim([min(rop) max(rop)]); + ylim([1e-5 0.3]); + grid on; + legend show; end -%% ================= helper ================= +%% -------------------- VIOLIN of ROP at FEC crossing -------------------- +% Build technique groups; each group is {N_distances x N_ch} cell entries +% You can pick a distance to show (e.g., eval_ptr=1 for 0 km, or last) +eval_ptr_violin = min(N_distances, max(1, N_distances)); % default: last available +% eval_ptr_violin = 1; + +% Choose which technique(s) to show +techNames = {'VNLE'}; +S_groups = {Svnle(eval_ptr_violin,:)}; + +% Only plot if violinplot exists +if exist('violinplot','file') == 2 + figure('Name', sprintf('FEC crossing violin @ %s', distLabel(distances, eval_ptr_violin))); hold on; + + colsTech = linspecer(numel(S_groups)); + + for i = 1:numel(S_groups) + S_cell = S_groups{i}; % 1 x N_ch cell, each cell is 1xK crossings + + % Pad to rectangular: rows = realizations (max K), cols = wavelengths (N_ch) + Kmax = max(cellfun(@numel, S_cell)); + if isempty(Kmax) || Kmax == 0 + continue; + end + S_mat = zeros(Kmax, N_ch); + for ch = 1:N_ch + k = numel(S_cell{ch}); + if k > 0 + S_mat(1:k, ch) = S_cell{ch}(:); + end + end + + catLabels = arrayfun(@(nm) sprintf('%d nm', nm), round(wavelengthplan), 'UniformOutput', false); + violinplot(S_mat, catLabels, ... + 'ViolinColor', colsTech(i,:), ... + 'ViolinAlpha', 0.10, ... + 'MarkerSize', 20, ... + 'ShowMedian', true, ... + 'EdgeColor', colsTech(i,:), ... + 'ShowWhiskers', false, ... + 'ShowData', true, ... + 'ShowBox', false, ... + 'Bandwidth', 0.05); + + % Add a legend proxy + plot(nan, nan, 'o', 'Color', colsTech(i,:), 'DisplayName', techNames{i}); + end + + ylim([-10,-6]) + ylabel('ROP at FEC crossing'); + title(sprintf('ROP to cross BER %.2e @ %s', fec, distLabel(distances, eval_ptr_violin))); + grid on; box on; + legend show; +else + fprintf('violinplot.m not found on path -> skipping violin plot.\n'); +end + + +%% ================= helpers ================= + +function cells2D = sliceCells4D(C4, ch, eval_ptr, N_rop) +% Returns a [N_rop x N_realiz] cell array (with N_realiz inferred) +% Works even if trailing dimensions are singleton. +% +% If C4 is smaller than expected, missing entries are returned as []. + +dims = size(C4); +if numel(dims) < 4, dims(end+1:4) = 1; end +N_realiz = dims(3); + +cells2D = cell(N_rop, N_realiz); + +% bounds (in case stored arrays are smaller) +chMax = dims(1); +ropMax = dims(2); +rMax = dims(3); +eMax = dims(4); + +if ch > chMax || eval_ptr > eMax + return; % all empty +end + +ropUse = min(N_rop, ropMax); +rUse = min(N_realiz, rMax); + +% Extract and place into a consistent sized cell matrix +tmp = squeeze(C4(ch, 1:ropUse, 1:rUse, eval_ptr)); +% squeeze may return vector/empty if rUse==1 etc. Normalize: +tmp = reshape(tmp, ropUse, rUse); + +cells2D(1:ropUse, 1:rUse) = tmp; +end + +function lbl = distLabel(distances, eval_ptr) +if eval_ptr <= numel(distances) + d = distances(eval_ptr); + if isfinite(d) + lbl = sprintf('%.0f km', d); + return; + end +end +lbl = sprintf('eval\\_%d', eval_ptr); +end + function plotBandMeanBL(x, Y, color, techLabel, qLow, qHigh, lineSpec, showLegend) % Y: (nPoints x nRealizations) -% Remove realizations that are entirely zero (like removeZeros behavior) +if isempty(Y), return; end + +% Remove realizations that are entirely zero badCols = all(Y == 0, 1); Y(:, badCols) = []; +if isempty(Y), return; end + +% Avoid log(0) Y(Y==0) = 1e-8; -% Stats across realizations -mu = mean(Y, 2, 'omitnan'); % mean line -lo = quantile(Y, qLow, 2); % lower bound -hi = quantile(Y, qHigh, 2); % upper bound -% Convert to asymmetric distances required by boundedline: -% b(:,1) = distance to lower side; b(:,2) = distance to upper side -b = [mu - lo, hi - mu]; +mu = mean(Y, 2, 'omitnan'); +lo = quantile(Y, qLow, 2); +hi = quantile(Y, qHigh, 2); +b = [mu - lo, hi - mu]; -% Call boundedline with alpha shading -[hl, hp] = boundedline(x(:), mu(:), b, lineSpec, 'alpha', ... - 'transparency', 0.18); -% Color styling -set(hl, 'Color', color, 'LineWidth', 1.4, 'MarkerSize', 4); -set(hp, 'FaceColor', color, 'HandleVisibility','off'); % patch hidden in legend +% boundedline optional; fallback to simple plot if missing +if 0%exist('boundedline','file') == 2 + [hl, hp] = boundedline(x(:), mu(:), b, lineSpec, 'alpha', 'transparency', 0.18); + set(hl, 'Color', color, 'LineWidth', 1.4, 'MarkerSize', 4); + set(hp, 'FaceColor', color, 'HandleVisibility','off'); +else + hl = plot(x(:), mu(:), lineSpec, 'LineWidth', 1.4, 'MarkerSize', 4); + set(hl, 'Color', color); +end -% Single legend entry per technique (use first wavelength only) if showLegend set(hl, 'DisplayName', techLabel); else set(hl, 'HandleVisibility','off'); end -% Optional: outline the bounds if outlinebounds is available -if exist('outlinebounds','file') == 2 +if 0%exist('outlinebounds','file') == 2 && exist('boundedline','file') == 2 ho = outlinebounds(hl, hp); - set(ho, 'linestyle', ':', 'color', color, 'linewidth', 1, ... - 'HandleVisibility','off'); + set(ho, 'linestyle', ':', 'color', color, 'linewidth', 1, 'HandleVisibility','off'); end end -function [S, noCrossingMask, Y_keep] = fecCrossings(rop, cells12xR, fec) -% cells12xR: 12xR cell array (one wavelength + scheme slice) -% each cell must be a struct with .metrics.BER -% rop: 12x1 numeric vector of ROP points -% fec: scalar FEC threshold (e.g., 3.8e-3) -% -% Outputs: -% S 1xK vector of crossing ROP per kept realization (NaN if none) -% noCrossingMask 1xK logical mask: true if no crossing for that realization -% Y_keep 12xK numeric BER matrix used for the crossing detection - -% 1) keep only complete realization columns -Y = extractCompleteBER(cells12xR); % -> 12 x K +function [S, noCrossingMask] = fecCrossings(rop, cellsNxR, fec) +% cellsNxR: N_rop x R cell array; each cell is struct with .metrics.BER +Y = extractCompleteBER(cellsNxR); % -> N_rop x K if isempty(Y) S = []; noCrossingMask = []; - Y_keep = Y; return; end -% 2) optionally drop realizations with mean BER > 0.1 ok = mean(Y,1,'omitnan') <= 0.1; Y = Y(:, ok); if isempty(Y) S = []; noCrossingMask = []; - Y_keep = Y; return; end -% 3) find crossings per realization nR = size(Y,2); S = nan(1,nR); noCrossingMask = true(1,nR); -rop = rop(:); % ensure column +rop = rop(:); for j = 1:nR y = Y(:,j); - - % sign change from >fec to <=fec (first time it drops below FEC) above = (y > fec); idx = find(above(1:end-1) & ~above(2:end), 1, 'first'); - if ~isempty(idx) - % linear interpolation between (x1,y1) and (x2,y2) x1 = rop(idx); y1 = y(idx); x2 = rop(idx+1); y2 = y(idx+1); - if isfinite(y1) && isfinite(y2) && y2 ~= y1 t = (fec - y1) / (y2 - y1); S(j) = x1 + t*(x2 - x1); @@ -209,37 +292,34 @@ for j = 1:nR end end end - -Y_keep = Y; end - - function Y = extractCompleteBER(cellSlice) -% cellSlice: 12xR cell array; each cell should be a struct with .metrics.BER -% Keep only those realization columns where ALL 12 ROP entries are valid. +% Keep only those realization columns where ALL N_rop entries are non-empty if isempty(cellSlice), Y = []; return; end nR = size(cellSlice,2); keep = false(1,nR); + for r = 1:nR col = cellSlice(:,r); - keep(r) = all(cellfun(@(c) ~isempty(c) , col)); + keep(r) = all(cellfun(@(c) ~isempty(c), col)); end + if ~any(keep), Y = []; return; end + Y = cellfun(@(c) c.metrics.BER, cellSlice(:,keep), 'UniformOutput', true); end function Y = extractCompleteAlphas(cellSlice) -% cellSlice: 12xR cell array; each cell should be a struct with .metrics.BER -% Keep only those realization columns where ALL 12 ROP entries are valid. if isempty(cellSlice), Y = []; return; end nR = size(cellSlice,2); keep = false(1,nR); + for r = 1:nR col = cellSlice(:,r); - keep(r) = all(cellfun(@(c) ~isempty(c) , col)); + keep(r) = all(cellfun(@(c) ~isempty(c), col)); end + if ~any(keep), Y = []; return; end Y = cellfun(@(c) c.metrics.Alpha, cellSlice(:,keep), 'UniformOutput', true); end - diff --git a/projects/WDM/WDM_model_10km.m b/projects/WDM/WDM_model_10km.m index d8977eb..779365c 100644 --- a/projects/WDM/WDM_model_10km.m +++ b/projects/WDM/WDM_model_10km.m @@ -67,8 +67,8 @@ fname = sprintf('WDM_%s_%s_%s_%dkm_%dch_%dghz_%s.mat', char(t), host, jobid, opt 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.1; -s.gamma = 0.0023; +s.pmd = 0;%0.1; +s.gamma = 0;%0.0023; s.M = 4; fsym = 112e9; @@ -108,7 +108,7 @@ Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",1 N = numel(s.wavelengthplan); f_plan = physconst('lightspeed')./(s.wavelengthplan.*1e-9); -margin = 20e12; % some THz left and right +margin = 2e12; % some THz left and right f_span = (max(f_plan)+margin)-(min(f_plan)-margin); f_nyq = f_span/2; @@ -124,13 +124,13 @@ signal_cell = {}; Symbols = {}; Tx_bits = {}; -s.rop = -12:0.75:-0.75; +s.rop = -12;%:0.75:-0; %% ---------- Intermediate evaluation points (distance dimension) ---------- % Evaluate BER at these intermediate distances (km), plus always include the final link_length if > 0. segment_length = 1; % km (must match fiber loop below) -eval_dist_km = [2 10]; +eval_dist_km = [10]; eval_dist_km = eval_dist_km(eval_dist_km <= link_length); % Always include final distance (if > 0) and avoid duplicates @@ -209,7 +209,7 @@ for realiz = 1:s.num_realiz % 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).process(El_sig); + "linewidth",laser_linewidth,"randomkey",s.random_key+l+realiz,"alpha",0.8).process(El_sig); signal_cell{l} = Polarization_Controller("mode","rot_power","desired_power",pol_rot(l)).process(Eml_out); end diff --git a/projects/WDM/WDM_model_10km_queue.m b/projects/WDM/WDM_model_10km_queue.m index 7e1832d..24b7e4e 100644 --- a/projects/WDM/WDM_model_10km_queue.m +++ b/projects/WDM/WDM_model_10km_queue.m @@ -1,18 +1,14 @@ function WDM_model_10km_queue(options) -%WDM_model_10km_queue -% Job-queue / pipeline version of your script: -% - Main thread runs fiber propagation segment-by-segment -% - At each eval distance: demux once and enqueue per-(l,ri) DSP jobs via parfeval -% - While continuing propagation: periodically harvest completed jobs (non-blocking) -% - End of realization: drain remaining jobs, then save -% -% Save this as a NEW file, e.g. "WDM_model_10km_queue.m" +%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 = 2; + options.fiber_length_km = 10; options.rand_key = 1; options.num_realiz = 1; options.fwm_mitigation_technique = "co"; @@ -65,12 +61,16 @@ 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.1; -s.gamma = 0.0023; +s.pmd = 0;%0.1; +s.gamma = 0;%0.0023; s.M = 4; fsym = 112e9; @@ -105,7 +105,7 @@ Pform = Pulseformer("fsym",fsym,"fdac",4*fsym,"pulse","rc","pulselength",1 N = numel(s.wavelengthplan); f_plan = physconst('lightspeed')./(s.wavelengthplan.*1e-9); -margin = 25e12; % some THz left and right +margin = 2e12; % some THz left and right f_span = (max(f_plan)+margin)-(min(f_plan)-margin); f_nyq = f_span/2; @@ -117,11 +117,11 @@ 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 = [2 10]; +eval_dist_km = [0]; eval_dist_km = eval_dist_km(eval_dist_km <= link_length); if link_length > 0 @@ -167,7 +167,7 @@ end for realiz = 1:s.num_realiz - % Per-realization TX storage + % 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); @@ -181,42 +181,62 @@ for realiz = 1:s.num_realiz for l = 1:N [Digi_sig,Symbols{l},Tx_bits{l}] = PAMsource( ... - "fsym",fsym,"M",s.M,"order",18,"useprbs",0, ... + "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(); + ).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).process(El_sig); + "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); @@ -225,20 +245,70 @@ for realiz = 1:s.num_realiz zdw = 1310; randomize_D = true; - Dvec = getDispersionVector(nSegments, d_local, zdw, randomize_D, s.random_key+realiz); + + % 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); - if 1 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); - end % -------- Evaluate at intermediate distance (enqueue jobs) -------- if eval_ptr <= nEval && seg == eval_seg(eval_ptr) @@ -246,39 +316,52 @@ for realiz = 1:s.num_realiz 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 - % Enqueue one job per (l,ri) + + 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}, ... % channel signal at this eval point - s.rop(ri), ... % rop [dB] - Symbols{l}, Tx_bits{l}, ...% reference + bits - s, l, realiz, ... % settings + indices + 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); + 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: store finished jobs (if any) while we continue propagation + % 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 - % Drain remaining jobs for this realization (blocking until queue is empty) 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; @@ -292,6 +375,9 @@ for realiz = 1:s.num_realiz 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 @@ -300,25 +386,39 @@ 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, realiz, ... - fdac, kover, fadc, fsym, len_tr, mu_dc, mu_ffe, mu_dfe, dfe_order, duob_mode) + 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, ... @@ -330,11 +430,28 @@ Scpe_sig = Scope("fsimu",fdac*kover,"fadc",fadc, ... "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], ... @@ -345,10 +462,11 @@ eq_ffe = EQ("Ne",ffe_order,"Nb",[0,0,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 -------------------- -dfe_order2 = [50, 0, 0]; -eq_dfe = EQ("Ne",dfe_order2,"Nb",[2,0,0], ... +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); @@ -356,6 +474,8 @@ eq_dfe = EQ("Ne",dfe_order2,"Nb",[2,0,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; @@ -371,11 +491,22 @@ mlse_ = MLSE("duobinary_output",0,'M',s.M,'trellis_states',PAMmapper(s.M,0).leve [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); - -dbt_results = duobinary_target(eq_v, mlse_db_, s.M, Rx_sig, Symbols_l, Tx_bits_l, ... +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 @@ -385,15 +516,12 @@ function [F, meta, output_ffe, output_dfe, output_vnle, output_mlse, output_dbt] if nargin < 9, block = false; end if isempty(F), return; end -% Non-blocking mode: try to fetch as many finished as possible with timeout=0 -% Blocking mode: wait for at least one result (timeout=Inf), then keep draining with timeout=0 if block timeout_first = Inf; else timeout_first = 0; end -% Fetch first (maybe blocking), then drain the rest non-blocking first = true; while ~isempty(F) if first @@ -406,12 +534,11 @@ while ~isempty(F) try [k, ffe_r, dfe_r, vnle_r, mlse_r, dbt_r] = fetchNext(F, timeout); catch - % fetchNext throws if timeout=0 and nothing finished in some releases; treat as "nothing ready" break; end if isempty(k) - break; % nothing ready + break; end m = meta(k); @@ -422,8 +549,177 @@ while ~isempty(F) 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; - % Remove completed future & meta entry (important: indices change afterwards) 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 + 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