function snakes_probability ( )

%*****************************************************************************80
%
%% SNAKES_PROBABILITY computes PDF and CDF for Snakes and Ladders.
%
%  Discussion:
%
%    This script determines the probability of finishing a one-player
%    game of chutes and ladders in exactly N moves (PDF) or in no more 
%    than N moves (CDF).
%
%    The game starts at square 0, and ends when square 100 is reached.
%
%    This program was originally named "chutes.m" and was written
%    by Desmond Higham and Nicholas Higham.
%
%  Modified:
%
%    15 July 2016
%
%  Author:
%
%    Desmond Higham, Nicholas Higham
%
%  Reference:
%
%    Desmond Higham, Nicholas Higham,
%    MATLAB Guide,
%    SIAM, 2005,
%    ISBN13: 9780898717891.
%
  timestamp ( )
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SNAKES_PROBABILITY\n' );
  fprintf ( 1, '  Determine the exact PDF and CDF for a version of the\n' );
  fprintf ( 1, '  game Snakes and Ladders.\n' );
%
%  N is the number of squares or "states" in the game.
%
  N = 100;
%
%  Define the transition matrix P.
%
  P = toeplitz(zeros(1,N+1),[0 ones(1,6) zeros(1,N-6)]);

  for k = N - 4 : N + 1
    P(k,k) = k - N + 5;
  end
  P = P / 6;
%
%  Define the chutes and ladders.
%  "+1" translates square (0 to 100) to state (1 to 101).
%
  top = [ 1  4  9 16 21 28 36 47 49 51 56 62 64 71 80  87 93 95 98] + 1;
  bot = [38 14 31  6 42 84 44 26 11 67 53 19 60 91 100 24 73 75 78] + 1;

  for k = 1 : length(top)
    r = top(k);
    s = bot(k);     % Chute or ladder from r to s.
    P(:,s) = P(:,s) + P(:,r);   % Add column r to column s.
  end
%
%  Remove starts of chutes and ladders.
%  Instead of 101 states, we now have 101-19 = 82 states.
%
  P(top,:) = [];
  P(:,top) = [];  
%
%  Figure 1: The transition matrix.
%
  figure ( 1 )
  spy ( P )
  title ( 'Snakes and Ladders Transition Matrix' );
  filename = 'snakes_transition.png';
  print ( '-dpng', filename );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  Transition matrix plot saved as "%s"\n', filename );
%
%  Determine the cumulative probabilities for game length.
%
  M = 200;
  cumprob = zeros ( M, 1 );
  cumprob(1) = P(1,end);
  v = P(1,:);
  for n = 2 : M,
    v = v * P;
    cumprob(n) = v(end);
  end
%
%  Figure 2: Game Length PDF.
%
  figure ( 2 )
  colormap ( [0.8,0.4,0.4] )
  bar ( diff ( [0;cumprob] ) )
  title ( 'Probability for Game Length', 'FontSize', 12, 'FontWeight', 'Bold' )
  grid on
  xlim ( [ 0 M ] )

  filename = 'snakes_gamelength_pdf.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Gamelength PDF plot saved as "%s"\n', filename );
%
%  Figure 3: Game Length CDF
%
  figure ( 3 )
  colormap ( [0.8,0.4,0.4] )
  bar ( cumprob )
  title ( 'Cumulative Probability for Game Length',...
    'FontSize', 12, 'FontWeight', 'Bold' )
  grid on
  xlim ( [0 M] )

  filename = 'snakes_gamelength_cdf.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Gamelength CDF plot saved as "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SNAKES_PROBABILITY:\n' );
  fprintf ( 1, '  Normal end of execution.\n' );
  fprintf ( 1, '\n' );
  timestamp ( )

  return
end

function timestamp ( )

%*****************************************************************************80
%
%% TIMESTAMP prints the current YMDHMS date as a timestamp.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    14 February 2003
%
%  Author:
%
%    John Burkardt
%
  t = now;
  c = datevec ( t );
  s = datestr ( c, 0 );
  fprintf ( 1, '%s\n', s );

  return
end
