% FUNCTION NAME:
%   desat
%
% Description:
%   Outsourced function that calculates the hypoxic burden index of an SpO2 signal
%
% AUTHOR:
%   Antoine Weihs <antoine.weihs@uni-greifswald.de>
% Version: 
%	1.0
%
% INPUT:
%   SpO2 - (double) vector containing the SpO2 signal sampled at 1Hz
%   ID - (string) currently analysed subject ID
%   Wave - (string) assesment wave
%   threshold - (double) threshold defining a desaturation event. Default:
%                        2%
%
% OUTPUT:
%   Output - (table) Output table of dimentions 1 x 10 with the
%                    following columns:
%                    * ID (string): subject ID
%                    * Wave (string): wave ID
%                    * AREA (double): area under the curve of desaturation
%                                     events (percentage * time in seconds)
%					 * Threshold (double): used threshold
%                    * TST (double): Assessed total sleep time - Artefacts
%                    * TotalAREA: (double) Total area under the cuve
%                    * HBI: (double) Hypoxic burden index (percentage
%                                    time in seconds / time in seconds)
%                    * numEvents (double): number of desaturation events
%                    * events_per_hour (double): Number of desaturation 
%                                                events per hour
%
% REFERENCES:
%   Code adapted from Esmaeili N et al., Hypoxic Burden Based on 
%   Automatically Identified Desaturations Is Associated with Adverse 
%   Health Outcomes. Ann Am Thorac Soc. 2023 Nov;20(11):1633-1641. 
%   doi: 10.1513/AnnalsATS.202303-248OC. PMID: 37531573; 
%   PMCID: PMC10632930.
% 
% License: This code is licensed under the MIT License.

function [Output] = HypoxicBurden(SpO2, ID, Wave, threshold)
    arguments
        SpO2 (1,:) double {mustBeNumeric,mustBeReal}
        ID (1,1) string
        Wave (1,1) string
        threshold (1,1) double {mustBeNumeric,mustBeReal, mustBePositive} = 2
    end

    Output = table('Size', [1, 9], ...
        'VariableTypes', {'string', 'string', 'double', 'double', 'double', 'double', 'double', 'double', 'double'}, ...
        'VariableNames', {'ID', 'Wave', 'Threshold', 'AREA', 'TST', 'TotalAREA', 'HBI', 'numEvents', 'events_per_hour'});
    

    desatEvents = [];
    l_index = 1;
    while l_index < length(SpO2) - 10
        temp = SpO2(l_index:(l_index + 9));
        if min(temp) < SpO2(l_index) - threshold
            desatEvents = [desatEvents, l_index+ find(temp == min(temp), 1, "first")];
            l_index = l_index +10;
        else
            l_index = l_index +1;
        end
    end
    clear l_index; clear temp;


    %calculate total area under the curve
    TotalAREA = SpO2(~isnan(SpO2));
    TotalAREA = trapz(TotalAREA);
    
    if isempty(desatEvents)
        Output.ID(1) = ID;
        Output.Wave(1) = Wave;
		Output.Threshold = threshold;
        Output.AREA(1) = 0;
        Output.TST(1) = length(find(~isnan(SpO2) == 1));
        Output.TotalAREA(1) = TotalAREA;
        Output.HBI(1) = 0;
        Output.numEvents(1) = 0;
        Output.events_per_hour(1) = 0;
        clear desatEvents;
    else
        %create desaturation curves
        desatWindows = nan(length(desatEvents), 201);
        for l_index = 1:length(desatEvents)
            if desatEvents(l_index)-100 < 0 || desatEvents(l_index)+100 > length(SpO2)
                continue
            end
            desatWindows(l_index, :) = transpose(SpO2((desatEvents(l_index)-100):(desatEvents(l_index)+100)));
        end
        clear l_index; clear temp;
        dirty = desatWindows;
        for l_index = 1:height(desatWindows)
            if anynan(desatWindows(l_index,:))
                desatWindows(l_index,1) = NaN;
            end
        end
        desatWindows = desatWindows(~isnan(desatWindows(:,1)),:);
        if isempty(desatWindows)
            errorStruct.message = 'No clean window to calculate desaturation window';
            errorStruct.identifier = 'MyScript:NoCleanWindow';
            error(errorStruct)
        end
    
        ENS_AVG = NaN(1, width(desatWindows));
        for l_index = 1:width(desatWindows)
            ENS_AVG(l_index) = mean(desatWindows(:, l_index), 'omitmissing');
        end
        clear l_index;
    
        [~, locs] = findpeaks(ENS_AVG(1:100));
        if length(locs) > 1
            peak1 = locs(find(locs < 90, 1, 'last'));
        else
            peak1 = locs;
        end
        [~, locs] = findpeaks(ENS_AVG(101:201));
        peak2 = locs(1) + 100;
        clear locs;

        %define search window
        padded1 = max(peak1 - round((peak2 - peak1)*0.1, TieBreaker="fromzero"), 1);
        padded2 = min(peak2 + round((peak2 - peak1)*0.1, TieBreaker="fromzero"), width(desatWindows));

        desatWindows_cropped = dirty(:, padded1:padded2);
        clear desatWindows; clear dirty; clear padded1; clear padded2;
    
        %calculate sum of the area under the curve -> Plain
        myArea = 0;
        for l_index = 1:height(desatWindows_cropped)
            temp = desatWindows_cropped(l_index, :);
            temp = temp(~isnan(temp));
            myArea = myArea + trapz(temp);
        end
        clear l_index; clear temp;

        Output.ID(1) = ID;
        Output.Wave(1) = Wave;
		Output.Threshold = threshold;
        Output.AREA(1) = myArea; %percentage * time in seconds
        Output.TST(1) = length(find(~isnan(SpO2) == 1)); %seconds
        Output.TotalAREA(1) = TotalAREA;
        Output.HBI(1) = Output.AREA(1)/Output.TST(1);
        Output.numEvents(1) = length(desatEvents); %count
        Output.events_per_hour(1) = length(desatEvents)/(length(find(~isnan(SpO2) == 1))/60/60); %count / time in hours
    
        clear desatEvents; clear myArea; clear peak1; clear peak2; clear ENS_AVG;
    end
end