% 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

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


% Form the state transition matrix 
state_transition_array = 0;
transition_probability_matrix = 0;
Q = 0;
for i = 1:N
    for j = 1:N
        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(state_array(i),state_array(i+1)) = state_transition_array(state_array(i),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:N
    for j = 1:N
        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

temp_matrix = 0
temp_matrix = transition_probability_matrix;
new_N = N;
for i = 1:N
    check_if_all_zero = 1;
    for j = 1:N
        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(N);


% % Optional plot for number of times in each state
% figure;
% 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:N-1
        for k=1:N
            next_state_probability(i,j,k)=0;
        end
    end
end

for k = 1:N
    for i = 1:1
        counter = 1;
        for j = 1:N-1
            if (j == k)
                counter = counter + 1;
            end
            next_state_probability(i,j,k) = counter;
            counter = counter + 1;
            if (temp_matrix(k,k)==2)
                next_state_probability(i+1,j,k) = 0;
            else
                next_state_probability(i+1,j,k) =  Q(k,next_state_probability(i,j,k))/(-1*Q(k,k));
            end
        end
    end
end


MMPP_state_array_size = 0;
MMPP_samples_time_array = 0;
MMPP_samples_state_array = 0;
initial_state = 0;
i = 1;
for j = 1:N
    initial_state(i,j)=j;
    initial_state(i+1,j)=j;
end
i = 1;
zero_initial_state_counter = 0;
for j = 1:N
    if (temp_matrix(j,j)==2)
        zero_initial_state_counter = zero_initial_state_counter  + 1;
    end
end
initial_state_prob = 0;
initial_state_prob = 1/(N - zero_initial_state_counter);
for j = 1:N
    if (temp_matrix(j,j)==2)
        initial_state(i+1,j)=0;
    else
        initial_state(i+1,j)=initial_state_prob;
    end
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(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 MayOriginalMMPPWorkspace;

% 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');