function [ G , F , R ] = SpAIN_Cal( W , X , G , F , Omega_G , Omega_F , Phi_G , Phi_F , D , k , lambda , Mean_F , mu , F_theo , N_iter )
% Author: Clément DORFFER
% Date: 26/11/2018
% @: clement.dorffer@ensta-bretagne.fr
% Goal: perform mobile sensor calibration using the Sparsity- and Average-Constrained Informed NMF calibraton technique. 
% If you use this code for research or educational purpose, please cite:
% 
% C. Dorffer, M. Puigt, G. Delmaire, G. Roussel, Informed Nonnegative Matrix Factorization Methods for 
% Mobile Sensor Network Calibration, IEEE Transactions on Signal and Information Processing over Networks, 
% Volume 4, Issue 4, pp. 667-682, December 2018.  
%

W2 = W.^2;
tilde_W2 = [W2,lambda*ones(size(W,1),1)];

Omega_Bar_G = ones(size(G))-Omega_G;
Omega_Bar_F = ones(size(F))-Omega_F;
Delta_G = G.*Omega_Bar_G;
Delta_F = F.*Omega_Bar_F;
m = size(F,2)-1;

R = zeros(N_iter+1,1);
R(1) = norm(F(2,1:end-1)-F_theo(2,1:end-1),2)/sqrt(size(F_theo,2)-1);

for i = 1 : N_iter
    
    x = OMP(G(:,2),D,k);
    tilde_X = [X,D*x];
    tilde_F = [F,[0;1]];
    
    % updating G
    Delta_G = Updt_Delta_G( tilde_W2 , tilde_X , Delta_G , Phi_G , tilde_F , Omega_Bar_G );
    G = Phi_G + Delta_G;
    
    % updating F
    Delta_F = Updt_Delta_F( W2 , X , G , Delta_F , Phi_F , Omega_Bar_F , Mean_F , mu , m );
    F = Phi_F + Delta_F;
    
    R(i+1) = norm(F(2,1:end-1)-F_theo(2,1:end-1),2)/sqrt(size(F_theo,2)-1);
    
end

end



function [ Delta_F ] = Updt_Delta_F( W , X , G , Delta_F , Phi_F , Omega_B_F , Mean_F , mu , m )
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Update rule for Delta_F %
%%%%%%%%%%%%%%%%%%%%%%%%%%%
Delta_F = Delta_F.*(Omega_B_F.*(G'*(W.*secu_plus(X-G*Phi_F)))+mu/m*diag(Mean_F)*Omega_B_F)./(Omega_B_F.*(G'*(W.*(G*Delta_F)))+mu/m^2*diag(Delta_F*ones(m+1,1))*Omega_B_F);
Delta_F(isnan(Delta_F))=0;
end



function [ Delta_G ] = Updt_Delta_G( W , X , Delta_G , Phi_G , F , Omega_B_G )
%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Update rule for Delta_G %
%%%%%%%%%%%%%%%%%%%%%%%%%%%
Delta_G = Delta_G.*(Omega_B_G).*((W.*secu_plus(X-Phi_G*F))*F')./((W.*(Delta_G*F)*F')); % mise à jour
Delta_G(isnan(Delta_G))=0;
end



function [toto] = secu_plus(tutu,s)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Goal: Security in the NMF procedure which project the negative data to
% epsilon (small user-defined threshold)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(nargin<2)
    s=1.e-12;
end
toto=max(tutu,s);
toto(isnan(tutu)) = 0;
end