MarketShare

function [] = demo()
    clc; clear all;
    %% Market share draft
    
    Lambdaiwant = 1.501;
    
    % Number of points
    npoints = 100;
    
    % Number of muN settings
    nmuN = 1;
    
    
    MSfast = zeros(nmuN,npoints);
    
    MSfastub = zeros(nmuN, npoints);
    MSfastlb = zeros(nmuN, npoints);
    Lambdavalues = zeros(nmuN, npoints);
    
    muNsettings = [1.1];
    
    settingnumber = 0;
    for muN=muNsettings
    
        settingnumber = settingnumber+1
    
        a = 0.5;
    
        low = 0.05;
        high = (1 + muN)*0.94;
    
        stepsize = (high - low)/(npoints - 1);
    
    
        count = 0;
    
        % Original Tmin
        Tmin = 10;
        for Lambda = low : stepsize : high
            count = count + 1
        
            [pis, iMax, jMax, Tmin] = StationaryProbabilities( muN, Lambda, a, 20);
            [lambdaslow, lambdafast] = FindArrivalRates( Lambda, iMax, jMax, muN, a, pis );
            Lambdavalues(settingnumber,count) = Lambda;
            MSfast(settingnumber,count) = lambdafast/Lambda;
    
        end;
    
    end;
end


function [ pis, iMax, jMax, T] = StationaryProbabilities( muN, Lambda, a, Tmin)
%% This script finds the stationary probabilities when both A and N announce. Without loss of generality, let \mu_A = 1 and \mu_N take a higher value
% T is the threshold that we move to attain an epsilon condition
% pis is a matrix of stationary probabilities imax X jmax. i indexes slow server A, j
% indexes fast server N.

% Direction = -1 implies muN is barely less than its value. +1 implies muN is
% barely more than its value.

if(Lambda/(1+muN) > 0.95)
    epsilon = 1e-3;
else
    epsilon = 1e-4  ;
end;

% What is the minimum I should allow T to be?
% Tmin = 20;

%Lambda = Lambda - 0.0001;
% Initialize
condition = 1;
T = Tmin-1;
% This would go into a loop
while(condition > epsilon)
    T = T+1;
    iMax = T;
    jMax = ceil(T * muN);
    % Solve Ax = b
    A = zeros((iMax+1) * (jMax+1), (iMax+1) * (jMax+1));
    b = zeros((iMax+1) * (jMax+1),1);

    % Put in the coefficients for A - one for each state (except 0,0) and one
    % for normalization
    eqcount = 0;
    for i=0:iMax
        for j=0:jMax
            % Ignore state (0,0)
            if(i==0 && j==0)
                continue;
            end;

            % Otherwise
            eqcount = eqcount + 1;
            b(eqcount) = 0;

            % Find the coefficients, which depend on the values of i and
            % j
            if(i==0)
                A(eqcount,Scalarize(i,j,jMax+1)) = Lambda + muN;
                A(eqcount,Scalarize(i+1,j,jMax+1)) = -1;
                if(j==1)
                    A(eqcount,Scalarize(0,0,jMax+1)) = -Lambda * a;
                end;
                if(j<muN*T)
                    A(eqcount,Scalarize(i,j+1,jMax+1)) = -muN;
                end;
                continue;
            end;

            if(j==0)
                A(eqcount,Scalarize(i,j,jMax+1)) = Lambda + 1;
                A(eqcount,Scalarize(i,j+1,jMax+1)) = -muN;
                if(i==1)
                    A(eqcount,Scalarize(0,0,jMax+1)) = -Lambda *(1-a);
                end;
                if(i<T)
                    A(eqcount,Scalarize(i+1,j,jMax+1)) = -1;
                end;            
                continue;
            end;

            if(i==T && j==muN*T)
                A(eqcount,Scalarize(i,j,jMax+1)) = 1 + muN;
                A(eqcount,Scalarize(i,j-1,jMax+1)) = -Lambda;
                A(eqcount,Scalarize(i-1,j,jMax+1)) = -Lambda;
                continue;
            end;

            if(j==muN*T)
                A(eqcount,Scalarize(i,j,jMax+1)) = Lambda + muN + 1;
                A(eqcount,Scalarize(i+1,j,jMax+1)) = -1;
                A(eqcount,Scalarize(i-1,j,jMax+1)) = -Lambda;
                
                if(j-1 == muN * i)
                    A(eqcount,Scalarize(i,j-1,jMax+1)) = -Lambda*a;
                elseif(j-1 < muN * i)
                    A(eqcount,Scalarize(i,j-1,jMax+1)) = -Lambda;
                end;
                
                %{
                if(direction == -1)
                    if(j-1 <= muN*i)
                        A(eqcount,Scalarize(i,j-1,jMax+1)) = -Lambda;
                    end;
                elseif(direction == 1)
                    if(j-1 < muN*i)
                        A(eqcount,Scalarize(i,j-1,jMax+1)) = -Lambda;
                    end;
                end;
               %}
            
                continue;
            end;

            if(i==T)
                A(eqcount,Scalarize(i,j,jMax+1)) = Lambda + muN + 1;
                A(eqcount,Scalarize(i,j+1,jMax+1)) = -muN;
                A(eqcount,Scalarize(i,j-1,jMax+1)) = -Lambda;
                
                if(j == muN * (i-1))
                    A(eqcount,Scalarize(i-1,j,jMax+1)) = -Lambda*(1-a);
                elseif(j > muN * (i-1))
                    A(eqcount,Scalarize(i-1,j,jMax+1)) = -Lambda;
                end;
                
                %{
                if(direction == -1)
                    if(j >= muN*(i-1))
                        A(eqcount,Scalarize(i-1,j,jMax+1)) = -Lambda;
                    end;
                elseif(direction == 1)
                    if(j > muN*(i-1))
                        A(eqcount,Scalarize(i-1,j,jMax+1)) = -Lambda;
                    end;
                end;
                %}
                continue;
            end;

            % Otherwise, it is a generic state (i,j)
            A(eqcount,Scalarize(i,j,jMax+1)) = Lambda + muN + 1;
            A(eqcount,Scalarize(i+1,j,jMax+1)) = -1;
            A(eqcount,Scalarize(i,j+1,jMax+1)) = -muN;
            
            if(j > muN*(i-1))
                    A(eqcount,Scalarize(i-1,j,jMax+1)) = -Lambda;
            elseif (j == muN*(i-1))
                    A(eqcount,Scalarize(i-1,j,jMax+1)) = -(1-a) * Lambda;
            end;
            
            if(j-1 < muN*i)
                    A(eqcount,Scalarize(i,j-1,jMax+1)) = -Lambda; 
            elseif(j-1 == muN * i)
                    A(eqcount,Scalarize(i,j-1,jMax+1)) = -a * Lambda;
            end;        
            
            %{
            if(direction==-1)
                if(j >= muN*(i-1))
                    A(eqcount,Scalarize(i-1,j,jMax+1)) = -Lambda;
                end;
                if(j-1 < muN*i)
                    A(eqcount,Scalarize(i,j-1,jMax+1)) = -Lambda; 
                end;        
            elseif(direction==1)
                if(j > muN*(i-1))
                    A(eqcount,Scalarize(i-1,j,jMax+1)) = -Lambda;
                end;
                if(j-1 <= muN*i)
                    A(eqcount,Scalarize(i,j-1,jMax+1)) = -Lambda; 
                end; 
            end;
            %}
        end;
    end;

    % Add the normalization equation
    eqcount = eqcount + 1;
    A(eqcount,:) = 1;
    b(eqcount) = 1;

    x = A\b;

    % Squarify
    pis = vec2mat(x,jMax+1);

    condition = sum(pis(iMax+1,:))+sum(pis(:,jMax+1));
    %condition = pis(iMax+1, jMax+1);
end;
%T
end