%% NMF by Vector-wise BCD, also called HALS
% Andersen Ang @ MARO, FPM, UMONS, Belgium
% Last update 2020-May-10
%%
clear, close all, clc
% load an data, for example the rice image
M = double(imread('rice.png'));
% get the size of X
[m,n] = size(M);
% set the rank r
r = 84; % just randomly chossen
%% Initialization
  % Initialize W and H
  Wini = rand(m,r);
  Hini = rand(r,n);
  % Set the epsilon 
  epsilon = 0;
  epsW = epsilon*ones(m,1);
  epsH = epsilon*ones(1,n);
  % inital error value
  err0 = norm(M-Wini*Hini,'fro')^2/norm(M,'fro')^2*100;  
  % max number of iteration
  maxiter = r*30;
  % pre-allocate space for error value
  errval = NaN*zeros(2,maxiter);
%% The main loop
W = Wini;
H = Hini;
for k = 1 : maxiter % or other stopping conditions
 % On varaible W
   % pick update index
   j = mod(k-1,r)+1;
   % pick column of W
   wj = W(:,j);
   % pick row of H
   hj = H(j,:);
   % compute Xi
   Mj = M - W*H + wj*hj;
   % update wi
   wj = max(Mj*hj'/(hj*hj'), epsW);
   % Put it back
   W(:,j) = wj;
   % Compute error value
   errval(1,k) = norm(M-W*H,'fro')^2/norm(M,'fro')^2*100;

 % On variable H
   % pick update index
   % i = mod(k-1,r)+1; % same so no need to do
   % pick column of W
   % wj = W(:,i);  % same so no need to do
   % pick row of H
   % hj = H(i,:); % same so no need to do
   % compute Xi
   Mj = M - W*H + wj*hj;
   % update hi
   hj = max((Mj'*wj)'/(wj'*wj), epsH);
   % Put it back
   H(j,:) = hj;
   % Compute error value
   errval(2,k) = norm(M-W*H,'fro')^2/norm(M,'fro')^2*100;
end
%% Plot the error curve
figure('units','normalized','outerposition',[0 0 1 1])
subplot(2,2,[1:2])
both_errval = errval(:);
semilogy(both_errval,'k','linewidth',1),hold on
semilogy(1:2:numel(both_errval), both_errval(1:2:end),'b.','markersize',13)
semilogy(2:2:numel(both_errval), both_errval(2:2:end),'r.','markersize',13)
grid on, axis tight
title('Relative Percentage Error verus iterations','interpreter','latex','fontsize',20)
xlabel('Iteration','interpreter','latex','fontsize',15)
ylabel('RPE = $\frac{\| X - WH \|_F^2}{\| X \|_F^2} \times 100\% $','interpreter','latex','fontsize',20)
lgd = legend('Error curve','After w is updated','After h is updated');
lgd.FontSize = 20;

subplot(223),imagesc(M),colormap gray
title('Original image','interpreter','latex','fontsize',20);
subplot(224),imagesc(W*H),colormap gray
title('Reconstrcted using HALS','interpreter','latex','fontsize',20);