%Comparsion of rSPF and LRAS (IRLS as Benchmark)


%% Parameters

dimX1 = 16;
dimX2 = 100;
numX = 30;

alpha = 0;
beta = 0.5;
rank = 1;
sparsity = 0.1;
N0 = 50;

measurements = rank*(40:5:115);

normX = 10;
tolerance = 0.1;
errorNrm = 0.0;
rec_prob = zeros(3,size(measurements,2));
meanErr = zeros(3,size(measurements,2));
t = [0,0,0];

%% Test
    
    for p = 1:size(measurements,2)
                
        numMeas = measurements(p);
        A = randn(numMeas,dimX1*dimX2)/(numMeas^0.5);
    
        for k = 1:numX     

            X = createSparseRandomLowRankInputCommonSupport(dimX1,dimX2,rank,sparsity);
            nrm = norm(X,'fro');
            X = X*normX/nrm;
            s = nnz(round(X(1,:),4));
            
            [p,k]
            
            %% LRAS
            
            y = getMeasurement(X,A);
            eps = rand(length(y),1);
            eps = eps*errorNrm*normX/norm(eps);
            y = y+eps;
            z = norm(y);
            [u0,v0] = StartValue(A*X(:),A,dimX1,dimX2,rank);

            tic
            [U,S,V] = LRAS(y,A,alpha,beta,u0,v0,N0,rank);
            t(1) = t(1)+toc;
            EX = norm(X-U*eye(dimX1,dimX2)*V','fro');

            if EX < tolerance*normX
                rec_prob(1,p) = rec_prob(1,p)+1;
                meanErr(1,p) = meanErr(1,p)+(EX/normX);
            end
            
            %% rSPF                
            
            tic
            y = A*X(:) + eps;
            [U,V] = rSPF(A,y,dimX1,dimX2,rank,1,sparsity,u0,v0,N0);
            t(2) = t(2)+toc;
            EX = norm(X-U*V','fro');
            
            if EX < tolerance*normX
                rec_prob(2,p) = rec_prob(2,p)+1;
                meanErr(2,p) = meanErr(2,p)+(EX/normX);
            end
            
            %% IRLS
            
            if numMeas/rank > 100
               
                tic
                Xe = HM_IRLS(dimX1,dimX2,numMeas,rank,N0,0.01,A,y);
                t(3) = t(3)+toc;
                EX = norm(X-Xe(:,:,size(Xe,3)),'fro');

                if EX < tolerance*normX
                    rec_prob(3,p) = rec_prob(3,p)+1;
                    meanErr(3,p) = meanErr(3,p)+(EX/normX);
                end 
                
            end
            
        end
        
    end
    
    
%% Check Results and Plot

meanErr = meanErr./rec_prob;
rec_prob = rec_prob/numX;
t = t/(numX*size(measurements,2));

plot(measurements,rec_prob(1,:),measurements,rec_prob(2,:),measurements,rec_prob(3,:));
axis([measurements(1) measurements(length(measurements)) 0 1]);
xlabel('#Measurements');
ylabel('Recovery Probability');
title('Recovery of randomly chosen 16x100 matrices with 10-sparse right-SVs and 15% error tolerance');
legend('LRAS','SPF','IRLS');
    
