## Copyright (C) 2015 Lachlan Andrew <lachlanbis@gmail.com>
##
## This program is free software; you can redistribute it and/or modify it under
## the terms of the GNU General Public License as published by the Free Software
## Foundation; either version 3 of the License, or (at your option) any later
## version.
##
## This program is distributed in the hope that it will be useful, but WITHOUT
## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
## details.
##
## You should have received a copy of the GNU General Public License along with
## this program; if not, see <http://www.gnu.org/licenses/>.

## -*- texinfo -*-
## @deftypefn  {Function File} randsample (@var{n}, @var{k}, [@var{Replacement}, [@var{weights}]])
## @deftypefnx {Function File} randsample (@var{population}, @var{k}, [@var{Replacement}, [@var{weights}]])
## Randomly generate a column vector of @var{k} samples from 1:@var{n} or from
## @var{population}.
##
## If @var{Replacement} is specified, it is a logical value (true or false),
## specifying whether the sampling is performd with or without replacement.
##
## If @var{weights} is specified, then it is a non-negative vector of length
## @var{n}, or of the same length as @var{population}, specifying the
## probability with which each item is chosen.
##
## Unlike the Matlab(tm) function of the same name, this function supports
## weighted sampling without replacement, but does not support non-default
## random streams.
##
## @seealso{randi, randperm}
## @end deftypefn

function sample = randsample(n, k, replacement = false, weight = [])
  if (isscalar (n))
    population = [];
  else
    population = n;
    n = length (population(:));
  endif

  if (isempty (weight))
    if (n > intmax)
      error ('randsample: n is %d but should not exceed intmax (%d)', ...
             n, intmax);
    endif
  else
    if (n != length (weight))
      error ('randsample: Length of weight does not match n.  %d != %d', ...
             length (weight), n);
    endif
  endif

  ## Select random indices into the population, with(out) weights/replacement
  if (replacement)
    if (isempty (weight))
      if (n > intmax)
        error ('randsample: n is %d but should not exceed intmax (%d)', n, intmax);
      endif
      sample = randi (n, k, 1);     # standard uniform random with replacement
    else
      sample = weighted (k, weight)';# defined below
    endif
  else
    if (k > n)
      error ('randsample: cannot select %d from %d without replacement', k, n);
    endif
    if (isempty (weight))
      sample = randperm (n, k)';    # standard uniform random without replace
    else
      weight = weight(:);
      sample = sort (weighted (k, weight));
      while 1
        Idup = find (diff (sample)==0);
        if length (Idup) == 0
          break
        else
          weight(sample(Idup)) = 0;
          sample(Idup) = weighted (length (Idup), weight);
          sample = sort (sample);
        endif
      endwhile
      sample = sample';
    endif
  endif

  ## extract the population member
  if (iscell (population))
    sample = population{sample};
  elseif (!isempty (population))
    sample = population(sample);
  endif
endfunction

# Weighted sampling of k values from 1 to length(w), with replacement
function sample = weighted (k, w)
  [~, sample] = histc (rand (1,k), cumsum ([0; (w(:) ./ sum (w))]));
endfunction
