m = 15000
l = 256
c = 8

erasures = 4
iterations = 4
tests = 1000

% Generate input messages
messages = random('unid',l,m,c);

% Convert into sparse messages
n = l * c;
sparsemessages = zeros(m,n);
for i=1:m
    for j=1:c
        sparsemessages(i,(j-1)*l+messages(i,j)) = 1;
    end
end

% Create network
network = zeros(n,n);
for i=1:m
    for j=1:c
        for k=1:c
            network((j-1)*l+messages(i,j),(k-1)*l+messages(i,k)) = 1;
        end
    end
end

% Compute density
density = (sum(sum(network))-sum(diag(network)))/(l*l*c*(c-1))

% Compute error rate
score = 0;
for t=1:tests
    init = transpose(sparsemessages(random('unid',m),:));
    indexes = randperm(c);
    input = init;
    for i=1:erasures
        input((indexes(i)-1)*l+1:indexes(i)*l) = zeros(l,1);
    end
    for iter=1:iterations
        prod = network*(input);
        out = zeros(0);
        for i=1:c
            maxi = max(prod((i-1)*l+1:i*l));        
            out = cat(1,out,prod((i-1)*l+1:i*l) == maxi);
        end
        input = out;
    end
    if isequal(init,input) 
        score = score + 1;
    end
end
errorrate = 1-score/tests
