% Program to compute RISK dice probabilities
% Written By - Aaron Clauset aaron at cs dot unm dot edu
% Hosted at  - http://cs.unm.edu/~aaron/
% 
% Disclaimer - This code is distributed AS IS and is neither supported
% nor intended for any commercial use; this code is provided for hobbyist
% purposes only. Enjoy!

defMax = 2;                                       % max number of defending dice
attMax = 3;                                       % max number of attacking dice
numericFlag      = 1;                             % 1/0 en/disables numeric simulation
combinatoricFlag = 1;                             % 1/0 en/disables combinatoric calculation

% The following section of code is a numeric simulation of throwing RISK
% dice to determine the losses of a battle. We simulate 100,000 dice throws
% to get reasonable accuracy, by default
if (numericFlag==1)
numThrows = 100000;                               % change if you'd like

fprintf('Begining numeric simulation with %i throws - this could take a while\n',numThrows);
for DefNum=1:defMax
     ndata{DefNum} = zeros(3,DefNum+1);
     for AttNum=1:attMax		
		fprintf('Simulating %i(d) vs %i(a)\n',AttNum,DefNum);
		for i=1:numThrows
             aDice = sort(ceil(6*rand(AttNum,1))); % attacker throws dice
             dDice = sort(ceil(6*rand(DefNum,1))); % defender throws dice
             dLost = 0;
             for j=1:min([AttNum DefNum])          % pair off high dice
                  if (dDice(end-j+1) < aDice(end-j+1))
                       dLost = dLost + 1;          % add up defender losses
                  end
             end
             ndata{DefNum}(AttNum,dLost+1) = ndata{DefNum}(AttNum,dLost+1) + 1;
		end;
		ndata{DefNum}./numThrows;
     end
     ndata{DefNum} = ndata{DefNum}./numThrows;
end
% Now display the results from the numeric simulation
fprintf('\nDefending with 1 die results (numeric):\n');
fprintf('   \t K_0  \t K_1\n');
for i=1:attMax
     fprintf('A_%i\t%6.4f\t%6.4f\n',i,ndata{1}(i,:));
end;
fprintf('\nDefending with 2 die results (numeric):\n');
fprintf('   \t K_0  \t K_1  \t K_2\n');
for i=1:attMax
     fprintf('A_%i\t%6.4f\t%6.4f\t%6.4f\n',i,ndata{2}(i,:));
end;
fprintf('\n\n');
end

% This section of the code simplifies the combinatorics calculations by
% simply doing the counting after enumerating all the possible combinations
% of dice values for attackers and defenders.

if (combinatoricFlag==1)
fprintf('Begining combinatoric tabulation - this should not take long\n',numThrows);
dicevals = [1 2 3 4 5 6]';
for DefNum=1:defMax
     cdata{DefNum} = zeros(3,DefNum+1);
     for AttNum=1:attMax		
		fprintf('Counting %i(d) vs %i(a)\n',AttNum,DefNum);

          % This loop enumerates all the dice combinations
          perms = dicevals;
          for i=1:DefNum+AttNum-1
               newperms = [];
               for j=1:size(dicevals,1)
                    newperms = [newperms; repmat(dicevals(j),size(perms,1),1)];
               end
               perms = [newperms repmat(perms,6,1)];
          end
          
          % Now tabulate the defender losses for each combination
		for i=1:size(perms,1)
               % pair off high dice and compute differences
               aDice = sort(perms(i,end-AttNum+1:end));
               dDice = sort(perms(i,1:DefNum));
               dLost = 0;
               for j=1:min([AttNum DefNum])          % pair off high dice
                    if (dDice(end-j+1) < aDice(end-j+1))
                         dLost = dLost + 1;          % add up defender losses
                    end
               end
               cdata{DefNum}(AttNum,dLost+1) = cdata{DefNum}(AttNum,dLost+1) + 1;
		end
          cdata{DefNum}(AttNum,:) = cdata{DefNum}(AttNum,:)./size(perms,1);
     end
end

% Now display the results from the numeric simulation
fprintf('\nDefending with 1 die results (combinatoric):\n');
fprintf('   \t K_0  \t K_1\n');
for i=1:attMax
     fprintf('A_%i\t%6.4f\t%6.4f\n',i,cdata{1}(i,:));
end;
fprintf('\nDefending with 2 die results (combinatoric):\n');
fprintf('   \t K_0  \t K_1  \t K_2\n');
for i=1:attMax
     fprintf('A_%i\t%6.4f\t%6.4f\t%6.4f\n',i,cdata{2}(i,:));
end;

meanerr = mean([mean(abs(ndata{1} - cdata{1})) mean(abs(ndata{2} - cdata{2}))]);
fprintf('\nEstimation of mean error in the simulation: %6.4f\n',meanerr);
end;
fprintf('\n\n');
