fast_in_situ_calibration/MATLAB/calib_meth/REMNeNMF/remnenmf.m

function [ T , RMSE ] = remnenmf( W , X , G , F , Omega_G, Omega_F, Phi_G, Phi_F , InnerMinIter , InnerMaxIter , Tmax , v, F_theo, delta_measure)

r = 0;
X0 = X;
Omega_G = (Omega_G == 1); % Logical mask is faster than indexing in matlab.
Omega_F = (Omega_F == 1); % Logical mask is faster than indexing in matlab.
nOmega_G = ~Omega_G; % Logical mask is faster than indexing in matlab.
nOmega_F = ~Omega_F; % Logical mask is faster than indexing in matlab.

[~, num_sensor] = size(F);
num_sensor = num_sensor-1;
em_iter_max = round(Tmax / delta_measure) ;
T = nan(1,em_iter_max);
RMSE = nan(2,em_iter_max);

X = G*F+W.*(X0-G*F);
[L,R]=RSI_compression(X,r);
% Compress left and right
X_L = L * X;
X_R = X * R;
G_L = L * G;
F_R = F * R;

GG = G_L' * G_L;
GX = G_L' * X_L;
GradF = GG * F - GX;

FF = F_R * F_R';
XF = X_R * F_R';
GradG = G * FF - XF;

d = Grad_P([GradG',GradF],[G',F]);
StoppingCritF = 1.e-3*d;
StoppingCritG = StoppingCritF;

tic
i = 1;
niter = 0;
RMSE(:,i) = vecnorm(F(:,1:end-1)- F_theo(:,1:end-1),2,2)/sqrt(num_sensor);
T(i) = toc;
while toc<Tmax
    
    % Estimation step
    X = G*F+W.*(X0-G*F);
    [L,R]=RSI_compression(X,r);
    % Compress left and right
    X_L = L * X;
    X_R = X * R;
    
    % Maximization step
    for j =1:v
        F_R = F * R;
        FF = F_R * F_R';
        XF = X_R * F_R' - Phi_G * FF;

        G(Omega_G) = 0; % Convert G to \Delta G
        [ G , iterG ] = MaJ_G_EM_NeNMF( FF , XF , G , InnerMinIter , InnerMaxIter , StoppingCritG , nOmega_G); % Update \Delta G
        G(Omega_G) = Phi_G(Omega_G); % Convert \Delta G to G
        niter = niter + iterG;
        if(iterG<=InnerMinIter)
            StoppingCritG = 1.e-1*StoppingCritG;
        end
        
        G_L = L * G;
        GG = G_L' * G_L;
        GX = G_L' * X_L - GG * Phi_F;

        F(Omega_F) = 0; % Convert F to \Delta F
        [ F , iterF ] = MaJ_F_EM_NeNMF( GG , GX , F , InnerMinIter , InnerMaxIter , StoppingCritF , nOmega_F); % Update \Delta F
        F(Omega_F) = Phi_F(Omega_F); % Convert \Delta F to F
        niter = niter + iterF;
        if(iterF<=InnerMinIter)
            StoppingCritF = 1.e-1*StoppingCritF;
        end

        if toc - i*delta_measure >= delta_measure
            i = i+1;
            if i > em_iter_max
                break
            end
            T(i) = toc;
            RMSE(:,i) = vecnorm(F(:,1:end-1) - F_theo(:,1:end-1),2,2)/sqrt(num_sensor);
        end
    end
    
end
niter
end

function [ L,R ] = RSI_compression(X,r)
%             Tepper, M., & Sapiro, G. (2016). Compressed nonnegative
%             matrix factorization is fast and accurate. IEEE Transactions
%             on Signal Processing, 64(9), 2269-2283.
%             see: https://arxiv.org/pdf/1505.04650
%             The corresponding code is originally created by the authors
%             Then, it is modified by F. Yahaya.
%             Date: 13/04/2018
%

compressionLevel=20;
[m,n]=size(X);

l = min(n, max(compressionLevel, r + 10));

OmegaL = randn(n,l);

Y = X * OmegaL;

for i=1:4
    
    [Y,~]=qr(Y,0);
    S=X'*Y;
    [Z,~]=qr(S,0);
    Y=X* Z;
end
[L,~]=qr(Y,0);
L=L';

OmegaR = randn(l, m);
Y = OmegaR * X;

for i=1:4
    [Y,~]=qr(Y',0);
    S=X*Y;
    [Z,~]=qr(S,0);
    
    Y=Z'*X;
end
Y=Y';
[R,~] = qr(Y,0);


end