function [W,A] = SVDstoch_icamp(X,K)

% [W,A] = SVDstoch_icamp(X,K)
% Code to perform stochastic matrix decomposition. 
% W and A are learned to approximate the non-zero entries in X.
% non-zero entries are picked randomly and gradient updates on A and W are
% performend.
% Input: X (matrix to be approximated), K (nr. of dimensions to project
% down to)
% Output: W and A such that X ~ W*A.
% In doing so, the zeros entries will be fileld in as well which produces
% predictions.

[Nattr,Ndata] = size(X);

showit = 1000;
T = 100000;
ErrR = zeros(1,T);

Mask = (X~=0);
Nv = sum(Mask(:));

alpha = 0;
beta = 0;
eta_W = 0.01;
eta_A = 0.01;
fac = 0.5;

W = 0.1*randn(Nattr,K);
A = 0.1*randn(K,Ndata);

ddW = zeros(Nattr,K);
ddA = zeros(K,Ndata);

t=0;
t_plot=0;
while t<T
    t=t+1;
    
    i = ceil(Nattr*rand);
    j = ceil(Ndata*rand);
    
    if X(i,j) ~= 0 % this should be made more efficient for large matrices. 
                   % We need to pick non-zero elements instead of 
                   % rejecting zero elements when we pick them.
        
        D = X(i,j) - W(i,:)*A(:,j);
        
        dW = D.*A(:,j)' - alpha*W(i,:);  %update for W
        ddW(i,:) = eta_W*dW + fac*ddW(i,:);
        W(i,:) = W(i,:) + ddW(i,:);
        
        dA = D.*W(i,:)' - beta*A(:,j);   %update for A
        ddA(:,j) = eta_A*dA + fac*ddA(:,j);
        A(:,j) = A(:,j) + ddA(:,j);
        
        
        %%%%%%%%%%%%%  plotting/printing  %%%%%%%%%%%%%%%%%%%
        
        
        
        if rem(t,showit) == 0
            t_plot=t_plot+1;
            
            fprintf('ITERATION: %.3d \n',t );
            
            R = W*A;
            ErrR(t_plot) = sqrt(sum(sum((Mask.*(X-R).^2))/(Nv)));
            
            figure(1);
            plot(ErrR(1:t_plot));
            title('Root Mean Square Error')
            xlabel('iteration')
            ylabel('RMSE');
            
            figure(2);
            colormap('hot')
            imagesc(X);
            title('Original Matrix M')
            
            figure(3);
            colormap('hot')
            imagesc(R)
            title('Approximated Matrix M_approx')
            
            figure(4);
            colormap('hot')
            imagesc(abs(Mask.*(X-R)));
            title('Difference: M-M_approx')
            
            figure(5);
            plot(W(:,1),W(:,2),'r.');
            title('2-d projection of movies genres')
            
            figure(6);
            plot(A(1,:),A(2,:),'b.');
            title('2-d projection of user attitudes')
            
            drawnow;
            
        end
        
    end
    
end%for t