function [rr, targets_SR_AF, targets_APB, pafBoundaries, pafEpisodeLength] = simPAF_gen_rr_intervals(rrLength, fibFreqz, APBrate, realRRon, stayInAF, AFburden)
%
% [] = simPAF_gen_rr_intervals() returns RR intervals during PAF. 
% The dominant length of PAF episodes, and the percentage of the total time AF is
% present (AF burden) is determined by a first-order two-state Markov chain.
%
% An atrioventricular node model in which the ventricles are assumed to be 
% activated by the arriving atrial impulses according to a Poisson process
% is used to generate realistic RR series during an episode of AF (Corino et al 2011)
%
% Ventricular rhythm during SR is simulated using RR interval generator 
% proposed by McSharry et al(2003) in which both the impact of 
% parasympathetic stimulation (respiratory sinus arrhythmia) and baroreflex 
% regulation (Mayer waves) is modeled by a bimodal power spectrum
%
% In case of real data the entire MIT-BIH Normal Sinus Rhythm database, 
% consisting of 18 long-term ECG recordings, was taken as a basis for the set
% of 18 RR interval series of sinus rhythm.
% 
% The Long Term Atrial Fibrillation database was used to compose a set
% of AF rhythm. In total, 69 RR interval series were extracted out of 
% 84 long-term ECG recordings. The remaining records were excluded since 
% they did not meet the criteria of <5000 RR intervals of AF.
%
% Copyright (C) 2017  Andrius Petrenas
% Biomedical Engineering Institute, Kaunas University of Technology
%
% Available under the GNU General Public License version 3
% - please see the accompanying file named "LICENSE"
%

% Decide which rhythm type to generate
if AFburden == 0
    rhythmType = 0; % SR
elseif AFburden == 1
    rhythmType = 1; % AF
elseif AFburden > 0 && AFburden < 1
    rhythmType = 2; % PAF
else
    error('AF burden must be a value between 0 and 1')
end

switch rhythmType  % 0 - sinus rhythm, 1 - AF, 2 - PAF
    case 0 % The entire rhythm is SR
        if realRRon == 1 % Use real RR series
            rr = simPAF_construct_real_RR(0, rrLength);
        else % Use simulated RR series
            rr = simPAF_gen_SR_RR_intervals(rrLength);
            rr = rr(1:rrLength);
        end
        targets_SR_AF = zeros(1,rrLength);
        pafBoundaries(1,1) = 0; 
        pafBoundaries(1,2) = 0;
        pafEpisodeLength = 0;
        
    case 1 % The entire rhythm is AF
        if realRRon == 1 % Use real RR series
            rr = simPAF_construct_real_RR(1, rrLength);
        else
            rr = simPAF_gen_AF_RR_intervals(fibFreqz,rrLength);
            rr = rr(1:rrLength); 
        end
        targets_SR_AF = ones(1,rrLength);
        pafBoundaries(1,1) = 1; 
        pafBoundaries(1,2) = rrLength;
        pafEpisodeLength = rrLength;
        
    case 2 % PAF
        % Generate alternating SR and AF episodes
        goToNS = 1-stayInAF;
        goToAF = (goToNS*AFburden)/(1-AFburden);
        stayInNS = 1 - goToAF;

        trans = [stayInNS,1-stayInNS;1-stayInAF, stayInAF];
        [~,targets_SR_AF] = hmmgenerate(rrLength,trans,[0; 1]);

        % Find the number of AF episodes and the length of each episode
        k = 1;
        pafEpisodeLength = [];
        for p = 1:length(targets_SR_AF)-1
            if targets_SR_AF(p) == 2
                if targets_SR_AF(p+1) == 2
                    k = k + 1;
                    if p == (length(targets_SR_AF)-1)
                        pafEpisodeLength = [pafEpisodeLength k]; 
                    end
                else
                pafEpisodeLength = [pafEpisodeLength k];
                k = 1;
                end         
            end
        end
        % Find boundaries of each PAF episode
        if rrLength == length(find(targets_SR_AF == 1))% The entire signal is SR   
            pafBoundaries(1,1) = 0; 
            pafBoundaries(1,2) = 0;
        elseif rrLength == length(find(targets_SR_AF == 2))% The entire signal is AF
            pafBoundaries(1,1) = 1; 
            pafBoundaries(1,2) = rrLength;
        elseif rrLength ~= length(find(targets_SR_AF == 2))% The signal with PAF
            diffTar = diff(targets_SR_AF);
            j = 1;
            k = 1;
            flag = 1;
            for i = 1:length(diffTar)
                if diffTar(i) == 1
                    pafBoundaries(j,1) = i + 1;
                    j = j + 1;
                    flag = 1;
                end
                if diffTar(i) == -1
                    pafBoundaries(k,2) = i;
                    k = k + 1;
                    flag = 2;
                end   
                if i == length(diffTar)
                    if flag == 1
                        pafBoundaries(k,2) = i+1;
                    end
                end
            end
        end
        
        % Generate RR intervals during AF
        rrLengthAF = length(find(targets_SR_AF == 2));
        if realRRon == 1 % Use real RR series
            rrAF = simPAF_construct_real_RR(1, rrLengthAF);
        else
            rrAF = simPAF_gen_AF_RR_intervals(fibFreqz,rrLengthAF);
            rrAF = rrAF(1:rrLengthAF); 
        end
        
        % Generate RR intervals during SR
        rrLengthSR = length(find(targets_SR_AF == 1));
        if realRRon == 1 % Use real RR series
            rrSR = simPAF_construct_real_RR(0, rrLengthSR);
        else % Use simulated RR series
            rrSR = simPAF_gen_SR_RR_intervals(rrLengthSR);
            rrSR = rrSR(1:rrLengthSR);
        end
        % Make sure that the average RR value during SR is larger than that in AF
        if mean(rrSR) < mean(rrAF)
            rrSR = rrSR + (mean(rrAF) - mean(rrSR));
        end
        
        % Construct PAF RR series
        j = 1;
        k = 1;
        for i = 1:rrLength
            if targets_SR_AF(i) == 1
                rr(i) = rrSR(j); 
                j = j+1;
            else
                rr(i) = rrAF(k); 
                k = k+1;
            end
        end
        targets_SR_AF = targets_SR_AF - 1;
end

%% Insert atrial premature beats
if APBrate > 0.5
    error('APB rate must be less or equal to 0.5')
elseif APBrate > 0
    targets_APB = zeros(1,rrLength);
    transAPB = [1-APBrate,APBrate;0.99999,0.00001];
    [~,statesAPB] = hmmgenerate(rrLength,transAPB, [0; 1]);
    % Remove APBs during AF
    statesAPB = statesAPB.*(abs(targets_SR_AF-1));
    rrAPB=rr;
    
    apbType = randi([0 3],1,1);

    switch apbType 
        case 0 % APBs with sinus reset
            for i=1:length(statesAPB)-1
                if statesAPB (i) == 2
                    rrAPB(i)=0.8*rrAPB(i);   
                    targets_APB(i) = 1;
                end
            end
        case 1 % Interpolated PACs
            for i=1:length(statesAPB)-1
                if statesAPB (i) == 2
                    rrAPBtemp = rrAPB(i);
                    rrAPB(i)=0.6*rrAPB(i);
                    rrAPB(i+1)=rrAPBtemp - rrAPB(i);   
                    targets_APB(i) = 1;
                end
            end   
        case 2 % PACs with delayed sinus reset
            for i=1:length(statesAPB)-1
                if statesAPB (i) == 2
                    rrAPBtemp = rrAPB(i);
                    rrAPB(i)= 0.8*rrAPB(i);
                    rrAPB(i+1)= 1.2*rrAPB(i+1);   
                    targets_APB(i) = 1;
                end
            end
            case 3 %PACs with full compensatory pause
                for i=1:length(statesAPB)-1
                    if statesAPB (i) == 2
                        rrAPBtemp = rrAPB(i);
                        rrAPB(i)= 0.8*rrAPB(i);
                        rrAPB(i+1)= 2*rrAPBtemp-rrAPB(i);   
                        targets_APB(i) = 1;
                    end
                end            
    end
    rr = rrAPB;
else
    targets_APB = zeros(1,rrLength);
end


% %% Switching between AF an SR
% goToNS = 1-stayInAF;
% goToAF = (goToNS*AFburden)/(1-AFburden);
% stayInNS = 1 - goToAF;
% 
% trans = [stayInNS,1-stayInNS;1-stayInAF, stayInAF];
% [~,states_SR_AF] = hmmgenerate(ecgLength,trans,[0; 1]);
% 
% %% Calculate the number of AF episodes generate by Markov chain
% k = 1;
% episodeLength = [];
% for p = 1:length(states_SR_AF)-1
%     if states_SR_AF(p) == 2
%         if states_SR_AF(p+1) == 2
%             k = k + 1;
%            if p == (length(states_SR_AF)-1)
%               episodeLength = [episodeLength k]; 
%            end
%         else
%             episodeLength = [episodeLength k]; 
%             k = 1;
%         end         
%     end
% end
% 
% %% Find boundaries of each PAF episode
% if ecgLength == length(find(states_SR_AF == 1))
%     % The entire signal is SR
%     beginOfAF(1,1) = 0; 
%     beginOfAF(1,2) = 0;
% elseif ecgLength == length(find(states_SR_AF == 2))
%     % The entire signal is AF
%     beginOfAF(1,1) = 1; 
%     beginOfAF(1,2) = ecgLength;
% elseif ecgLength ~= length(find(states_SR_AF == 2))
%     diffTar = diff(states_SR_AF);
%     j = 1;
%     k = 1;
%     flag = 1;
%     for i = 1:length(diffTar)
%         if diffTar(i) == 1
%             beginOfAF(j,1) = i + 1;
%             j = j + 1;
%             flag = 1;
%         end
%         if diffTar(i) == -1
%             beginOfAF(k,2) = i;
%             k = k + 1;
%             flag = 2;
%         end
%     
%         if i == length(diffTar)
%             if flag == 1
%             beginOfAF(k,2) = i+1;
%             end
%         end
%     end
% end
% 
% %% 
%    
% 
end