function K = kmer_kernel(S,k)

%function K = kmer_kernel(S,k)
%
% Computes the k-mer kernel between the strings stored in S
%
%INPUTS
% S = a cell array containing all strings
% k = an integer indicating the size of the k-mer
%
%OUTPUTS
% K = the kmer kernel evaluated between all strings
%
%
%For more info, see www.kernel-methods.net
%
%Author: Tijl De Bie, march 2004.

% Bug fixed october 2004 (thanks to Neil Lawrence).
% Added a test for speed-up in the iterates function, july 2005.

nS=length(S);

alphabet=[];
for j=1:nS
    for i=1:length(S{j})
        alphabet=union(alphabet,S{j}(i));
    end
end

na=length(alphabet);

K=zeros(nS,nS);
for i=1:na
    i
    present=[];
    clear newS newind
    ct=0;
    for j=1:nS
        newind_store=[];
        newind_store=find(S{j}(1:end-k+1)==alphabet(i));
        if ~isempty(newind_store)
            ct=ct+1;
            present=[present ; j];
            newS{ct}=S{j};
            newind{ct}=newind_store;
        end
    end
    if ~isempty(present)
        K(present,present)=K(present,present)+iterates(newS,newind,newind,1,k,alphabet(i),alphabet);
    end
end


function K=iterates(S,indstart,indend,it,maxit,a,alphabet)

%function K=iterates(S,indstart,indend,it,maxit,a,alphabet)
%
% A help function for kmer_kernel.m
%
%
%Author: Tijl De Bie, march 2004.

na=length(alphabet);
nS=length(S);

if it==maxit
    for i=1:nS
        number(i,1)=length(indend{i});
    end
    K=number*number';
    return
elseif length(indstart)==1 && length(indstart{1})==1
    K=1;
%    disp(['done the speed-up at depth ' num2str(it)])
    return
end

K=zeros(nS,nS);
for i=1:na
    ct=0;
    present=[];
    clear newindstart newindend newS
    for j=1:nS
        newindstart_store=[];
        newindend_store=[];
        for k=1:length(indstart{j})
            if S{j}(indend{j}(k)+1)==alphabet(i)
                newindstart_store=[newindstart_store ; indstart{j}(k)];
                newindend_store=[newindend_store ; indend{j}(k)+1];
            end
        end
        if ~isempty(newindend_store)
            ct=ct+1;
            newindstart{ct}=newindstart_store;
            newindend{ct}=newindend_store;
            present=[present ; j];
            newS{ct}=S{j};
        end
    end
    if ct~=0
        K(present,present)=K(present,present)...
            +iterates(newS,newindstart,newindend,it+1,maxit,alphabet(i),alphabet);
    end
end