% 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(y,u,iy,p,yval,uval,T,beta);

if nargin<8,
Tau=[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;

dataval=iddata(yval,uval);

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

esigma=[0:0.3:10];
esigma=[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>1) & (beta_old>beta(ibeta-1)))
        cont=0;
    end
    nt=0;
    update=0;
    while ((nt<ny+nu)&cont)
        nt=nt+1;

        if ((nt_old < nt-1) & update)
            nt=ny+nu;
        else
            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_v2(y,u,iy,p,nt,beta(ibeta),sigma(isigma),T);
                    yvalh=predict(Model_new,dataval);
                    yvalh=yvalh.y;
                    eval=yval-yvalh;
                    eval = eval(T+1:end);
                    std_new=norm(eval)/sqrt(length(eval));
                    if std_new<std_old,
                        update_sigma=1;
                        update=1;
                        
                        if ((nt_old<nt2)|((abs(std2-std_old)/std_old)>0.03))
                        Model2=Model_old;
                        beta2=beta_old;
                        sigma2=sigma_old;
                        std2=std_old;
                        nt2=nt_old;
                        end;
                        Model_old=Model_new;
                        beta_old=beta(ibeta);
                        sigma_old=sigma(isigma);
                        std_old=std_new;
                        nt_old=nt;
                    end;
                end;

            end;
        end;

    end; %cont=1;
end;
Model=Model_old;
beta=beta_old;
sigma=sigma_old;
nt=nt_old;


