function polyomino_monohedral_tiling_print ( r, p, x, label )

%*****************************************************************************80
%
%% POLYOMINO_MONOHEDRAL_TILING_PRINT prints monohedral polyomino tilings.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    08 June 2018
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer R(MR,NR), a matrix of 0's and 1's representing the 
%    region.
%
%    Input, integer P(MP,NP), a matrix of 0's and 1's representing the 
%    polyomino.
%
%    Input, integer X(NVAR,X_NUM), a binary vector describing one or more
%    solutions of the tiling problem.  Every possible placement of a copy of 
%    P has been assigned an index.  If X(I,J) is 1, then the I-th placement 
%    is included in solution J.
%
%    Input, string LABEL, a label for the printout, such as 
%    'Solution of the Reid Example'.
%

%
%  Get the sizes of R, P, and X.
%
  [ mp, np ] = size ( p );
  [ mr, nr ] = size ( r );
  [ mx, nx ] = size ( x );
%
%  Construct all possible variants (rotations and reflections) of P.
%
  [ mn_num, mn_v, nm_num, nm_v ] = polyomino_monohedral_variants ( mp, np, p );
%
%  Process each solution.
%
  for x_index = 1 : nx
%
%  The array S is our representation of the tiling.
%  S starts as a zeroed out version of R.
%
    s = zeros ( mr, nr );
%
%  Start the variable count at 0.
%
    var = 0;
%
%  Handle the variants which have shape (MP,NP) like P.
%  Store the actual representation of the variant in Q.
%
    for k = 1 : mn_num

      mq = mp;
      nq = np;
      q = mn_v(1:mq,1:nq,k);
%
%  Determine the number of ways, and the corresponding offsets, by which
%  Q can be placed on R.
%
      number = polyomino_embed_number ( mr, nr, r, mq, nq, q );

      list = polyomino_embed_list ( mr, nr, r, mq, nq, q, number );
%
%  For each placement VAR, check whether it is part of the solution.
%  If so, we can mark the corresponding entries of S with VAR.
%
      for l = 1 : number
        var = var + 1;
        if ( x(var,x_index) == 1 )
          ioff = list(l,1);
          joff = list(l,2);
          s(1+ioff:mq+ioff,1+joff:nq+joff) = ...
            s(1+ioff:mq+ioff,1+joff:nq+joff) + var * q(1:mq,1:nq);
        end
      end

    end
%
%  Similar procedure for the variants of P that have shape (N,M).
%
    for k = 1 : nm_num

      mq = np;
      nq = mp;
      q = nm_v(1:mq,1:nq,k);

      number = polyomino_embed_number ( mr, nr, r, mq, nq, q );
 
      list = polyomino_embed_list ( mr, nr, r, mq, nq, q, number );

      for l = 1 : number
        var = var + 1;
        if ( x(var,x_index) == 1 )
          ioff = list(l,1);
          joff = list(l,2);
          s(1+ioff:mq+ioff,1+joff:nq+joff) = ...
            s(1+ioff:mq+ioff,1+joff:nq+joff) + var * q(1:mq,1:nq);
        end
      end

    end
%
%  How big is the largest integer we have to print?
%  Create appropriate formats to integers (and blanks) line up.
%
    smax = max ( max ( s ) );
    smax = max ( smax, 1 );

    ssize = floor ( log10 ( smax ) ) + 1;
    blanks = '';
    for i = 0 : ssize
      blanks = [ blanks, ' ' ];
    end
    ssize_string = num2str ( ssize );
    format1 = [ ' ', '%', ssize_string, 'd' ];
%
%  The array S should contain integer values only.
%  However, it was computed with real arithmetic.
%  To ensure that the printout below doesn't accidentally start using
%  scientific notation, convert S to integer format.
%
    s_int = int32 ( s );
%
%  Print the tiling.
%
    fprintf ( 1, '\n' );
    fprintf ( 1, '  Solution %d\n', x_index );
    fprintf ( 1, '%s\n', label );
    fprintf ( 1, '\n' );

    for i = 1 : mr
      fprintf ( 1, ' ' );
      for j = 1 : nr
        if ( s_int(i,j) == 0 )
          fprintf ( 1, '%s', blanks );
        else
          fprintf ( 1, format1, s_int(i,j) );
        end
      end
      fprintf ( 1, '\n' );
    end

  end

  return
end
