function [result, K] = anova(s,t,p)
%ANOVA
%        -ANOVA decomposition kernel using a dynamic alignment algorithm for two vectors
%         of equal length, s and t, and character length, p.
%        -Note that s and t are both row vectors of the same length.
%
%        -The following algorithm is used:
%         K[0](x,y) = 1 for all x,y
%         K[p](s,'empty string') = 0 for all p > 0, all s
%         K[p](s[1:m],t[1:m]) = K[p](s[1:m-1],t[1:m-1]) + (s[m]*t[m]) * K[p-1](s[1:m-1],t[1:m-1])
%
%        -Simply prompting the function will return the value K(s,t), however
%         using the function as [result,K] = K(s,t) will also return the matrix K[p].
%
%        -Example: anova([5 4],[3 2], 2) returns a value of 120.
%            (Note that anova([5 4],[3 2],2)=anova([3 2],[5 4],2) since K(s,t,p) = K(t,s,p) ).
%        -Example: anova([5 4],[3 2], 1) returns a value of 23.
%         
%
%
%USAGE:   scalar = anova([vector 1],[vector 2], p);    (where p is the length of the subsequence)
%
%         [scalar, matrix] = anova([vector 1],[vector 2], p);
%
%

%For more information, visit Http://www.kernel-methods.net/
%Written and tested in Matlab 6.0, Release 12.
%Copyright 2003, Manju M. Pai 4/2003
%manju@kernel-methods.net
%------------------------------------------------------------------------------------------

%Obtain lengths of vectors
[num_rows_s, n] = size(s);
[num_rows_t, m] = size(t);


%Error checking statements:
  %Both vectors need to have the same number of components
  if n ~= m
      error('Error: s and t vectors need to be of the same dimension.');
  end;

  %Make sure input vectors are horizontal.
  if (num_rows_s ~= 1 | num_rows_t ~= 1)  
     error('Error: s and t must be horizontal vectors.');
  end;
  
  %Make sure input vectors consist only of numbers
  if (ischar(s) | ischar(t) )
      error('Error: vectors must consist only of numbers.');
  end;

  %If p is less than zero, program should quit due to faulty variable input.
  if p <= 0
      error('Error: p needs to be greater than 0.');
  end;
%End of error checking


%Initially set every matrix index to -1 to show value has not yet been found
ans = repmat(-1, [n, m, p]);

%Fill in the matrix using the function anova_kernel(s,t,K,p)
for h=1:p
    for i=1:n
        for j=1:m      
            ans(i,j,h) = anova_kernel(s(1:i), t(1:j), ans, h);
        end;
    end;
end;

result = ans(n,m,p);
K = ans(:,:,p);


%------------------------------------------------------------------------------------------

function ans = anova_kernel(sa, ta, K, p)
%This function is called by anova(s,t,p).
%Type 'help anova' for a description of the program.
%

%------------------------------------------------------------------------------------------


%Obtain lengths of both strings
n = length(sa);
m = length(ta);

%truncate last character of string
s = sa(1:n-1);
t = ta(1:m-1);
    


%Start algorithm:

  % 1) Split main algorithm into two parts:
    % a) K[p](s,t)
  
      if (length(s) == 0) | (length(t) == 0)
        %This is a base case where 0 is returned if either vector has 0 components
        ans = 0;
      elseif( K( length(s), length(t) , p) == -1 )
        % Value has not yet been calculated
        ans = anova_kernel(s,t,K,p);
      else
        % Value has already been calculated
        ans = K( length(s), length(t), p);
      end;

    % b) K[p-1](s,t) * k(sa,ta)

      %First calculate k(sa,ta).
      little_k = sa(n) * ta(m);
        
      %Now calculate K[p-1](s,t)
      if (p-1) == 0
        %This is a base case where 1 should always be returned if p = 0;
        result = 1;
      elseif (length(s) == 0) | (length(t) == 0)
        %This is a base case where 0 is returned if either vector has 0 components
        result = 0;
      elseif( K( length(s) , length(t) , p-1) == -1)
        % Value has not yet been calculated
        result = anova_kernel( s , t, K, p-1 );
      else
        % Value has already been calculated.
        result = K( length(s), length(t), p-1 );
      end;
        
      ans = ans + (little_k * result);

    return
% End of algorithm