From 8c4edf2490344a97edbeebfba41510128975674e Mon Sep 17 00:00:00 2001 From: Silas Oettinghaus Date: Tue, 14 Jul 2026 11:50:22 +0200 Subject: [PATCH] Final MPI Calculations here --- .../Equalizer/MPI_reduction/FFE_DCTracking.m | 235 +++++++++++++- Functions/EQ_recipes/mpi_recipe_dev.m | 84 ++++- Functions/Job_Processing/dsp_runid.m | 15 +- .../PLOT_mpi_reduction_db_ber_vs_sir.m | 254 +++++++++++++++ .../db/appendMpiReductionDspOutput.m | 99 ++++++ .../db/appendMpiReductionWarehouse.m | 307 ++++++++++++++++++ .../MPI_revisit/investigate_mpi_algorithms.m | 21 +- workerError.mat | Bin 1173 -> 1193 bytes 8 files changed, 1003 insertions(+), 12 deletions(-) create mode 100644 projects/Diss/MPI_revisit/algorithms/PLOT_mpi_reduction_db_ber_vs_sir.m create mode 100644 projects/Diss/MPI_revisit/db/appendMpiReductionDspOutput.m create mode 100644 projects/Diss/MPI_revisit/db/appendMpiReductionWarehouse.m diff --git a/Classes/04_DSP/Equalizer/MPI_reduction/FFE_DCTracking.m b/Classes/04_DSP/Equalizer/MPI_reduction/FFE_DCTracking.m index b730b54..24e2bde 100644 --- a/Classes/04_DSP/Equalizer/MPI_reduction/FFE_DCTracking.m +++ b/Classes/04_DSP/Equalizer/MPI_reduction/FFE_DCTracking.m @@ -11,14 +11,22 @@ classdef FFE_DCTracking < FFE_plain dc_tracking_persistence_gain dc_tracking_power_exponent dc_tracking_buffer_len + dc_tracking_mu_eff_max e_dc + + optimize_dc_tracking_params = 0 + dc_tracking_optimization + dc_tracking_optimization_iter = 0 + dc_tracking_optimization_len + dc_tracking_optimization_max_evals + dc_tracking_optimization_delay_weight + dc_tracking_optimization_smoothing_len end properties (Access = private) dc_tracking_mu_min = 1e-6 dc_tracking_mu_max = 3e-1 dc_tracking_mu_eff_min = 0 - dc_tracking_mu_eff_max = inf end methods @@ -42,6 +50,7 @@ classdef FFE_DCTracking < FFE_plain options.dc_tracking_persistence_gain = 0 options.dc_tracking_power_exponent = 2 options.dc_tracking_buffer_len = 1 + options.dc_tracking_mu_eff_max = inf options.decide = false @@ -50,6 +59,11 @@ classdef FFE_DCTracking < FFE_plain options.mu_optimization_len = 2^15 options.plot_mu_optimization = 0 options.mu_optimization_fignum = 3010 + options.optimize_dc_tracking_params = 0 + options.dc_tracking_optimization_len = 2^15 + options.dc_tracking_optimization_max_evals = 20 + options.dc_tracking_optimization_delay_weight = 1e-3 + options.dc_tracking_optimization_smoothing_len = 501 end obj@FFE_plain( ... @@ -75,8 +89,18 @@ classdef FFE_DCTracking < FFE_plain obj.dc_tracking_persistence_gain = options.dc_tracking_persistence_gain; obj.dc_tracking_power_exponent = options.dc_tracking_power_exponent; obj.dc_tracking_buffer_len = floor(options.dc_tracking_buffer_len); + obj.dc_tracking_mu_eff_max = options.dc_tracking_mu_eff_max; obj.e_dc = 0; + obj.optimize_dc_tracking_params = options.optimize_dc_tracking_params; + obj.dc_tracking_optimization_len = options.dc_tracking_optimization_len; + obj.dc_tracking_optimization_max_evals = max(1, ... + floor(options.dc_tracking_optimization_max_evals)); + obj.dc_tracking_optimization_delay_weight = ... + max(0,options.dc_tracking_optimization_delay_weight); + obj.dc_tracking_optimization_smoothing_len = max(1, ... + floor(options.dc_tracking_optimization_smoothing_len)); + assert(obj.dc_tracking_buffer_len >= 0); end @@ -90,6 +114,11 @@ classdef FFE_DCTracking < FFE_plain obj.resetTrackingState(); end + if obj.optimize_dc_tracking_params + obj.optimizeDcTrackingParams(X.signal,D.signal); + obj.resetTrackingState(); + end + training = true; showviz = false; obj.equalize(X.signal,D.signal,obj.mu_tr,obj.epochs_tr,obj.len_tr,training,showviz); @@ -134,7 +163,8 @@ classdef FFE_DCTracking < FFE_plain if obj.dd_mode vars = [vars, optimizableVariable("mu_dd",mu_range,"Transform","log")]; end - optimize_dc_tracking_mu = obj.dc_tracking_mu ~= 0; + optimize_dc_tracking_mu = obj.dc_tracking_mu ~= 0 && ... + ~obj.optimize_dc_tracking_params; if optimize_dc_tracking_mu vars = [vars, optimizableVariable("dc_tracking_mu",[1e-5,1e-1],"Transform","log")]; end @@ -175,6 +205,61 @@ classdef FFE_DCTracking < FFE_plain end end + function optimizeDcTrackingParams(obj,x,d) + [x_opt,d_opt,N_opt] = obj.optimizationSignals( ... + x,d,obj.dc_tracking_optimization_len); + + vars = optimizableVariable( ... + "dc_tracking_mu",[1e-5,1e-1],"Transform","log"); + initial_values = min(max(obj.dc_tracking_mu,1e-5),1e-1); + initial_names = "dc_tracking_mu"; + if obj.dc_tracking_adaptive_enabled + vars = [vars, ... + optimizableVariable("dc_tracking_persistence_gain",[0,2]), ... + optimizableVariable("dc_tracking_mu_eff_max",[1e-3,3e-1], ... + "Transform","log")]; + initial_values = [initial_values, ... + min(max(obj.dc_tracking_persistence_gain,0),2), ... + min(max(obj.dc_tracking_mu_eff_max,1e-3),3e-1)]; + initial_names = [initial_names, ... + "dc_tracking_persistence_gain", "dc_tracking_mu_eff_max"]; + end + initial_x = array2table(initial_values, ... + "VariableNames",cellstr(initial_names)); + + obj.dc_tracking_optimization_iter = 0; + fprintf("FFE_DCTracking DC opt uses fixed mu_tr=%9.3e, mu_dd=%9.3e on %d samples / %d symbols\n", ... + obj.mu_tr,obj.mu_dd,N_opt,numel(d_opt)); + + old_rng = rng; + cleanup_rng = onCleanup(@()rng(old_rng)); + rng(42,"twister"); + obj.dc_tracking_optimization = bayesopt( ... + @(p)obj.dcTrackingObjective(p,x_opt,d_opt),vars, ... + "MaxObjectiveEvaluations",obj.dc_tracking_optimization_max_evals, ... + "InitialX",initial_x, ... + "AcquisitionFunctionName","expected-improvement-plus", ... + "IsObjectiveDeterministic",true, ... + "Verbose",0, ... + "PlotFcn",[]); + clear cleanup_rng + + best = obj.dc_tracking_optimization.XAtMinObjective; + obj.dc_tracking_mu = best.dc_tracking_mu; + if obj.dc_tracking_adaptive_enabled + obj.dc_tracking_persistence_gain = ... + best.dc_tracking_persistence_gain; + obj.dc_tracking_mu_eff_max = best.dc_tracking_mu_eff_max; + fprintf("\nFFE_DCTracking DC opt done: dc_tracking_mu=%9.3e, persistence_gain=%6.3f, mu_eff_max=%9.3e, objective=%9.3e\n", ... + obj.dc_tracking_mu,obj.dc_tracking_persistence_gain, ... + obj.dc_tracking_mu_eff_max, ... + obj.dc_tracking_optimization.MinObjective); + else + fprintf("\nFFE_DCTracking DC opt done: dc_tracking_mu=%9.3e, objective=%9.3e\n", ... + obj.dc_tracking_mu,obj.dc_tracking_optimization.MinObjective); + end + end + function objective = muObjective(obj,params,x,d) old_state = obj.captureTrackingObjectiveState(); cleanup = onCleanup(@()obj.restoreTrackingObjectiveState(old_state)); @@ -218,6 +303,49 @@ classdef FFE_DCTracking < FFE_plain clear cleanup end + function objective = dcTrackingObjective(obj,params,x,d) + old_state = obj.captureTrackingObjectiveState(); + cleanup = onCleanup(@()obj.restoreTrackingObjectiveState(old_state)); + + obj.applyDcTrackingObjectiveParams(params); + obj.save_debug = 1; + obj.resetTrackingState(); + + N_tr = min(obj.len_tr,numel(x)); + obj.equalize(x,d,obj.mu_tr,obj.epochs_tr,N_tr,true,false); + if obj.dd_mode + [signal,~] = obj.equalize( ... + x,d,obj.mu_dd,obj.epochs_dd,numel(x),false,false); + else + [signal,~] = obj.equalize(x,d,0,1,numel(x),false,false); + end + + [ber,errors] = obj.berObjective(signal,d); + [delay_symbols,delay_corr] = obj.dcTrackingDelayObjective(signal,d); + delay_penalty = obj.dc_tracking_optimization_delay_weight * ... + abs(delay_symbols) / max(numel(d),1); + objective = ber + delay_penalty; + if ~isfinite(objective) + objective = inf; + end + + obj.dc_tracking_optimization_iter = ... + obj.dc_tracking_optimization_iter + 1; + if obj.dc_tracking_adaptive_enabled + fprintf("\rFFE_DCTracking DC opt %02d: dc_tracking_mu=%9.3e, persistence_gain=%6.3f, mu_eff_max=%9.3e, BER=%9.3e, delay=%7.0f, corr=%6.3f, obj=%9.3e, errors=%d", ... + obj.dc_tracking_optimization_iter,params.dc_tracking_mu, ... + params.dc_tracking_persistence_gain, ... + params.dc_tracking_mu_eff_max,ber,delay_symbols,delay_corr, ... + objective,errors); + else + fprintf("\rFFE_DCTracking DC opt %02d: dc_tracking_mu=%9.3e, BER=%9.3e, delay=%7.0f, corr=%6.3f, obj=%9.3e, errors=%d", ... + obj.dc_tracking_optimization_iter,params.dc_tracking_mu, ... + ber,delay_symbols,delay_corr,objective,errors); + end + + clear cleanup + end + function [y,d_hat] = equalize(obj,x,d,mu,epochs,N,training,showviz) arguments obj @@ -392,6 +520,9 @@ classdef FFE_DCTracking < FFE_plain state.save_debug = obj.save_debug; state.debug_struct = obj.debug_struct; state.dc_tracking_mu = obj.dc_tracking_mu; + state.dc_tracking_persistence_gain = ... + obj.dc_tracking_persistence_gain; + state.dc_tracking_mu_eff_max = obj.dc_tracking_mu_eff_max; end function restoreTrackingObjectiveState(obj,state) @@ -403,6 +534,106 @@ classdef FFE_DCTracking < FFE_plain obj.save_debug = state.save_debug; obj.debug_struct = state.debug_struct; obj.dc_tracking_mu = state.dc_tracking_mu; + obj.dc_tracking_persistence_gain = ... + state.dc_tracking_persistence_gain; + obj.dc_tracking_mu_eff_max = state.dc_tracking_mu_eff_max; + end + + function applyDcTrackingObjectiveParams(obj,params) + var_names = string(params.Properties.VariableNames); + if any(var_names == "dc_tracking_mu") + obj.dc_tracking_mu = params.dc_tracking_mu; + end + if any(var_names == "dc_tracking_persistence_gain") + obj.dc_tracking_persistence_gain = ... + params.dc_tracking_persistence_gain; + end + if any(var_names == "dc_tracking_mu_eff_max") + obj.dc_tracking_mu_eff_max = params.dc_tracking_mu_eff_max; + end + end + + function [delay_symbols,delay_corr] = ... + dcTrackingDelayObjective(obj,signal,d) + delay_symbols = 0; + delay_corr = 0; + if ~isfield(obj.debug_struct,"dc_tracking_est") || ... + isempty(obj.debug_struct.dc_tracking_est) + return + end + + smooth_len = obj.dc_tracking_optimization_smoothing_len; + dc_tracking_est_s = movmean( ... + obj.debug_struct.dc_tracking_est(:),smooth_len,"omitnan"); + avg_lvl_dc = obj.averageLevelTrace(signal,d,smooth_len); + + xcorr_len = min(numel(dc_tracking_est_s),numel(avg_lvl_dc)); + if xcorr_len < 2 + return + end + + inv_dc_xcorr = -dc_tracking_est_s(1:xcorr_len); + avg_lvl_xcorr = avg_lvl_dc(1:xcorr_len); + inv_dc_xcorr = fillmissing( ... + inv_dc_xcorr,"linear","EndValues","nearest"); + avg_lvl_xcorr = fillmissing( ... + avg_lvl_xcorr,"linear","EndValues","nearest"); + inv_dc_xcorr = inv_dc_xcorr - mean(inv_dc_xcorr,"omitnan"); + avg_lvl_xcorr = avg_lvl_xcorr - mean(avg_lvl_xcorr,"omitnan"); + + if rms(inv_dc_xcorr) <= eps || rms(avg_lvl_xcorr) <= eps + return + end + + [dc_level_xcorr,dc_level_lags] = xcorr( ... + inv_dc_xcorr,avg_lvl_xcorr,"coeff"); + [delay_corr,delay_idx] = max(dc_level_xcorr); + delay_symbols = dc_level_lags(delay_idx); + if ~isfinite(delay_corr) + delay_corr = 0; + delay_symbols = 0; + end + end + + function avg_lvl_dc = averageLevelTrace(obj,signal,d,smooth_len) + signal = signal(:); + d = d(:); + n_symbols = min(numel(signal),numel(d)); + signal = signal(1:n_symbols); + d = d(1:n_symbols); + levels = unique(d); + avg_for_lvl = NaN(numel(levels),n_symbols); + + for level_idx = 1:numel(levels) + level_mask = d == levels(level_idx); + level_samples = signal(level_mask); + if isempty(level_samples) + continue + end + + smooth_window = min(smooth_len,numel(level_samples)); + avg_for_lvl(level_idx,level_mask) = movmean( ... + level_samples,smooth_window,"omitnan","Endpoints","shrink"); + avg_for_lvl(level_idx,:) = obj.interpolateMissingAverage( ... + avg_for_lvl(level_idx,:)); + end + + avg_lvl_dc = mean(avg_for_lvl,1,"omitnan").'; + end + + function level_average = interpolateMissingAverage(~,level_average) + valid_samples = isfinite(level_average); + if nnz(valid_samples) == 0 + return + elseif nnz(valid_samples) == 1 + level_average(:) = level_average(valid_samples); + return + end + + t = 1:numel(level_average); + level_average(~valid_samples) = interp1( ... + t(valid_samples),level_average(valid_samples), ... + t(~valid_samples),"linear","extrap"); end function N = validSampleLength(obj,N,x,d) diff --git a/Functions/EQ_recipes/mpi_recipe_dev.m b/Functions/EQ_recipes/mpi_recipe_dev.m index 8f4022a..bd7413a 100644 --- a/Functions/EQ_recipes/mpi_recipe_dev.m +++ b/Functions/EQ_recipes/mpi_recipe_dev.m @@ -66,9 +66,12 @@ if run_conventional_ffe "plot_mu_optimization", options.debug_plots, ... "save_debug", eq_save_debug); + storageName = "conventional_ffe"; [ffe_results,equalized_signal] = runFfe(eq_ffe, "Conventional FFE", ... Scpe_sig, Symbols, Tx_bits, options); - output.conventional_ffe = ffe_results; + ffe_results = attachMpiReductionConfig(ffe_results, eq_ffe, storageName, ... + "plain_ffe", "baseline", block_update); + output.(char(storageName)) = ffe_results; if plot_output_signals plotEqSignals(equalized_signal,Symbols,options,400,-1); @@ -85,7 +88,7 @@ if run_a2_tracked_levels "len_tr", eq_len_tr, ... "epochs_tr", eq_epochs_tr, ... "mu_tr", eq_mu_tr, ... - "dd_mode", false, ... + "dd_mode", true, ... "epochs_dd", eq_epochs_dd, ... "mu_dd", 0.012, ... "dc_smoothing_a2", 1, ... @@ -94,9 +97,12 @@ if run_a2_tracked_levels "dc_level_weights_a2", [1], ... "save_debug", eq_save_debug); + storageName = "a2_adaptive_levels"; [ffe_results,equalized_signal] = runFfe(eq_ffe, "A2 tracked levels", ... Scpe_sig, Symbols, Tx_bits, options); - output.a2_adaptive_levels = ffe_results; + ffe_results = attachMpiReductionConfig(ffe_results, eq_ffe, storageName, ... + "a2_tracked_levels", "tracked_levels", block_update); + output.(char(storageName)) = ffe_results; if plot_output_signals plotEqSignals(equalized_signal,Symbols,options,410,-1); @@ -122,9 +128,12 @@ if run_a2_residual "dc_level_weights_a2", [0.6], ... "save_debug", eq_save_debug, "optmize_mus",0); + storageName = "a2_residual"; [ffe_results,equalized_signal] = runFfe(eq_ffe, "A2 residual", ... Scpe_sig, Symbols, Tx_bits, options); - output.a2_residual = ffe_results; + ffe_results = attachMpiReductionConfig(ffe_results, eq_ffe, storageName, ... + "a2_residual", "residual_correction", block_update); + output.(char(storageName)) = ffe_results; if plot_output_signals plotEqSignals(equalized_signal,Symbols,options,420,-1); @@ -149,9 +158,12 @@ if run_a1 "dc_avg_update_blocklength_a1", block_update, ... "save_debug", eq_save_debug); + storageName = "a1_ff_dc_avg"; [ffe_results,equalized_signal] = runFfe(eq_ffe, "A1", ... Scpe_sig, Symbols, Tx_bits, options); - output.a1_ff_dc_avg = ffe_results; + ffe_results = attachMpiReductionConfig(ffe_results, eq_ffe, storageName, ... + "a1_moving_average", "ff_dc_avg", block_update); + output.(char(storageName)) = ffe_results; if plot_output_signals plotEqSignals(equalized_signal,Symbols,options,430,-1); @@ -176,12 +188,18 @@ if run_tracking_adaptive "dc_tracking_persistence_gain", 0, ... "dc_tracking_buffer_len", block_update, ... "optmize_mus", false, ... + "optimize_dc_tracking_params", true, ... + "dc_tracking_optimization_len", 2^15, ... + "dc_tracking_optimization_max_evals", 20, ... "plot_mu_optimization", options.debug_plots, ... "save_debug", eq_save_debug); + storageName = "dc_tracking"; [ffe_results,equalized_signal] = runFfe(eq_ffe, "DC tracking adaptive", ... Scpe_sig, Symbols, Tx_bits, options); - output.dc_tracking = ffe_results; + ffe_results = attachMpiReductionConfig(ffe_results, eq_ffe, storageName, ... + "dc_tracking", "adaptive", block_update); + output.(char(storageName)) = ffe_results; if plot_output_signals plotEqSignals(equalized_signal,Symbols,options,450,-1); @@ -200,6 +218,60 @@ function [ffe_results,equalized_signal] = runFfe(eq_ffe,description,Scpe_sig,Sym ffe_results.metrics.print("description",resultDescription(description,options)); end +function ffe_results = attachMpiReductionConfig(ffe_results, eq_ffe, storageName, ... + algorithm, algorithmVariant, block_update) + ffe_results.mpi_reduction_config = struct( ... + "storage_name", char(storageName), ... + "algorithm", char(algorithm), ... + "algorithm_variant", char(algorithmVariant), ... + "eq_class", class(eq_ffe), ... + "params", collectMpiReductionParams(eq_ffe, block_update)); +end + +function params = collectMpiReductionParams(eq_ffe, block_update) + params = struct(); + params.block_update = block_update; + + whitelistedProps = [ ... + "sps", ... + "order", ... + "decide", ... + "adaption_technique", ... + "len_tr", ... + "epochs_tr", ... + "mu_tr", ... + "dd_mode", ... + "epochs_dd", ... + "mu_dd", ... + "optmize_mus", ... + "plot_mu_optimization", ... + "save_debug", ... + "dc_smoothing_a1", ... + "dc_avg_bufferlength_a1", ... + "dc_avg_update_blocklength_a1", ... + "dc_smoothing_a2", ... + "dc_level_avg_bufferlength_a2", ... + "dc_level_update_blocklength_a2", ... + "dc_level_weights_a2", ... + "dc_tracking_mu", ... + "dc_tracking_adaptive_enabled", ... + "dc_tracking_persistence_gain", ... + "dc_tracking_buffer_len", ... + "dc_tracking_mu_eff_max", ... + "optimize_dc_tracking_params", ... + "dc_tracking_optimization_len", ... + "dc_tracking_optimization_max_evals", ... + "dc_tracking_optimization_delay_weight", ... + "dc_tracking_optimization_smoothing_len"]; + + for propIdx = 1:numel(whitelistedProps) + propName = char(whitelistedProps(propIdx)); + if isprop(eq_ffe, propName) + params.(propName) = eq_ffe.(propName); + end + end +end + function description = resultDescription(prefix,options) sir = options.dataTable.sir; if numel(sir) > 1 diff --git a/Functions/Job_Processing/dsp_runid.m b/Functions/Job_Processing/dsp_runid.m index 498eaed..e101297 100644 --- a/Functions/Job_Processing/dsp_runid.m +++ b/Functions/Job_Processing/dsp_runid.m @@ -3,6 +3,9 @@ function output = dsp_runid(run_id, options) arguments run_id options.append_to_db = 0; + options.append_mpi_reduction_db (1,1) logical = false; + options.mpi_reduction_study_name string = "mpi_reduction_v1"; + options.mpi_reduction_writer_path string = ""; options.max_occurences = 4; options.start_occurence = 1; options.userParameters = struct(); @@ -24,7 +27,11 @@ try database = []; inputSource = normalizeDspInputSource(options.mode); - if inputSource == "run_id" || options.append_to_db + if options.append_mpi_reduction_db && strlength(options.mpi_reduction_writer_path) > 0 + addpath(options.mpi_reduction_writer_path); + end + + if inputSource == "run_id" || options.append_to_db || options.append_mpi_reduction_db database = DBHandler("dataBase", [options.dataBase], "type", options.database_type, ... "user", options.user, "password", options.password, ... "server", options.server, "port", options.port); @@ -60,6 +67,12 @@ try if options.append_to_db appendDspOutputToDatabase(database, run_id, dspOutput); end + if options.append_mpi_reduction_db + occurrence_idx = options.start_occurence + r - 1; + appendMpiReductionDspOutput(database, run_id, occurrence_idx, dspOutput, options, ... + "study_name", options.mpi_reduction_study_name, ... + "verbose", false); + end end catch ME diff --git a/projects/Diss/MPI_revisit/algorithms/PLOT_mpi_reduction_db_ber_vs_sir.m b/projects/Diss/MPI_revisit/algorithms/PLOT_mpi_reduction_db_ber_vs_sir.m new file mode 100644 index 0000000..81daa33 --- /dev/null +++ b/projects/Diss/MPI_revisit/algorithms/PLOT_mpi_reduction_db_ber_vs_sir.m @@ -0,0 +1,254 @@ +%% BER over SIR from MpiReductionResults +% 1) gather data from the database +% 2) clean data and remove per-curve outliers +% 3) plot BER over SIR grouped by algorithm + +clear; clc; + +%% 1) Gather data + +studyName = "pam4_112_greater_3153"; +pathLengthToPlot = 0; +selectedBlockUpdate = 1; % set [] to pool all block_update values +selectedPamLevels = []; % set [] to use all PAM levels in the query result +selectedAlgorithms = []; % set [] to use all algorithms in the query result + +useBoundedLines = true; +usePolyfit = true; +polyfitOrderMax = 4; +maxBerForPlot = 0.1; + +db = DBHandler( ... + "dataBase", "labor", ... + "type", "mysql", ... + "server", "192.168.178.192", ... + "user", "silas", ... + "password", "silas"); +db.refresh(); + +fp = QueryFilter(); +fp.where('Runs', 'interference_path_length', 'EQUALS', pathLengthToPlot); +fp.where('Runs', 'fiber_length', 'EQUALS', 0); +fp.where('Runs', 'db_mode', 'EQUALS', '"no_db"'); +fp.where('MpiReductionResults', 'study_name', 'EQUALS', char(studyName)); +if ~isempty(selectedBlockUpdate) + fp.where('MpiReductionResults', 'block_update', 'EQUALS', selectedBlockUpdate); +end + +selectedFields = { ... + 'Runs.run_id', ... + 'Runs.sir', ... + 'Runs.pam_level', ... + 'Runs.symbolrate', ... + 'Runs.interference_path_length', ... + 'Runs.power_mpi_interference', ... + 'MpiReductionResults.occurrence_idx', ... + 'MpiReductionResults.storage_name', ... + 'MpiReductionResults.algorithm', ... + 'MpiReductionResults.algorithm_variant', ... + 'MpiReductionResults.eq_class', ... + 'MpiReductionResults.block_update', ... + 'MpiReductionResults.BER', ... + 'MpiReductionResults.BER_precoded', ... + 'MpiReductionResults.SNR', ... + 'MpiReductionResults.GMI', ... + 'MpiReductionResults.AIR'}; +selectedFields = selectedFields(:); + +[rawData, query] = db.queryDB(fp, selectedFields); +disp(query); +fprintf("Fetched %d MPI reduction rows.\n", height(rawData)); + +%% 2) Clean data + +data = rawData; +numericFields = ["run_id", "sir", "pam_level", "symbolrate", ... + "interference_path_length", "occurrence_idx", "block_update", ... + "power_mpi_interference", "BER", "BER_precoded", "SNR", "GMI", "AIR"]; +for fieldIdx = 1:numel(numericFields) + fieldName = numericFields(fieldIdx); + if ismember(fieldName, string(data.Properties.VariableNames)) + data.(char(fieldName)) = numericColumn(data.(char(fieldName))); + end +end + +stringFields = ["storage_name", "algorithm", "algorithm_variant", "eq_class"]; +for fieldIdx = 1:numel(stringFields) + fieldName = stringFields(fieldIdx); + if ismember(fieldName, string(data.Properties.VariableNames)) + data.(char(fieldName)) = stringColumn(data.(char(fieldName))); + end +end + +data = data(isfinite(data.BER) & data.BER > 0 & data.BER < maxBerForPlot, :); +data.sir_exact = -7 - data.power_mpi_interference; +if isempty(data) + warning("PLOT_mpi_reduction_db_ber_vs_sir:NoRows", ... + "No rows remain after initial BER/path/study/block filtering."); + return +end + +if isempty(selectedPamLevels) + selectedPamLevels = unique(data.pam_level(isfinite(data.pam_level))).'; +else + data = data(ismember(data.pam_level, selectedPamLevels), :); +end + +if isempty(selectedAlgorithms) + selectedAlgorithms = unique(data.algorithm, "stable").'; +else + selectedAlgorithms = string(selectedAlgorithms); + data = data(ismember(data.algorithm, selectedAlgorithms), :); +end + +if isempty(data) + warning("PLOT_mpi_reduction_db_ber_vs_sir:NoSelectedRows", ... + "No rows remain after selectedPamLevels/selectedAlgorithms filtering."); + return +end + +data.clean_keep = true(height(data), 1); +[groupId, groupPam, groupAlgorithm, groupSir] = findgroups(data.pam_level, data.algorithm, data.sir); +for groupIdx = 1:max(groupId) + rowMask = groupId == groupIdx; + berValues = data.BER(rowMask); + if nnz(rowMask) > 3 + data.clean_keep(rowMask) = ~isoutlier(berValues); + end +end +cleanData = data(data.clean_keep, :); + +summaryTable = groupsummary(cleanData, ["pam_level", "algorithm", "sir"], ... + {"mean", "min", "max"}, "BER"); +sirExactTable = groupsummary(cleanData, ["pam_level", "algorithm", "sir"], ... + "median", "sir_exact"); +summaryTable.sir_exact = sirExactTable.median_sir_exact; +summaryTable = sortrows(summaryTable, ["pam_level", "algorithm", "sir"]); + +fprintf("Cleaned to %d rows across %d PAM/algorithm/SIR groups.\n", ... + height(cleanData), height(summaryTable)); +disp(groupcounts(cleanData, ["pam_level", "algorithm"])); + +%% 3) Plot data + +if exist("linspecer", "file") + algColors = linspecer(numel(selectedAlgorithms)); +else + algColors = lines(numel(selectedAlgorithms)); +end + +figure(); clf; +tiledlayout(numel(selectedPamLevels), 1, "TileSpacing", "compact"); + +for pamIdx = 1:numel(selectedPamLevels) + pamLevel = selectedPamLevels(pamIdx); + nexttile; hold on; + + for algIdx = 1:numel(selectedAlgorithms) + algorithmName = selectedAlgorithms(algIdx); + algColor = algColors(algIdx, :); + rawMask = cleanData.pam_level == pamLevel & cleanData.algorithm == algorithmName; + curveMask = summaryTable.pam_level == pamLevel & summaryTable.algorithm == algorithmName; + + if ~any(curveMask) + continue + end + scatter(cleanData.sir_exact(rawMask), cleanData.BER(rawMask), ... + 30, ... + "Marker", ".", ... + "MarkerEdgeColor", algColor, ... + "HandleVisibility", "off"); + + sirValues = summaryTable.sir_exact(curveMask).'; + meanBer = summaryTable.mean_BER(curveMask).'; + minBer = summaryTable.min_BER(curveMask).'; + maxBer = summaryTable.max_BER(curveMask).'; + valid = isfinite(sirValues) & isfinite(meanBer) & meanBer > 0; + + if useBoundedLines && exist("boundedline", "file") && any(valid) + yLower = max(meanBer - minBer, 0); + yUpper = max(maxBer - meanBer, 0); + yBounds = [yLower(:), yUpper(:)]; + + [hl, hp] = boundedline(sirValues(valid).', meanBer(valid).', yBounds(valid, :), ... + 'alpha', 'transparency', 0.08, ... + 'cmap', algColor, ... + 'nan', 'fill', ... + 'orientation', 'vert'); + set(hl, "LineStyle", "none", "Marker", "none", "HandleVisibility", "off"); + set(hp, "LineStyle", "none", "HandleVisibility", "off"); + end + + plot(sirValues(valid), meanBer(valid), ... + "LineStyle", "none", ... + "Marker", "o", ... + "MarkerSize", 5, ... + "LineWidth", 1.2, ... + "Color", algColor, ... + "MarkerFaceColor", algColor, ... + "DisplayName", char(algorithmName)); + + 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(3.8e-3, "LineWidth", 1, "LineStyle", "--", "HandleVisibility", "off"); + yline(2e-2, "LineWidth", 1, "LineStyle", "--", "HandleVisibility", "off"); + + title(sprintf("PAM %.0f, path %.0f m, block update %s", ... + pamLevel, pathLengthToPlot, blockUpdateLabel(selectedBlockUpdate))); + xlabel("SIR (dB)"); + ylabel("BER"); + set(gca, "YScale", "log"); + ylim([9e-5, maxBerForPlot]); + grid on; + box on; + legend("Location", "best", "Interpreter", "none"); +end + +if exist("beautifyBERplot", "file") + beautifyBERplot("logscale", true, "setcolors", false, "setmarkers", false); +end + +%% Local helpers + +function values = numericColumn(values) + if iscell(values) + values = string(values); + end + if isstring(values) || ischar(values) + values = str2double(values); + end + values = double(values); +end + +function values = stringColumn(values) + if iscell(values) + values = string(values); + elseif ischar(values) + values = string(values); + end + values = strip(string(values)); +end + +function label = blockUpdateLabel(selectedBlockUpdate) + if isempty(selectedBlockUpdate) + label = "pooled"; + else + label = string(selectedBlockUpdate); + end +end diff --git a/projects/Diss/MPI_revisit/db/appendMpiReductionDspOutput.m b/projects/Diss/MPI_revisit/db/appendMpiReductionDspOutput.m new file mode 100644 index 0000000..50b003c --- /dev/null +++ b/projects/Diss/MPI_revisit/db/appendMpiReductionDspOutput.m @@ -0,0 +1,99 @@ +function [summary, rows] = appendMpiReductionDspOutput(db, run_id, occurrence_idx, dspOutput, dsp_options, options) +%appendMpiReductionDspOutput Store one DSP occurrence in MpiReductionResults. + +arguments + db + run_id + occurrence_idx + dspOutput struct + dsp_options struct = struct() + options.study_name string = "mpi_reduction_v1" + options.dry_run (1,1) logical = false + options.verbose (1,1) logical = false +end + +summary = emptySummary(); +rows = {}; + +packageNames = fieldnames(dspOutput); +if isempty(packageNames) + return +end + +block_update = resolveBlockUpdate(dspOutput, dsp_options); +wh = DataStorage(struct( ... + "run_id", double(run_id), ... + "block_update", double(block_update))); + +hasPackages = false; +for packageIdx = 1:numel(packageNames) + packageName = packageNames{packageIdx}; + package = dspOutput.(packageName); + if isempty(package) + continue + end + + if iscell(package) + packageCell = package; + else + packageCell = {package}; + end + + wh.addStorage(packageName); + wh.addValueToStorageByLinIdx(packageCell, packageName, 1); + hasPackages = true; +end + +if ~hasPackages + return +end + +dsp_options_local = dsp_options; +dsp_options_local.start_occurence = occurrence_idx; +dsp_options_local.userParameters.block_update = block_update; + +[summary, rows] = appendMpiReductionWarehouse(db, wh, run_id, dsp_options_local, ... + "study_name", options.study_name, ... + "dry_run", options.dry_run, ... + "verbose", options.verbose); + +end + +function summary = emptySummary() + summary = struct( ... + "rows_considered", 0, ... + "rows_inserted", 0, ... + "rows_skipped_duplicate", 0, ... + "rows_skipped_existing_key", 0, ... + "rows_skipped_empty", 0); +end + +function block_update = resolveBlockUpdate(dspOutput, dsp_options) + if isfield(dsp_options, "userParameters") && ... + isfield(dsp_options.userParameters, "block_update") && ... + isscalar(dsp_options.userParameters.block_update) + block_update = dsp_options.userParameters.block_update; + return + end + + packageNames = fieldnames(dspOutput); + for packageIdx = 1:numel(packageNames) + package = dspOutput.(packageNames{packageIdx}); + if isempty(package) + continue + end + if iscell(package) + package = package{1}; + end + if isstruct(package) && ... + isfield(package, "mpi_reduction_config") && ... + isfield(package.mpi_reduction_config, "params") && ... + isfield(package.mpi_reduction_config.params, "block_update") + block_update = package.mpi_reduction_config.params.block_update; + return + end + end + + error("appendMpiReductionDspOutput:MissingBlockUpdate", ... + "Could not resolve scalar block_update from dsp_options or package metadata."); +end diff --git a/projects/Diss/MPI_revisit/db/appendMpiReductionWarehouse.m b/projects/Diss/MPI_revisit/db/appendMpiReductionWarehouse.m new file mode 100644 index 0000000..5d48d9f --- /dev/null +++ b/projects/Diss/MPI_revisit/db/appendMpiReductionWarehouse.m @@ -0,0 +1,307 @@ +function [summary, rows] = appendMpiReductionWarehouse(db, wh, run_ids, dsp_options, options) +%appendMpiReductionWarehouse Store MPI reduction warehouse entries in MySQL. + +arguments + db + wh + run_ids + dsp_options struct + options.study_name string = "mpi_reduction_v1" + options.dry_run (1,1) logical = false + options.verbose (1,1) logical = true +end + +ensureTableVisible(db); + +storageNames = fieldnames(wh.sto); +summary = struct( ... + "rows_considered", 0, ... + "rows_inserted", 0, ... + "rows_skipped_duplicate", 0, ... + "rows_skipped_existing_key", 0, ... + "rows_skipped_empty", 0); +rows = {}; + +for storageIdx = 1:numel(storageNames) + storageName = storageNames{storageIdx}; + storageArray = wh.sto.(storageName); + + for linIdx = 1:numel(storageArray) + packageCell = storageArray{linIdx}; + if isempty(packageCell) + summary.rows_skipped_empty = summary.rows_skipped_empty + 1; + continue + end + + phys = physicalCoordinates(wh, linIdx); + run_id = resolveRunId(phys, run_ids); + block_update = resolveBlockUpdate(phys, dsp_options); + packageCell = normalizePackageCell(packageCell); + + for packageIdx = 1:numel(packageCell) + package = packageCell{packageIdx}; + if isempty(package) || ~isstruct(package) || ~isfield(package, "metrics") + summary.rows_skipped_empty = summary.rows_skipped_empty + 1; + continue + end + + occurrence_idx = startOccurrence(dsp_options) + packageIdx - 1; + row = buildRow(db, package, storageName, run_id, occurrence_idx, ... + block_update, options.study_name); + summary.rows_considered = summary.rows_considered + 1; + rows{end+1,1} = row; %#ok + + if options.dry_run + continue + end + + if resultHashExists(db, row.result_hash) + summary.rows_skipped_duplicate = summary.rows_skipped_duplicate + 1; + if options.verbose + fprintf("MPI reduction DB: skipped duplicate %s, run_id=%d, occurrence=%d, block_update=%d\n", ... + row.storage_name, row.run_id, row.occurrence_idx, row.block_update); + end + continue + end + + if uniqueResultKeyExists(db, row) + summary.rows_skipped_existing_key = summary.rows_skipped_existing_key + 1; + if options.verbose + fprintf("MPI reduction DB: skipped existing key %s, run_id=%d, occurrence=%d, block_update=%d\n", ... + row.storage_name, row.run_id, row.occurrence_idx, row.block_update); + end + continue + end + + db.appendToTable("MpiReductionResults", row); + summary.rows_inserted = summary.rows_inserted + 1; + end + end +end + +if options.verbose + fprintf("MPI reduction DB: considered %d rows, inserted %d, duplicates %d, existing keys %d, empty %d\n", ... + summary.rows_considered, summary.rows_inserted, ... + summary.rows_skipped_duplicate, summary.rows_skipped_existing_key, ... + summary.rows_skipped_empty); +end + +end + +function ensureTableVisible(db) + if ~isfield(db.tables, "MpiReductionResults") + db.refresh(); + end + + if ~isfield(db.tables, "MpiReductionResults") + error("appendMpiReductionWarehouse:MissingTable", ... + "Table MpiReductionResults is not visible to DBHandler. Create the table and call db.refresh()."); + end +end + +function phys = physicalCoordinates(wh, linIdx) + phys = struct(); + [physValues, physNames] = wh.getPhysIndicesByLinIndex(linIdx); + + for physIdx = 1:numel(physNames) + name = char(physNames{physIdx}); + phys.(name) = physValues{physIdx}; + end +end + +function run_id = resolveRunId(phys, run_ids) + if isfield(phys, "run_id") + run_id = phys.run_id; + return + end + + if isscalar(run_ids) + run_id = run_ids; + return + end + + error("appendMpiReductionWarehouse:MissingRunId", ... + "Warehouse has no run_id axis, but run_ids is not scalar."); +end + +function block_update = resolveBlockUpdate(phys, dsp_options) + if isfield(phys, "block_update") + block_update = phys.block_update; + return + end + + if isfield(dsp_options, "userParameters") && ... + isfield(dsp_options.userParameters, "block_update") && ... + isscalar(dsp_options.userParameters.block_update) + block_update = dsp_options.userParameters.block_update; + return + end + + error("appendMpiReductionWarehouse:MissingBlockUpdate", ... + "Could not resolve block_update from warehouse coordinates or scalar userParameters."); +end + +function packageCell = normalizePackageCell(packageCell) + if ~iscell(packageCell) + packageCell = {packageCell}; + end +end + +function occurrence = startOccurrence(dsp_options) + occurrence = 1; + if isfield(dsp_options, "start_occurence") + occurrence = dsp_options.start_occurence; + end +end + +function row = buildRow(db, package, storageName, run_id, occurrence_idx, block_update, study_name) + config = fieldOr(package, "mpi_reduction_config", fallbackConfig(storageName, block_update)); + params = fieldOr(config, "params", struct()); + if isempty(params) + params = struct(); + end + + metrics = package.metrics; + if isobject(metrics) && ismethod(metrics, "toStruct") + metrics = metrics.toStruct(); + end + + param_json = jsonencode(params); + if strlength(string(param_json)) == 0 + param_json = "{}"; + end + + row = struct(); + row.run_id = double(run_id); + row.occurrence_idx = double(occurrence_idx); + row.study_name = char(study_name); + row.storage_name = char(fieldOr(config, "storage_name", storageName)); + row.algorithm = char(fieldOr(config, "algorithm", "unknown")); + row.algorithm_variant = char(fieldOr(config, "algorithm_variant", "")); + row.eq_class = char(fieldOr(config, "eq_class", "unknown")); + row.block_update = double(block_update); + + row.sps = scalarParam(params, "sps"); + row.eq_order = scalarParam(params, "order"); + row.len_tr = scalarParam(params, "len_tr"); + row.epochs_tr = scalarParam(params, "epochs_tr"); + row.mu_tr = scalarParam(params, "mu_tr"); + row.dd_mode = logicalParam(params, "dd_mode"); + row.epochs_dd = scalarParam(params, "epochs_dd"); + row.mu_dd = scalarParam(params, "mu_dd"); + + row.numBits = metricScalar(metrics, "numBits", 0); + row.numBitErr = metricScalar(metrics, "numBitErr", 0); + row.BER = metricScalar(metrics, "BER", NaN); + row.numBitErr_precoded = metricScalar(metrics, "numBitErr_precoded", NaN); + row.BER_precoded = metricScalar(metrics, "BER_precoded", NaN); + row.SNR = metricScalar(metrics, "SNR", NaN); + row.SNR_level = metricJson(metrics, "SNR_level"); + row.GMI = metricScalar(metrics, "GMI", NaN); + row.AIR = metricScalar(metrics, "AIR", NaN); + row.EVM = metricScalar(metrics, "EVM", NaN); + row.EVM_level = metricJson(metrics, "EVM_level"); + row.STD = metricScalar(metrics, "STD", NaN); + row.STD_level = metricJson(metrics, "STD_level"); + row.STDrx = metricScalar(metrics, "STDrx", NaN); + row.STDrx_level = metricJson(metrics, "STDrx_level"); + row.Alpha = firstFiniteScalar(fieldOr(metrics, "Alpha", NaN)); + + row.param_json = char(param_json); + row.param_hash = db.calcHash(params); + row.result_hash = db.calcHash(row); +end + +function config = fallbackConfig(storageName, block_update) + params = struct("block_update", block_update); + config = struct( ... + "storage_name", storageName, ... + "algorithm", "unknown", ... + "algorithm_variant", "", ... + "eq_class", "unknown", ... + "params", params); +end + +function value = fieldOr(s, fieldName, defaultValue) + fieldName = char(fieldName); + if isstruct(s) && isfield(s, fieldName) && ~isempty(s.(fieldName)) + value = s.(fieldName); + else + value = defaultValue; + end +end + +function value = scalarParam(params, fieldName) + value = firstFiniteScalar(fieldOr(params, fieldName, NaN)); +end + +function value = logicalParam(params, fieldName) + value = firstFiniteScalar(fieldOr(params, fieldName, false)); + if isnan(value) + value = false; + else + value = logical(value); + end +end + +function value = metricScalar(metrics, fieldName, defaultValue) + value = firstFiniteScalar(fieldOr(metrics, fieldName, defaultValue)); +end + +function value = metricJson(metrics, fieldName) + metricValue = fieldOr(metrics, fieldName, []); + value = char(jsonencode(metricValue)); +end + +function value = firstFiniteScalar(valueIn) + if islogical(valueIn) + valueIn = double(valueIn); + end + + if isempty(valueIn) + value = NaN; + return + end + + if isnumeric(valueIn) + finiteValues = valueIn(isfinite(valueIn)); + if isempty(finiteValues) + value = NaN; + else + value = double(finiteValues(1)); + end + return + end + + value = str2double(string(valueIn)); + if isnan(value) + value = NaN; + end +end + +function exists = resultHashExists(db, resultHash) + query = sprintf("SELECT mpi_result_id FROM MpiReductionResults WHERE result_hash = '%s' LIMIT 1", ... + char(resultHash)); + existing = db.fetch(query); + exists = ~isempty(existing); +end + +function exists = uniqueResultKeyExists(db, row) + query = sprintf( ... + "SELECT mpi_result_id FROM MpiReductionResults " + ... + "WHERE run_id = %d " + ... + "AND occurrence_idx = %d " + ... + "AND study_name = '%s' " + ... + "AND storage_name = '%s' " + ... + "AND block_update = %d " + ... + "AND param_hash = '%s' LIMIT 1", ... + row.run_id, row.occurrence_idx, ... + sqlString(row.study_name), sqlString(row.storage_name), ... + row.block_update, sqlString(row.param_hash)); + existing = db.fetch(query); + exists = ~isempty(existing); +end + +function value = sqlString(value) + value = strrep(char(value), '''', ''''''); +end diff --git a/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m b/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m index 7670ffd..7dc53ac 100644 --- a/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m +++ b/projects/Diss/MPI_revisit/investigate_mpi_algorithms.m @@ -7,6 +7,12 @@ dsp_options.start_occurence = 1; dsp_options.max_occurences = 15; dsp_options.debug_plots = false; +write_mpi_reduction_db = 1; +stream_mpi_reduction_db = 1; +mpi_reduction_study_name = "block_update_sweep"; +mpi_reduction_writer_path = fullfile(fileparts(mfilename('fullpath')),"db"); +addpath(mpi_reduction_writer_path); + dsp_options.database_type = "mysql"; dsp_options.dataBase = "labor"; @@ -16,6 +22,9 @@ dsp_options.server = "192.168.178.192"; dsp_options.port = 3306; dsp_options.user = "silas"; dsp_options.password = "silas"; +dsp_options.append_mpi_reduction_db = write_mpi_reduction_db && stream_mpi_reduction_db; +dsp_options.mpi_reduction_study_name = mpi_reduction_study_name; +dsp_options.mpi_reduction_writer_path = mpi_reduction_writer_path; db = DBHandler("dataBase", [dsp_options.dataBase],... "type", dsp_options.database_type,... @@ -45,7 +54,6 @@ for i = 1 [dataTable,~] = db.queryDB(fp, db.getTableFieldNames('Runs')); - [~, sortIdx] = sort(dataTable.sir, 'descend'); dataTable = dataTable(sortIdx, :); @@ -55,11 +63,10 @@ for i = 1 %% === Warehouse setup === dsp_options.userParameters = struct(); - dsp_options.userParameters.block_update = 1;%[1,2,4,8,16,32,64,112,224,448,512,1024,2048,2048*2,2048*4];%%logspace(-3.8,-1,22);%[linspace(2,4096,22)]; + dsp_options.userParameters.block_update = [1,2,4,8,16,32,64,112,224,448,448*2,1024,2048,4096,8192,16384];%%logspace(-3.8,-1,22);%[linspace(2,4096,22)]; % dsp_options.userParameters.block_update = linspace(1,224,22); wh = DataStorage(dsp_options.userParameters); - %% n_realizations = (dsp_options.max_occurences - dsp_options.start_occurence + 1); n_userparams = prod(wh.dim); @@ -76,6 +83,14 @@ for i = 1 [results, wh] = submitJobs(run_ids, dsp_options, processingMode.parallel, ... "wh", wh, ... "waitbar", true); + + %% + + if write_mpi_reduction_db && ~stream_mpi_reduction_db + db.refresh(); + appendMpiReductionWarehouse(db, wh, run_ids, dsp_options, ... + "study_name", mpi_reduction_study_name); + end % savepath = fullfile("C:","Users","Silas","Documents","MATLAB","imdd_simulation","projects","Diss","MPI_revisit","algorithms",sprintf("pam_%d_results.mat",M)); % save(wh,savepath) diff --git a/workerError.mat b/workerError.mat index b33c61c71ed4dd97732ddc29026a029402df862b..fe91e0e0f34c99a2cc479db5c37601e94b47f8e7 100644 GIT binary patch delta 1009 zcmV0C=43SV2$QFcfyP z!YE^GJAkHXhl&%23AI(x5SJAS6-?_u6aid#a^z?1XY42J zC(KS9a7{{}1qRxl^y)mv@AbX!IWO@W2LS9$0x-n^L6pKa0v6hXp&ifwn)?8e4kRYBVhQ@~f6!qH6x3bcM`nOs0-DPQCi0?7 z6NVIlR4SxTcTp}ymgh#m<@hxrKauX`#i$AFlKaEb{iom}--+GfYxw6$4vkCW570Z$ zi7#@V=cM(+DeL7p3t~Q|-ygVsNz%`!=;tNMG*!;AxIL}SY zn{H@Zf1ll813`=L_avUv47*VR34EUbo?{ z1p}xr*VpRB7PcYPk?ER7rBs>K9SDR4rH8Otn&aZ>%pgjjFL>l*<)k z7YB4|M#BGRzIS~{uuH2CFBo%20f7zzHYO;xe|M4*{}&~GF4yPvEsA`q$6+nGX8d&N zG2l9OEX$zS32n$9FTD|x=#iuPM|am)_g#}cvY51Wo|y$=wGN575t0^e6KK;Jhx@gX z;n#%X7nf7bfdTjJQKR-Ca9X0?D)GXg;kyotZq*tLr%UDr#R>!I-9~_rjJHqAlNbhZa*58`=# zBz3RwJ2h?_YoJPV@1vzF+#!$Lv7)xNp7w zzSYwj zxSU57&xB+=6aVG&s2Xo_7VmW)PwN+fE<5=<)gR07SB;ef4{DQ(VZr zB1L_ubL!cUVwCixV(M-$b@i!~MwR3z=#5K4a$NF$D`;28rHQM*QKxPi*H!se*WYx1 fUq@elucq~LcYus*EX(~|`FnAYelDj!$&WS_Me4Ouf=FA*YB~mhx_QG4?Pe0h9T^M zWIuXetgk9K`1rFzcHAe0@$s7IE3)efYd=~0E!scUn(W!JXg$_~uO0$DNOuq-br{;c z0ALG|0CzF+frBW8HUbvv!muCE0GeL_B0We=`rCnoZu zTN8#9fmA7{caKRSMVI@Q!TI<#q2H076~$E(+$GOT((}7ul<&lEDF;4J0#uQ%ryx7e z2^Kleuch&&l<{($WzjFB#_gX;_C>{B#>?e_^xpjmc6|JrV#oEy_Tjtjovk)uQv(xZ z?S&*te?9kiJjLHr#UIz(gv}k-Ciok&s$Jo%>iX3I^fB!8QFXnGZ4x@`n?u}xm$KJ$ zlH;`Xe4n*0TbF=92Z26=FY3b-+-(In*Xvt)nAxf8(0*YS3eXbnIA`L9r9| zA&*E;c%u652ZM!-Jw$cK!{Xks&|uRM@zti_f2k;VxIXS`-g$LU%fri|g5z*5Uf#ML zqv^XQD?OOlxf!T{`ITueF zo@E6O*Aq8$@$u(VEcI%)3pv!dYP`p6AP<@I9grxEe;YYL98f}=YF ze=nZzerB658;L0!)I4BSV1OwObW8?_2H1uHGI(`>ecSbM&~qZ9IKJi-e2*1;T+jBp z2;N!|Z>Qz#`d(7(d3{gv`2vXg9!XT>BZ(_L7t`&B4&is~BPKwcSunzVyAP!A1GC$| zoa_7KcK=D=e`e3et4qztP22Zo&tvUUfAert_uJX+pU;(#+kJkveg3-Y*j3xj@qfP% z`)T=qqTuFwcK%m}etVMNv~vy_SDt!Z0^>Ni9{>EG=U?R^5t)bNy=N0}CTB5OIMdDx z!I+=orP{Af!K=m})g$d{+*U2{)G3Z&9|)c Ld@lb1RuxptrOE&z