% Read trace file and initialize variables
pps_fname = 'packsPerSec';
time = 0;
pkts = 0;
peak = 0;
[pkts] = textread(pps_fname, '%u');
time = [1:length(pkts)]';
peak = max (pkts);
datamin = min (pkts);

% ************PART 1************************
% Algorithm Lambda from reference [1]
% Determine rates and number of states
lambda_array = 0;
j = 0;
a = 2;
w = min(pkts);
N = 1;
lambda_array(1) = ((sqrt (1 + peak)) - 1)^2;
r = lambda_array(1) - (a * sqrt(lambda_array(1)));
while (r>w)
    j = N + 1;
    lambda_array(j) = (sqrt(lambda_array(j-1)) - a)^2;
    r = lambda_array(j) - (a * sqrt(lambda_array(j)));
    N = j;
end

% % Optional code to show the states on the same plot as data
% ocplotthis = 0;
% ocsizeof = 0;
% ocsizeof = size (time);
% ocY = ones(ocsizeof,1);
% ocplotthis= ocY*(lambda_array(1));
% for ocj = 2:N
%     ocplotthis= [ocplotthis ocY*(lambda_array(ocj))];
%     for oci=1:ocsizeof
%         ocplotthis(oci,ocj) = lambda_array(ocj);
%     end
% end
% ocplotthis = [pkts ocplotthis];
% figure ;
% plot (time', ocplotthis,'.-');
% grid on;

% ************PART 2************************
% Assign each sample to a state
state_array = 0;
datasize = 0;
[datasize, trash]= size(pkts);
for i = 1:datasize
    for j = 1:N
        if (((lambda_array(j) - (a* sqrt(lambda_array(j)))) < pkts(i)) && ((lambda_array(j) + (a* sqrt(lambda_array(j)))) >= pkts(i)));
            break;
        end
    end
    state_array(i) = j;
end

m = 0; w = 0; x = 0;
m = 1000; w = 1000; x = [1:N]; 
U = 0; Un = 0; S = 0; n_SS=1; S(1)=0;
D = 0; V = 0; junk1=0; junk2=0;
current_index = m;
DS = cell(1,1); VS = cell(1,1); Uindex = cell(1,1); 
[D junk1 V junk2] = find_dominant(1, m, state_array, N);
DS{1,n_SS} = D; VS{1,n_SS} = V;  
startstopindices = 0; corresponding_state_set = 0;

markhere = 0;
for i = 1:datasize
    markhere(i) = 0;
end

first_SS_separator_found = 0;
comparable_SS_pointer = 1;
while ((current_index+w)<datasize)
    diffset1=0; diffset2=0;
    [D junk1 V junk2] = find_dominant(current_index, (current_index + w), state_array, N);
    diffset1 = setdiff(DS{1,comparable_SS_pointer}, D);
    diffset2 = setdiff(D, DS{1,comparable_SS_pointer});
    if (((size(diffset1,2)) + (size(diffset2,2)))>2)
        markhere(current_index)=500;
        if (first_SS_separator_found==0)
            first_SS_separator_found = 1;
            startstopindices(1)=current_index;
            corresponding_state_set(1)=1;
        else
            startstopindices((size(startstopindices,2))+1) = current_index;
            corresponding_state_set((size(corresponding_state_set,2))+1) = comparable_SS_pointer;
        end
        
        same_SS_found = 0;
        for nss_counter=1:n_SS
            if ((size(setdiff(DS{1,nss_counter},D),2) +  size(setdiff(D,DS{1,nss_counter}),2))<=2)
                same_SS_found = 1;
                comparable_SS_pointer = nss_counter;
                VS{1,comparable_SS_pointer} = union(VS{1,comparable_SS_pointer}, V);
                break;
            end
        end
        if (same_SS_found == 0)
            n_SS = n_SS + 1;
            DS{1,n_SS} = D; VS{1,n_SS} = V;
            comparable_SS_pointer = n_SS;
        end
    else
        VS{1,comparable_SS_pointer} = union(VS{1,comparable_SS_pointer}, V);
    end
    current_index = current_index + (w);
end
startstopindices((size(startstopindices,2))+1) = datasize+1;
corresponding_state_set((size(corresponding_state_set,2))+1) = comparable_SS_pointer;

for i = 1:n_SS
    U((Un+1):(Un+size(VS{1,i},2))) = VS{1,i};
    Uindex{1,i} = [(Un+1):(Un+size(VS{1,i},2))];
    Un = (Un+size(VS{1,i},2));
end

new_state_array = 0; current_SS_pointer_index = 1; state_in_SS_pointer = 0;
for i = 1 : datasize
    if (i==startstopindices(current_SS_pointer_index))
        current_SS_pointer_index = current_SS_pointer_index + 1;
    end
    state_in_SS_pointer = find((VS{1,corresponding_state_set(current_SS_pointer_index)})==state_array(i));
    new_state_array(i)=Uindex{1,corresponding_state_set(current_SS_pointer_index)}(state_in_SS_pointer);
end

% save octwkspace;


% % Optional plot of states versus time
% figure 
% plot (time, new_state_array')


% Form the state transition matrix 
state_transition_array = 0;
transition_probability_matrix = 0;
Q = 0;
for i = 1:Un
    for j = 1:Un
        state_transition_array(i,j) = 0;
        transition_probability_matrix(i,j) = 0;
        Q(i,j) = 0;
    end
end
for i = 1:(datasize-1)
    state_transition_array(new_state_array(i),new_state_array(i+1)) = state_transition_array(new_state_array(i),new_state_array(i+1)) + 1;
end

% ************PART 3************************
% Form the state transition probability matrix from state transition matrix
total_time_in_state = 0;
total_time_in_state = sum (state_transition_array');
for i = 1:Un
    for j = 1:Un
        if (total_time_in_state (i) ~= 0)
            transition_probability_matrix(i,j) = (state_transition_array(i,j) / total_time_in_state (i));
        else
            transition_probability_matrix(i,j) = 0;
        end
    end
end

% With the new method, there will be no states that are never visited
% 
% temp_matrix = 0
% temp_matrix = transition_probability_matrix;
% new_N = Un;
% for i = 1:Un
%     check_if_all_zero = 1;
%     for j = 1:Un
%         if ((temp_matrix(i,j) ~= 0)&&(temp_matrix(i,j) ~= 2))
%             check_if_all_zero = 0;
%             break;
%         end
%     end
%     if (check_if_all_zero == 1)
%         for j = 1:N
%             temp_matrix(i,j) = 2;
%             temp_matrix(j,i) = 2;
%         end
%     end
% end
% 
%         
% Form infinitismal generator matrix Q from state transition probablility matrix
Q = transition_probability_matrix - eye(Un);
% 
% 
% % % Optional plot for number of times in each state
% % figure (4);
% % plot (total_time_in_state, '.-');
% % grid on;
% % xlabel ('state');
% % ylabel ('time in state');
% 

% ************PART 4************************

% % Generate a MMPP sample from the Q matrix
% % Also count number of packets per second generated by MMPP and form the matrix MMPP_pkts
next_state_probability = 0;
for i=1:2
    for j=1:Un-1
        for k=1:Un
            next_state_probability(i,j,k)=0;
        end
    end
end

for k = 1:Un
    for i = 1:1
        counter = 1;
        for j = 1:Un-1
            if (j == k)
                counter = counter + 1;
            end
            next_state_probability(i,j,k) = counter;
            counter = counter + 1;
            next_state_probability(i+1,j,k) =  Q(k,next_state_probability(i,j,k))/(-1*Q(k,k));
        end
    end
end


MMPP_state_array_size = 0;
MMPP_samples_time_array = 0;
MMPP_samples_state_array = 0;
initial_state_prob = 0;
initial_state_prob = 1/Un;
initial_state = 0;
i = 1;
for j = 1:Un
    initial_state(i,j)=j;
    initial_state(i+1,j)=initial_state_prob;
end

        
current_state = randsrc (1,1,initial_state);
current_time = 0;
next_interrupt_time = 0;
i = 1;
MMPP_samples_time_array (i) = current_time;
MMPP_samples_state_array (i) = current_state;
next_interrupt_time = current_time + exprnd(1/(-Q(current_state, current_state)));
while (next_interrupt_time <= time (datasize))
    i = i + 1;
    current_time = next_interrupt_time;
    current_state = randsrc (1,1,next_state_probability(:,:,current_state));
    MMPP_samples_time_array (i) = current_time;
    MMPP_samples_state_array (i) = current_state;
    next_interrupt_time = current_time + exprnd(1/(-Q(current_state, current_state)));
end
MMPP_state_array_size = i;

% 
% % % Optional plot
% % figure ;
% % plot (MMPP_samples_time_array, MMPP_samples_state_array);
% 
% 

% ************PART 5************************
% Now that the states and times are known, generate the packets
current_time = 0;
packet_counter = 0;
next_interrupt_time_delta = 0;
MMPP_pkts = 0;
% MMPP_pkt_arvl_time_array = 0;
i = 1;
k = 1;
% j = 1;
while (i<= MMPP_state_array_size)
    next_interrupt_time_delta = exprnd(1/(lambda_array(U(MMPP_samples_state_array(i)))));
    if (i == MMPP_state_array_size)
        current_time = min((current_time + next_interrupt_time_delta), datasize);
        if (current_time == datasize)
            i = i + 1;
        end
    else
        current_time = min((current_time + next_interrupt_time_delta), MMPP_samples_time_array(i+1));
        if (current_time == MMPP_samples_time_array(i+1))
            i = i + 1;
        end
    end
    if (current_time > k)
        MMPP_pkts(k) = packet_counter;
        k = k + 1;
        packet_counter = 0;
    end
    packet_counter = packet_counter + 1;
%     MMPP_pkt_arvl_time_array(j) = current_time;
%     j = j + 1;
end
MMPP_pkts(k) = packet_counter;
save OriginalMMPPWorkspace;
% 
% figure ;
% plot (time, pkts);
% grid on;
% xlabel('Time (Seconds)');
% ylabel('Arrivals (Packets/Second)');
% saveas(gcf,'raw_data.fig');
% 
% figure ;
% plot (time, MMPP_pkts);
% grid on;
% xlabel('Time (Seconds)');
% ylabel('Arrivals (Packets/Second)');
% saveas(gcf,'unimproved_MMPP_data.fig');
% 
% figure;
% qqplot(pkts,MMPP_pkts);
% xlabel('Trace Quantiles');
% ylabel('Quantiles of fitted process');
% grid on;
% saveas(gcf,'quantile_plot_original.fig');