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      = config.Bound_phen;
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;
var_n           = config.var_n;

%% Scene simulation 

n_pic = 15;
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([ 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);
    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_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 = (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));
    list_nosen = 1:N_sousCpt;
    list_nosen(sen) = [];
    maxPhen_nosen = norm(Bound_phen(2*list_nosen));
    for sor = list_nosen
        f_theo_nosen = rand(1,N_Cpt).*C/(sqrt(N_sousCpt)*maxPhen_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);

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);

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
    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);
    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;

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);
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);
    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
    end
end
% Finit = F_theo;

end