Final update

MATLAB/calib_meth/EMNeNMF/emnenmf.m

@@ -1,5 +1,7 @@
 function [ T , RMSE , MERs ] = emnenmf( X, X_theo, W, F_theo, Omega_G, Omega_F, Phi_G, Phi_F, G, F, config)
 
+disp(['RMSE 2ème colonne de Ginit : ',num2str(norm(G(:,2) - X_theo(:,end),2)/sqrt(config.sceneWidth*config.sceneLength))])
+
 %% loading the config parameters 
 Tmax = config.Tmax;
 delta_measure = config.delta_measure;
@@ -50,17 +52,7 @@ while toc<Tmax
     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;
 
@@ -71,6 +63,17 @@ while toc<Tmax
         if(iterF<=InnerMinIter)
             StoppingCritF = 1.e-1*StoppingCritF;
         end
+        
+        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
 
         if toc - i*delta_measure >= delta_measure
             i = i+1;
@@ -80,7 +83,7 @@ while toc<Tmax
             [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);
+            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);
         

MATLAB/calib_meth/INCAL/IN_Cal.m

@@ -6,8 +6,7 @@ delta_measure = config.delta_measure;
 
 %% 
 
-[~, num_sensor] = size(F);
-num_sensor = num_sensor-1;
+num_sensor = config.numSensor;
 
 W2 = W.^2;
 Omega_Bar_G = ones(size(G))-Omega_G;
@@ -23,7 +22,7 @@ MERs = nan(1+config.numSubSensor,iter_max);
 tic
 i = 1;
 T(i) = toc;
-RMSE(:,i) = vecnorm(F(:,1:end-1)-F_theo(:,1:end-1),2,2)/sqrt(num_sensor);
+RMSE(:,i) = vecnorm(F(:,1:num_sensor)- F_theo(:,1:num_sensor),2,2)/sqrt(num_sensor);
 
 while toc<Tmax
     
@@ -43,7 +42,7 @@ while toc<Tmax
         [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);
+        RMSE(:,i) = vecnorm(F(:,1:num_sensor)- F_theo(:,1:num_sensor),2,2)/sqrt(num_sensor);
     end
     
 end

MATLAB/calib_meth/REMNeNMF/remnenmf.m

@@ -1,4 +1,4 @@
-function [ T , RMSE ] = remnenmf( X, X_theo, W, F_theo, Omega_G, Omega_F, Phi_G, Phi_F, G, F, config)
+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;
@@ -17,11 +17,11 @@ 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;
+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);
@@ -45,14 +45,14 @@ T_M = [];
 tic
 i = 1;
 niter = 0;
-RMSE(:,i) = vecnorm(F(:,1:end-1)- F_theo(:,1:end-1),2,2)/sqrt(num_sensor);
+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 i>1
+%     if mod(i-1,5) ~= 0
 %         [L,R]=RSI_compression(X,r,L,R);
 %     else
 %         [L,R]=RSI_compression(X,r);
@@ -101,8 +101,10 @@ while toc<Tmax
             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);
+            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);
         

MATLAB/calib_meth/test/remnenmf_full.m

@@ -1,160 +0,0 @@
-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

@@ -17,8 +17,7 @@ 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;
+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);
@@ -46,7 +45,7 @@ T_M = [];
 tic
 i = 1;
 niter = 0;
-RMSE(:,i) = vecnorm(F(:,1:end-1)- F_theo(:,1:end-1),2,2)/sqrt(num_sensor);
+RMSE(:,i) = vecnorm(F(:,1:num_sensor)- F_theo(:,1:num_sensor),2,2)/sqrt(num_sensor);
 T(i) = toc;
 while toc<Tmax
     
@@ -105,7 +104,7 @@ while toc<Tmax
             [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);
+            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);
         

MATLAB/common/data_gen.m

@@ -7,7 +7,7 @@ s_length        = config.sceneLength;
 N_Ref           = config.numRef;
 N_Cpt           = config.numSensor;
 N_sousCpt       = config.numSubSensor;
-Bound_phen      = config.Bound_phen;
+Bound_phen      = reshape(config.Bound_phen(1:2*N_sousCpt),[2 N_sousCpt])';
 Mu_beta         = config.Mu_beta;
 Mu_alpha        = config.Mu_alpha;
 Bound_beta      = config.Bound_beta;
@@ -15,11 +15,11 @@ Bound_alpha     = config.Bound_alpha;
 gamma           = config.gamma;
 MV              = config.mvR;
 RV              = config.rdvR;
-var_n           = config.var_n;
+noise           = config.noise;
 
 %% Scene simulation 
 
-n_pic = 15;
+n_pic = 5;
 s_n = s_width*s_length; % Total number of areas in the scene
 [xx,yy] = meshgrid((-1:2/(s_width-1):1),(-1:2/(s_length-1):1));
 xxyy = cat(2,xx(:),yy(:));
@@ -28,12 +28,10 @@ for sensor = 1:N_sousCpt
     g = zeros(s_n,1);
     for pic = 1:n_pic
         mu = 2*(rand(1,2)-0.5);
-        sig = diag([ Bound_phen((sensor-1)*2+1) + (Bound_phen((sensor-1)*2+2) - Bound_phen((sensor-1)*2+1))*abs(randn()+0.5) , ...
-            Bound_phen((sensor-1)*2+1) + (Bound_phen((sensor-1)*2+2) - Bound_phen((sensor-1)*2+1))*abs(randn()+0.5) ]);
+        sig = diag([0.2,0.45]);
         g = g + mvnpdf(xxyy,mu,sig);
     end
-    g = (Bound_phen((sensor-1)*2+2) - Bound_phen((sensor-1)*2+1))/(max(g)-min(g))*(g-min(g)) + Bound_phen((sensor-1)*2+1);
-    % .5*(g/max(g))+1e-5
+    g = (Bound_phen(sensor,2) - Bound_phen(sensor,1))/(max(g)-min(g))*(g-min(g)) + Bound_phen(sensor,1);
     G_theo = cat(2, G_theo, g);
 end
 
@@ -54,12 +52,12 @@ for sen = 1:N_sousCpt
         Mu_alpha(sen)+(Bound_alpha((sen-1)*2+2)-Bound_alpha((sen-1)*2+1))*0.25*randn(1,N_Cpt))), 1);
     F_theo(sen+1, (sen-1)*(N_Cpt+1)+1:sen*(N_Cpt+1)) = f_theo;
     
-    C = (1-gamma)/gamma*(F_theo(1,(sen-1)*(N_Cpt+1)+1:sen*(N_Cpt+1)-1)+Bound_phen((sen-1)*2+1)*f_theo(1:end-1));
+    C = 10^(-gamma/20)*mean(Bound_phen(sen,:),2)*f_theo(1:end-1);
     list_nosen = 1:N_sousCpt;
     list_nosen(sen) = [];
-    maxPhen_nosen = norm(Bound_phen(2*list_nosen));
+    meanPhen_nosen = norm(mean(Bound_phen(list_nosen,:)));
     for sor = list_nosen
-        f_theo_nosen = rand(1,N_Cpt).*C/(sqrt(N_sousCpt)*maxPhen_nosen);
+        f_theo_nosen = rand(1,N_Cpt).*C/(sqrt(N_sousCpt)*meanPhen_nosen);
         other_f_theo = cat(2, f_theo_nosen, 0);
         F_theo(sor+1, (sen-1)*(N_Cpt+1)+1:sen*(N_Cpt+1)) = other_f_theo;
     end
@@ -116,10 +114,18 @@ end
 
 W = repmat(W, 1, N_sousCpt);
 
-N = var_n*randn(s_n,N_sousCpt*(N_Cpt+1)); % Noise simulation
-N(:,(N_Cpt+1)*(1:N_sousCpt)) = 0;
-N = max(N,-X_theo);
-
+if isinf(noise)
+    N = zeros(s_n,N_sousCpt*(N_Cpt+1));
+else
+    N = randn(s_n,N_sousCpt*(N_Cpt+1)); % Noise simulation
+    N(:,(N_Cpt+1)*(1:N_sousCpt)) = 0; % No noise for the references
+    for i = 1:N_sousCpt
+        noise_i = N(:,(i-1)*(N_Cpt+1)+1:i*(N_Cpt+1)-1);
+        X_i = X_theo(:,(i-1)*(N_Cpt+1)+1:i*(N_Cpt+1)-1);
+        N(:,(i-1)*(N_Cpt+1)+1:i*(N_Cpt+1)-1) = (max(X_i(:))-min(X_i(:)))/(max(noise_i(:))-min(noise_i(:)))*10^(-noise/20)*noise_i;
+    end
+    N = max(N,-X_theo);
+end
 X = W.*(X_theo+N); % Data matrix X
 
 Omega_G = [ones(s_n,1),W(:,(N_Cpt+1)*(1:N_sousCpt))]; % Mask on known values in G (see eq.(14) of [1])
@@ -132,17 +138,19 @@ Phi_F = F_theo .* Omega_F; % Known values in F (see eq.(15) of [1])
 
 Ginit = ones(s_n,1);
 for sensor = 1:N_sousCpt
-    g = zeros(s_n,1);
-    for pic = 1:n_pic
-        mu = 2*(rand(1,2)-0.5);
-        sig = diag([ Bound_phen((sensor-1)*2+1) + (Bound_phen((sensor-1)*2+2) - Bound_phen((sensor-1)*2+1))*abs(randn()+0.5) , ...
-            Bound_phen((sensor-1)*2+1) + (Bound_phen((sensor-1)*2+2) - Bound_phen((sensor-1)*2+1))*abs(randn()+0.5) ]);
-        g = g + mvnpdf(xxyy,mu,sig);
+    if ~isempty(find(sensor==config.CalibratedSensor,1))
+        g = X(:, (sensor-1)*(N_Cpt+1)+1:sensor*(N_Cpt+1)-1);
+        g(~W(:,(sensor-1)*(N_Cpt+1)+1:sensor*(N_Cpt+1)-1)) = NaN;
+        g = 1/Mu_alpha(sensor)*(nanmean(g,2)-Mu_beta(sensor));
+        g(isnan(g))=0;
+    else
+        g = Phi_G(:,sensor+1);
+        g(g==0) = mean(g(g~=0));
+        g(Phi_G(:,sensor+1)==0) = abs(1+0.05*randn(size(find(Phi_G(:,sensor+1)==0)))).*g(Phi_G(:,sensor+1)==0);
     end
-    g = (Bound_phen((sensor-1)*2+2) - Bound_phen((sensor-1)*2+1))/(max(g)-min(g))*(g-min(g)) + Bound_phen((sensor-1)*2+1);
-    % .5*(g/max(g))+1e-5
     Ginit = cat(2, Ginit, g);
 end
+
 Ginit = (1-Omega_G).*Ginit+Phi_G;
 % Ginit = G_theo;
 
@@ -151,28 +159,26 @@ Finit = zeros(N_sousCpt+1, N_sousCpt*(N_Cpt+1));
 for sensor = 1:N_sousCpt
     Finit(1,(sensor-1)*(N_Cpt+1)+1:sensor*(N_Cpt+1)) = ...
         cat(2, max(Bound_beta((sensor-1)*2+1), min(Bound_beta((sensor-1)*2+2), ...
-        Mu_beta(sensor)+(Bound_beta((sensor-1)*2+2)-Bound_beta((sensor-1)*2+1))*0.25*randn(1,N_Cpt))), 0);
+        Mu_beta(sensor)+(Bound_beta((sensor-1)*2+2)-Bound_beta((sensor-1)*2+1))*0.55*randn(1,N_Cpt)))+eps, 0);
 end
 
 for sen = 1:N_sousCpt
     finit = cat(2, max(Bound_alpha((sen-1)*2+1),...
         min( Bound_alpha((sen-1)*2+2), ...
-        Mu_alpha(sen)+(Bound_alpha((sen-1)*2+2)-Bound_alpha((sen-1)*2+1))*0.25*randn(1,N_Cpt))), 1);
+        Mu_alpha(sen)+(Bound_alpha((sen-1)*2+2)-Bound_alpha((sen-1)*2+1))*0.25*randn(1,N_Cpt)))+eps, 1);
     Finit(sen+1, (sen-1)*(N_Cpt+1)+1:sen*(N_Cpt+1)) = finit;
-    
-%     C = (1-gamma)/gamma*(Finit(1,(sen-1)*(N_Cpt+1)+1:sen*(N_Cpt+1)-1)+Bound_phen((sen-1)*2+1)*finit(1:end-1));
-%     list_nosen = 1:N_sousCpt;
-%     list_nosen(sen) = [];
-%     maxPhen_nosen = norm(Bound_phen(2*list_nosen));
-    for sor = 1:N_sousCpt
-        if sen~=sor
-%             finit_nosen = rand(1,N_Cpt).*C/(sqrt(N_sousCpt)*maxPhen_nosen);
-%             other_finit = cat(2, finit_nosen, 0);
-%             Finit(sor+1, (sen-1)*(N_Cpt+1)+1:sen*(N_Cpt+1)) = other_finit;
-            Finit(sor+1, (sen-1)*(N_Cpt+1)+1:sen*(N_Cpt+1)) = zeros(1,N_Cpt+1); % initialization of other dependencies at zero
-        end
+    C = 10^(-0/20)*mean(Bound_phen(sen,:),2)*finit(1:end-1);
+    list_nosen = 1:N_sousCpt;
+    list_nosen(sen) = [];
+    meanPhen_nosen = norm(mean(Bound_phen(list_nosen,:)));
+    for sor = list_nosen
+        finit_nosen = rand(1,N_Cpt).*C/(sqrt(N_sousCpt)*meanPhen_nosen)+eps;
+        other_finit = cat(2, finit_nosen, 0);
+        Finit(sor+1, (sen-1)*(N_Cpt+1)+1:sen*(N_Cpt+1)) = other_finit;
     end
 end
+
+
 % Finit = F_theo;
 
 end

MATLAB/config.json

@@ -1,22 +1,22 @@
 {
-    "sceneWidth" : 10,
-    "sceneLength" : 10,
+    "sceneWidth" : 20,
+    "sceneLength" : 20,
     "numSensor" : 100,
     "numSubSensor" : 1,
     "CalibratedSensor" : [1],
-    "numRef" : 4,
+    "numRef" : 10,
     "rdvR" : 0.3,
     "mvR" : 0.5,
     "calibrationMethods" : ["IN_Cal"],
     "numRuns" : 1,
-    "Bound_phen" : [0.05,0.35 , 10,20],
-    "Mu_beta" : [2 , 5],
-    "Mu_alpha" : [30 , 0.5],
-    "Bound_beta" : [1,3 , 4.5,5.5],
-    "Bound_alpha" : [20,40 , 0.4,0.6],
-    "gamma" : 0.8,
-    "var_n" : 0,
-    "Tmax" : 30,
+    "Bound_phen" : [0,0.5 , 0,0.5, 0,0.5],
+    "Mu_beta" : [0.9 , 0.9, 0.9],
+    "Mu_alpha" : [5 , 5, 5],
+    "Bound_beta" : [0,1.5 , 0,1.5, 0,1.5],
+    "Bound_alpha" : [3.5,6.5 , 3.5,6.5, 3.5,6.5],
+    "gamma" : 15,
+    "noise" : Inf,
+    "Tmax" : 60,
     "nu" : 12,
     "save2dat" : 1,
     "display" : 1,
@@ -24,4 +24,4 @@
     "M_loop" : 50,
     "InnerMinIter" : 5,
     "InnerMaxIter" : 20
-}
+}

MATLAB/config2.json

@@ -3,19 +3,19 @@
     "sceneLength" : 10,
     "numSensor" : 100,
     "numSubSensor" : 2,
-    "CalibratedSensor" : [1],
+    "CalibratedSensor" : [1,2],
     "numRef" : 4,
     "rdvR" : 0.3,
     "mvR" : 0.5,
     "calibrationMethods" : ["emnenmf"],
     "numRuns" : 1,
-    "Bound_phen" : [0.05,0.35 , 10,20],
-    "Mu_beta" : [2 , 5],
-    "Mu_alpha" : [30 , 0.5],
-    "Bound_beta" : [1,3 , 4.5,5.5],
-    "Bound_alpha" : [20,40 , 0.4,0.6],
-    "gamma" : 0.8,
-    "var_n" : 0,
+    "Bound_phen" : [0.05,0.35 , 0.05,0.35],
+    "Mu_beta" : [2 , 2],
+    "Mu_alpha" : [30 , 30],
+    "Bound_beta" : [1,3 , 1,3],
+    "Bound_alpha" : [20,40 , 20,40],
+    "gamma" : 15,
+    "noise" : Inf,
     "Tmax" : 30,
     "nu" : 12,
     "save2dat" : 1,

MATLAB/main.m

@@ -4,5 +4,17 @@ clear all
 
 addpath(genpath(pwd))
 
-run_config("config2.json")
-% run_config_forget_other_sen("config2.json")
+%% run config with maag et al values of the phenomenon
+
+% run_config("config3.json")
+% res = run_config("config.json");
+% run_config("config4.json") % Not working for IN_CAL ???
+% run_config("config2.json")
+% run_config_forget_other_sen("config3.json")
+
+% run_config("SNR_100dB.json")
+% run_config("SNR_infdB.json")
+
+% res = run_config("config10x10.json");
+res = run_config("test_rapport.json");
+% res = run_config("30x30farincal.json");

MATLAB/run_config.m

@@ -1,4 +1,4 @@
-function [res] = run_config(varargin)
+function [RMSEs] = run_config(varargin)
     if nargin == 0
         file_name = "default_config.json";
     elseif nargin == 1
@@ -18,14 +18,15 @@ function [res] = run_config(varargin)
         [X, X_theo, W, F_theo, Omega_G, Omega_F, Phi_G, Phi_F, Ginit, Finit] = data_gen(config, run);
         seed = floor(rand*10000);
         
-        % We want to only calibrate the sensor of interest
-        X = X(:,1:config.numSensor+1);
-        X_theo = X_theo(:,1:config.numSensor+1);
-        W = W(:,1:config.numSensor+1);
-        F_theo = F_theo(:,1:config.numSensor+1);
-        Omega_F = Omega_F(:,1:config.numSensor+1);
-        Phi_F = Phi_F(:,1:config.numSensor+1);
-        Finit = Finit(:,1:config.numSensor+1);
+        % We want to only calibrate the sensor(s) of interest according to config.CalibratedSensor
+        idx_to_calibrate = kron(config.CalibratedSensor'-1,(config.numSensor+1)*ones(1,config.numSensor+1))+kron(ones(1,numel(config.CalibratedSensor)),1:config.numSensor+1);
+        X = X(:,idx_to_calibrate);
+        X_theo = X_theo(:,idx_to_calibrate);
+        W = W(:,idx_to_calibrate);
+        F_theo = F_theo(:,idx_to_calibrate);
+        Omega_F = Omega_F(:,idx_to_calibrate);
+        Phi_F = Phi_F(:,idx_to_calibrate);
+        Finit = Finit(:,idx_to_calibrate);
 
 %         X = X(:,1:config.numSensor+1);
 %         X_theo = X_theo(:,1:config.numSensor+1);
@@ -58,16 +59,77 @@ function [res] = run_config(varargin)
     end
     if config.display
         figure
+        colors = ['b','r','g','k','y','c','m'];
         for it_meth = 1:length(methods)
-            for measure = 1:(config.numSubSensor+1)
-                subplot(length(methods),config.numSubSensor+1,(it_meth-1)*(config.numSubSensor+1)+measure)
-                for k = 1:size(RMSEs{it_meth,measure},1)
-                    semilogy(RMSEs{it_meth,measure}(k,:))
-                    hold on
-                end
-                hold off
+            for measure = 2 % Only the RMSE for the gain
+%                 subplot(1,config.numSubSensor+1,measure) % Only the RMSE for the gain
+                enveloppe_min = min(RMSEs{it_meth,measure},[],1);
+                enveloppe_max = max(RMSEs{it_meth,measure},[],1);
+                semilogy(enveloppe_max,colors(it_meth),'HandleVisibility','off');
+                hold on
+                semilogy(enveloppe_min,colors(it_meth),'DisplayName',"Enveloppe "+strrep(methods{it_meth},'_','\_'));
+                legend show
             end
         end
+        for it_meth = 1:length(methods)
+            for measure = 2 % Only the RMSE for the gain
+                med = median(RMSEs{it_meth,measure});
+                semilogy(med,"--"+colors(it_meth),'LineWidth',3,'DisplayName',"Médiane "+strrep(methods{it_meth},'_','\_'));
+                hold on
+                legend show
+            end
+        end 
+        xlabel("Time (s)")
+        ylabel("RMSE")
+        grid on
+        set(gca,'FontUnits','points','FontSize',18,'FontName','Times')
+        fig = gcf;
+        legend('AutoUpdate','off')
+        Lines = findall(gcf,'Type','line');
+        nenv = numel(Lines)/3;
+        for k = 1:nenv
+            line1 = Lines(3*nenv-(k-1)*2);
+            line2 = Lines(3*nenv-(k-1)*2-1);
+            fill([line1.XData flip(line2.XData)],[line1.YData flip(line2.YData)],colors(k),'LineStyle','none')
+        end
+
+        alpha(0.25)
+%         for it_meth = 1:length(methods)
+%             for measure = 2 % Only the RMSE for the gain
+% %                 subplot(1,config.numSubSensor+1,measure) % Only the RMSE for the gain
+%                 for k = 1:size(RMSEs{it_meth,measure},1)
+%                     if k == size(RMSEs{it_meth,measure},1)
+%                         semilogy(RMSEs{it_meth,measure}(k,:),colors(it_meth),'DisplayName',"Enveloppe "+strrep(methods{it_meth},'_','\_'));
+%                     else
+%                         semilogy(RMSEs{it_meth,measure}(k,:),colors(it_meth),'HandleVisibility','off');
+%                     end
+%                     hold on
+%                 end
+%                 legend show
+%                 set(gca,'FontUnits','points','FontSize',12,'FontName','Times')
+%             end
+%         end
     end
+
+%     if config.display
+%         figure
+%         for it_meth = 1:length(methods)
+%             for measure = 1:(config.numSubSensor+1)
+%                 subplot(length(methods),config.numSubSensor+1,(it_meth-1)*(config.numSubSensor+1)+measure)
+%                 for k = 1:size(RMSEs{it_meth,measure},1)
+%                     semilogy(RMSEs{it_meth,measure}(k,:))
+%                     hold on
+%                 end
+%                 hold off
+%             end
+%         end
+%     end
 end
 
+% set(gcf,'OuterPosition',[200 200 400 440]);
+% axis square
+% pos = get(gcf,'paperposition');
+% set(gcf,'paperposition',[pos(1),pos(2), 4, 4]);
+
+% print -depsc epsFig
+

MATLAB/run_config_forget_other_sen.m

@@ -1,4 +1,4 @@
-function [res] = run_config(varargin)
+function [res] = run_config_forget_other_sen(varargin)
     if nargin == 0
         file_name = "default_config.json";
     elseif nargin == 1