clear; close all;

% K-SVD dictionary update Demo.
% 
% main_KSVD generates random overcomplete dictionary, random sparse
% coeficients to produce the training data. Then a trained dictionary is
% learned by K-SVD (includes sparse coding stage and
% dictionary update stage)
% 
% 
% References:
% "The K-SVD: An Algorithm for Designing of Overcomplete Dictionaries for 
% Sparse Representation", written by M. Aharon, M. Elad, and A. M. Bruckstein
% and appeared in the IEEE Trans. On Signal Processing, Vol. 54, no. 11, 
% pp. 4311-4322, November 2006
%
%
% See also main_regSimCO
%
% Wei Dai, Tao Xu, Wenwu Wang
% Imperial College London, University of Surrey
% wei.dai1@imperial.ac.uk  t.xu@surrey.ac.uk  w.wang@surrey.ac.uk
% October 2011

%% initialization
m = 16; %the dimension of the training signal
d = 32; %the number of atoms in training dictionary
S = 4; %the number of nonzero elements in coefficient vector
Nv = [64:32:512]; %the set of different number of training signals
NvLen = length(Nv); %the number of tested value for Nv
n_o_realization = 100; %the number of random experiments
itN = 50; %the number of iterations of KSVD algorithm
IPara.itN = 1; %the number of iterations in dictionary update stage
IPara.gmin = 1e-5; %the minimum value of gradient
IPara.dispN = 20; %the display interval
IPara.DebugFlag = 0; %the flag for debugging ill-conditioned problem
IPara.mu = 0; %K-SVD without regularized item, so IPara.mu is set to 0
f_save = zeros(NvLen,n_o_realization); %save the cost function value of each realization
d_save = zeros(NvLen,n_o_realization); %save the average distance between original dictionary and trained dictionary
time_save = zeros(1,NvLen); %save the computational cost of different number of training signals

%test different number Nv(i) of training signals
for nn = 1:NvLen
    randn('state',0);
    rand('state',0);
    n = Nv(nn); 
    
    time_start = clock; 
    for rn = 1:n_o_realization
        %generate a dictionary and normalize it
        D_true = randn(m,d);
        DCol2 = sqrt(sum(D_true.*D_true,1));
        D_true = D_true .* repmat(1./DCol2,m,1);

        %generate a sparse coefficients matrix X_true
        X_true = randn(d,n);
        for cn = 1:n
            NZIndex = randsample(d,S);
            X_true( setdiff(1:d,NZIndex),cn ) = 0;
        end
        
        %generate Y implicates the sparsity pattern
        Y = D_true*X_true;
        Yn2 = sqrt(sum(Y.*Y,1));
        Y = Y.*repmat(1./Yn2,m,1);
        
        %generate initial dictionary for the test
        D=randn(m,d);
        DCol2 = sqrt(sum(D.*D,1));
        D = D .* repmat(1./DCol2,m,1); 
        
        %dictionary learning step
        fprintf('Dictionary Learning: KSVD: n=%d realization=%d\n',n,rn);
        
        for itn = 1:itN
            X = OMP(D,Y,4); %sparse coding stage
            
            [D,X,OPara] = K_DictUpdate02(Y,D,X,IPara);
        end
        
        Yr = Y - D*X;
        ftt = sum(sum(Yr.*Yr)); 
        fprintf('** f=%3.3e n=%d realization=%d\n\n',ftt,n,rn);
        
        %save data
        f_save(nn,rn) = ftt;
        d_save(nn,rn) = DictDist01(D,D_true);
    end
    
    time_end = clock;
    time_save(nn) = etime(time_end,time_start);
    
    save maindata_KSVD Nv n_o_realization m d S f_save dc_save time_save;
end

