% Probability matrix
% Modes N= Normal, T=Tachycardia, B=Bradycardia
tic;
p_NN=0.9;
p_NT=0.05;
p_NB=0.05;

p_BN=0.9;
p_BT=0.05;
p_BB=0.05;

p_TN=0.9;
p_TT=0.05;
p_TB=0.05;

Prob = [p_NN,p_NT,p_NB;p_TN,p_TT,p_TB;p_BN,p_BT,p_BB];
time=0;

% Matrix of n-combinations (n_comb) of 3 modes
n_comb=8;
[COMB]=combinations(n_comb);

% Initial probability distribution per mode [N T B]
Init = [0.4, 0.3, 0.3];
pr_final = 0;
pr_final_after = 0;
exp_energy = 0;

fprintf(1,'Total number of iterations %d\n', size(COMB,1));
for i=1:size(COMB,1),
    if COMB(i,1) == 0,
            TH1=TH_N;
            N_STATES1=N_STATES_N;
    elseif COMB(i,1) == 1,
            TH1=TH_T;
            N_STATES1=N_STATES_T;
    else
            TH1=TH_B;
            N_STATES1=N_STATES_B;
    end
        
    for j=1:n_comb-1,
        if COMB(i,j+1) == 0,
            TH2=TH_N;
            N_STATES2=N_STATES_N;
        elseif COMB(i,j+1) == 1,
            TH2=TH_T;
            N_STATES2=N_STATES_T;
        else
            TH2=TH_B;
            N_STATES2=N_STATES_B;
        end
        
        % Combine the all transition matrices
        [TH, N_STATES ] = combinetransitions( N_STATES1, TH1, N_STATES2, TH2);
        TH1=TH;
        N_STATES1=N_STATES;
    end
    
    c_state_HEART=TH1(1,1);
    
    % Take the composition
    TH_COMPOSITION = composition(N_STATES1, TH1, TP, IM, c_state_PACEMAKER, c_state_HEART );
    
    
    % Check the correctness of the path
    correct=monitor(TH1,60,120);
    
    if correct==0,
        % Computing probabilities 
        [prob_val] = prob_computation(Prob,COMB(i,:));
        pr_final=pr_final+prob_val*Init(COMB(i,1)+1);
    end
    
        % Check the correctness of the path
    correct=monitor(TH_COMPOSITION,60,120);
    
     % Computing probabilities 
    [prob_val_after] = prob_computation(Prob,COMB(i,:));
    if correct==0,
            pr_final_after=pr_final_after+prob_val_after*Init(COMB(i,1)+1);
    end
    
   exp_energy=exp_energy+count_all_beats(TH_COMPOSITION)*prob_val_after*Init(COMB(i,1)+1);
    
fprintf(1,'Iteration number %d of %d\n',i, size(COMB,1));    
end

time=toc;
fprintf(1,'Running time of the experiment %d\n',time);