%-----------------------------------------------|
%--SINGLE REGIME - FIRST ORDER - FINITE BUFFER--|
%--------------ZERO RATE INCLUDED---------------|
%-----------------------------------------------|

%---------Output: Statewise CDF Plot------------|
%Author: H. Emre KANKAYA
ctime=cputime;
%-----------------PARAMETERS--------------------|
B=1;%Buffer Length                             %|
slottime=1/24;
Q=[-3,1,1,1;1,-3,1,1;1,1,-2,0;1,1,1,-3];       %|

rates=[-1,0,0,3];                              %|
noofsamples=1000;%For Plotting                 %|
%-----------------------------------------------|

%-------------------PREPARE---------------------|
n=length(rates);%number of states              %|
%-   -   -   -   -   -   -   -   -   -   -   - %|
%A)Find Zero Rate States:                      %|
%------------------------                      %|
noofzero=0;                                    %|
for i=1:n                                      %|
    if rates(i)==0                             %| 
        noofzero=noofzero+1;                   %| 
    end                                        %| 
end                                            %| 
%-   -   -   -   -   -   -   -   -   -   -   - %|
%B)Pi(steady-state) Distribution:              %|
%---------------------------------             %|
Qnew=[Q ones(n,1)];                            %|
colnew=[zeros(1,n), 1];                        %|
sol=colnew/Qnew;%steady state distribution     %|
sol2=zeros(1,n-noofzero);                      %|
ratesnz=zeros(1,n-noofzero);                   %|
j=1;                                           %|
for i=1:n                                      %|
    if ~(rates(i)==0)                          %|
        ratesnz(j)=rates(i);                   %|
        sol2(j)=sol(i);                        %|
        j=j+1;                                 %|
    end                                        %|
end                                            %|
sol=sol2;                                      %|
%-   -   -   -   -   -   -   -   -   -   -   - %|
%C)New Q & New R:                              %|
%------------------                            %|
R=diag(ratesnz);                               %| 
Qb=zeros(noofzero);                            %|
Qa=zeros(n-noofzero,noofzero);                 %|
Q0=zeros(noofzero,n-noofzero);                 %|
Qnz=zeros(n-noofzero,n-noofzero);              %|
ib=1;                                          %|
ia=1;                                          %|
i0=1;                                          %|
inz=1;                                         %|
for i=1:n                                      %|
    jb=1;                                      %|
    ja=1;                                      %|
    j0=1;                                      %|
    jnz=1;                                     %|
    for j=1:n                                  %|   
        if rates(i)==0&rates(j)==0             %|
            Qb(ib,jb)=Q(i,j);                  %|
            jb=jb+1;                           %|
        elseif ~(rates(i)==0)&rates(j)==0      %|
            Qa(ia,ja)=Q(i,j);                  %|
            ja=ja+1;                           %|
        elseif rates(i)==0&~(rates(j)==0)      %|
            Q0(i0,j0)=Q(i,j);                  %|
            j0=j0+1;                           %|
        else                                   %|
            Qnz(inz,jnz)=Q(i,j);               %|
            jnz=jnz+1;                         %|
        end                                    %|
    end                                        %|
    if rates(i)==0                             %|
        ib=ib+1;                               %|
        i0=i0+1;                               %|
    else                                       %|
        ia=ia+1;                               %|
        inz=inz+1;                             %|
    end                                        %|
end                                            %|
Q=Qnz-Qa*(Qb)^-1*Q0;                           %|
n=n-noofzero;                                  %|
rates=ratesnz;                                 %|
%D)Find Eigenvalues and Eigenvectors:          %|
%------------------------------------          %|
[v,d]=eig(Q',R');                              %|
%# of 0 eigenvales=1                           %|
%# of - eigenvalues=noofL                      %|
%# of + eigenvalues=noofU-1                    %|                                                   
%-   -   -   -   -   -   -   -   -   -   -   - %|
%B)Reorder the eigenvalues and eigenvectors:   %|
%-------------------------------------------   %|
[y,index]=sort(diag(d));                       %|
dnew=zeros(n,n);%ordered eigenvalues           %|
vnew=zeros(n,n);%ordered eigenvectors          %|
for i=1:n                                      %|
    dnew(i,i)=d(index(i),index(i));            %|
    vnew(i,:)=(v(:,index(i)))';                %|
end                                            %|
%-----------------------------------------------|



%-------------------SOLUTION--------------------|
%A) Find # of + and - States:                  %|
%----------------------------                  %|
noofL=0;%# of states with - rate               %|
for i=1:n                                      %|
    if rates(i)<0                              %|
        noofL=noofL+1;                         %|
    end                                        %|
end                                            %| 
Lindex=zeros(1,noofL);                         %|
z=1;                                           %|
for i=1:n                                      %|
    if rates(i)<0                              %|
        Lindex(z)=i;                           %|
        z=z+1;                                 %|
    end                                        %|
end                                            %|
noofU=n-noofL;%# of states with + rate         %|
Uindex=zeros(1,noofU);                         %|
z=1;                                           %|
for i=1:n                                      %| 
    if rates(i)>0                              %|
        Uindex(z)=i;                           %|
        z=z+1;                                 %|
    end                                        %|
end                                            %|
%-   -   -   -   -   -   -   -   -   -   -   - %|
%B) Boundary Conditions:                       %|
%-----------------------                       %|
%Boundary Conditions at zero,                  %|
%for states with + rate                        %|
vu=zeros(n,noofU);                             %|
for i=1:noofU                                  %|
    vu(:,i)=vnew(:,Uindex(i));                 %|
end                                            %|
%Boundary Conditions at B                      %|
%for states with - rate                        %|
vLex1=zeros(n);                                %|
for i=1:n                                      %|
    vLex1(i,:)=vnew(i,1:n)*exp(dnew(i,i)*B);   %|
end                                            %|
vLex2=zeros(n,noofL);                          %|
for i=1:noofL                                  %|
    vLex2(:,i)=vLex1(:,Lindex(i));             %|
end                                            %|
%-   -   -   -   -   -   -   -   -   -   -   - %|
%C)Solution:                                   %|
%-----------                                   %|
totalmatrix=[vLex2 vu];                        %|
totalvec=[ sol(1:noofL) zeros(1,noofU)];       %|
sol2=totalvec/totalmatrix;                     %|
%-----------------------------------------------|
%-------------------PLOTS-----------------------|
no=noofsamples-2;                              %|
deltay=B/no;                                   %|
Fg=zeros(noofsamples+1,n+noofzero);            %|
for i=1:no+1                                   %|
    for j=1:n                                  %|
        Fg(i,1:n)=Fg(i,1:n)+sol2(j)*vnew(j,:)*exp(dnew(j,j)*(i-1)*deltay);
    end                                        %|
end                                            %|
for i=no+1:noofsamples+1                       %|
        Fg(i,1:n)=sol;                         %|
end                                            %|
if noofzero>0                                  %|
    Fg(:,n+1:n+noofzero)=-Fg(:,1:n)*Qa*Qb^-1;  %|
end                                            %|
i=0:deltay:noofsamples*deltay;                 %|
figure(1);                                     %|
hold on;                                       %|
plot(i,Fg)                                     %|
str='';                                        %|
hold on;                                       %|
for k=1:n                                      %|
    str1=strcat('state',num2str(k));           %|
    str=strvcat(str,str1);                     %|
end                                            %|
legend(str,0);                                 %|
figure(2);                                     %|
plot(i,Fg*ones(n+noofzero,1))                  %|
%-----------------------------------------------|

%-------------------OUTPUTS---------------------|    
disp('--------------OUTPUTS---------------')%|
disp('# OF STATES WITH POSITIVE RATE:')        %|
disp(noofL)                                    %|
disp('# OF STATES WITH NEGATIVE RATE:')        %|
disp(noofU)                                    %|
disp('COEFFICIENTS:')                          %|
disp(sol2)                                     %|
disp('ACCUMULATIONS AT 0:')                    %|
disp(Fg(1,:))                                  %|
disp('ACCUMULATIONS AT B:')                    %|
disp(Fg(no+1,:)-Fg(no,:))                      %|
disp('SYSTEM IS SOLVED IN(sec.)')              %|
disp(-ctime+cputime)                            %|
disp('---------------------------------------')%|
%-----------------------------------------------|