function score = getAlignmentScore(s1, s2, SM, GS)

%GETALIGNMENTSCORE  Compute score of a sequence alignment
%   score = getAlignmentScore(s1, s2, SM, GS), where
%     s1 and s2 are aligned sequences (must have the same length)
%     SM is the scoring matrix
%     GS is the gap score

if (length(s1) ~= length(s2))
   error('Aligned sequences must have the same length');
end

score = 0;
for i=1:length(s1)
   if (or(s1(i) == '-', s2(i) == '-'))
      score = score + GS;
   else
      score = score + lookupScore(s1(i), s2(i), SM);
   end
end



function score = lookupScore(a1, a2, SM)

% Returns the score of two amino acids.
% Assumes that the scoring matrix SM is indexed in the standard 
% amino acid order: that is, ARNDCQEGHILKMFPSTWYV

aa_abbrev = 'ARNDCQEGHILKMFPSTWYV';
i = findstr(aa_abbrev, upper(a1));
j = findstr(aa_abbrev, upper(a2));
score = SM(i,j);