function [obsDiff, pVal] = nestetPermutationTest(sample1, sample2, ...
                                                 groupsCell, reorderedFlyIndsCell, ...
                                                 permutations, sidedness, plotresult)

allobservations = [sample1, sample2];
observeddifference = nanmean(groupsCell{1}) - ...
                     nanmean(groupsCell{2});
effectsize = observeddifference / ...
             nanmean([nanstd(groupsCell{1}), nanstd(groupsCell{2})]);

exact = 0;
showprogress = permutations / 10;

for n = 1:permutations
    if showprogress && mod(n,showprogress) == 0, display(sprintf('Permutation %d of %d', n, permutations)); end

    % selecting either next combination, or random permutation
    if exact, permutation = [allcombinations(n,:), setdiff(1:numel(allobservations), allcombinations(n,:))];
    else, permutation = randperm(length(allobservations)); end

    % dividing into two samples
    randomSample1 = cell(1, numel(sample1));
    randomSample2 = cell(1, numel(sample2));
    jCount = 1;
    for iElem = permutation
        % First fill sample1, then check if the sample coes from first or
        % second genotype.
        if jCount <= length(sample1)
            if iElem > length(sample1)
                sampleInds = reorderedFlyIndsCell{2} == iElem - length(sample1);
                randomSample1{jCount} = groupsCell{2}(sampleInds);
            else
                sampleInds = reorderedFlyIndsCell{1} == iElem;
                randomSample1{jCount} = groupsCell{1}(sampleInds);
            end
        else
            if iElem > length(sample1)
                sampleInds = reorderedFlyIndsCell{2} == iElem - length(sample1);
                randomSample2{jCount - numel(sample1)} = groupsCell{2}(sampleInds);
            else
                sampleInds = reorderedFlyIndsCell{1} == iElem;
                randomSample2{jCount - numel(sample1)} = groupsCell{1}(sampleInds);
            end
        end
        jCount = jCount + 1;
    end

    % saving differences between the two samples
    randomdifferences(n) = nanmean(cat(2, randomSample1{:})) - ...
                           nanmean(cat(2, randomSample2{:}));
end


% getting probability of finding observed difference from random permutations
if strcmp(sidedness, 'both')
    p = (length(find(abs(randomdifferences) > abs(observeddifference)))+1) / (permutations+1);
elseif strcmp(sidedness, 'smaller')
    p = (length(find(randomdifferences < observeddifference))+1) / (permutations+1);
elseif strcmp(sidedness, 'larger')
    p = (length(find(randomdifferences > observeddifference))+1) / (permutations+1);
end

% plotting result
if plotresult
    figure;
    hist(randomdifferences, 20);
    hold on;
    xlabel('Random differences');
    ylabel('Count')
    od = plot(observeddifference, 0, '*r', 'DisplayName', sprintf('Observed difference.\nEffect size: %.2f,\np = %f', effectsize, p));
    legend(od);
end

obsDiff = observeddifference;
pVal = p;

end