fast_in_situ_calibration/MATLAB/calib_meth/EMNeNMF/emnenmf.m

function [ T , RMSE , MERs ] = emnenmf( 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;

%% 

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);

X = G*F+W.*(X0-G*F);
GG = G'*G;
GX = G'*X;
GradF = GG*F-GX;

FF = F*F';
XF = X*F';
GradG = G*FF-XF;

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

tic
i = 1;
T(i) = toc;
RMSE(:,i) = vecnorm(F(:,1:num_sensor)- F_theo(:,1:num_sensor),2,2)/sqrt(num_sensor);
niter = 0;
T_E = [];
T_M = [];
while toc<Tmax
    
    t_e = toc;
    X = G*F+W.*(X0-G*F);
    T_E = cat(1,T_E,toc - t_e);
    
    for j =1:M_loop
        
        t_m = toc;
        FF = F*F';
        XF = X*F' - 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

        GG = G'*G;
        GX = G'*X-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
            [MER,~]=bss_eval_mix(F_theo',F');
            MERs(:,i) = MER;
            T(i) = toc;
            RMSE(:,i) = vecnorm(F(:,1:end-1) - F_theo(:,1:end-1),2,2)/sqrt(num_sensor);
        end
        T_M = cat(1,T_M,toc - t_m);
        
    end
  
end
niter  
disp(['em E step : ',num2str(mean(T_E))])
disp(['em M step : ',num2str(mean(T_M))])
disp(['RMSE 2ème colonne de G : ',num2str(norm(G(:,2) - X_theo(:,end),2)/sqrt(config.sceneWidth*config.sceneLength))])
end