% RMUSIC-DANSER Toy Example

clear;
clc;
close all

%% parameter settings
load Lib_demo 
normA  = sqrt(sum(Lib.^2));
Dict_E = min(normA);
normB = sqrt(sum(Lib.^2));
alpha = .8;
delta = Dict_E *(1-alpha)/(1+alpha);
 
Lib_o = Lib; % save the orginal library without modeling error
K_o = size(Lib,2);
N = 6;       % number of endmembers

[M, n_materiais]=size(Lib);
sel_mat = randperm(n_materiais);%
rand_index = sel_mat;
sel_mat = sel_mat(1:N);  % Randomly pick up N signatures
Bands = 224;
% -------- random varation --------------
Variab = randn(Bands,n_materiais); %
max_norm = max(sqrt(sum(Variab.^2)));
for i=1:n_materiais
    Variab(:,i)=Variab(:,i)/max_norm;
end

%-----------------------------------------------
A_o = Lib_o(1:Bands,sel_mat);        %A_o is the true endmembers 
clear A;
A = A_o;                             %Ground truth mixing matrix without modeling error



%% create HSI

SIGNATURES_TYPE = 1;        % Uniform in [0,1]
p = N;                      % number of endmembers
L = 5000;                   % number of pixels
SNR = 35;                   % signal-to-noise ratio (E ||y||^2/E ||n||^2) in dBs
SHAPE_PARAMETER  = 1   ;    % uniform over the simplex
MAX_PURIRY = .7;           % there are pure pixels in the data set
OUTLIERS  = 0;              % Number of outliers in the data set

[X,s,noise] = spectMixGen(A,L,'Source_pdf', 'Diri_id','pdf_pars',SHAPE_PARAMETER,...
    'max_purity',MAX_PURIRY*ones(1,p),'no_outliers',OUTLIERS, ...
    'violation_extremes',[1,1.2],'snr', SNR, ...
    'noise_shape','gaussian','noise_pars', [2,50,20] ,'pure_pixels','no');

%% projection matrix for MUSIC and RMUSIC
[w,Rw]=estNoise(X,'off');
[kf,Ek,E,delta_p]=hysime(X,w,Rw,'off');
E = E(:,1:kf);
Pyort = (eye(M)-E*E');


%% create library with modeling error
DMER = 20;
sigma = Dict_E / 10^(DMER/20);
Lib= Lib_o +  sigma*Variab;            % Suppose that our library contains some errors


%% ROBUST MUSIC 
Mp = 40; % specify the number of library members we want to keep 
th = delta;

[ind_rmusic,dddd,err_rmusic] = RMUSIC( Lib,th,Pyort,Mp,sel_mat );
Dict_rmusic  = Lib(:,ind_rmusic);

figure(1)
subplot(2,1,1)
semilogy([1:n_materiais],err_rmusic,'b.');
hold on
semilogy(sel_mat,err_rmusic(sel_mat),'ro')
title('RMUSIC','fontsize',14)
xlabel('index of materials (k)','fontsize',14)
ylabel('residual','fontsize',14)
set(gca,'fontsize',14)
xlim([1 n_materiais])

%% Ordinary MUSIC
[ ind_music, detect_MUSIC, err_music] = MUSIC_endm_detect( Lib,Pyort,Mp,sel_mat );
Dict_music = Lib(:, ind_music );

 
figure(1) 
subplot(2,1,2)
semilogy([1:n_materiais],err_music,'b.');
hold on
semilogy(sel_mat,err_music(sel_mat),'ro')
title('MUSIC','fontsize',14)
xlabel('index of materials (k)','fontsize',14)
ylabel('residual','fontsize',14)
set(gca,'fontsize',14)
xlim([1 n_materiais]) 

%% Unmixing Part 

%----------------------------------------
%                     RMUSIC-CSR
%----------------------------------------

[X_hat_l21] = clsunsal(Dict_rmusic,X,'POSITIVITY','yes','VERBOSE','no','ADDONE','yes', ...
    'lambda',1,'AL_ITERS',2000, 'TOL', 1e-8);
S = zeros(K_o,L);
S(ind_rmusic(1:Mp),:) = X_hat_l21;
S1=S;
S_norm = sum(S1,2);

S_real = zeros(1,n_materiais);
S_real(sel_mat)= max(S_norm)+100;

h=figure(2);
subplot(2,1,1)
plot([1:n_materiais],S_norm);hold on
stem([1:n_materiais],S_real,'--k','linewidth',0.5,'marker','none')
title('RMUSIC-CSR','fontsize',14)
xlabel('index of materials (k)','fontsize',14)
ylabel('ck','fontsize',14)
set(gca,'fontsize',14)
ylim([0 max(S_norm)+100])
xlim([1 n_materiais])

%----------------------------------------
%                     RMUSIC-DANSER
%----------------------------------------
epsilon=delta;
p= .75; lambda = .5; gamma= 10000;
maxiter =10000;
[S_rmusic] = DNANSER_mod(Dict_rmusic, Dict_rmusic, X_hat_l21, p, X, lambda, gamma, delta, maxiter );

S = zeros(K_o,L);
S(ind_rmusic(1:Mp),:) = S_rmusic;

S_norm = sum(S,2);
S_real(sel_mat)= max(S_norm)+100;
figure(2)
subplot(2,1,2)
plot([1:n_materiais],S_norm);hold on
stem([1:n_materiais],S_real,'--k','linewidth',0.5,'marker','none')
title('RMUSIC-DANSER','fontsize',14)
xlabel('index of materials (k)','fontsize',14)
ylabel('ck','fontsize',14)
set(gca,'fontsize',14)
ylim([0 max(S_norm)+100])
xlim([1 n_materiais])



     
       



