clear
M = 20; % Signal length
N = 40; % Coefficient Space Dimension
L = 32*N; % Number of Training Samples
R = 3; % Sparsity
IT = 1000; % Number of alternating sparse approximation and dictionary update
map = 1; % Debiasing. 0 = No, 1 = Yes
maxIT = 1000; % Inner-loop maximum iteration number.
lambda = 2*.2; % Lagrangian multiplier.
epsx = 10^-7; % Stopping criterion for iterative softthresholding
epsd = 10^-7; % Stopping criterion for MM dictionary update
cvset = 0; % Dictionary constraint. 0 = Non convex ||d|| = 1, 1 = Convex ||d||<=1
Tre = .99; % Threshold for accepting too atoms identical
%%%% Generative Dictionaries
Do = randn(M,N); % Generative Dictionary
Do = Do*(diag(sum((Do'*Do).*eye(length(Do))).^-.5)); % Normalization
%%%% Sparse signal generation %%%%% 
Xo = zeros(N,L); % Original Sparse Coefficients
for l = 1:L
    r = 1;
    while r<=R
        ind = fix(rand(1)*N)+ones(1);
        a = rand(1);
        if Xo(ind)==0
            Xo(ind,l) = (.8*rand(1)+.2)*((a>=.5)-(a<.5));
            r = r+1;
        end
    end
end
Y = Do*Xo; % Sparse Signals
%%%% Algorithm initialization
D = randn(M,N); % Initial Dictionary
D = D*(diag(sum((D'*D).*eye(length(D))).^-.5)); % Normalization
X = ones(size(Xo)); % Initial coefficients
for it = 1:IT,
    it
    to = .1+svds(D,1); 
    [X,cost(it)] = mm1(D,Y,X,to,lambda,maxIT,epsx,map);   
    plot(cost);
    [D,X] = dict_update_REG_cn(D,Y,X,maxIT,epsd,cvset);    
end
%%%
success = sum(max(abs((Do'*D)))>=Tre);
display(['  ------------------'])
display(['  ',num2str(success),'% of the atoms successfully recovered after ',num2str(IT),' iterations.']);