% estimate predictors using LARS in Reproducing Kernel Hilbert Spaces,
% estimate hyperparametrs using validation data set

% y= output datal
% u = input data;
% yval = validation output datal
% uval = validation input data;
% iy = index of the component for which the predictor is being built
% p = number of eigenfunctions
% T = length of eigenfucntions
% beta = decay rate of eigenfunctions (need to find a smart rule)
% (optional, estimated if not given)


function [Model,beta,sigma,nt,Model2,beta2,sigma2,nt2] = RKHS_Pred_LARS_Validation_all_k(y,u,iy,p,yval,uval,T,beta);

if nargin<8,
    %Tau=[1 2 5 8 12 16 20 25];
    %Tau=[1  3  5  7  9 1 12 15 18 22 26];
    Tau=[0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 6 7 8 9 10 12 15 18 22 26 30];
    beta=1./Tau;
end;


[N,ny]=size(y);
[N,nu]=size(u);

dataval=iddata(yval(:,iy),[yval(:,1:iy-1) yval(:,iy+1:ny) uval]);


esigma=[0:0.3:10];
esigma=[1 log10(50) 2 3 4 7 10 25 50 75 100];
sigma = 10.^esigma;

Model_old=[];
std_old =1000000000000;
beta_old=1;
nt_old=1;
ibeta=0;
cont=1;
nt2=ny+nu;
std2=+inf;
sigma_old=100000;
beta_old=1000000;

while ((ibeta<length(beta)) & cont)
    ibeta=ibeta+1;

    if ((ibeta>2) & (beta_old>beta(ibeta-1)))
        cont=0;
    end


    [X,Phi,Q] = RKHS_Pred_LARS_all_k_computeRegressors_fixed_betaT(y,u,iy,p,beta(ibeta),T);
    isigma=0;
    update_sigma=0;
    cont_sigma=1;
    while ((isigma<length(sigma)) & cont_sigma)
        isigma=isigma+1;
        if ((isigma>1) & (sigma_old < sigma(isigma-1)) & update_sigma)
            cont_sigma=0;
        else
            Model_new = RKHS_Pred_LARS_all_k_fixed_betaT(y,u,X,iy,Phi,Q,sigma(isigma),beta(ibeta));
            %Model_new = RKHS_Pred_LARS_all_k(y,u,iy,p,beta(ibeta),sigma(isigma),T);

            nt=0;
            std_new_prev=+inf;
            cont_nt=1;
            while (cont_nt & nt<ny+nu),
            %for nt=1:nu+ny,
                nt=nt+1;
                yvalh=predict(Model_new{nt},dataval);
                yvalh=yvalh.y;
                eval=yval(:,iy)-yvalh;
                eval = eval(T+1:end);
                std_new=norm(eval)/sqrt(length(eval));
                if std_new > std_new_prev
                    cont_nt=0;
                else
                    std_new_prev=std_new;
                    if ((std_new<std_old)|(((abs(std_old-std_new)/std_new)<0.02)&(nt<nt2))),
                        update_sigma=1;
                        if ((nt<=nt2)|((abs(std2-std_new)/std_new)>0.02))
                            Model2=Model_new{nt};
                            beta2=beta(ibeta);
                            sigma2=sigma(isigma)';
                            std2=std_new;
                            %nt2=nt
                        end;
                        if (std_new<std_old)
                            Model_old=Model_new{nt};
                            beta_old=beta(ibeta);
                            sigma_old=sigma(isigma)';
                            std_old=std_new;
                            %nt_old=nt
                        end;

                    end;
                end;
            end;
            %std_old
%             nt_old
%             sigma_old
%             beta_old
        end;

    end;
end;


Model=Model_old;
beta=beta_old;
sigma=sigma_old
nt=nt_old;


