plot_recovered_spectra.m #1

clearvars;
if ~exist('first','var')
    run('custom_settings')
    load([getenv('datadir') 'simulated'],'gendata');
    voigtp = load([getenv('datadir') 'simulated_voigt']);
    nmf = load([getenv('datadir') 'simulated_nmf']);
    first = 1;
end
opt.fig.size = 90*[opt.fig.col2width 2.5];

N = length(gendata.A);
W = length(gendata.B);
Mest = 1000;

V_Bhat = mean(voigtp.Bm(end-Mest:end,:));
V_bhat = mean(voigtp.bm(end-Mest:end,:));
V_BBhat = V_bhat'*V_Bhat;
V_Ahat = mean(voigtp.Am(:,:,end-Mest:end),3);
V_Vhat = mean(voigtp.Vm(:,:,end-Mest:end),3);
V_tauhat = mean(voigtp.taum(end-Mest:end));

N_Bhat = mean(nmf.Bm(:,:,end-Mest:end),3);
N_Ahat = mean(nmf.Am(:,:,end-Mest:end),3);
% normalize
factors = max(N_Ahat);
N_Ahat=N_Ahat./repmat(factors,W,1);
N_Bhat=N_Bhat.*repmat(factors',1,N);



%% plot recovered weights
figure(1)
clf
[vals inx] = sort(gendata.A,'descend');
plot(1:N,V_Ahat(inx))
xlim([1 N])
hold on
plot(1:N,N_Bhat(1,inx)')
ylabel(['Intensity'])
xlabel(['Measurement index'])
prepfig('',opt.fig)
plot(1:N,vals,'--','linewidth',2)
print( gcf,'-dpdf',[strFIGdir 'simulated_loadings.pdf']);

%% plot recovered background
figure(2)
clf
plot(1:W,V_Bhat)
xlim([1 W])
hold on
plot(1:W,N_Ahat(:,2))
ylabel(['Intensity'])
xlabel(['Wavenumber'])
prepfig('',opt.fig)
plot(1:W,gendata.B,'--','linewidth',2)
print( gcf,'-dpdf',[strFIGdir 'simulated_background.pdf']);
%% plot recovered spectrum
figure(3)
clf
plot(voigt(V_Vhat,1:W))

hold on
plot(1:W,N_Ahat(:,1))
ylabel(['Intensity'])
xlabel(['Wavenumber'])
prepfig('',opt.fig)
plot(voigt(gendata.theta,1:W),'--','linewidth',2)
print( gcf,'-dpdf',[strFIGdir 'simulated_recovered_spectrum.pdf']);