%---------------------------------------------------|
%-------------~ANALYTICAL SOLUTION~-----------------|
%------FIRST ORDER FEEDBACK FLUID QUEUES------------|
%-----------WITH MULTIPLE THRESHOLDS----------------|
%---------------FINITE BUFFER-----------------------|
%---------(WITHOUT SCHUR DECOMPOSITION)-------------|
%---------------------------------------------------|
%Author:H. Emre KANKAYA
clear                                              %|
ctime=cputime;                                     %|
%---------Output: Statewise PDF Plot----------------|
%---------Output: Total PDF Plot--------------------|

disp('---------------')                            %|
disp('PROGRAM STARTED')                            %|

%------------------PARAMETERS-----------------------|                                                   
%Q{i}=Q at region i                                %|
%QT{i}=Q at threshold i-1                          %|
%R{i}=rates at region i                            %|
%ratesT{i}=rates at threshold i-1                  %|
%EXAMPLE-1                                         %|
Q{1}=[-3,1,1,1;1,-3,1,1;1,1,-3,1;1,1,1,-3];        %|
Q{2}=Q{1};                                         %|
QT{1}=Q{1};                                        %|
QT{2}=Q{2};                                        %|
QT{3}=Q{2};                                        %|
R{1}=[-1,2,3,1];                                   %|
R{2}=[-3,-1,2,1];                                  %|
thresholds=[0,1,3];                                %|
noofsamples=2000;%For Plotting                     %|
%---------------------------------------------------|

disp('PREPARATION FOR THE SOLUTION')
%-------------------PREPARE-------------------------|
s=length(R{1});%number of states                   %|
r=length(R);%number of regions                     %|
for i=1:r                                          %|
    %-   -   -   -   -   -   -   -   -   -   -   - %|
    %A)Find Eigenvalues and Eigenvectors:          %|
    %------------------------------------          %|
    [v,d]=eig(Q{i}',diag(R{i})');                  %|
    %# of 0 eigenvales=1                           %|
    %# of - eigenvalues=noofL                      %|
    %# of + eigenvalues=noofU-1                    %|                                                   
    %-   -   -   -   -   -   -   -   -   -   -   - %|
    %B)Reorder the eigenvalues and eigenvectors:   %|
    %-------------------------------------------   %|
    [y,index]=sort(diag(d));                       %|
    dnew=zeros(s,s);%ordered eigenvalues           %|
    vnew=zeros(s,s);%ordered eigenvectors          %|
    for j=1:s                                      %|
        dnew(j,j)=d(index(j),index(j));            %|
        vnew(j,:)=(v(:,index(j)))';                %|
    end                                            %|
    dnewcell{i}=dnew;                              %|
    vnewcell{i}=vnew;                              %|
end                                                %|
accsstate=zeros(1,s);%takes values in 1->Nstates   %|
accsbound=zeros(1,s);%takes values in 1->Nregions  %|
%if accstate(i)=1, accsbound(i)=2                  %|
%then there is an accumulation atstate i boundary 2%| 
%-   -   -   -   -   -   -   -   -   -   -   -   - %|
%C)Find Accumulations:                             %|
%---------------------                             %|
k=1;                                               %|
%Accumulations at zero                             %|
for i=1:s                                          %|
    if R{1}(i)<0                                   %|
        accsstate(k)=i;                            %|
        accsbound(k)=0;                            %|
        k=k+1;                                     %|
    end                                            %|
end                                                %|
%Accumulation at Intermediate Boundaries           %|
for j=1:r-1                                        %|
    for i=1:s                                      %|
        if R{j}(i)>0&R{j+1}(i)<0                   %|
            accsstate(k)=i;                        %|
            accsbound(k)=j;                        %|
            k=k+1;                                 %|
        end                                        %|
    end                                            %|
end                                                %|
%Accumulations at Buffer Length                    %|
for i=1:s                                          %|
    if R{r}(i)>0                                   %|
        accsstate(k)=i;                            %|
        accsbound(k)=r;                            %|
        k=k+1;                                     %|
    end                                            %|
end                                                %|
%---------------------------------------------------|


disp('SOLVING THE SYSTEM')
%---------------------SOLUTION----------------------|
totalmatrix=zeros(r*s+s,r*s+s+1);                  %|
%A) Boundary Conditions:                           %|
%-----------------------                           %|
%Boundary C. at Zero                               %|
mat1=vnewcell{1}*diag(R{1});                       %|
totalmatrix(1:s,1:s)=mat1;                         %|
%Boundary C. at thresholds                         %| 
Ex1=zeros(s);                                      %|
Ex2=zeros(s);                                      %|
for i=1:r-1                                        %|
    for j=1:s                                      %|
        Ex1(j,j)=exp(dnewcell{i}(j,j)*thresholds(i+1));
        Ex2(j,j)=exp(dnewcell{i+1}(j,j)*thresholds(i+1));
    end                                            %|
    mat1=Ex1*vnewcell{i}*diag(R{i});               %|
    mat2=Ex2*vnewcell{i+1}*diag(R{i+1});           %|
    totalmatrix((i-1)*s+1:i*s,i*s+1:(i+1)*s)=-mat1;%|
    totalmatrix(i*s+1:(i+1)*s,i*s+1:(i+1)*s)=mat2; %|
end                                                %|
%Boundary C. at Buffer Length                      %| 
for j=1:s                                          %|
    Ex1(j,j)=exp(dnewcell{r}(j,j)*thresholds(r+1));%|
end                                                %|
    mat1=Ex1*vnewcell{r}*diag(R{r});               %|
    totalmatrix(s*(r-1)+1:s*r,s*r+1:s*(r+1))=mat1; %|
%Normalization C.                                  %|
for i=1:r                                          %|
    for j=1:s                                      %|
        if dnewcell{i}(j,j)==0                     %|
            Ex1(j,j)=thresholds(i);                %|
            Ex2(j,j)=thresholds(i+1);              %|
        else                                       %|
            Ex1(j,j)=exp(dnewcell{i}(j,j)*thresholds(i))/dnewcell{i}(j,j);
            Ex2(j,j)=exp(dnewcell{i}(j,j)*thresholds(i+1))/dnewcell{i}(j,j);
        end                                        %|
    end                                            %|
    mat3=(Ex2-Ex1)*vnewcell{i}*ones(s,1);          %|
    totalmatrix((i-1)*s+1:i*s,r*s+s+1)=mat3;       %|
end                                                %|
totalmatrix(r*s+1:(r+1)*s,r*s+s+1)=ones(s,1);      %|
for i=1:s                                          %|
    if accsbound(i)==0                             %|
        totalmatrix(r*s+i,1:s)=-QT{1}(accsstate(i),:);
    elseif accsbound(i)==r                         %|
        totalmatrix(r*s+i,r*s+1:(r+1)*s)=QT{r+1}(accsstate(i),:);
    else                                           %|
        totalmatrix(r*s+i,accsbound(i)*s+1:(accsbound(i)+1)*s)=-QT{accsbound(i)+1}(accsstate(i),:);
    end                                            %|
end                                                %|
%-   -   -   -   -   -   -   -   -   -   -   -   - %|
%B)Solution:                                       %|
%-----------                                       %|
totalvec=[zeros(1,s*r+s) 1];                       %|
solution=totalvec/totalmatrix;                     %|
%---------------------------------------------------|
disp('PLOTTING')
%-------------------PLOTS---------------------------|
samples=noofsamples+1;                             %|
deltax=thresholds(r+1)/noofsamples;                %|
pdfs=zeros(s,samples);                             %|
totalpdfs=zeros(1,samples);                        %|
for i=1:r                                          %|
    a{i}=solution((i-1)*s+1:i*s);                  %|
end                                                %|
k=1;                                               %|
for j=1:samples                                    %|
    if (j-1)*deltax<=thresholds(k+1)               %|
        pdfs(:,j)=(a{k}*expm(dnewcell{k}*(j-1)*deltax)*vnewcell{k})';
    else                                           %|
        k=k+1;                                     %|
        j=j-1;                                     %|
    end                                            %|
end                                                %|
t=0:deltax:deltax*(samples-1);                     %|
figure(1);                                         %|
hold on                                            %|
for i=1:s                                          %|
    plot(t,pdfs(i,:),'r');                         %|
end                                                %|
totalpdfs=ones(1,s)*pdfs;                          %|
t=0:deltax:deltax*(samples-1);                     %|
figure(2);                                         %|
hold on                                            %|
plot(t,totalpdfs,'r');                             %|
%---------------------------------------------------|
disp('TIME')                                       %|
disp(cputime-ctime)                                %|
disp('PROGRAM TERMINATED')                         %|
disp('---------------')                            %|