function K=K_def(method,par1,par2,par3)
% function K=K_def(method,par1,par2,par3)
% creates a kernel
% K- kernel (structure):
% K.type  - type of the kernel:
%      'opt' - optimal kernel (default)
%      'fun' - external function
%      'str' - string (variable denoted as 'x')
%      'pol' - polynomial kernel
% K.name - string with kernel expression or function name,
%               ignored for optimal and polynomial kernel
% K.coef - coefficients of the optimal or polynomial kernel
% K.support - support of the kernel, default [-1,1]
% K.ny, K.k, K.mi - parameters of the kernel (order etc.)
% K.var  - variance of the kernel
% K.beta - k-th moment of the kernel
%
% pre-defined types: 'epan','quart','rect','trian','gauss'
% call e.g. K=K_def('epan');
%

tol=1000*eps;
default=0;
K=struct('type',char(0),'name',char(0),'coef',[],'support',[-1,1],'ny',NaN,'k',NaN,'mi',NaN,'var',NaN,'beta',NaN);
if nargin==0 default=1;
elseif strncmpi(method,'defa',4) | strncmpi(method,'epan',4)
    default=1;
end
if default
    K.type='opt';
    K.coef=[-0.75,0,0.75];
    K.ny=0;
    K.k=2;
    K.mi=0;
    p1=K.coef;
    p2=conv(p1,p1);
    p3=polyint(p2);
    K.var=polyval(p3,1)-polyval(p3,-1);
    p2=conv(p1,[1 0 0]);
    p3=polyint(p2);
    K.beta=polyval(p3,1)-polyval(p3,-1);
elseif strncmpi(method,'rect',4)
    K.type='opt';
    K.coef=0.5;
    K.ny=0;
    K.k=2;
    K.mi=-1;
    p1=K.coef;
    p2=conv(p1,p1);
    p3=polyint(p2);
    K.var=polyval(p3,1)-polyval(p3,-1);
    p2=conv(p1,[1 0 0]);
    p3=polyint(p2);
    K.beta=polyval(p3,1)-polyval(p3,-1);
elseif strncmpi(method,'tria',4)
    K.type='str';
    K.name='1-abs(x)';
    K.ny=0;
    K.k=2;
    K.mi=0;
    K.var=2/3;
    K.beta=1/6;
elseif strncmpi(method,'quar',4)
    K.type='opt';
    K.coef=conv([-1,0,1],[-1,0,1])*15/16;
    K.ny=0;
    K.k=2;
    K.mi=1;
    p1=K.coef;
    p2=conv(p1,p1);
    p3=polyint(p2);
    K.var=polyval(p3,1)-polyval(p3,-1);
    p2=conv(p1,[1 0 0]);
    p3=polyint(p2);
    K.beta=polyval(p3,1)-polyval(p3,-1);
elseif strcmpi(method,'opt')
    K.type='opt';
    ny=par1;
    k=par2;
    mi=par3;
    K.ny=ny;
    K.k=k;
    K.mi=mi;
    K.coef=K_opt(ny,k,mi);
    p1=K.coef;
    p2=conv(p1,p1);
    p3=polyint(p2);
    K.var=polyval(p3,1)-polyval(p3,-1);
    p2=conv(p1,[1 zeros(1,k)]);
    p3=polyint(p2);
    K.beta=polyval(p3,1)-polyval(p3,-1);
elseif strcmpi(method,'pol')
    K.type='pol';
    P=par1;
    a=par2;
    b=par3;
    ny=0;
    P0=P;
    P1=polyint(P0);
    I=polyval(P1,b)-polyval(P1,a);
    while abs(I)<tol
        ny=ny+1;
        P0=conv(P0,[1,0]);
        P1=polyint(P0);
        I=polyval(P1,b)-polyval(P1,a);
    end
    P=(-1)^ny*prod(1:ny)/I*P;
    K.ny=ny;
    k=ny+1;
    P0=conv(P,[1,zeros(1,k)]);
    P1=polyint(P0);
    I=polyval(P1,b)-polyval(P1,a);
    while abs(I)<tol
        k=k+1;
        P0=conv(P0,[1,0]);
        P1=polyint(P0);
        I=polyval(P1,b)-polyval(P1,a);
    end
    K.k=k;
    K.beta=I;
    P0=conv(P,P);
    P1=polyint(P0);
    I=polyval(P1,b)-polyval(P1,a);
    K.var=I;
    mi=-1;
    P0=P;
    while abs(polyval(P0,a)) < tol & abs(polyval(P0,b)) < tol
        mi=mi+1;
        P0=polyder(P0);
    end
    K.mi=mi;
    K.coef=P;
    K.support=[a,b];
elseif strcmpi(method,'str')
    K.type='str';
    fs=strrep(par1,'.*','*');
    fs=strrep(fs,'./','/');
    fs=strrep(fs,'.^','^');
    K.name=fs;
    K.support=[par2,par3];
    if exist('sym')
        a=par2;
        b=par3;
        ny=0;
        f=sym(fs);
        f0=f;
        I=int(f0,a,b);
        while abs(eval(I))<tol
            ny=ny+1;
            f0=f0*sym('x');
            I=int(f0,a,b);
        end
        f=(-1)^ny*prod(1:ny)/I*f;
	fs=char(f);
        K.name=fs;
        K.ny=ny;
        k=ny+1;
        f0=f*sym(['x^',num2str(k)]);
        I=int(f0,a,b);
        while abs(eval(I))<tol
            k=k+1;
            f0=f0*sym('x');
            I=int(f0,a,b);
        end
        K.k=k;
        K.beta=eval(I);
        f0=f*f;
        I=int(f0,a,b);
        K.var=eval(I);
        mi=-1;
        f0=f;
        fa=limit(f0,a);
        fb=limit(f0,b);
        while (abs(eval(fa)) < tol & abs(eval(fb)) < tol) & mi<20
            mi=mi+1;
            f0=diff(f0);
            fa=limit(f0,a);
            fb=limit(f0,b);
        end
	if mi>=20 mi=inf; end
        K.mi=mi;
    end
elseif strcmpi(method,'fun')
    K.type='fun';
    K.name=par1;
    K.support=[par2,par3];
elseif strncmpi(method,'gaus',4) % Gauss kernel
    K.type='str';
    K.name='1/sqrt(2*pi)*exp(-x^2/2)';
    K.support=[-inf,inf];
    K.ny=0;
    K.k=2;
    K.mi=inf;
    K.var=1/(2*sqrt(pi));
    K.beta=1;
else
    error('Unrecognized parameter')
end
