% 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 [Y_calc, error, best_nb] = local_polynomial( X, Y, X_eval, Y_eval , bw_list, param_samplerate)
    Include;
    %           SAMPLE SPLITTING
    % 
    % Anzahl der Datenpunkte für Schritt 1:
    % Wie viele Datensätze sind zum Lernen? aufrunden
    [n d]=size(X);
    n_l=ceil(n*param_samplerate);
    n_t=n-n_l;
    %learning
    X_l=X(1:n_l,:);
    Y_l=Y(1:n_l);
    %testing
    X_t=X(n_l+1:n,:);
    Y_t=Y(n_l+1:n);
    
    % und uebferfluessige Variablen loeschen
    clear X Y;
    
    deg_poly=d;%min(2*d,10)+1;
    
    if d>1
       matrix_monomials=allVL1(d,deg_poly,'<=');
    else
       matrix_monomials=(0:(deg_poly-1))';
    end
    
    number_poly=nchoosek(d+deg_poly , d);
    
    %evaluate X_l,X_t
    X_l_evaluated=Evaluate_monomials(X_l,n_l,matrix_monomials,number_poly,deg_poly,d);
    X_t_evaluated=Evaluate_monomials(X_t,n_t,matrix_monomials,number_poly,deg_poly,d);
    
    
    number_bws=length(bw_list); 
    errors=zeros(deg_poly+1,number_bws);
    %betas=cell(deg_poly+1,number_bws);
    
    %compute distances
    distances=cell(number_bws,1);
    distance=zeros(n_t,n_l);

    for 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
    for nd=0:deg_poly
        number_poly_tmp=nchoosek(d+nd , d);
        B_Matrix=sparse(X_l_evaluated(:,1:number_poly_tmp));
        for nb=1:number_bws
            display(bw_list(nb))
            parfor train_data=1:n_t
                
                K_matrix=diag(sparse(distances{nb}(train_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(train_data)=X_t_evaluated(train_data,1:number_poly_tmp)*beta;

            end
            
            errors(nd+1,nb)=sum(abs(Y_calc-Y_t));       
            display(errors(nd+1,nb))
        end
    
    end
    [dummy_er, dummy_id]=min(errors,[],1);
    
    [~,best_nb]=min(dummy_er);
    best_nd=dummy_id(best_nb)-1;
    
    number_poly=nchoosek(d+best_nd,d);
    n_eval=size(X_eval,1);
    X_eval_evaluated=Evaluate_monomials(X_eval,n_eval,matrix_monomials(1:number_poly,:),number_poly,best_nd,d);
    B_Matrix=sparse(X_l_evaluated(:,1:number_poly));
    
    distance=zeros(n_eval,n_l);
    for 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));
    for eval_data=1:n_eval
                
        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)=X_eval_evaluated(eval_data,:)*beta;

    end

    error=sum(abs(Y_calc-Y_eval));  
    best_nb=bw_list(best_nb);
    
    
end

