% Rohde's Host Specificity Index program

% array dimenstions
maxJ = 10;                 % maximum number of host species
maxI = 10;                  % maximum number of species

% initialize arrays to zero
x     = zeros(maxI, maxJ); % number of parasite species of i'th species in j'th host species
n     = zeros(maxI, maxJ); % number of host individuals
S     = zeros(maxI, 1);    % host specicifity index results
Snorm = zeros(maxI, 1);    % normalised HSI results
r     = zeros(maxI, maxJ); % rank
ratio = zeros(maxI, maxJ); % x(i,j) / n(i,j)
jMax  = zeros(maxI);       % maximum j for each i

% x(i,j) data - (No. of parasited found)%
%x(1,1:10) = [14 12 4379 3 1 2 16 32 6 1]; % A
%x(2,1:7)  = [1 18 4 1 651 13 93]; % B
%x(3,1:5)  = [48 4000 11 2 1]; % C
%x(4,1:2)  = [6 1]; % D
%x(5,1:1)  = [2]; % E
%x(6,1:7)  = [300 205 302 400 650 500 920]; % F
%x(7,1:6)  = [480 802 110 295 105 3]; % G
%x(8,1:2)  = [605 901]; % H
%x(9,1:2)  = [9 78]; % I
%x(10,1:1) = [200]; % J

% x(i,j) data - (No. of hosts infected)
x(1,1:10) = [10 5 11 1 1 2 2 7 6 1]; % A
x(2,1:7)  = [1 1 1 1 10 6 3]; % B
x(3,1:5)  = [32 71 3 2 1]; % C
x(4,1:2)  = [5 1]; % D
x(5,1:1)  = [1]; % E
x(6,1:7)  = [12 1 79 30 8 14 21]; % F
x(7,1:6)  = [41 83 38 68 105 3]; % G
x(8,1:2)  = [12 617]; % H
x(9,1:2)  = [3 7]; % I
x(10,1:1) = [1]; % J

% n(j) data
n(1,1:10) = [21 95 15 7 83 32 84 21 82 64]; % A
n(2,1:7)  = [15 1 83 38 83 143 21]; % B
n(3,1:5)  = [41 83 38 68 143]; % C
n(4,1:2)  = [12 117]; % D
n(5,1:1)  = [3]; % E
n(6,1:7)  = [15 1 83 38 8 14 21]; % F
n(7,1:6)  = [41 83 38 68 143 51]; % G
n(8,1:2)  = [12 617]; % H
n(9,1:2)  = [3 7]; % I
n(10,1:1) = [3]; % J

% maximum j for each i
jMax(1)  = 10;
jMax(2)  = 7;
jMax(3)  = 5;
jMax(4)  = 2;
jMax(5)  = 1;
jMax(6)  = 7;
jMax(7)  = 6;
jMax(8)  = 2;
jMax(9)  = 2;
jMax(10) = 1;

% calculate ratios, x(i,j) / n(i,j)
for i = 1:maxI
    for j = 1:maxJ
        if n(i,j) ~= 0
            ratio(i,j) = x(i,j) / n(i,j);
        end
    end
end
        
% calculate ranks
for i = 1:maxI
    [A, B] = sortrows(-ratio(i,:)');
    [C, D] = sort(B);
    r(i,:) = D';
end

% calculate S(i)

% sum over all i
for i = 1:maxI
    sumNum = 0;
    sumDen = 0;
    % sum over all j
    for j = 1:jMax(i)
        sumNum = sumNum + (x(i,j) / (n(i,j) * r(i,j)));
        sumDen = sumDen + (x(i,j) / n(i,j));
    end
    S(i) = sumNum / sumDen;
end

% display results
S

% calculate normalised results
for i = 1:maxI
    % determine the normalisation factor
    normFact = (1 / jMax(i)) * sum(1 ./ (1:jMax(i)));
    
    % if j is not 1
    if jMax(i) ~= 1
        % normalise S
        Snorm(i) = (S(i) - normFact) / (1 - normFact);
    % otherwise    
    else
        % use the original S(i)
        Snorm(i) = S(i);
    end
end

% display normalised results
Snorm