Push B4 het

MATLAB/calib_meth/EMNeNMF/emnenmf.m

@@ -30,10 +30,17 @@ i = 1;
 T(i) = toc;
 RMSE(:,i) = vecnorm(F(:,1:end-1)- F_theo(:,1:end-1),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:v
+        
+        t_m = toc;
         FF = F*F';
         XF = X*F' - Phi_G*FF;
 
@@ -64,10 +71,14 @@ while toc<Tmax
             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(median(T_E))])
+disp(['em M step : ',num2str(median(T_M))])
 end
 
 

MATLAB/calib_meth/REMNeNMF/remnenmf.m

@@ -1,6 +1,6 @@
-function [ T , RMSE ] = remnenmf( W , X , G , F , Omega_G, Omega_F, Phi_G, Phi_F , InnerMinIter , InnerMaxIter , Tmax , v, F_theo, delta_measure)
+function [ T , RMSE ] = remnenmf( W , X , G , F , Omega_G, Omega_F, Phi_G, Phi_F , InnerMinIter , InnerMaxIter , Tmax , v, F_theo, delta_measure, nu)
 
-r = 0;
+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.
@@ -13,25 +13,24 @@ 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;
+nW = (1-W);
+% X = G*F+W.*(X0-G*F);
+X = X0 + nW.*(G*F);
 
-GG = G_L' * G_L;
-GX = G_L' * X_L;
+
+GG = G' * G;
+GX = G' * X ;
 GradF = GG * F - GX;
 
-FF = F_R * F_R';
-XF = X_R * F_R';
-GradG = G * FF - XF;
+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 = StoppingCritF;
+StoppingCritG = 1.e-3*d;
+T_E = [];
+T_M = [];
 
 tic
 i = 1;
@@ -40,34 +39,48 @@ RMSE(:,i) = vecnorm(F(:,1:end-1)- F_theo(:,1:end-1),2,2)/sqrt(num_sensor);
 T(i) = toc;
 while toc<Tmax
     
+    t_e = toc;
     % Estimation step
-    X = G*F+W.*(X0-G*F);
+    X = X0 + nW.*(G*F);
+%     if i>1
+%         [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: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
+        
+        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;
-        GX = G_L' * X_L - GG * Phi_F;
+%         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;
@@ -81,13 +94,17 @@ while toc<Tmax
             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(['rem E step : ',num2str(median(T_E))])
+disp(['rem M step : ',num2str(median(T_M))])
 end
 
-function [ L,R ] = RSI_compression(X,r)
+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.
@@ -97,16 +114,25 @@ function [ L,R ] = RSI_compression(X,r)
 %             Date: 13/04/2018
 %
 
-compressionLevel=20;
+compressionLevel=2;
 [m,n]=size(X);
 
-l = min(n, max(compressionLevel, r + 10));
-
-OmegaL = randn(n,l);
+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:4
+for i=1:q
     
     [Y,~]=qr(Y,0);
     S=X'*Y;
@@ -116,10 +142,10 @@ end
 [L,~]=qr(Y,0);
 L=L';
 
-OmegaR = randn(l, m);
+
 Y = OmegaR * X;
 
-for i=1:4
+for i=1:q
     [Y,~]=qr(Y',0);
     S=X*Y;
     [Z,~]=qr(S,0);

MATLAB/calib_meth/common_nesterov/MaJ_F_EM_NeNMF.m

@@ -21,11 +21,11 @@ for i = 1 : InnerMaxIter
     F1 = F2;
     alpha(1) = alpha(2);
     
-    if mod(i , InnerMinIter) == 0
-        if(Grad_P(Grad_y , Y)<=StoppingCritF)
-            break
-        end
-    end
+%     if mod(i , InnerMinIter) == 0
+%         if(Grad_P(Grad_y , Y)<=StoppingCritF)
+%             break
+%         end
+%     end
     
 end
 end

MATLAB/calib_meth/common_nesterov/MaJ_G_EM_NeNMF.m

@@ -21,11 +21,11 @@ for i = 1 : InnerMaxIter
     G1 = G2;
     alpha(1) = alpha(2);
     
-    if mod(i , InnerMinIter) == 0
-        if(Grad_P(Grad_y , Y)<=StoppingCritG)
-            break
-        end
-    end
+%     if mod(i , InnerMinIter) == 0
+%         if(Grad_P(Grad_y , Y)<=StoppingCritG)
+%             break
+%         end
+%     end
     
 end
 end

MATLAB/calib_meth/test/remnenmf_full.m

@@ -0,0 +1,160 @@
+function [ T , RMSE ] = remnenmf_full( W , X , G , F , Omega_G, Omega_F, Phi_G, Phi_F , InnerMinIter , InnerMaxIter , Tmax , v, F_theo, delta_measure, 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] = 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);
+
+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:end-1)- F_theo(:,1:end-1),2,2)/sqrt(num_sensor);
+T(i) = toc;
+while toc<Tmax
+    
+    t_e = toc;
+    % Estimation step
+    X = X0 + nW.*(G*F);
+%     if i>1
+%         [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:v
+        
+        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
+            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(['rem E step : ',num2str(median(T_E))])
+disp(['rem M step : ',num2str(median(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
+

MATLAB/calib_meth/test/remnenmf_not_full.m

@@ -0,0 +1,160 @@
+function [ T , RMSE ] = remnenmf_not_full( W , X , G , F , Omega_G, Omega_F, Phi_G, Phi_F , InnerMinIter , InnerMaxIter , Tmax , v, F_theo, delta_measure, 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] = 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);
+
+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:end-1)- F_theo(:,1:end-1),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:v
+        
+        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
+            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(['rem E step : ',num2str(median(T_E))])
+disp(['rem M step : ',num2str(median(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
+

MATLAB/main_emnenmf_vs_remnenmf.m

@@ -1,6 +1,5 @@
-clear
-close all
-clc
+clear all
+
 
 %% Adding every files in the path
 
@@ -11,7 +10,7 @@ addpath(genpath(pwd))
 s_width = 20;  % Scene width
 s_length = 20; % Scene length
 N_Ref = 4; % Nb. of reference measurements
-N_Cpt = 25; % Nb. of mobile sensors
+N_Cpt = 100; % Nb. of mobile sensors
 Mu_beta = .9; % Mean sensors gain
 Mu_alpha = 5; % Mean sensors offset
 Bound_beta = [.01;1.5]; % Gain boundaries
@@ -25,9 +24,12 @@ runs = 1; % Total number of runs
 
 %% Nesterov parameters
 
-InnerMaxIter = 300;
-InnerMinIter = floor(InnerMaxIter/10);
-Tmax = 30;
+InnerMinIter = 5;
+InnerMaxIter = 20;
+Tmax = 60;
+
+%% Compression parameters
+nu = 10;
 
 %% 
 delta_measure = 1;
@@ -45,11 +47,14 @@ RMSE_gain_emnenmf = nan(runs, iter_max);
 for run = 1:runs
     % data generation
     [X, X_theo, W, F_theo, Omega_G, Omega_F, Phi_G, Phi_F, Ginit, Finit] = data_gen(s_width, s_length, run, N_Ref, N_Cpt, Mu_beta, Mu_alpha, Bound_beta, Bound_alpha, MV, RV, var_n);
+    seed = floor(rand*10000);
     % remnenmf
-    [T_remnenmf, RMSE] = remnenmf( W , X , Ginit , Finit, Omega_G, Omega_F, Phi_G, Phi_F , InnerMinIter , InnerMaxIter , Tmax , M_loop, F_theo, delta_measure);
+    rng(seed)
+    [T_remnenmf, RMSE] = remnenmf_not_full( W , X , Ginit , Finit, Omega_G, Omega_F, Phi_G, Phi_F , InnerMinIter , InnerMaxIter , Tmax , M_loop, F_theo, delta_measure, nu);
     RMSE_offset_remnenmf(run,:) = RMSE(1,:);
     RMSE_gain_remnenmf(run,:) = RMSE(2,:);
     % emnenmf
+    rng(seed)
     [T_emnenmf, RMSE] = emnenmf( W , X , Ginit , Finit, Omega_G, Omega_F, Phi_G, Phi_F , InnerMinIter , InnerMaxIter , Tmax , M_loop, F_theo, delta_measure);
     RMSE_offset_emnenmf(run,:) = RMSE(1,:);
     RMSE_gain_emnenmf(run,:) = RMSE(2,:);    
@@ -65,6 +70,7 @@ min_gain_remnenmf = min(RMSE_gain_remnenmf,[],1,'omitnan');
 med_gain_remnenmf = median(RMSE_gain_remnenmf,1,'omitnan');
 max_gain_remnenmf = max(RMSE_gain_remnenmf,[],1,'omitnan');
 
+figure
 subplot(121)
 semilogy(T_remnenmf,min_offset_remnenmf,'b')
 hold on

MATLAB/main_full_vs_not_full_E_step.m

@@ -0,0 +1,119 @@
+clear all
+
+
+%% Adding every files in the path
+
+addpath(genpath(pwd))
+
+%% Simulation parameters
+
+s_width = 20;  % Scene width
+s_length = 20; % Scene length
+N_Ref = 4; % Nb. of reference measurements
+N_Cpt = 100; % Nb. of mobile sensors
+Mu_beta = .9; % Mean sensors gain
+Mu_alpha = 5; % Mean sensors offset
+Bound_beta = [.01;1.5]; % Gain boundaries
+Bound_alpha = [3.5;6.5]; % Offset boundaries
+
+MV = .9; % Missing Value prop.
+RV = 0.3; % RendezVous prop.
+var_n = 0; % Noise variance
+M_loop = 25; % Number of M loop per E step
+runs = 1; % Total number of runs
+
+%% Nesterov parameters
+
+InnerMinIter = 5;
+InnerMaxIter = 20;
+Tmax = 60;
+
+%% Compression parameters
+nu = 10;
+
+%% 
+delta_measure = 1;
+iter_max = round(Tmax / delta_measure);
+
+
+%% Allocation for the RMSE values
+% remnenmf
+RMSE_offset_remnenmf = nan(runs, iter_max);
+RMSE_gain_remnenmf = nan(runs, iter_max);
+% emnenmf
+RMSE_offset_emnenmf = nan(runs, iter_max);
+RMSE_gain_emnenmf = nan(runs, iter_max);
+
+for run = 1:runs
+    % data generation
+    [X, X_theo, W, F_theo, Omega_G, Omega_F, Phi_G, Phi_F, Ginit, Finit] = data_gen(s_width, s_length, run, N_Ref, N_Cpt, Mu_beta, Mu_alpha, Bound_beta, Bound_alpha, MV, RV, var_n);
+    seed = floor(rand*10000);
+    % remnenmf
+    rng(seed)
+    [T_remnenmf, RMSE] = remnenmf_full( W , X , Ginit , Finit, Omega_G, Omega_F, Phi_G, Phi_F , InnerMinIter , InnerMaxIter , Tmax , M_loop, F_theo, delta_measure, nu);
+    RMSE_offset_remnenmf(run,:) = RMSE(1,:);
+    RMSE_gain_remnenmf(run,:) = RMSE(2,:);
+    % emnenmf
+    rng(seed)
+    [T_emnenmf, RMSE] = remnenmf_not_full( W , X , Ginit , Finit, Omega_G, Omega_F, Phi_G, Phi_F , InnerMinIter , InnerMaxIter , Tmax , M_loop, F_theo, delta_measure, nu);
+    RMSE_offset_emnenmf(run,:) = RMSE(1,:);
+    RMSE_gain_emnenmf(run,:) = RMSE(2,:);    
+    
+end
+
+% remnenmf
+min_offset_remnenmf = min(RMSE_offset_remnenmf,[],1,'omitnan');
+med_offset_remnenmf = median(RMSE_offset_remnenmf,1,'omitnan');
+max_offset_remnenmf = max(RMSE_offset_remnenmf,[],1,'omitnan');
+
+min_gain_remnenmf = min(RMSE_gain_remnenmf,[],1,'omitnan');
+med_gain_remnenmf = median(RMSE_gain_remnenmf,1,'omitnan');
+max_gain_remnenmf = max(RMSE_gain_remnenmf,[],1,'omitnan');
+
+figure
+subplot(121)
+semilogy(T_remnenmf,min_offset_remnenmf,'b')
+hold on
+semilogy(T_remnenmf,med_offset_remnenmf,'b')
+o_e = semilogy(T_remnenmf,max_offset_remnenmf,'b');
+hold off
+
+subplot(122)
+semilogy(T_remnenmf,min_gain_remnenmf,'b')
+hold on
+semilogy(T_remnenmf,med_gain_remnenmf,'b')
+g_e = semilogy(T_remnenmf,max_gain_remnenmf,'b');
+hold off
+
+% emnenmf
+min_offset_emnenmf = min(RMSE_offset_emnenmf,[],1,'omitnan');
+med_offset_emnenmf = median(RMSE_offset_emnenmf,1,'omitnan');
+max_offset_emnenmf = max(RMSE_offset_emnenmf,[],1,'omitnan');
+
+min_gain_emnenmf = min(RMSE_gain_emnenmf,[],1,'omitnan');
+med_gain_emnenmf = median(RMSE_gain_emnenmf,1,'omitnan');
+max_gain_emnenmf = max(RMSE_gain_emnenmf,[],1,'omitnan');
+
+subplot(121)
+hold on
+semilogy(T_emnenmf,min_offset_emnenmf,'r')
+semilogy(T_emnenmf,med_offset_emnenmf,'r')
+o_i = semilogy(T_emnenmf,max_offset_emnenmf,'r');
+hold off
+
+subplot(122)
+hold on
+semilogy(T_emnenmf,min_gain_emnenmf,'r')
+semilogy(T_emnenmf,med_gain_emnenmf,'r')
+g_i = semilogy(T_emnenmf,max_gain_emnenmf,'r');
+hold off
+
+% adding title and labels
+subplot(121)
+title('RMSE offset')
+legend([o_e o_i],'new L and R','old L and R')
+
+subplot(122)
+title('RMSE gain')
+legend([g_e g_i],'new L and R','old L and R')
+

MATLAB/main_incal_vs_emnenmf.m

@@ -11,7 +11,7 @@ addpath(genpath(pwd))
 s_width = 50;  % Scene width
 s_length = 50; % Scene length
 N_Ref = 4; % Nb. of reference measurements
-N_Cpt = 25; % Nb. of mobile sensors
+N_Cpt = 100; % Nb. of mobile sensors
 Mu_beta = .9; % Mean sensors gain
 Mu_alpha = 5; % Mean sensors offset
 Bound_beta = [.01;1.5]; % Gain boundaries
@@ -20,12 +20,12 @@ Bound_alpha = [3.5;6.5]; % Offset boundaries
 MV = .5; % Missing Value prop.
 RV = 0.3; % RendezVous prop.
 var_n = 0; % Noise variance
-M_loop = 5; % Number of M loop per E step
+M_loop = 50; % Number of M loop per E step
 runs = 1; % Total number of runs
 
 %% Nesterov parameters
-InnerMinIter = 10;
-InnerMaxIter = 100;
+InnerMinIter = 5;
+InnerMaxIter = 50;
 Tmax = 30;
 
 %% 

MATLAB/main_incal_vs_remnenmf.m

@@ -11,13 +11,13 @@ addpath(genpath(pwd))
 s_width = 50;  % Scene width
 s_length = 50; % Scene length
 N_Ref = 4; % Nb. of reference measurements
-N_Cpt = 25; % Nb. of mobile sensors
+N_Cpt = 100; % Nb. of mobile sensors
 Mu_beta = .9; % Mean sensors gain
 Mu_alpha = 5; % Mean sensors offset
 Bound_beta = [.01;1.5]; % Gain boundaries
 Bound_alpha = [3.5;6.5]; % Offset boundaries
 
-MV = .5; % Missing Value prop.
+MV = .9; % Missing Value prop.
 RV = 0.3; % RendezVous prop.
 var_n = 0; % Noise variance
 M_loop = 50; % Number of M loop per E step
@@ -25,8 +25,11 @@ runs = 1; % Total number of runs
 
 %% Nesterov parameters
 InnerMinIter = 5;
-InnerMaxIter = 100;
-Tmax = 30;
+InnerMaxIter = 20;
+Tmax = 300;
+
+%% Compression parameters
+nu = 10;
 
 %% 
 delta_measure = 1;
@@ -44,11 +47,14 @@ RMSE_gain_incal = nan(runs, iter_max);
 for run = 1:runs
     % data generation
     [X, X_theo, W, F_theo, Omega_G, Omega_F, Phi_G, Phi_F, Ginit, Finit] = data_gen(s_width, s_length, run, N_Ref, N_Cpt, Mu_beta, Mu_alpha, Bound_beta, Bound_alpha, MV, RV, var_n);
+    seed = floor(rand*10000);
     % remnenmf
-    [T_remnenmf, RMSE] = remnenmf( W , X , Ginit , Finit, Omega_G, Omega_F, Phi_G, Phi_F , InnerMinIter , InnerMaxIter , Tmax , M_loop, F_theo, delta_measure);
+    rng(seed)
+    [T_remnenmf, RMSE] = remnenmf_not_full( W , X , Ginit , Finit, Omega_G, Omega_F, Phi_G, Phi_F , InnerMinIter , InnerMaxIter , Tmax , M_loop, F_theo, delta_measure, nu);
     RMSE_offset_remnenmf(run,:) = RMSE(1,:);
     RMSE_gain_remnenmf(run,:) = RMSE(2,:);
     % incal
+    rng(seed)
     [T_incal, RMSE] = IN_Cal( W , X , Ginit , Finit , Omega_G , Omega_F , Phi_G , Phi_F , F_theo , Tmax, delta_measure );
     RMSE_offset_incal(run,:) = RMSE(1,:);
     RMSE_gain_incal(run,:) = RMSE(2,:);