fast_in_situ_calibration/MATLAB/common/data_gen.m

function [X, X_theo, W, F_theo, Omega_G, Omega_F, Phi_G, Phi_F, Ginit, Finit] = data_gen(config, run) 

rng(run+1306) % Random seed

s_width         = config.sceneWidth;
s_length        = config.sceneLength;
N_Ref           = config.numRef;
N_Cpt           = config.numSensor;
N_sousCpt       = config.numSubSensor;
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;
Bound_alpha     = config.Bound_alpha;
gamma           = config.gamma;
MV              = config.mvR;
RV              = config.rdvR;
noise           = config.noise;

%% Scene simulation 

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(:));
G_theo = 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([0.2,0.45]);
        g = g + mvnpdf(xxyy,mu,sig);
    end
    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


%% Sensors simulation

F_theo = zeros(N_sousCpt+1, N_sousCpt*(N_Cpt+1));

for sensor = 1:N_sousCpt
    F_theo(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.55*randn(1,N_Cpt))), 0);
end

for sen = 1:N_sousCpt
    f_theo = 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);
    F_theo(sen+1, (sen-1)*(N_Cpt+1)+1:sen*(N_Cpt+1)) = f_theo;
    
    C = 10^(-gamma/20)*mean(Bound_phen(sen,:),2)*f_theo(1:end-1);
    list_nosen = 1:N_sousCpt;
    list_nosen(sen) = [];
    meanPhen_nosen = norm(mean(Bound_phen(list_nosen,:)));
    for sor = list_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
end

% F_theo = [];
% 
% for sensor = 1:N_sousCpt
%     F_theo = cat(2, F_theo, cat(2, max(Bound_beta(1), min(Bound_beta(2), Mu_beta+.5*randn(1,N_Cpt))), 0));
% end
% 
% 
% for sen = 1:N_sousCpt
%     f_theo = [];
%     for sor = 1:N_sousCpt
%         if sen==sor
%             f_theo = cat(2, f_theo, cat(2, max(Bound_alpha(1), min(Bound_alpha(2), Mu_alpha+.5*randn(1,N_Cpt))), sen==sor));
%         else
%             f_theo = cat(2, f_theo, cat(2, 0*max(Bound_alpha(1), min(Bound_alpha(2), Mu_alpha+.5*randn(1,N_Cpt))), sen==sor));
%         end
%     end
%     F_theo = cat(1, F_theo, f_theo);
% end

%% Data simulation

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

W = zeros(s_n,N_Cpt+1);

idx_Ref = randperm(s_n);
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)*(s_n-N_Ref)); % Nb. of measurements in data matrix X
xCpt = 1 : s_n;
xCpt(idx_Ref) = []; % Reference free locations
[xx,yy] = meshgrid(xCpt,1:N_Cpt); % Possibly sensed locations

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

W = repmat(W, 1, N_sousCpt);

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])
Omega_F = zeros(N_sousCpt+1,N_sousCpt*(N_Cpt+1)); 
Omega_F(:,(N_Cpt+1)*(1:N_sousCpt)) = 1; % Mask on known values in F (see eq.(15) of [1])
Phi_G = [ones(s_n,1),X(:,(N_Cpt+1)*(1:N_sousCpt))]; % Known values in G (see eq.(14) of [1])
Phi_F = F_theo .* Omega_F; % Known values in F (see eq.(15) of [1])

%% Initialization

Ginit = ones(s_n,1);
for sensor = 1:N_sousCpt
    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
    Ginit = cat(2, Ginit, g);
end

Ginit = (1-Omega_G).*Ginit+Phi_G;
% Ginit = G_theo;

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.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)))+eps, 1);
    Finit(sen+1, (sen-1)*(N_Cpt+1)+1:sen*(N_Cpt+1)) = finit;
    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