% a trial run of the barrel model
% this is a very simple program that generates the diagrams appearing the first part. If you have MATLAB installed on your
% computer and interested in the topic, you can explore the model changing its parameters. You can also change the model
% itself and play with it to test its assumptions, conclusions, and general validity. I will try to answer all the reasonable
% questions. I can only repeat that the algorithm really is extremely simple.
clear
clc
depth = 19; % the depth of the trees
p = 1/3; % a node rho
epoch = 240;
rlimit = 0.5;
perfs = zeros (depth, epoch); % performances
heap = zeros (depth, epoch); % accrues reserves
t = 2 ^ (depth + 1) - 1; % tree size
reserves = zeros (1, t); % reserves per node
n = 2 ^ depth; % leaf row size
rho = makeRho (depth, p); % random binary rho
fuse = rho2Leaf (rho); % fusion leaf
perm = randperm (n); % permutate the links
demand = fuse(perm); % initial demand leave
lap = 15; % check lap to estimate epoch duration
t1 = datetime ('now');
for i = 1 : epoch
if i == lap % display the time estimate
t2 = datetime ('now');
dura = (t2 - t1) * epoch / lap;
disp (dura) % epoch estimate hours:minutes:seconds
% epoch parameters
fprintf ('Epoch \t\t\t\t\t%u\n', epoch)
fprintf ('Depth \t\t\t\t\t%u\n', depth)
fprintf ('Reserve limit \t\t\t%.1f\n', rlimit)
end
supply = demand; % copy supply from past demand
x = supply; % propagating supply row
perm = randperm (n); % permutated links per cycle
demand = fuse(perm); % next cycle demand
y = demand; % propagating demand row
for d = depth : -1 : 1
m = length (x); % current row length
r = reserves(m:2*m-1); % current row reserves
u = x + r; % total available row supply
ratio = u ./ y; % row supply / row demand
perf = min (ratio); % current level performance
perfs(d,i) = perf; % historic performance
reserve = u - perf * y; % unused supply = reserve
reserve(reserve < 0) = 0; % beautify
% calculate heap by cutoff row reserve
[val,ind] = sort (reserve); % large to trim
cum(2:m) = cumsum (val(1:end-1)); % large to trim
q = m : -1 : 1; % in reverse order
w = cum(1:m) + val .* q; % accumulated reserve
k = find (w > rlimit, 1, 'first'); % cutoff ref
if ~isempty (k) % limit the reserve
t = (rlimit - cum(k)) / q(k);
val(k:end) = t; % trim surplus reserve
reserve = val(ind); % adjusted reserve
end
%
reserves(m:2*m-1) = reserve; % adjusted reserves
heap(d,i) = sum (reserve); % historic reserve
x = upRow (x); % next supply level
y = upRow (y); % next demand level
end
end
reserve = zeros (depth, 1); % final reserves per level
for d = 1 : depth
s = 2 ^ d;
reserve(d) = sum (reserves(s:2*s-1));
end
% mean performance, reserve per level
% in the second half of the epoch
meanPerf = mean (perfs(:,round(end/2):end), 2);
meanPerf(meanPerf > 1) = 1;
meanHeap = mean (heap(:,round(end/2):end), 2);
formatSpec = '\t%2u)\t%6.3f\t\t%6.2f\n';
fprintf ('\nPerformance, Reserves:\n')
for i = 1 : depth
fprintf (formatSpec, i, meanPerf(i), meanHeap(i))
end
figure (1)
bar (meanPerf)
figure (2)
bar (meanHeap)
function rho = makeRho (depth, ratio)
% create and randomly assign nodes' rho given depth, ratio
n = 2 ^ depth - 1; % node count
rho = randi (2, 1, n); % uniform 1 and 2
rho(rho==1) = ratio;
rho(rho==2) = 1 - ratio;
function leaf = rho2Leaf (rho)
% calculate the leaf row from the nodes' rho
m = length (rho) + 1; % leaf length
leaf = zeros (1, m);
leaf(1) = 1;
s = 1; % start
e = 1; % end
while s < m
rho(s:e) = rho(s:e) .* leaf(1:s); % left leaves
leaf(2:2:2*s) = leaf(1:s) - rho(s:e);
leaf(1:2:2*s) = rho(s:e);
s = e + 1; % new start
e = s + e; % new end
end
function up = upRow (row)
% calculate the nodes one level up
m = length (row);
assert (m > 1)
up = row(1:2:m); % odd nodes
up = up + row(2:2:m); % even nodes
Comments