fast_in_situ_calibration/MATLAB/calib_meth/REMNeNMF/remnenmf.m

function [ T , RMSE , MERs ] = remnenmf( X, X_theo, W, F_theo, Omega_G, Omega_F, Phi_G, Phi_F, G, F, config)

%% loading the config parameters 
Tmax = config.Tmax;
delta_measure = config.delta_measure;
InnerMinIter = config.InnerMinIter;
InnerMaxIter = config.InnerMaxIter;
M_loop = config.M_loop;
nu = config.nu;

%% 

r = nu;
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 = config.numSensor;
em_iter_max = round(Tmax / delta_measure) ;
T = nan(1,em_iter_max);
RMSE = nan(1+config.numSubSensor,em_iter_max);
MERs = nan(1+config.numSubSensor,em_iter_max);

nW = (1-W);
% X = G*F+W.*(X0-G*F);
X = X0 + nW.*(G*F);


GG = G' * G;
GX = G' * X ;
GradF = GG * F - GX;

FF = F * F';
XF = X * F' ;
GradG = nOmega_G.*(G * FF - XF);

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

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

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

        F(Omega_F) = 0; % Convert F to \Delta F
%         F = F - Phi_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
%         F = F + Phi_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
            [MER,~]=bss_eval_mix(F_theo',F');
            MERs(:,i) = MER;
            T(i) = toc;
            RMSE(:,i) = vecnorm(F(:,1:num_sensor) - F_theo(:,1:num_sensor),2,2)/sqrt(num_sensor);
        end
        T_M = cat(1,T_M,toc - t_m);
        
    end
    
end
niter
disp(['rem E step : ',num2str(mean(T_E))])
disp(['rem M step : ',num2str(mean(T_M))])
end

function [ L,R ] = RSI_compression(X,r,varargin)
%             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=2;
[m,n]=size(X);

l = min(min(n,m), max(compressionLevel, r ));

switch nargin
    case 2
        OmegaL = randn(n,l);
        OmegaR = randn(l, m);
        q = 4;
    case 4
        OmegaL = varargin{2};
        OmegaR = varargin{1};
        q = 1;
end

Y = X * OmegaL;

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


Y = OmegaR * X;

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


end