% estimate predictors using LARS in Reproducing Kernel Hilbert Spaces

% y= output datal
% u = input data;
% iy = index of the component for which the predictor is being built
% p = number of eigenfunctions
% nt = number of nonzero components
% beta = decay rate of eigenfunctions (need to find a smart rule)
% sigma = scalar regularization parameter (RKHS norm, sigma large = little
% regularization)
% T= length of eigenfunctions

function Model = RKHS_Pred_LARS_v2(y,u,iy,p,nt,beta,sigma,T);

%load Eigenfunctions05

dT=1;  % sampling time
if isempty(T),
    T = 50; % length of eigenfunctions (not a bias issue..) 
end;
[Phi lambda] = stable_spline(p,beta,T,dT);

Phi = Phi(:,1:p);
lambda = lambda(1:p);

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

% compute regression matrix;
for j=1:ny,
    for i=1:p,
        X(:,(j-1)*p+i)=filter(Phi(:,i),1,y(:,j)); % assume one delay from all outputs
    end;
end;

for j=1:nu,
    for i=1:p,
        X(:,p*ny+(j-1)*p+i)=filter(Phi(:,i),1,u(:,j)); % assume one delay from all inputs
    end;
end;

%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_fix_k(yGLARS,XGLARS,Nx,nt,Q);
ThetaY = Theta(:,:,1:ny);
ThetaU = Theta(:,:,ny+1:ny+nu);


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 = idpoly([1 -squeeze(ThetaY(:,:,iy))'*Phi'],[Fu],1,1,ones(ny-1+nu,1),[],1);
else
Model = 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,Data);
yh=yh.y;
e=y-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 = idpoly([1 -squeeze(ThetaY(:,:,iy))'*Phi'],[Fu],1,1,ones(ny-1+nu,1),nv,1);
else
Model = idpoly([1 -squeeze(ThetaY(:,:,iy))'*Phi'],[Fy; Fu],1,1,ones(ny-1+nu,1),nv,1);
end;


