% run code for system realization

m = 10; 
r = 10;
i = 2*r+1;
n = m; p = m;
T = 999;
freq = 10;
maxiter = inf;

repeats = 10;
fid = fopen('sys.txt','w');
fclose(fid);

fid = fopen('sys2.txt','w');
fclose(fid);

for j = [100 200 300]
  
  q = zeros(m,n*(i+j-1));
  q(:,1:n*j) = ones(m,n*j);
  
  for mu = [0.05 0.1 0.5]
    
    PPG = [];
    ADMM = [];
    PROX = [];
    
    for repeat = 1:repeats
      
      A = randn(r,r); A = A/norm(A); % normalize A
      B = randn(r,p); B = B/norm(B); % normalize B
      C = randn(m,r); C = C/norm(C);
      
      x = randn(r,1);
      ynew = zeros(n,T+1); 
      for k = 1:T+1  
        e = randn(n,1);
        xnew = A*x + B*e;
        ynew(:,k) = C*x + e;
        x = xnew;  
      end
      % Add noise
      by = ynew + 5e-2*randn(n,T+1);
      
      bx = zeros(m,n*(i+j-1));
      for k = 1:j
        bx(:,(k-1)*n+1:k*n) = zeros(n,n);
        for t = 1:T+2-k
          bx(:,(k-1)*n+1:k*n) = bx(:,(k-1)*n+1:k*n) + by(:,t+k-1)*by(:,t)';
        end
      end
      bx = bx/(T+1);
      
      if mu == 0.05
        [zbest,fbest,dval,time,iter,dfeas] = rcompletePPG2(q,bx,m,n,i,j,mu,1e-4,freq,maxiter);
      else
        [zbest,fbest,dval,time,iter,dfeas] = rcompletePPG(q,bx,m,n,i,j,mu,1e-4,freq,maxiter);
      end
      PPG = [PPG;iter time fbest -dval dfeas];
      
      [x,fval,dval,t,iter,dfeas] = rcompleteADM2(q,bx,m,n,i,j,mu,1e-4,1,freq,maxiter);
      ADMM = [ADMM;iter t fval dval dfeas];
      
      [zbest,fbest,dval,time,iter,dfeas] = rcompletemFBS_uv(q,bx,m,n,i,j,mu,1e-4,freq,maxiter);
      PROX = [PROX;iter time fbest -dval dfeas];
      
    end
    PPG = mean(PPG,1);
    ADMM = mean(ADMM,1);
    PROX = mean(PROX,1);
    
    fid = fopen('sys.txt','a');
    fprintf(fid,'%3d & %3.2f & %0.0f & %0.1f & %0.3e/%0.3e/%0.1e & %0.0f & %0.1f & %0.3e/%0.3e/%0.1e\n', j, mu, PPG, ADMM);
    fclose(fid);
    
    fid = fopen('sys2.txt','a');
    fprintf(fid,'%3d & %3.2f & %0.0f & %0.1f & %0.3e/%0.3e/%0.1e\n', j, mu, PROX);
    fclose(fid);
    
    
  end
end