% 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_na] = b_splines( X, Y, X_eval, Y_eval , alpha_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);   
    
    number_alphas=length(alpha_list);
    dim_list=d*ones(1,number_alphas);
    Lb_list=zeros(1,number_alphas);
    Ub_list=ones(1,number_alphas);
    M_list=2*(ones(1,number_alphas));

    [Splines Combinations]=Make_splines(dim_list,alpha_list,M_list,Lb_list,Ub_list);

    % und uebferfluessige Variablen loeschen
    clear X Y;
    
    errors=zeros(number_alphas,1);
    betas=cell(number_alphas,1);
    
    %best choice
    for na=1:number_alphas
        display(alpha_list(na))
        
        B_Matrix=sparse(Make_B_Matrix([X_l;X_t],dim_list(na),Splines{na},Combinations{na}));
        
        warning('off','MATLAB:rankDeficientMatrix');
        betas{na}=spqr_solve (B_Matrix(1:n_l,:)'*B_Matrix(1:n_l,:),B_Matrix(1:n_l,:)'*Y_l,struct ('solution','min2norm')); 
        warning('on','MATLAB:rankDeficientMatrix');
        
        Y_calc=B_Matrix((n_l+1):n,:)*betas{na};
       
        errors(na)=sum(abs(Y_calc-Y_t));       
        display(errors(na))
    end
    
    [~, best_na]=min(errors);
    B_calc=sparse(Make_B_Matrix(X_eval,dim_list(best_na),Splines{best_na},Combinations{best_na}));

    Y_calc=B_calc*betas{best_na};%max(B_calc*beta(best_na),0);
    error=sum(abs(Y_calc-Y_eval));  
    best_na=alpha_list(best_na);
    
    
end

