function [K,Y1,Y2] = chebkernel(a,b,c,p)
% P is chebychev polynomial order
% for the kernel value K(u,v) u and v are vectors 
% a is scalar where a= <u,u> or xsup*xsup'
% b is scalar where b= <v,v> or x*x'
% c is scalar where c= <u,v>  or x*xsup'
% K is the output value for the kernel: K(u,v)=K
% This function also returns Y1 and Y2 values (the values before the
% multiplication
%
%
% If you would like to use Second Kind of Chebyshev polynomials in kernel
% functions instead of first kind polynomials, then simply change the code
% on this file: 
% t=chevpoly(p);  <==================== Change this code on this page
%
% t=chevpolyu(p);  <=================== with this one
%
%
% Written by Sedat OZER
% NOTE: The code is not optimized to work faster on MATLAB. 
%
%
% Usage: Use this code (copy and paste) into your svmkernel.m kernel
% function file for generalized Chebyshev kernel function as:
% (where 'eps' is a small constant to avoid division by zero)
%
%     case 'gencheb' 
%     eps=0.0001;
%     nn=length(xsup(1,:));
%     for i=1:n1
%         for j=1:n
%             u=xsup(j,:);
%             v=x(i,:);
%             a=u*u';
%             b=v*v';
%             c=u*v';
%             K(i,j)=(chebkernel(a,b,c,kerneloption))/((sqrt(nn-c)+eps));             
%        end
%     end      
%
%
% Example code piece for "Modified Chebyshev Kernel Function" to be used in svmkernel.m function 
% 
%     case 'modcheb' 
%     gamma=0.01;
%     for i=1:n1
%         for j=1:n
%             u=xsup(j,:);
%             v=x(i,:);
%             a=u*u';
%             b=v*v';
%             c=u*v';
%             K(i,j)=(chebkernel(a,b,c,kerneloption))/(exp((norm(u-v)^2)/gamma));                 
%        end
%     end      
% 
% 
%
%
%NOTE: This function (chebkernel.m) gives you only the numerator part of the kernel.
%Therefore output of this file has to be divided with corresponding
%weighting function :(for gen cheb kernel, you need to use
%((sqrt(nn-c)+eps)) and
%similarly for Modified Cheb kernel you need to use the function
%(exp((norm(u-v)^2))) as the weighting function.
%
% Thank you 
% Sedat OZER
%


if p==0
    K=1;
elseif p==1
    K=c;
else
t=chevpoly(p);
Y1=[];
Y2=[];
for j=1:p+1
    if (t(j,1)~=0)
        r=0;
        e=0;
        for k=1:p+1
            r=r+t(j,k)*a^((k-1)/2);
            e=e+t(j,k)*b^((k-1)/2);
        end
        Y1(j)=r;
        Y2(j)=e;
    else
        r=0;
        e=0;
        for k=1:p+1
            r=r+t(j,k)*a^((k-2)/2);
            e=e+t(j,k)*b^((k-2)/2);
        end
        Y1(j)=r;
        Y2(j)=e*c;
    end
end

% for k=1:p+1
%     if (mod(k,2)==0)
%         Y2(k)=Y2(k)*c;
%     end
% end
Y1=Y1';
Y2=Y2';
K=Y1'*Y2;
end




















