% Preprocessing für den Partitionsschätzer. Hier werden die Daten
% dimensionsweise auf das Intervall [0,1] gelegt. Nur die Extrema an beiden
% Rändern (in 1/100) werden aus dem Intervall gefiltert.

function [result] = local_polynomial_kernel(X_l, Y_l, X_t, Y_t, X_eval, Y_eval, bw_list)
    
%     X_l=X_l+1;
%     X_t=X_t+1;
%     X_eval=X_eval+1;
    
    [n_t d]=size(X_t);
    n_l=size(X_l,1);
        
    max_deg=1;
    number_poly=(max_deg+1)^d;
    matrix_monomials=combinator(max_deg+1,d,'p','r')-1;
    
    %evaluate X_l,X_t
    X_l_evaluated=Evaluate_monomials(X_l,n_l,matrix_monomials,number_poly,max_deg,d);
    X_t_evaluated=Evaluate_monomials(X_t,n_t,matrix_monomials,number_poly,max_deg,d);
    
    
    number_bws=length(bw_list); 
    errors=zeros(1,number_bws);
    %betas=cell(1,number_bws);
    
    %compute distances
    %distances=cell(number_bws,1);
    distance=zeros(n_t,n_l);

    parfor j=1:n_l    
        for t=1:d
            distance(:,j)=distance(:,j)+(X_t(:,t)-X_l(j,t)).^2;
        end   
    end 
    
%     for nb=1:number_bws
%         distances{nb}=exp(-distance/bw_list(nb));
%     end
    Y_calc=zeros(n_t,1);
    %best choice
    B_Matrix=sparse(X_l_evaluated);
    for nb=1:number_bws
        distance_nb=exp(-distance/bw_list(nb));
        %display(nb)
        parfor train_data=1:n_t
            %display(train_data)
            K_matrix=diag(sparse(distance_nb(train_data,:)));
            warning('off','MATLAB:rankDeficientMatrix');
            betas=spqr_solve (B_Matrix'*K_matrix*B_Matrix,B_Matrix'*K_matrix*Y_l,struct ('solution','min2norm'));
            warning('on','MATLAB:rankDeficientMatrix');
            Y_calc(train_data)=max(X_t_evaluated(train_data,:)*betas,0);

        end

        errors(nb)=sum((Y_calc-Y_t).^2);       
        %display(errors(nb))
    end
    
    [~,best_nb]=min(errors);
    
    n_eval=size(X_eval,1);
    X_eval_evaluated=Evaluate_monomials(X_eval,n_eval,matrix_monomials,number_poly,max_deg,d);
    B_Matrix=sparse(X_l_evaluated);
    Y_calc=zeros(n_eval,1);
    
    distance=zeros(n_eval,n_l);
    parfor j=1:n_l    
        for t=1:d
            distance(:,j)=distance(:,j)+(X_eval(:,t)-X_l(j,t)).^2;
        end   
    end 
    distance=exp(-distance/bw_list(best_nb));
    
    parfor eval_data=1:n_eval
        %display(eval_data)        
        K_matrix=diag(sparse(distance(eval_data,:)));
        warning('off','MATLAB:rankDeficientMatrix');
        beta=spqr_solve (B_Matrix'*K_matrix*B_Matrix,B_Matrix'*K_matrix*Y_l,struct ('solution','min2norm'));
        warning('on','MATLAB:rankDeficientMatrix');
        Y_calc(eval_data)=max(X_eval_evaluated(eval_data,:)*beta,0);

    end

    result.mse=mean((Y_calc-Y_eval).^2);
    result.res=Y_calc-Y_eval;
    result.param=bw_list(best_nb);
    
    
end

