function K = kmer_wildcard_kernel(S,k,maxw,lambda)

%function K = kmer_wildcard_kernel(S,k,maxw,lambda)
%
% Computes the k-mer wildcard kernel between the strings stored in S;
% w wildcards are allowed. The features are weighted by an
% exponential lambda^-m, where m is the actual number of wildcards for
% that feature.
%
%INPUTS
% S = a cell array containing all strings
% k = an integer indicating the size of the k-mer
% maxw = an integer <k indicating the number of wildcards allowed
% lambda = the forget factor
%
%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).

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_wildcard{j}=[1:length(S{j})-k+1];
        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,0,maxw,lambda,alphabet);
    end
end
if maxw>0
    K=K+lambda*iterates(S,newind_wildcard,newind_wildcard,1,k,1,maxw,lambda,alphabet);
end


function K=iterates(S,indstart,indend,it,maxit,nw,maxw,lambda,alphabet)

%function K=iterates(S,indstart,indend,it,maxit,nw,maxw,lambda,alphabet)
%
% A help function for kmer_wildcard_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
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,nw,maxw,lambda,alphabet);

    end
end
if nw<maxw
    clear indend_wildcard
    for j=1:length(indend)
        indend_wildcard{j}=indend{j}+1;
    end
    K=K+lambda*iterates(S,indstart,indend_wildcard,it+1,maxit,nw+1,maxw,lambda,alphabet);
end