function [result] = kernel(X_l, Y_l, X_t, Y_t, X_eval, Y_eval , bandwiths_list)
   
    number_bandwiths=length(bandwiths_list);
    
    errors=zeros(number_bandwiths,1);
    
    %best choice
    parfor bw=1:number_bandwiths
        
        %calclulate nearest neighbour estimate
        Y_calc=kernel_estimate(X_t,X_l,Y_l,bandwiths_list(bw));
        errors(bw)=sum((Y_calc-Y_t).^2);       
        %display(errors(bw))
    end
    
    [~, best_bw]=min(errors);
    Y_calc=kernel_estimate(X_eval,X_l,Y_l,bandwiths_list(best_bw));
    
    result.mse=mean((Y_calc-Y_eval).^2);
    result.res=Y_calc-Y_eval;
    result.param=bandwiths_list(best_bw);
    
end

function [Y_calc]=kernel_estimate(X_t,X_l,Y_l,bw)
    
    [n,d] = size(X_t);
    m=size(X_l,1);
    Y_calc=zeros(n,1);
    %go through testing data A and compute distances to all learning data B
    for l=1:n
        
        distance=zeros(m,1);
        for t=1:d
            distance=distance+(X_l(:,t)-X_t(l,t)).^2;
        end
        
        weights=exp(-distance/bw);
        denominator=sum(weights);
        if denominator>0
           Y_calc(l)=sum(weights.*Y_l)/denominator;
        end
       
    end
    
end





