terror={};
H={};
error=[];
ly2sd=[];
ly2sv=[];

us  =   0.08;
ys1 =   0.005;
ys21=   .05;
ys22=   2;
ds  =   .05;
ns  =  0.5;

nd = 1;
%load binary_gain
load('Gdataopen_Mary_colas_LV_tobeLin','G1','Gy','G1d','Gyd')
Gy_all = Gy;
Gyd_all = Gyd;
G2  = [1 0; Gy(20,:)];
G2d = [zeros(1,nd); Gyd(20,1:nd)];


% number of runs
mlist =[200, 32];
   
    
for series = 1:150
    
    series
    m = mlist(1);
    lu1d = us*randn(2,m);
    ly1sd = ys1*randn(2,m);
    ly2sd(1,:) = ys21*randn(1,m);
    ly2sd(2,:) = ys22*randn(1,m);
    
    ldd = ds*randn(nd,m);
    lnydd = ns*randn(41,m);
    
    lu1v = us*randn(2,30);
    ly1sv = ys1*randn(2,30);
    ly2sv(1,:) = ys21*randn(1,30);
    ly2sv(2,:) = ys22*randn(1,30);
    ldv = ds*randn(nd,30);
    lnyvv = ns*randn(41,30);
    
for sampletype = 1:2
    m=mlist(sampletype);

    u1d = lu1d(:,1:m);
    y1sd = ly1sd(:,1:m); 
    y2sd = ly2sd(:,1:m); 
    dd = ldd(:,1:m); 
    nydd = lnydd(:,1:m);
    
    u1v = lu1v(:,1:30);
    y1sv = ly1sv(:,1:30); 
    y2sv = ly2sv(:,1:30);
    dv = ldv(:,1:30);
    nyvv = lnyvv(:,1:30);
    
runs = 4:41;
for run = runs
% number of variables    
    ydata_d = [25,15,5,35,20,22,17,27,34,4,14,24,3,37,19,23,6,36,16,26,2,38,21,18,28,12,39,9,29,40,10,30,7,11,31,8,32,13,33,41,1];
        
    %ydata_d = [21 20 22 19 23 18 24 17 25 16 26 15 27 14 28 13 29 12 30 11 31 10 32 9 33 8 34 7 35 6 36 5 37 4 38 3 39 2 40 1 41];
    ydata = ydata_d(1:run);
    Gy = Gy_all(ydata,:);
    Gyd = Gyd_all(ydata,1:nd);

    % model estimators 
    n = size(Gy,1);
    wny  = diag(ns*ones(1,n));
    wu   = diag(us*ones(1,2));
    wd   = diag(ds*ones(1,nd));
    wy1s = diag(ys1*ones(1,2));
    wy2s = diag([ys21, ys22]);


    F = Gyd - Gy*inv(G1)*G1d;
    Y = [G1*inv(G1)*wy1s zeros(2,size(F*wd,2)) zeros(2,size(wny,2))];
    X = [Gy*inv(G1)*wy1s F*wd wny];
    H{1} = Y*pinv(X,0);

    % closed-loop
    F = Gyd - Gy*inv(G1)*G1d;
    Ftilde = [F*wd wny];
    Htemp = pinv(Ftilde*Ftilde')*Gy*inv(Gy'*pinv(Ftilde*Ftilde')*Gy)*G1';
    H{2} = Htemp';

    nyd=nydd(ydata,:);


    % data generation for estimators 

    u = inv(G1)*(y1sd-G1d*dd);
    y1=y1sd;
    y  = Gy*u+Gyd*dd;
    ym = y + nyd;
    y2 = G2*u+G2d*dd;

    % data estimators
    % LS
    H{3} = lsn(y1,ym);

    % H data
    [u,s,v] = svd(y1);
    y1t = y1*v;
    ymt = ym*v;

    p = size(y1t,1);
    G1t = y1t(:,1:p);
    Gyt = ymt(:,1:p);
    Ftilde = ymt(:,p+1:end);
    invF = inv(Ftilde*Ftilde');
    Htemp = invF*Gyt*inv(Gyt'*invF*Gyt)*G1t';
    H{4} = Htemp';

    nchs = 3;
    if run <= nchs
        nch = run;
    else
        nch = nchs;
    end

     Htemp = hdata(y1,pinvn(ym,nch));
    H{5} = Htemp;

%% PCR

A = svd(ym);
pctvar = cumsum(A)./sum(A);
ncpcr = find(pctvar >= .95);
if run <= 4 
    ncpcr = run;
 elseif ncpcr(1) < 4
     ncpcr = 4;
 end
%disp(['nc - pcr: ',num2str(ncpcr(1))])
B=pcrn(y1,ym,ncpcr(1));
H{6} = B;

% pls - analytisk
ncplss = 3;
if run <= ncplss
    ncpls = run;
else
    ncpls = ncplss;
end

[XL,YL,XS,YS,B,PCTVAR,MSE] = plsregress(ym',y1',3);
B=B';
H{7}=B(:,2:end);

% validation of estimators
nH = size(H,2);
nyv=nyvv(ydata,:);

for k = 1:nH
    
    if k == 10
        u = -inv(H{k}*[ones(1,size(Gy,2));Gy])*H{k}*([ones(1,size(dv,2));Gyd*dv+nyv])+inv(H{k}*[ones(1,size(Gy,2));Gy])*y1sv;
    else
        u = -inv(H{k}*Gy)*H{k}*(Gyd*dv+nyv)+inv(H{k}*Gy)*y1sv;
    end
    y1 = G1*u+G1d*dv;
    y  = Gy*u+Gyd*dv;
    ym = y + nyv;
    y2 = G2*u+G2d*dv;
    if k == 10
        yhat{k} = H{k}*[ones(1,size(ym,2));ym];
    else
    yhat{k} = H{k}*ym;
    end
    ycors{k} = y1;
end

for k = 1:nH
    error(k,run)=norm(ycors{k}-yhat{k});
end

end
terror{sampletype}(:,:,series) = error;

end
end
%%
for plotn= 1:2
results = median(terror{plotn},3);
figure(plotn)
pl = plot(runs,results([1,2,3,4,6,7,5],runs));

 set(pl(1), 'LineStyle','--','LineWidth',1.5)
 set(pl(2), 'LineStyle','--','LineWidth',1.5)
 set(pl(3), 'LineWidth',1.5)
 set(pl(4), 'LineWidth',1.5)
 set(pl(5), 'LineWidth',1.5)
 set(pl(6), 'LineWidth',1.5)
 set(pl(7), 'LineWidth',1.5)
 

xlab = xlabel('Number of measurements','FontSize',19);
ylab = ylabel ('Median of 150 prediction error norms','FontSize',19);
axis([0 41 0.004 0.024])
set(gca,'FontSize',14)
grid on



set(ylab,'FontSize',14,'interpreter','tex')
set(xlab,'FontSize',14,'interpreter','tex')

leg = legend('$\mathbf{H}_{2}$','$\mathbf{H}_{4}$','$\mathbf{B}_{\rm{LS}}$','$\mathbf{B}_{\rm{CL}}$','$\mathbf{B}_{\rm{PCR}}$','$\mathbf{B}_{\rm{PLS}}$','$\mathbf{B}^\dagger_{\rm{CL}}$','location','best');
set(leg,'interpreter','latex')
end

