% plot recovered spectrum (figure 2)
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_Vhat = mean(voigtp.Vm(:,:,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);

figure
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)
opt.fig.setProps = false;
prepfig([strFIGdir 'simulated_recovered_spectrum'],opt.fig)