function [ml_mlse_results] = ml_mlse(eq_, M, rx_signal, tx_symbols, tx_bits, options)
% 
%
% Inputs:
%   eq_ - Equalizer object
%   M - Modulation order
%   rx_signal - Received signal
%   tx_symbols - Transmitted symbols
%   tx_bits - Transmitted bits
%   options - Optional parameters
%
% Outputs:
%   ffe_results - Results from FFE processing

arguments
    eq_
    M
    rx_signal
    tx_symbols
    tx_bits
    options.precode_mode db_mode
    options.eth_style_symbol_mapping = 0;
    options.postFFE = [];

end

%% Process signals through equalizer

[eq_signal_hd,y_ref] = eq_.process(rx_signal,tx_symbols);

%% Calculate BER based on precoding mode
[bits, errors, ber, error_pos, errors_precoded, ber_precoded] = calculateBER(eq_signal_hd, tx_symbols, tx_bits, options.precode_mode, M, options.eth_style_symbol_mapping);


% Create FFE results structure
ml_mlse_results = struct();
try
    eq_.e = [];
    eq_.e2 = [];
    eq_.e3 = [];
    eq_.b = [];
    eq_.b2 = [];
    eq_.b3 = [];
end

ml_mlse_results.config = Equalizerstruct();

eq_ = strip_eq(eq_, 10);

ml_mlse_results.config.eq = jsonencode(eq_);
ml_mlse_results.config.equalizer_structure = int32(equalizer_structure.ml_mlse);
ml_mlse_results.config.comment = 'function: ML-based MLSE';

ml_mlse_results.metrics = Metricstruct;
% ml_mlse_results.metrics.result_id = NaN;
% ml_mlse_results.metrics.run_id = NaN;
% ml_mlse_results.metrics.eqParam_id = NaN;
ml_mlse_results.metrics.date_of_processing = datetime('now');
ml_mlse_results.metrics.BER = ber;
ml_mlse_results.metrics.numBits = bits;
ml_mlse_results.metrics.numBitErr = errors;
ml_mlse_results.metrics.BER_precoded = ber_precoded;
ml_mlse_results.metrics.numBitErr_precoded = errors_precoded;
% ml_mlse_results.metrics.SNR = NaN;
% ml_mlse_results.metrics.SNR_level = NaN;
% ml_mlse_results.metrics.STD = NaN;
% ml_mlse_results.metrics.STD_level = NaN;
% ml_mlse_results.metrics.STDrx = NaN;
% ml_mlse_results.metrics.STDrx_level = NaN;
% ml_mlse_results.metrics.GMI = NaN;
% ml_mlse_results.metrics.AIR = NaN;
% ml_mlse_results.metrics.EVM = NaN;
% ml_mlse_results.metrics.EVM_level = NaN;
% ml_mlse_results.metrics.Alpha = NaN;


end

%% Helper Functions
function [bits, errors, ber, error_pos, errors_precoded, ber_precoded] = calculateBER(eq_signal_hd, tx_symbols, tx_bits, precode_mode, M, eth_style)
% Calculate BER based on precoding mode
mapper = PAMmapper(M, 0, "eth_style", eth_style);

switch precode_mode
    case db_mode.no_db
        % TX Data is not precoded
        % A) Emulate diff precoding
        eq_signal_hd_precoded = Duobinary().encode(eq_signal_hd, "M", M);
        eq_signal_hd_precoded = Duobinary().decode(eq_signal_hd_precoded, "M", M);
        
        tx_symbols_precoded = Duobinary().encode(tx_symbols);
        tx_symbols_precoded = Duobinary().decode(tx_symbols_precoded);
        
        tx_bits_precoded = mapper.demap(tx_symbols_precoded);
        
        rx_bits = mapper.demap(eq_signal_hd_precoded);
        [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits.signal, tx_bits_precoded.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1);
        
        % B) Just determine BER
        rx_bits = mapper.demap(eq_signal_hd);
        [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1);
        
    case db_mode.db_precoded
        % Data is precoded on TX side
        % A) Decode at Rx if no DB targeting was applied
        eq_signal_hd_decoded = Duobinary().encode(eq_signal_hd, "M", M);
        eq_signal_hd_decoded = Duobinary().decode(eq_signal_hd_decoded, "M", M);
        rx_bits_decoded = mapper.demap(eq_signal_hd_decoded);
        [~, errors_precoded, ber_precoded, ~] = calc_ber(rx_bits_decoded.signal, tx_bits.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1);
        
        % B) Omit the Coding by comparing with demapped TX symbol sequence
        tx_bits_demapped = mapper.demap(tx_symbols);
        rx_bits = mapper.demap(eq_signal_hd);
        [bits, errors, ber, error_pos] = calc_ber(rx_bits.signal, tx_bits_demapped.signal, "skip_front", 30000, "skip_end", 150, "returnErrorLocation", 1);
end
end


function eq_out = strip_eq(eq_, max_elems)
    % strip_eq removes all large fields from the ML_MLSE object
    % eq_out = strip_eq(eq_, max_elems)
    % max_elems ... maximum number of elements to keep (default = 10)

    if nargin < 2
        max_elems = 10; % default threshold
    end

    props = properties(eq_);
    for i = 1:numel(props)
        val = eq_.(props{i});
        if ~isempty(val)
            % Count total number of elements
            if numel(val) > max_elems
                eq_.(props{i}) = [];
            end
        end
        if issparse(val)
            eq_.(props{i}) = find(eq_.(props{i}));
        end
    end
    eq_out = eq_;
end