%-----------------------------------------------------|
%---------------FLUID SIMULATION-2--------------------|
%------------Steady State Distribution----------------|              
%-----------------------------------------------------|

%DESCRIPTION
%------------------------------------------------------

%This program simulates first order feedback fluid queues
%with multiple thresholds and plots the steady state pdf 
%and cdf of the system concerned.

%Time epochs, when the system touches to the thresholds or  
%when the state changes, are set as the steps of the simul.
%Look at fluidsimul1.m for another (slower) simulation.

%All sorts of boundaries can be simulated.
%------------------------------------------------------

clear;
a=cputime;
disp('----------------------')
disp('SIMULATION-2 STARTED...')
disp('INITIALIZATIONS');
%PARAMETERS FOR THE SIMULATION
%------------------------------------------------------
simtime=30000;     % simulate system for simtime secs.
deltax=0.003;       %histogram intervals
deltat=0.001;      %time intervals


%PARAMETERS FOR THE SYSTEM
%------------------------------------------------------
%EXAMPLE 1
thresholds=[0,2,3]; %terminal and intermediate boundary points
Q{1}=[-4,1,1,1,1;1,-4,1,1,1;1,1,-4,1,1;1,1,1,-4,1;1,1,1,1,-4];%infinitesimal generator in region 1          
Q{2}=[-4,1,1,1,1;1,-4,1,1,1;1,1,-4,1,1;1,1,1,-4,1;1,1,1,1,-4];%infinitesimal generator in region 2                    
QT{1}=Q{1};
QT{2}=Q{1};
QT{3}=Q{2};
R{1}=[1,-2,2,-4,0];%rates in region 1: states 1 and 2 respectively     
R{2}=[-1,1,3,-1,2];
RT{1}=[1,0,2,0,0];
RT{2}=[0,0,2,-4,0];
RT{3}=[-1,0,0,-1,0];

%SYSTEM IDENTIFICATION
%------------------------------------------------------
N=length(thresholds)-1;%number of regions
S=length(R{1});          %number of states


%INITIALIZATIONS&DEFINITIONS
%------------------------------------------------------
t=0;                %time parameter
cstate=1;           %current state 1->s
cregime=1;          %current regime 1->N if inregion=1;currentregime=1->N+1 if inregion=0
coccupancy=0;       %current occupancy THRESHOLDS(1)->THRESHOLDS(3)
change=0;           %the change of queue level up to the next epoch
schange=0;          %1 if state changes in the next epoch,0 o.w.
touch=0;            %1 if queue level touches to a threshold,0 o.w.
touchtime=0;        %time of that touch
noofinter=ceil(thresholds(N+1)/deltax)+1;%# of space intervals 
ss=zeros(S,noofinter);% vector holds histogram data
teski=0;            %time of the last epoch
oeski=0;            %occupancy at the last epoch
reski=1;            %regime at the last epoch
statemat=zeros(2,S-1);
inregion=0;%inregion=0 if the occupancy is at a threshold, 1 otherwise
for k=1:S    
    for j=1:N
        cQ=Q{j};
        statematI{j}{k}=zeros(2,S-1);
        for i=1:S
            if i<k
                statematI{j}{k}(1,i)=-cQ(k,i)/cQ(k,k);
                statematI{j}{k}(2,i)=i;
            elseif i>k
                statematI{j}{k}(1,i-1)=-cQ(k,i)/cQ(k,k);
                statematI{j}{k}(2,i-1)=i;
            end
        end
        for i=2:S-1
            statematI{j}{k}(1,i)=statematI{j}{k}(1,i)+statematI{j}{k}(1,i-1);
        end    
    end
end
for k=1:S
    for j=1:N+1
        cQ=QT{j};
        statematT{j}{k}=zeros(2,S-1);
        for i=1:S
            if i<k
                statematT{j}{k}(1,i)=-cQ(k,i)/cQ(k,k);
                statematT{j}{k}(2,i)=i;
            elseif i>k
                statematT{j}{k}(1,i-1)=-cQ(k,i)/cQ(k,k);
                statematT{j}{k}(2,i-1)=i;
            end
        end
        for i=2:S-1
            statematT{j}{k}(1,i)=statematT{j}{k}(1,i)+statematT{j}{k}(1,i-1);
        end
    end
end
%------------------------------------------------------
disp('SIMULATING THE SYSTEM')
while t<simtime
    %UPDATE PARAMETERS    
    teski=t;
    oeski=coccupancy;    
    reski=cregime; 
    %FIND THE RATE OF THE SYSTEM IN THIS INTERVAL
    if inregion==1%intermediate points        
        rate=R{cregime}(cstate);
        cQ=Q{cregime};
        statemat=statematI{cregime}{cstate};
        r1=exprnd(-1/cQ(cstate,cstate));%time of the next state change
        change=r1*rate;
            %A) IF THE QUEUE STAYS IN THE SAME REGIME
        if change+coccupancy<thresholds(1+cregime)&change+coccupancy>thresholds(cregime)
            coccupancy=coccupancy+change;        
            schange=1;
            %B) IF THE QUEUE TOUCHES TO THE LOWER THRESHOLD          
        elseif change+coccupancy<=thresholds(cregime)         
        %i)if current regime is 1 
        %     a)queue level>0, it touches to 0,and stays there until a state change        
            if RT{cregime}(cstate)==0
                touch=1;
                touchtime=(thresholds(cregime)-coccupancy)/rate;                
                coccupancy=thresholds(cregime);
                schange=1;
                inregion=0;
            else                                 
                r1=(thresholds(cregime)-coccupancy)/rate;
                coccupancy=thresholds(cregime);
                inregion=0;
            end
            %C) IF THE QUEUE TOUCHES TO THE UPPER THRESHOLD          
        elseif change+coccupancy>=thresholds(cregime+1)            
            if RT{cregime+1}(cstate)==0
                touch=1;
                touchtime=(thresholds(cregime+1)-coccupancy)/rate;                
                coccupancy=thresholds(cregime+1);
                schange=1;
                inregion=0;
                cregime=cregime+1;
            else                                 
                r1=(thresholds(cregime+1)-coccupancy)/rate;
                coccupancy=thresholds(cregime+1);
                inregion=0;
                cregime=cregime+1;
            end
        end       
                                    
    else%inregion==0 boundary points        
        rate=RT{cregime}(cstate);
        cQ=QT{cregime};
        statemat=statematT{cregime}{cstate};
        if rate>0
           r1=deltat;           
           coccupancy=coccupancy+rate*deltat;           
           inregion=1;    
        elseif rate<0
           r1=deltat;
           cregime=cregime-1;
           coccupancy=coccupancy+rate*deltat;
           inregion=1;
        else%rate=0            
           schange=1;
           r1=exprnd(-1/cQ(cstate,cstate));%time of the next state change
        end            
    end
     t=t+r1;      
ostate=cstate;
    %FIND THE NEXT STATE OF THE SYSTEM:
    if schange==1        
        ostate=cstate;
        r2=rand;%uniform random number used to find the next state 
        k=1;%temp variable, next state
        while k<S
            prob=statemat(1,k);        
            if r2<prob;
                    cstate=statemat(2,k);
                    k=S+1;
            end            
            k=k+1;
        end
        schange=0;
    end       
    %CONSTRUCT THE HISTOGRAM BY SAMPLING THE QUEUE CONTENT CHANGE
    %A)If the queue level does not touch to terminal boundary points
    if touch==0
        m=oeski;    
        while teski<=t    
            m=m+rate*deltat;
            if m>thresholds(N+1)
                m=thresholds(N+1);
            elseif m<0
                m=0;
            end            
            c=ceil(m/deltax)+1;           
            ss(ostate,c)=ss(ostate,c)+1;
            teski=teski+deltat;
        end
        
    %A)If the queue level touches to terminal boundary points    
    else
        t2=teski+touchtime;    
        m=oeski;
        while teski<t2           
            m=m+rate*deltat;
            if m>thresholds(N+1)
                m=thresholds(N+1);
            elseif m<0
                m=0;
            end
            c=ceil(m/deltax)+1;                        
            ss(ostate,c)=ss(ostate,c)+1;
            teski=teski+deltat;
        end
        m=coccupancy;     
        c=ceil(m/deltax)+1;
        d=ceil((t-teski)/deltat)+1;        
        ss(ostate,c)=ss(ostate,c)+d;        
        touch=0;
    end
    
end

%PLOTS&OUTPUTS
%------------------------------------------------------
ss2=ss;
for i=2:noofinter
    ss2(:,i)=ss2(:,i-1)+ss2(:,i);
end
ss2=ss2/sum(ss2(:,noofinter));
tt=0:deltax:(noofinter-1)*deltax;
figure(4);
hold on
plot(tt,ss2);
figure(3);
hold on
plot(tt,ones(1,S)*ss2);

figure(2);
ss3=ss./sum(sum(ss))/deltax;
plot(tt,ones(1,S)*ss3)
figure(1);
hold on;
for i=1:S
    plot(tt,ss3(i,:));
end
disp('SIMULATION TIME:')
disp(cputime-a);
disp('SIMULATION-2 ENDED...')
disp('----------------------')