function [SaFre] = SA(Fre,n)
    % Generates a random sample of the frequency vector. 
    if numel(find(Fre==1)) == numel(Fre)
        SaFre = zeros(size(Fre));
        IX    = randN(n,sum(Fre));  % Generate random index vector
        SaFre(IX) = 1;
    else
        [EX,U]= expand(Fre);        % Create an expanded multiset.
        IX    = randN(n,sum(Fre));  % Generate random index vector
        sFre  = EX(IX);             % Create random sample
        SaFre = fold(U,Fre,sFre);   % fold back to frequency form 
    end
end

function [out,residual] = randN(n,N)
% Generate random n integers from the set of N integers; n<N

    t2=0; out=[]; N1=1:1:N;

    while t2<n
        % generate random indices
        t1 = floor(random('unif',1, numel(N1)+1-1e-5, 1, n-t2));

        tU = unique(t1);       % discard repeating elements
        out=[out N1(tU)];      % combine with previous elem.    
        t2 = numel(out);       % count current elements    
        N1=setdiff(N1,N1(tU)); % remove the found elements
    end
    residual=N1;
end

function [EXP, U] = expand(Fre)
% Expand multiset from a frequency form 
    EXP = zeros(1, sum(Fre));  % blank expanded multiset    
    U = 1:1:numel(Fre);        % unique identifiers 

    h = waitbar(0,'');
    binE=0;
    for i=1:numel(Fre)    
        binE = binE+Fre(i);     % bin end
        binS = binE-Fre(i)+1;   % bin start    
        EXP(binS:binE)=U(i); % fill each bin
        if mod( i, round(numel(Fre)/100) )==0
            waitbar(i/numel(Fre),h,'converting multiset');
        end
    end
    delete(h);
end

function [SaFre] = fold(U,Fre,sFre)
% collapse the expanded multised to frequency form.
    temp1 = sort(sFre, 'ascend')';
    [found, temp2, ~] = unique(temp1, 'first');

    temp3     = [temp2(2:end); numel(temp1)+1];
    Occur = temp3-temp2; 
    clear temp*

    [~,notfound]    = setdiff(U,found);
    SaFre           = Fre;
    SaFre(notfound) = 0;
    SaFre(found)    = Occur;
end
