function [ y, a ] = walker_build ( n, x )

%*****************************************************************************80
%
%% WALKER_BUILD sets up the data for a Walker sampler.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    21 February 2016
%
%  Author:
%
%    Original C version by Warren Smith.
%    This MATLAB version by John Burkardt.
%
%  Reference:
%
%    Warren Smith,
%    How to sample from a probability distribution,
%    April 18, 2002.
%
%  Parameters:
%
%    Input, integer N, indicates the size of X.
%
%    Input, real X(0:N+1), contains in X(1) through X(N) the
%    probabilities of outcomes 1 through N.  
%
%    Output, real Y(0:N+1), the Walker threshold vector.
%
%    Output, integer A(0:N+1), the Walker index vector.
%

%
%  Initialize A.
%
  a = zeros ( n + 2, 1 );
  for i = 0 : n + 1
    a(i+1) = i;
  end
%
%  Initialize B to the "stay here" value, and set sentinel values at the ends.
%
  b = zeros ( n + 2, 1 );
  for i = 0 : n + 1
    b(i+1) = i;
  end
%
%  Copy Y from X.
%  Scale the probability vector and set sentinel values at the ends.
%
  y = zeros ( n + 2, 1 );
  y(0+1) = 0.0;
  y(2:n+1) = x(2:n+1) * n;
  y(n+1+1) = 2.0;

  i = 0;
  j = n + 1;

  while ( 1 )
%
%  Find i so Y(B(i)) needs more.
%
    while ( 1 )
      i = i + 1;
      if ( 1.0 <= y(b(i+1)+1) )
        break;
      end
    end
%
%  Find j so Y(B(j)) wants less.
%
    while ( 1 )
      j = j - 1;
      if ( y(b(j+1)+1) < 1.0 )
        break;
      end
    end

    if ( j <= i )
      break;
    end
%
%  Swap B(i) and B(j).
%
    k      = b(i+1);
    b(i+1) = b(j+1);
    b(j+1) = k;

  end

  i = j;
  j = j + 1;

  while ( 0 < i )
%
%  Find J such that Y(B(j)) needs more.
%
    while ( y(b(j+1)+1) <= 1.0 )
      j = j + 1;
    end
%
%  Meanwhile, Y(B(i)) wants less.
%
    if ( n < j )
      break;
    end
%
%  B(i) will donate to B(j) to fix up.
%
    y(b(j+1)+1) = y(b(j+1)+1) - ( 1.0 - y(b(i+1)+1) );
    a(b(i+1)+1) = b(j+1);
% 
%  Y(B(j)) now wants less so readjust ordering.
%
    if ( y(b(j+1)+1) < 1.0 )

      k      = b(i+1);
      b(i+1) = b(j+1);
      b(j+1) = k;
      j = j + 1;

    else

      i = i - 1;

    end

  end

  return
end

