clear

% here, output is p and input is m

ordervec = [3;3;4;5;8;3;4;6];
% fid = fopen(' hihi.txt','w');

for ii = 1:8
  
  if ii == 1
    load CD_player_arm.dat;
    U = CD_player_arm(:,1:2);
    Y = CD_player_arm(:,3:4);
    m = 2; p = 2;
    N = 200; N1 = 600;		% Size of training set used by Liu & Vandenberghe
    u = U(1:N,:)';  uv = U(1:N1,:)';
    bx = Y(1:N,:)'; bxv = Y(1:N1,:)';
  end
  if ii == 2
    load cstr.dat;
    U = cstr(:,2);  
    Y = cstr(:,3:4);
    m = 1;  p = 2;
    N = 250;	N1 = 700;	% Size of training set used by Liu & Vandenberghe
    u = U(1:N,:)';  uv = U(1:N1,:)';
    bx = Y(1:N,:)'; bxv = Y(1:N1,:)';
  end
  if ii == 3
    load dryer.dat;
    U = dryer(:,1);  
    Y = dryer(:,2);
    m = 1;  p = 1; 
    N = 150; N1 = 400;	% Size of training set used by Liu & Vandenberghe
    u = U(1:N,:)';  uv = U(1:N1,:)';
    bx = Y(1:N,:)'; bxv = Y(1:N1,:)';
  end
  if ii == 4
    load exchanger.dat;
    U = exchanger(:,2);
    Y = exchanger(:,3);
    m = 1;  p = 1; 
    N = 400;	N1 = 1000;	% Size of training set used by Liu & Vandenberghe
    u = U(1:N,:)';  uv = U(1:N1,:)';
    bx = Y(1:N,:)'; bxv = Y(1:N1,:)';
  end
  if ii == 5
    load thermic_res_wall.dat;
    U = thermic_res_wall(:,1:2);
    Y = thermic_res_wall(:,3);
    m = 2;  p = 1; 
    N = 400;	N1 = 1000;	% Size of training set used by Liu & Vandenberghe
    u = U(1:N,:)';  uv = U(1:N1,:)';
    bx = Y(1:N,:)'; bxv = Y(1:N1,:)';
  end
  if ii == 6
    load winding.dat;
    U = winding(:,1:5);
    Y = winding(:,6:7);
    m = 5;  p = 2; 
    N = 250;	 N1 = 600;	% Size of training set used by Liu & Vandenberghe
    u = U(1:N,:)';  uv = U(1:N1,:)';
    bx = Y(1:N,:)'; bxv = Y(1:N1,:)';
  end
  if ii == 7
    load glassfurnace.dat;
    U = glassfurnace(:,2:4);
    Y = glassfurnace(:,5:10);
    m = 3;  p = 6; 
    N = 150;	N1 = 400;%Size of training set used by Liu & Vandenberghe
    u = U(1:N,:)'; uv = U(1:N1,:)'; 
    bx = Y(1:N,:)';bxv = Y(1:N1,:)'; 
  end
  if ii == 8
    load dryer2.dat;
    U = dryer2(:,2:4);
    Y = dryer2(:,5:7);
    m = 3;  p = 3; 
    N = 250;	N1 = 600;	%Size of training set used by Liu & Vandenberghe
    u = U(1:N,:)';  uv = U(1:N1,:)';
    bx = Y(1:N,:)'; bxv = Y(1:N1,:)';
  end
  
  fprintf(' Problem %d\n',ii);
  
  data = iddata(bx',u'); % add Ts?
  
  % default order

  output = n4sid(data,'best','nk',zeros(1,m),'InitialState','Estimate','CovarianceModel','None','Focus','Stability','DisturbanceModel','None');
  A = output.A;
  B = output.B;
  C = output.C;
  D = output.D;
  x0 = output.x;
  
  x = x0;
  y = zeros(p,N1);
  for i = 1:N1
    xnew = A*x + B*uv(:,i);
    y(:,i) = C*x + D*uv(:,i);
    x = xnew;
  end
  
 %% identification error
  ymeas = bx;
  ybarI = (ymeas*ones(N,1))/N;
  ybarI = repmat(ybarI,N,1);
  ierror = norm(reshape(y(:,1:N) - ymeas,N*p,1),2)/norm(reshape(ymeas,N*p,1) - ybarI,2);

%% Validation error
  ymeas = bxv;
  ybarI = (ymeas*ones(N1,1))/N1;
  ybarI = repmat(ybarI,N1,1);
  verror = norm(reshape(y - ymeas,N1*p,1),2)/norm(reshape(ymeas,N1*p,1)- ybarI,2);
  
  fprintf(' order = %d ierror = %g, verror = %g\n', size(A,1), ierror, verror);
  

  % nuclear norm input order
  
  output = n4sid(data,ordervec(ii),'nk',zeros(1,m),'InitialState','Estimate','CovarianceModel','None','Focus','Stability','DisturbanceModel','None');
  A = output.A;
  B = output.B;
  C = output.C;
  D = output.D;
  x0 = output.x;
  
  x = x0;
  y = zeros(p,N1);
  for i = 1:N1
    xnew = A*x + B*uv(:,i);
    y(:,i) = C*x + D*uv(:,i);
    x = xnew;
  end
  
 %% identification error
  ymeas = bx;
  ybarI = (ymeas*ones(N,1))/N;
  ybarI = repmat(ybarI,N,1);
  ierror = norm(reshape(y(:,1:N) - ymeas,N*p,1),2)/norm(reshape(ymeas,N*p,1) - ybarI,2);

%% Validation error
  ymeas = bxv;
  ybarI = (ymeas*ones(N1,1))/N1;
  ybarI = repmat(ybarI,N1,1);
  verror = norm(reshape(y - ymeas,N1*p,1),2)/norm(reshape(ymeas,N1*p,1)- ybarI,2);
  
  fprintf(' order = %d ierror = %g, verror = %g\n', size(A,1), ierror, verror);
  fprintf('\n');
end
% fclose(fid)