 % Author: Khaled Nakhleh1
%{
Implementation of the Variance-Weighted Deficit (VWD) 
scheduling policy, and the index policy 
%}

%format long

clc, clear all 

%% PARAMETERS %%%%%%%%%%%%
global numClients deficits means clientVars p q deliveriesCounter 
global totDeliveries activations AOIVectors selectedIndices Bn BnSquared
global clientArrivals An runningAOICounter lambdas aoiWeights numberOfTimesZeroActivations
global activationsVectorsChecks

numClients = 10;

% for 5 clients 

%p = [0.8898, 0.328, 0.3951, 0.3334, 0.2964];
%q = [0.7731, 0.6062, 0.3825, 0.543, 0.7491];
%lambdas = [0.086, 0.1113, 0.1373, 0.0552, 0.1505];

% for 10 clients

p = [0.0952, 0.5897, 0.2638, 0.6574, 0.8943, 0.7186, 0.4419, 0.0864, 0.0761, 0.6505];
q = [0.1276, 0.5019, 0.3984, 0.6473, 0.5765, 0.8358, 0.7736, 0.4821, 0.2352, 0.5572];
lambdas = [0.0698, 0.0529, 0.0466, 0.0574, 0.0208, 0.0982, 0.0212, 0.0463, 0.0195, 0.0464];

% for 20 clients
%p = [0.1485, 0.5255, 0.7394, 0.411, 0.5312, 0.7332, 0.8441, 0.4989, 0.6409, 0.3974, 0.7899, 0.1507, 0.5752, 0.8472, 0.4723, 0.6425, 0.2328, 0.6795, 0.8541, 0.7496]
%q = [0.6441, 0.9489, 0.4097, 0.5701, 0.6243, 0.6413, 0.2532, 0.6749, 0.1473, 0.5634, 0.8448, 0.4533, 0.8036, 0.442, 0.6971, 0.2481, 0.2652, 0.9278, 0.1075, 0.2036]
%lambdas = [0.037, 0.0255, 0.0359, 0.0243, 0.0455, 0.0053, 0.0287, 0.0129, 0.0469, 0.0083, 0.0375, 0.027, 0.0249, 0.0452, 0.0238, 0.026, 0.0283, 0.0374, 0.0376, 0.0443] 


% vwd = 1, indexpolicy = 2, 3 = stationary randomized policy  
% 4 = age-based max-weight

selectedSchedulingPolicy = 1
timeslots = 5000;
RUNS = 1000;

seed = 256 % seed value to make sure optimizer gives the same mean,var
rng(seed)

aoiWeights = ones(1, numClients) % equal weights
%aoiWeights = (5-1).*rand(1, numClients) + 1
[ms, var_s, means, clientVars] = optimizer(numClients, p, q, lambdas, aoiWeights);

means;
clientVars;

seeds = linspace(10, 2^32 - 10, RUNS); % seeds for different runs

assert(length(p) == numClients);
assert(length(q) == numClients);
assert(length(lambdas) == numClients);

AOIRuns = zeros(numClients, RUNS); % includes the summed AoI of each run
deliveriesCounter = zeros(numClients, timeslots, RUNS);
runningAOICounter = zeros(numClients, timeslots, RUNS);
numberOfTimesZeroActivations = zeros(1, RUNS);
activationsVectorsChecks = zeros(1, RUNS);
Bn = zeros(numClients, RUNS); % store the summed Bn for each client and run
BnSquared = zeros(numClients, RUNS); % stores Bn^2
averageThroughput = zeros(numClients, RUNS);

%% Start scheduling performance 

for run = 1 : RUNS

rng(seeds(run))
fprintf("current run %d\n", run)
activations = zeros(numClients, timeslots);
clientArrivals = zeros(numClients, timeslots);
deficits = zeros(numClients, timeslots);
totDeliveries = zeros(1, numClients);
AOIVectors = zeros(numClients, timeslots);
An = zeros(numClients, timeslots);
selectedIndices = [];

for i = 1 : numClients
    
    mc = createMC(p(i), q(i));
    walkIns = simulate(mc, timeslots);
    walkIns(walkIns == 2) = 0;
    activations(i, :) = walkIns(1:timeslots);
    clientArrivals(i, :) = setArrivals(lambdas(i), timeslots);
    AOIVectors(i,1) = 1;
end
activationsVectorsChecks(run) = sum(sum(activations, 1) == 0);
for time = 1 : timeslots  % main loop
    
    if selectedSchedulingPolicy == 1
    %disp('vwd selected')
    selectedClient = vwd(time, run);
    updateAn(selectedClient, time);
    updateAOIVals(time, selectedClient);
    elseif selectedSchedulingPolicy == 2
    %disp('index policy selected')
    selectedClient = indexPolicy(time, run);
    updateAn(selectedClient, time);
    updateAOIVals(time, selectedClient);
    
    elseif selectedSchedulingPolicy == 3
    selectedClient = stationaryRandomizedPolicy(time, run);
    updateAn(selectedClient, time);
    updateAOIVals(time, selectedClient);
        
    elseif selectedSchedulingPolicy == 4
    selectedClient = maxWeightPolicy(time, run);
    updateAn(selectedClient, time);
    updateAOIVals(time, selectedClient);
            
    end
end

calculateBn(run);

deficits;
activations;
clientArrivals;
selectedIndices;
AOIVectors = AOIVectors(:, 1:timeslots);
totDeliveries;
An;
Bn;
BnSquared;
averageThroughput(:, run) = totDeliveries;

for x = 1 : numClients
runningAOICounter(x, 1:timeslots, run) = AOIVectors(x,:);
AOIRuns(x, run) = sum(AOIVectors(x,:)); 
end

end

AOIRuns;
runningAOICounter;
aoiWeights;
numberOfTimesZeroActivations
activationsVectorsChecks
runningAOICounter = aoiWeights' .* sum(runningAOICounter, 3);

deliveriesCounter;
averageThroughput;
empiricalAOIs = AOIRuns./timeslots;
averageThroughput = averageThroughput ./ timeslots;

% average AOI from all runs for each client
Bn = sum(Bn, 2)./ RUNS;
BnSquared = sum(BnSquared, 2)./RUNS;

averageThroughput = sum(averageThroughput, 2) ./RUNS;
empiricalAOIs = sum(empiricalAOIs, 2)./ RUNS;


%% storing variance values for all clients
finalVar = zeros(timeslots, numClients); % variance values
finalMean = zeros(timeslots, numClients); % mean values

for client = 1 : numClients
vals = deliveriesCounter(client,:,:);
summedDeliveries = zeros(1, timeslots);
for counter = 1 : RUNS
    summedDeliveries = summedDeliveries + vals(:,:, counter);
end
summedDeliveries;
empiricalMu = summedDeliveries / (RUNS);

for run = 1 : RUNS
for time = 1 : timeslots

    empiricalVar(run, time) = (deliveriesCounter(client,time,run) - empiricalMu(1,time))^2 / time;
    

end
end
empiricalVar;
finalMean(:, client) = empiricalMu;
finalVar(:, client) = sum(empiricalVar, 1) / (RUNS - 1);

end

finalVar;
finalMean;
thepVals = p';
theqVals = q';
thelambdaVals = lambdas';

if selectedSchedulingPolicy == 1 
fileNameVar = sprintf('var_mean/var_vwd_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);
fileNameMean = sprintf('var_mean/mean_vwd_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);
fileNameValues = sprintf('var_mean/values_vwd_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);


elseif selectedSchedulingPolicy == 2
fileNameVar = sprintf('var_mean/var_index_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);
fileNameMean = sprintf('var_mean/mean_index_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);
fileNameValues = sprintf('var_mean/values_index_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);


elseif selectedSchedulingPolicy == 3
fileNameVar = sprintf('var_mean/var_stationary_random_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);
fileNameMean = sprintf('var_mean/mean_stationary_random_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);
fileNameValues = sprintf('var_mean/values_stationary_random_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);


elseif selectedSchedulingPolicy == 4
    
fileNameVar = sprintf('var_mean/var_max_weight_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);
fileNameMean = sprintf('var_mean/mean_max_weight_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);
fileNameValues = sprintf('var_mean/values_max_weight_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);

end
meansVals = means';
clientVarsVals = clientVars';
varTable = table(finalVar);
meanTable = table(finalMean);
valuesTable = table(thepVals, theqVals, thelambdaVals, meansVals, clientVarsVals);
writetable(valuesTable, fileNameValues);
writetable(varTable, fileNameVar);
writetable(meanTable, fileNameMean);


%% Storing output 

averageThroughput
empiricalAOIs
empiricalAOIs = empiricalAOIs .* aoiWeights'
BnSquared = BnSquared';
Bn = Bn';
BnAOI = aoiWeights.*(((BnSquared ./ (2.*Bn)) + 0.5) + (1./lambdas - 1));
BnAOI = BnAOI'
sumOfAvgerageAOIs = (repelem(sum(empiricalAOIs), numClients))';

means = means';
clientVars = clientVars';
ms =  (repelem(ms, numClients))';
varChannel = (repelem(var_s, numClients))';
p = p';
q = q';
sumOfAvgerageAOIs;
means = means'
clientVars = clientVars'
theoreticalAOI = aoiWeights.*((0.5 .* ((clientVars./means.^2) + (1 ./ means) + 1)) + (1./lambdas - 1));
theoreticalAOI = theoreticalAOI'
aoiRatio = repelem(sum(empiricalAOIs) / sum(theoreticalAOI), numClients)'
BnSquared = BnSquared';
Bn = Bn';
means = means';
clientVars = clientVars';
lambdas = lambdas';
results = table(Bn, BnSquared, BnAOI, averageThroughput, empiricalAOIs, theoreticalAOI, p, q, sumOfAvgerageAOIs, lambdas, aoiRatio, means, clientVars, ms, varChannel);

runningAOI = runningAOICounter';
AOIsRunningTable = table(runningAOI);

if selectedSchedulingPolicy == 1 
fileName = sprintf('results/vmd_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);
aoiFileName = sprintf('running_aoi_results/vmd_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);

elseif selectedSchedulingPolicy == 2
fileName = sprintf('results/whittle_index_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);
aoiFileName = sprintf('running_aoi_results/whittle_index_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);

elseif selectedSchedulingPolicy == 3 
fileName = sprintf('results/stationary_random_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);
aoiFileName = sprintf('running_aoi_results/stationary_random_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);

elseif selectedSchedulingPolicy == 4 
fileName = sprintf('results/max_weight_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);
aoiFileName = sprintf('running_aoi_results/max_weight_numClients_%d_timesteps_%d_runs_%d.csv', numClients, timeslots, RUNS);

end
writetable(results, fileName);
writetable(AOIsRunningTable, aoiFileName);

%% storing output for single client


if numClients == 1
    
   singleClientFileName =  sprintf('single_client_results/p_%.2f_q_%.2f_lambda_%.2f_runs_%d.csv', p, q, lambdas, RUNS);
   singleClientTable = table(empiricalAOIs, theoreticalAOI, thepVals, theqVals, thelambdaVals, means, clientVars);
   writetable(singleClientTable, singleClientFileName);
end



%% Functions 

function arrivals = setArrivals(lambda, timeslots)

% set a sequence of packet arrvials based on arrival rate lambda
arrivals = randsrc(1, timeslots, [1 0; lambda 1-lambda]);

end

function calculateBn(currentRun)
global An Bn BnSquared numClients timeslots

BVals = zeros(numClients, timeslots);

for i = 1 : numClients
    AnVals = nonzeros(An(i, :));
    BVals(i, 1) = AnVals(1) - 0;
    for j = 1 : (length(AnVals) - 1)
    BVals(i, j+1) = AnVals(j+1) - AnVals(j);
    end
    Bn(i, currentRun) = sum(BVals(i,:)) / length(nonzeros(BVals(i,:))); 
    BnSquared(i, currentRun) = sum(BVals(i,:).^2) / length(nonzeros(BVals(i,:))); 
end


end

function updateAn(selectedClient, time)

global An numClients

for i = 1 : numClients
if i == selectedClient 
An(selectedClient, time) = time;
else
    An(i, time) = 0;
end
end

end

function mc = createMC(P, G)

TRANS = [1-P, P ;
         G, 1-G];

mc =  dtmc(TRANS, 'StateNames', ["Good", "Bad"]);

end



function selectedClient = vwd(currentTimeSlot, run)

global numClients deficits means clientVars deliveriesCounter
global totDeliveries activations selectedIndices numberOfTimesZeroActivations

gamma = zeros(1, numClients) + -99999999999999999999 ;


for x = 1 : numClients
    if currentTimeSlot == 1
        deficits(x, currentTimeSlot) = 0;   
    else
        deficits(x, currentTimeSlot) = currentTimeSlot * means(x) - totDeliveries(x);
    end

    if activations(x, currentTimeSlot) == 1
    if currentTimeSlot == 1 
        gamma(x) = deficits(x, currentTimeSlot) / sqrt(clientVars(x));
    else 
        gamma(x) = deficits(x, currentTimeSlot - 1) / sqrt(clientVars(x));
    end
    end
end

[maxVal, index] = max(gamma);

if gamma(index) == -99999999999999999999    
    selectedClient = -1;
    numberOfTimesZeroActivations(run) = numberOfTimesZeroActivations(run) + 1;
else
    selectedClient = index;
    totDeliveries(selectedClient) = totDeliveries(selectedClient) + 1;
end


for x = 1 : numClients
    if currentTimeSlot == 1
        if x == selectedClient
        deliveriesCounter(x, currentTimeSlot, run) = 1;
        else
        deliveriesCounter(x, currentTimeSlot, run) = 0;
        end
    else
        if x == selectedClient
        deliveriesCounter(x, currentTimeSlot, run) = deliveriesCounter(x, currentTimeSlot - 1, run) + 1;
        else
        deliveriesCounter(x, currentTimeSlot, run) = deliveriesCounter(x, currentTimeSlot - 1, run);
        end        
    end
end

selectedIndices(currentTimeSlot) = selectedClient;


end

function selectedClient = indexPolicy(currentTimeSlot, run)

global numClients p q selectedIndices deliveriesCounter
global activations AOIVectors totDeliveries numberOfTimesZeroActivations
gamma = zeros(1, numClients) + -99999999999999999999 ; % this is the values we select the client from 

for x = 1 : numClients
    aoi = AOIVectors(x, currentTimeSlot);
    pTrans = 1 - (p(x) / (p(x) + q(x)));
    index = (aoi^2) / 2 - (aoi/2) + aoi/pTrans;
    
    if activations(x, currentTimeSlot) == 1
        gamma(x) = index;
    end
end

[maxVal, index] = max(gamma);

if gamma(index) == -99999999999999999999
    selectedClient = -1;
    numberOfTimesZeroActivations(run) = numberOfTimesZeroActivations(run) + 1;
else
    selectedClient = index;
    totDeliveries(selectedClient) = totDeliveries(selectedClient) + 1;
end


selectedIndices(currentTimeSlot) = selectedClient;


for x = 1 : numClients
    if currentTimeSlot == 1
        if x == selectedClient
        deliveriesCounter(x, currentTimeSlot, run) = 1;
        else
        deliveriesCounter(x, currentTimeSlot, run) = 0;
        end
    else
        if x == selectedClient
        deliveriesCounter(x, currentTimeSlot, run) = deliveriesCounter(x, currentTimeSlot - 1, run) + 1;
        else
        deliveriesCounter(x, currentTimeSlot, run) = deliveriesCounter(x, currentTimeSlot - 1, run);
        end        
    end
end

end

function selectedClient = stationaryRandomizedPolicy(currentTimeSlot, run)

global numClients activations means selectedIndices totDeliveries
global deliveriesCounter numberOfTimesZeroActivations

gamma = zeros(1, numClients);

for x = 1 : numClients
    if activations(x, currentTimeSlot) == 1
        gamma(x) = means(x);
    end
end 


gamma = gamma ./ sum(gamma);


try
selectedClient = randsrc(1,1, [1:1:numClients; gamma]);
totDeliveries(selectedClient) = totDeliveries(selectedClient) + 1;
catch
    selectedClient = -1;
    numberOfTimesZeroActivations(run) = numberOfTimesZeroActivations(run) + 1;
end

selectedIndices(currentTimeSlot) = selectedClient;

for x = 1 : numClients
    if currentTimeSlot == 1
        if x == selectedClient
        deliveriesCounter(x, currentTimeSlot, run) = 1;
        else
        deliveriesCounter(x, currentTimeSlot, run) = 0;
        end
    else
        if x == selectedClient
        deliveriesCounter(x, currentTimeSlot, run) = deliveriesCounter(x, currentTimeSlot - 1, run) + 1;
        else
        deliveriesCounter(x, currentTimeSlot, run) = deliveriesCounter(x, currentTimeSlot - 1, run);
        end        
    end
end
end

function selectedClient = maxWeightPolicy(currentTimeSlot, run)

global numClients activations means deliveriesCounter numberOfTimesZeroActivations
global selectedIndices totDeliveries AOIVectors lambdas

gamma = zeros(1, numClients) + -99999999999999999999;
for x = 1 : numClients
    if activations(x, currentTimeSlot) == 1
        gamma(x) = (1/means(x))*(AOIVectors(x, currentTimeSlot) - (1/lambdas(x)));
    end
end

gamma;

[maxVal, index] = max(gamma);

if gamma(index) == -99999999999999999999
    selectedClient = -1;
    selectedIndices(currentTimeSlot) = selectedClient; 
    numberOfTimesZeroActivations(run) = numberOfTimesZeroActivations(run) + 1;
else
selectedClient = index; 
totDeliveries(selectedClient) = totDeliveries(selectedClient) + 1;
selectedIndices(currentTimeSlot) = selectedClient; 
end


for x = 1 : numClients
    if currentTimeSlot == 1
        if x == selectedClient
        deliveriesCounter(x, currentTimeSlot, run) = 1;
        else
        deliveriesCounter(x, currentTimeSlot, run) = 0;
        end
    else
        if x == selectedClient
        deliveriesCounter(x, currentTimeSlot, run) = deliveriesCounter(x, currentTimeSlot - 1, run) + 1;
        else
        deliveriesCounter(x, currentTimeSlot, run) = deliveriesCounter(x, currentTimeSlot - 1, run);
        end        
    end
end

end

function updateAOIVals(currentTimeSlot, selectedClient)

global numClients AOIVectors clientArrivals


for x = 1 : numClients
    if (x == selectedClient) && (clientArrivals(x, currentTimeSlot) == 1)

        AOIVectors(x, currentTimeSlot+1) = 1;
    
    elseif (x == selectedClient) && (clientArrivals(x, currentTimeSlot) == 0)
        
        timeline = find(clientArrivals(x, (1:currentTimeSlot)) == 1);
        
        try 
            recentPacketTime = currentTimeSlot+1 - timeline(end);
        catch 
           
           recentPacketTime = currentTimeSlot + 1;
        end
            
        AOIVectors(x, currentTimeSlot+1) = recentPacketTime;
    else
        
        AOIVectors(x, currentTimeSlot+1) =  AOIVectors(x, currentTimeSlot) + 1;
    end
    

end


end






