function [id, Slength] = findSeqSimilarity(S,level);
% findSeqSimilarity is used to find the similarity between all sequences
% pairs stored in the sequence alignment S.
% If nessesary the amino acid must be converted from a character to an 
% integer representation. Gaps in both seqences are not counted.
% Program by Steinar Thorvaldsen, april 2004. Latest modified 4. Dec. 2005.

% Input
% S: character or integer array containing the sequence from M genes
% level: vecor containing 2D or 3D information about the sequence

% Output
% id: M x M x levels matrix with the identity ratio of each of the M sequences.
% Slength: M x levels matrix with the number of aa in each of the M sequences.
% When the input level is omitted, data for the overall sequences are found

% If nessesary, convert alignment from a character to an integer representation:
if ischar (S)
    S = aa2int(S); 
end
% To determine the number of sequences M, and length N:
[M,N] = size (S);

% The function is made flexible conserning input parameters:
if nargin < 2
    level=ones(N,1,'int8'); % no levels
    D3=1;
else
    D3=max(level);
    %The number of levels may be defined globally in the main script 
    %DeltaProt.m to stay constant throught all seq analysis:
    %global maxlevel;
    %D3=maxlevel;
end

if N > length(level)
   error('Sequences and level vector must be same length!')
end;

id = zeros (M, M, D3);
Slength = zeros (M, D3);

match = zeros (M, M, D3);
totals = zeros (M, M, D3);    
for n=1:N %loop through all sequence sites
  d = level(n);
  for m1=1:M
    aaNo1 = S(m1,n);
    if aaNo1 >0 && aaNo1 <21
        Slength(m1,d)=Slength(m1,d)+1; %one aa found in seq
    end
    for m2=1:M   
       aaNo2 = S(m2,n);
       if aaNo1 >0 && aaNo2 >0 && (aaNo1 <21 | aaNo2 <21) %at least one aa
       %if aaNo1 >0 && aaNo2 >0 && aaNo1 <21 && aaNo2 <21 %two aa (no gaps)
           totals(m1,m2,d) = totals(m1,m2,d) + 1;
           if aaNo1==aaNo2 %match
               match(m1,m2,d) = match(m1,m2,d) + 1;
           end
       end        
    end
  end
end

id = match./totals;