% This matlab script provides a demonstration of the informed and
% constrained NMF algorithms proposed in the article:

% [1] C. DORFFER, M. Puigt, G. Delmaire and G. Roussel, "Informed Nonnegative
% Matrix Factorization Methods for Mobile Sensor Network Calibration", in
% IEEE Transactions on Signal and Information Processing over Networks,
% vol. 4, no. 4, pp. 667-682, Dec. 2018.
%
% Author: Clément DORFFER
% Date: 26/11/2018
% @: clement.dorffer@ensta-bretagne.fr



%%

clear
close all
clc

addpath('functions/')
rng(3)


%% Simulation parameters

N_Ref = 3; % Nb. of reference measurements
N_Cpt = 25; % 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; % RendezVous prop.
var_n = 0; % Noise variance


%% Scene simulation

[xx,yy] = meshgrid((-1:2/9:1),(-1:2/9:1));
sig = [1;1];
y=@(xx,yy) exp(-[xx;yy]'*diag(sig)*[xx;yy]);
for i = 1:100
    g(i) = y(xx(i),yy(i));
end
g = g-min(g);
g = .5*(g/max(g))+1e-5;
G_theo = [ones(100,1),g']; % Theoretical matrix G (see eq.(3) of [1])

figure
subplot(221)
imagesc(reshape(g,10,10))


%% Sensors simulation

F_theo = [max(Bound_beta(1),min(Bound_beta(2),Mu_beta+.5*randn(1,N_Cpt)));...
     max(Bound_alpha(1),min(Bound_alpha(2),Mu_alpha+.5*randn(1,N_Cpt)))];
F_theo = [F_theo,[0;1]]; % Theoretical matrix F (see eq.(3) of [1])


%% Data simulation

X_theo = G_theo*F_theo; %  Theoretical matrix X (see eq.(2) of [1])

W = zeros(100,N_Cpt+1);

idx_Ref = randperm(100);
idx_Ref = idx_Ref(1:N_Ref); % Reference measurement locations
W(idx_Ref,end) = 1;

N_RV = round(N_Cpt*RV); % Nb. of sensors having a RendezVous
idx_CptRV = randperm(N_Cpt);
idx_CptRV = idx_CptRV(1:N_RV); % Selection of sensors having a RendezVous
idx_RefRV = randi(N_Ref,1,N_Cpt);
idx_RefRV = idx_Ref(idx_RefRV(1:N_RV)); % Selection of the references for each RendezVous

for i = 1 : N_RV
    W(idx_RefRV(i),idx_CptRV(i)) = 1;
end

N_data = round((1-MV)*(N_Cpt)*(100-N_Ref)); % Nb. of measurements in data matrix X
xCpt = 1 : 100;
xCpt(idx_Ref) = []; % Reference free locations
[xx,yy] = meshgrid(xCpt,1:N_Cpt); % Possibly sensed locations

idx_data = randperm((100-N_Ref)*N_Cpt);
for i = 1 : N_data
    W(xx(idx_data(i)),yy(idx_data(i))) = 1; % Sensor measurement placement
end

N = var_n*randn(100,N_Cpt+1); % Noise simulation
N(:,end) = 0;
N = max(N,-X_theo);
SNR = snr(X_theo(W~=0),N(W~=0));

X = W.*(X_theo+N); % Data matrix X
subplot(222)
imagesc(X)


%% Calibration parameters

% % Common parameters
N_iter = 1.e5; % Maximum nb. of iterations
Omega_G = [ones(100,1),W(:,end)]; % Mask on known values in G (see eq.(14) of [1])
Omega_F = [zeros(2,N_Cpt),[1;1]]; % Mask on known values in F (see eq.(15) of [1])
Phi_G = [ones(100,1),X(:,end)]; % Known values in G (see eq.(14) of [1])
Phi_F = [zeros(2,N_Cpt),[0;1]]; % Known values in F (see eq.(15) of [1])
Ginit = abs(randn(100,2)); % Initial matrix G
Ginit = (1-Omega_G).*Ginit+Phi_G;
Finit = [max(Bound_beta(1),min(Bound_beta(2),Mu_beta+.5*randn(1,N_Cpt)));...
     max(Bound_alpha(1),min(Bound_alpha(2),Mu_alpha+.5*randn(1,N_Cpt)))];
Finit = [Finit,[0;1]]; % Initial matrix F
Finit = (1-Omega_F).*Finit+Phi_F;

% % Parameters for the "average constrained" approach (called ACIN-Cal in [1])
Mean_F = mean(F_theo(:,1:N_Cpt),2);
mu = 10; % Regularization weight

% % Parameters for the Sparsity based regularization (called SpIN-Cal in [1])
% Dictionary construction
D = real(ifft(diag(ones(100,1))));
D = [D(:,1:15),g',D(:,1:15),g'];
D(51:end,1:16) = 0;
D(1:50,17:end) = 0;
D = D*diag(1./sqrt(sum(D.^2))); % Atoms normalisation

k = 2; % Nb. of atoms to be choosen by the OMP
lambda = 10; % Regularization weight


%% Calibration

% % IN-Cal
[ G_IN_Cal , F_IN_Cal , RMSE_IN_Cal ] = IN_Cal( W , X , Ginit , Finit , Omega_G , Omega_F , Phi_G , Phi_F , F_theo , N_iter );

% % ACIN-Cal
[ G_ACIN_Cal , F_ACIN_Cal , RMSE_ACIN_Cal ] = ACIN_Cal( W , X , Ginit , Finit , Omega_G , Omega_F , Phi_G , Phi_F , Mean_F , mu , F_theo , N_iter );

% % SpIN-Cal
[ G_SpIN_Cal , F_SpIN_Cal , RMSE_SpIN_Cal ] = SpIN_Cal( W , X , Ginit , Finit , Omega_G , Omega_F , Phi_G , Phi_F , D , k , lambda , F_theo , N_iter );

% % SpAIN-Cal
[ G_SpAIN_Cal , F_SpAIN_Cal , RMSE_SpAIN_Cal ] = SpAIN_Cal( W , X , Ginit , Finit , Omega_G , Omega_F , Phi_G , Phi_F , D , k , lambda , Mean_F , mu , F_theo , N_iter );


subplot(223)
semilogy(RMSE_IN_Cal)
axis([0 N_iter 1.e-16 1])
hold all
semilogy(RMSE_ACIN_Cal)
semilogy(RMSE_SpIN_Cal)
semilogy(RMSE_SpAIN_Cal)