% estimate predictors using LARS in Reproducing Kernel Hilbert Spaces


% use pre-computed regressor matrices for fixed beta, T
% y= output data
% u = X regressors;
% iy = index of the component for which the predictor is being built
% nu = number of inputs
% sigma = scalar regularization parameter (RKHS norm, sigma large = little
% regularization)

function Model = RKHS_Pred_LARS_all_k_fixed_betaT(y,u,X,iy,Phi,Q,sigma,beta);


Q = 1/sigma*Q; % for computing the norm in the RKHS
Qnorm=1*Q;
p=size(Phi,2);
T=size(Phi,1);
[N,ny]=size(y);
[N,nu]=size(u);
Nx = p*ones(1,nu+ny);

%keyboard

trans = min(T,floor(1/beta*log(100))); % neglect transient 

yGLARS = y(trans+2:end,iy);
XGLARS = X(trans+1:end-1,:);
lambdaRR = zeros(sum(Nx),sum(Nx));
for i=1:length(Nx),
    lambdaRR((i-1)*p+1:i*p,(i-1)*p+1:i*p)=Qnorm;
end;

[yGLARS XGLARS]=prepare_data_LARS_RR(yGLARS,XGLARS,Nx,lambdaRR);

Theta = Group_Lars_all_k(yGLARS,XGLARS,Nx,Q);

for nk=1:nu+ny, % models with increasing number of non-zero components

    ThetaY = Theta(:,:,1:ny,nk);
    ThetaU = Theta(:,:,ny+1:ny+nu,nk);


    for i=1:iy-1,
        Fy(i,:) = [0 squeeze(ThetaY(:,:,i))'*Phi'];
    end;

    for i=iy+1:ny,
        Fy(i-1,:) = [0 squeeze(ThetaY(:,:,i))'*Phi'];
    end;

    for i=1:nu,
        Fu(i,:) = [0 squeeze(ThetaU(:,:,i))'*Phi'];
    end;





    if ny==1
        Model{nk} = idpoly([1 -squeeze(ThetaY(:,:,iy))'*Phi'],[Fu],1,1,ones(ny-1+nu,1),[],1);
    else
        Model{nk} = idpoly([1 -squeeze(ThetaY(:,:,iy))'*Phi'],[Fy; Fu],1,1,ones(ny-1+nu,1),[],1);
    end;

    Data = iddata(y(:,iy),[y(:,1:iy-1) y(:,iy+1:end) u]);

    yh=predict(Model{nk},Data);
    yh=yh.y;
   
    e=y(:,iy)-yh;
    nv=std(e(T+1:end))^2;

    % figure
    % plot(y(T+1:end)),hold  on, plot(yh(T+1:end),':r')
    % pause

    if ny==1
        Model{nk} = idpoly([1 -squeeze(ThetaY(:,:,iy))'*Phi'],[Fu],1,1,ones(ny-1+nu,1),nv,1);
    else
        Model{nk} = idpoly([1 -squeeze(ThetaY(:,:,iy))'*Phi'],[Fy; Fu],1,1,ones(ny-1+nu,1),nv,1);
    end;

end;

