random_prune_idx_amd_fb.m

function prunemask=random_prune_idx_amd_fb(niter, LDmat, defvec, repeats)
% 02.07 added warning for improper use of LDmat

if ~exist( 'repeats', 'var' )
    repeats = 'maxout';
end

nsnp = size(LDmat, 1);
prunemask = false(nsnp, niter);

%tic
%LDmat(1:nsnp+1:end) = 0; % to prevent "self-pruning" % AMD: what is this??

fprintf('\n   Generating random prune indices... ')
niterstrlen = length( num2str( niter ) );
fprintf( '%*d/%*d', niterstrlen, 0, niterstrlen, niter )
for k=1:niter
    fprintf('\b');
    for n=1:niterstrlen, fprintf('\b\b'); end
    fprintf( '%*d/%*d', niterstrlen, k, niterstrlen, niter )
    tmplogp = unifrnd(0, 1, [nsnp 1]);

    tmplogp(~defvec) = NaN; % AMD: remove SNPs not defined in all traits 
    
    tmplogp_p = FastPrune(tmplogp, LDmat);

    prunemask(:, k) = isfinite(tmplogp_p);
end

if strcmp( repeats, 'none' )
    % 'none' allows *no* repetitions
    % find repeated snps ..
    repeated_snps = sum(prunemask, 2) > 1;
    % .. and create a prunemask slice out of them only
    prunemask_slice = prunemask(repeated_snps, :);
    % get the indexes of the slice ..
    [a, b] = ind2sub( [ sum(repeated_snps), niter ], find(prunemask_slice) );
    % .. and permute them so we don't pick always the same
    myorder = randperm(length(a));
    a = a(myorder);
    b = b(myorder);
    % pick the first unique indexes after permutation
    [~, myindex, ~] = unique(a);
    a = a(myindex);
    b = b(myindex);
    % re-initialize the slice ..
    prunemask_slice(:, :) = false;
    % .. set the unique indexes to 'true' ..
    prunemask_slice(sub2ind(size(prunemask_slice), a, b)) = true;
    % .. and re-insert the slice in the whole
    prunemask(repeated_snps, :) = prunemask_slice;
elseif strcmp( repeats, 'maxout' )
    % 'maxout' allows a maximum number of repetitions equal to the minimum number of repetitions
    % present in the set generated by random selection with replacement, as to avoid the peak of
    % singleton over-representation
    % find repeated snps ..
    repeated_snps = sum(prunemask, 2) == niter;
    repeated_cnt = sum(repeated_snps);
    % .. and create a prunemask slice out of them only
    prunemask_slice = prunemask(repeated_snps, :);
    % get the minimum iteration presence count in 'prunemask'
    mincnt = max(1, min(histc(sum(prunemask, 2), 1:niter)) - 1);
    % re-initialize the slice ..
    prunemask_slice(:, :) = false;
    % pick 'mincnt' random indexes from the slice
    prunemask_slice(randperm(repeated_cnt, mincnt), :) = true;
    % .. and re-insert the slice in the whole
    prunemask(repeated_snps, :) = prunemask_slice;
end

fprintf('\n');

%toc
end