% plot recovered weights (figure 4)
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_Ahat = mean(voigtp.Am(:,:,end-Mest:end),3);
N_Bhat = mean(nmf.Bm(:,:,end-Mest:end),3);
N_Ahat = mean(nmf.Am(:,:,end-Mest:end),3);
% normalize
factors = max(N_Ahat);
N_Bhat=N_Bhat.*repmat(factors',1,N);

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