%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%   Generation of Figure 5 of the Read et al. manuscript (under review)          %%%%
%%   Code edited 14 Sep 2015							  %%%%
%%   contact: qread@utk.edu or M.B.Eppinga@uu.nl for questions                    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% INITIALIZATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Clear the Matlab workspace
clear all

% Prepare the genotype dataset
AllGenotypes = [-1 -1 -1; 1 -1 -1; -1 1 -1; 1 1 -1; -1 -1 1; 1 -1 1; -1 1 1; 1 1 1];  %% All theoretically possible combinations (max 3 markers)
NoM          =   3;                                                                   %% Number of Markers
Genotypes    = AllGenotypes(1:(2^NoM),1:NoM);                                         %% Select relevant number of genotypes (rows) and markers (columns)

% Define details of the simulation
SampleSizes = [10 100 1000 10000 100000];                                             %% Select the sample sizes for the simulations (i.e. the x coordinates)
sig_all     = [0.1 0.25 0.5 1 2 5 10];                                                %% Levels of genotypic variance (equal for both spp.)
NoReps      = 100;                                                                    %% Number of replicate simulations for each sample size
H2          = 0.5;                                                                    %% Heritability parameter
mu1         = 0;                                                                      %% Mean trait value species 1
mu2         = 0;                                                                      %% Mean trait value species 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% START SIMULATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Create a waitbar
h = waitbar(0,'Simulation running');

% Start a loop repeating the analysis for all levels of genotypic variance
[X1 X2]      = size(sig_all);
[X3 X4]      = size(SampleSizes);
DeltamuTot   = zeros(X4,X2);
DeltamuError = zeros(X4,X2);
for sig_i    = 1:X2
    sig      = sig_all(sig_i); 
    
    % Start a loop to calculate species differences for each sample size
    for x_j        = 1:X4                                  
        SampleSize = SampleSizes(x_j);
        S1         = zeros(1,SampleSize);
        S2         = zeros(1,SampleSize);
        
        % Starting a loop within, which is going to compare mean trait values between the two populations, for all replicates 
        Deltamu     = zeros(1, NoReps);
        for Rep     = 1:NoReps
            
            %% Calculate G1 and G2  
            G1 = Genotypes*(sig.*randn(NoM,1));
            G2 = Genotypes*(sig.*randn(NoM,1));

            %% Calculate V1 and V2 for all individuals
            V1 = mu1 + (1 - sqrt(H2))./sig .* randn(1, SampleSize);
            V2 = mu2 + (1 - sqrt(H2))./sig .* randn(1, SampleSize);
         
            %% Draw individuals at random from the total genotype pool, and store them in two population vectors, S1 and S2:
            for i = 1:SampleSize
                Y1      = randperm(2.^NoM);
                Y2      = randperm(2.^NoM);
                S1(1,i) = G1(Y1(1),1);
                S2(1,i) = G2(Y2(1),1);
            end

            %% Calculate the trait values for both populations, stored as Y1 and Y2
            Z1 = S1+V1;
            Z2 = S2+V2;

            %% Calculate deviations in observed trait differences as compared to projected trait differences, store result in vector Deltamu
            Deltamu(Rep) = abs((mean(Z1)-mu1)-(mean(Z2)-mu2));
        end

        %% Calculate mean and standard error of the observed trait deviations of all NoRep replicates (and then clear the vector)
        DeltamuTot(x_j  , sig_i) = mean(Deltamu);
        DeltamuError(x_j, sig_i) = std(Deltamu)./sqrt(NoReps);
        waitbar((sig_i-1)./X2+x_j./(X2.*X4),h)
    end
    waitbar(sig_i/X2,h)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% EXPORT RESULTS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%% Save the results as CSV files
csvwrite('deltamu.csv', DeltamuTot)
csvwrite('deltamuerror.csv', DeltamuError)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%