clear

M = 3; % this is the number of different memories
K1 = 5;
pe1 = 0.25;
lambda1= 0.3;

K2 = 10;
pe2 = 0.25;
lambda2 = 0.4;

K3 = 5;
pe3 = 0.5;
lambda3 = 0.3;

B1 = [ 1 zeros(1,K1-2)];
S01 = [ zeros(K1-2,1); pe1 ]
S1 = zeros(K1-1,K1-1);

for i=1:K1-1
    S1(i,i) = 1-pe1;
    if (i < K1-1)        
        S1(i,i+1) = pe1;
    end
end

B2 = [ 1 zeros(1,K2-2)];
S02 = [ zeros(K2-2,1); pe2 ]
S2 = zeros(K2-1,K2-1);

for i=1:K2-1
    S2(i,i) = 1-pe2;
    if (i < K2-1)        
        S2(i,i+1) = pe2;
    end
end

B3 = [ 1 zeros(1,K3-2)];
S03 = [ zeros(K3-2,1); pe3 ]
S3 = zeros(K3-1,K3-1);

for i=1:K3-1
    S3(i,i) = 1-pe3;
    if (i < K3-1)        
        S3(i,i+1) = pe3;
    end
end

%now build the combined service process
K = K1 + K2 + K3;

S = [ S1 zeros(K1-1, (K2-1)+(K3-1)) ;
      zeros(K2-1, K1-1) S2 zeros(K2-1, K3-1) ;
      zeros(K3-1, (K1-1)+(K2-1)) S3 ]

S0 = [S01; 
      S02;
      S03]
  
B = [ lambda1 zeros(1, K1-2) lambda2 zeros(1, K2-2) lambda3 zeros(1, K3-2)]

S_prime = [ 0 B ;
            S0 S ];
        
Max_index = 120;

Counts = zeros(1, Max_index);

iter = 100000;
        
for i = 1:iter
    curr_count = 0;    
    prob_gone = 0;
    p = 1;
    q = 1;
    done = 0;
    while (done == 0)
        curr_rand = rand(1);        
        if (curr_rand <= (S_prime(p, q) / (1 - prob_gone)))
            %S_prime_val = S_prime(p, q)
            %blah = (S_prime(p, q) / (1 - prob_gone))
            %curr_rand
            %prob_gone
            %p
            %q
            p = q;
            q = 1;
            prob_gone = 0;
            curr_count = curr_count + 1;
            if (p == 1)
                done = 1;
            end
        else
            prob_gone = prob_gone + S_prime(p, q);
            q = q + 1;            
        end
    end
    Counts(curr_count) = Counts(curr_count) + 1;
    if (mod(i,1000) == 0)
        i
    end
        
end

total = 0;
i = 1;
while (total < 0.999)
    Dist(i) = Counts(i) / iter;
    Num(i) = i;
    total = total + Dist(i);
    i = i + 1;    
end        

bar(Num,Dist)

figure;

total = 0;
i = K1;
go = 1;
while ((total < 0.999) && (i < 500) && (go == 1))
    calc = choose(i-1,K1-1)*(1-pe1)^(i-K1)*pe1^(K1);
            
    if ((calc ~= NaN) && (calc ~= Inf))
        total = total + calc;
        go = 1;
        Dist1(i-K1+1) = calc;
        Num1(i-K1+1) = i;
        if ((i-K1+1) == 1)
            Cumlative1(i-K1+1) = 0;
        else
            Cumlative1(i-K1+1) = Cumlative1(i-K1) + calc;
        end
    else
        go = 0;        
        i        
    end;
    i = i+1;
end
bar(Num1,Dist1)

figure;

total = 0;
i = K2;
go = 1;
while ((total < 0.999) && (i < 500) && (go == 1))
    calc = choose(i-1,K2-1)*(1-pe2)^(i-K2)*pe2^(K2);
            
    if ((calc ~= NaN) && (calc ~= Inf))
        total = total + calc;
        go = 1;
        Dist2(i-K2+1) = calc;
        Num2(i-K2+1) = i;
        if ((i-K2+1) == 1)
            Cumlative2(i-K2+1) = 0;
        else
            Cumlative2(i-K2+1) = Cumlative2(i-K2) + calc;
        end
    else
        go = 0;        
        i        
    end;
    i = i+1;
end
bar(Num2,Dist2)

figure;

total = 0;
i = K3;
go = 1;
while ((total < 0.999) && (i < 500) && (go == 1))
    calc = choose(i-1,K3-1)*(1-pe3)^(i-K3)*pe3^(K3);
            
    if ((calc ~= NaN) && (calc ~= Inf))
        total = total + calc;
        go = 1;
        Dist3(i-K3+1) = calc;
        Num3(i-K3+1) = i;
        if ((i-K3+1) == 1)
            Cumlative3(i-K3+1) = 0;
        else
            Cumlative3(i-K3+1) = Cumlative3(i-K3) + calc;
        end
    else
        go = 0;        
        i        
    end;
    i = i+1;
end
bar(Num3,Dist3)
