n = str2num(input('How many treatments are in the set (n)? ', 's'));
k = str2num(input('How many treatments are included in each block (k)? ', 's'));
lambda = str2num(input('How many observations of each pair do you want (lambda)? ', 's'));
m = str2num(input('How many designs do you want the algorithm to generate before picking the best one (m)? ','s'));

tic
filename = 'ExperimentDesign';
i = 1;
while isfile([filename num2str(i) '.csv'])
    i = i + 1;
end
filename = [filename num2str(i) '.csv'];
f = fopen(filename,'wt');
fprintf(f,'Block,\t');
for i = 1:k
    fprintf(f,'Treat%d,\t',i);
end
fprintf(f,'\n');

for i = 1:m
    comparisons = lambda.* (tril(ones(n,n),-1) + triu(ones(n,n),1)); 
    current_design = zeros(1,k);
    while sum(sum(comparisons)) > 0
        current_block = zeros(1,k);
        for stim = 1:k
            degree = sum(comparisons);
            links = zeros(n,1);
            for s = 1:n
                for s1 = 1:stim-1
                    links(s) = links(s) + comparisons(s,current_block(s1));
                end
            end
            links(current_block(1:stim-1)) = NaN;
            idx = find(links == max(links)); 
            if length(idx) == 1
                 current_block(stim) = idx;
            else 
                idx2 = find(degree == max(degree(idx))); 
                idx3 = intersect(idx,idx2); 
                current_block(stim) = idx3(randi(length(idx3)));
            end    
            for s = 1:stim
                comparisons(current_block(stim),current_block(s)) = max([comparisons(current_block(stim),current_block(s))-1,0]);
                comparisons(current_block(s),current_block(stim)) = max([comparisons(current_block(s),current_block(stim))-1,0]);
            end
        end
        current_design = [current_design; current_block];
    end
    current_design = current_design(2:end,:);
    if i == 1
        best_design = current_design;
    else
        [current_rows, ~] = size(current_design);
        [best_rows,~] = size(best_design);
        if current_rows < best_rows
            best_design = current_design;
        end
    end
end
[best_rows,~] = size(best_design);
for block = 1:best_rows
    fprintf(f,'%d,\t',block);
    for stimulus = 1:k-1
        fprintf(f,'%d,\t',best_design(block,stimulus));
    end
    fprintf(f,'%d\n',best_design(block,k));    
end
fclose(f);
disp(['Saved as ' filename '.']);
toc1 = toc;
disp(['Executed in ' num2str(toc1) ' seconds.']);