% 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] = additive_b_splines(X_l, Y_l, X_t, Y_t, X_eval, Y_eval, alpha_list)
    
    X_l=X_l+1;
    X_t=X_t+1;
    X_eval=X_eval+1;
    
    [n_l d]=size(X_l);
    n=n_l+size(X_t,1);  
    
    number_alphas=length(alpha_list);
    deg_list=0:3;
    
    number_deg=length(deg_list);
    
    [Splines]=Make_additive_splines(alpha_list,deg_list);
    
    errors=zeros(number_alphas,1);
    best_deg=zeros(number_alphas,1);
    betas=cell(number_alphas,number_deg);
    
    %best choice
    for na=1:number_alphas
        
        error_dummy=zeros(number_deg,1);
        
        parfor ndeg=1:number_deg
        %display(ndeg)
        B_Matrix=sparse(Make_additive_B_Matrix([X_l;X_t],d,Splines{na,ndeg}));
        
        warning('off','MATLAB:rankDeficientMatrix');
        betas{na,ndeg}=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,ndeg};
        
        error_dummy(ndeg)=sum((Y_calc-Y_t).^2);       
        
        end
        
        [errors(na) best_deg(na)]=min(error_dummy);
        
    end
    
    [~, best_na]=min(errors);
    B_calc=sparse(Make_additive_B_Matrix(X_eval,d,Splines{best_na,best_deg(best_na)}));

    Y_calc=max(B_calc*betas{best_na,best_deg(best_na)},0);%max(B_calc*beta(best_na),0);
        
    result.mse=mean((Y_calc-Y_eval).^2);
    result.res=Y_calc-Y_eval;
    result.param=[alpha_list(best_na),deg_list(best_deg(best_na))];
    
    
end

