% This code calls sysHankelgp_cv.m to solve the nuclear norm minimization problem
% using accelerated gradient projection algorithm, and then computes
% identification and verification errors for various parameters.

%% Generate input

clear

% Read U and bx from benchmark problem

ii = 4; % change this number to specify one of the 8 problems

if ii == 1 % eyeball
  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,:)';
%   sigmatol = 3e-2;
end
if ii == 2 % eyeball
  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,:)';
%   sigmatol = 1e-3;
end
if ii == 3 % eyeball
  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,:)';
%   sigmatol = 1e-3;
end
if ii == 4 % tradeoff not clear, hard to eyeball. Use threshold
  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,:)';
  sigmatol = 1e-4; 
end
if ii == 5 % eyeball
  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,:)';
%   sigmatol = 1e-3;
end
if ii == 6 % eyeball
  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,:)';
%   sigmatol = 1e-2;
end
if ii == 7 % tradeoff not clear, hard to eyeball. Use threshold
  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,:)'; 
  sigmatol = 1e-2;
end
if ii == 8 % eyeball
  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,:)';
%   sigmatol = 1e-3;
end

%% parameters for the algorithm

tol = 1e-4;
freq = 10;
mu = sort(logspace(1,-4,100)');
lamb = 1./mu; 
r = floor((N+1)/(p+m+1)); % r in H_r, Used in Vandenberghe

U = Hrx(u,r); % block hankel input matrix
U = null(U); % Uperp            
Ebx = Hrx(bx,r)*U;   
c = svd(Ebx)/max(svd(Ebx)); % ratio of singular values of Ebx

l1 = size(Ebx,1);
l2 = size(Ebx,2);

YUperpM = zeros(l1*l2,size(lamb,1)); % array that stores optimal YUperp,
Yopt = zeros(p*N,size(lamb,1));      % array that stores optimal Y
Y0 = zeros(p*(r+1),size(U,2));

ran = zeros(1,size(lamb,1)); % records the rank of solutions

%% DAGP for different values of lambdas
 
for in2 = 1:size(lamb,1)
  lambda = lamb(in2,1);
  fprintf(' mu = %g\n', 1/lambda)

  [y,Y0] = sysHankelgp_cv(U,bx,r,1/lambda,tol,freq,inf,Y0);
                             
  Y = Hrx(y,r)*U;
  YUperpM(:,in2) = vec(Y);    %vec(H_r(Y)Uperp) for all delta
  Yopt(:,in2) = vec(y);   %vec(Y) for all delta
end

save Data_raw YUperpM Yopt bx bxv uv r m p lamb

%% Record the ranks of solution matrices
for in2 = 1:size(lamb,1)
  
  d1 = reshape(YUperpM(:,in2),l1,l2);
  e1 = svd(d1);
  enew = e1(1)*sigmatol;
  n = 2;
  while e1(n) >= enew
    n = n + 1;
    if n == length(e1) + 1
      break
    end
  end
  ran(1,in2) = min(r,n - 1);
%% uncomment the next three lines for user-input ranks
%   semilogy(e1(1:20)/e1(1),'o')
%   [xcoor,~] = ginput(1);
%   ran(1,in2) = round(xcoor);
  
end

%% Compute identification and validation errors
nvec = ran';
Nid = N; Nv = N1;
ymeasi = bx; ymeasv = bxv;

save Datas3 YUperpM Yopt bx bxv uv nvec r Nid Nv m p lamb
Totalresult = Inputiderror_acc(YUperpM,Yopt,ymeasi,ymeasv,uv,nvec,r,Nid,Nv,m,p,lamb);
Totalresult(:,[1 2 4 5])

% validation_plot