function [ A, b ] = splitting_multihedral_matrix ( r_shape, p_num, p_shapes, d, C )

%*****************************************************************************80
%
%% splitting_multihedral_matrix() constructs a multihedral tiling linear system.
%
%  Discusssion:
%
%    This function constructs the linear system Ax = b for a pariomino 
%    splitting problem associated with a multihedral polyomino tiling problem.
%
%    In a multihedral tiling problem, the tiling is to be carried out using
%    copies of several polyomino shapes.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    16 February 2023
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer r_shape(rowsR,colsR): a rectangular matrix whose entries are 0 or 1.
%    The 1 values define the shape of the tiling region.
%
%    integer p_num: the number of polyomino shapes.
%
%    integer p_shapes(rows,cols,p_num): a rectangular matrix whose entries are 0 or 1.
%    The 1 values define the shape of the polyomino used for tiling.
%
%    integer C: the number of colors.
%
%  Output:
%  
%    integer A(C+p_num,C*p_num): the matrix of constraint coefficients.
%
%    integer b(C+p_num): the right hand sides of the constraint equations.
%

%
%  Allocate A and b.
%  We could set m = C + p_num - 1 because one of the constraints is redundant.
%
  n = C * p_num;
  m = C + p_num;
  b = zeros ( m, 1 );
  a = zeros ( m, n ); 
%
%  Construct the RHS vector b first. 
%
%  If r_shape does NOT represent a rectangle, use the first of the C 
%  variants associated with the r_shape.
%  If r_shape DOES represent a rectangle, use the variant of r_shape 
%  that has a '1' in the (1,1) position to be consistent
%  with a lemma in the paper.
%
  [ rowsR, colsR ] = size ( r_shape );
  
  if ( any ( r_shape(:) ==0 ) )
    for j = 1 : colsR 
      for i = 1 : rowsR
        if ( r_shape(i,j) == 1 )
          r_shape_pari(i,j) = mod ( i + j, C ) + 1;
        end
      end
    end
  else 
    for page = 1 : C
      variantsR(:,:,page) = zeros ( rowsR, colsR);
    end
    for j = 1 : colsR
      for i = 1 : rowsR
        if ( r_shape(i,j) == 1 )
          for page = 1 : C
            variantsR(i,j,page) = mod ( i + j + page - 1, C ) + 1;
          end
        end
      end
    end
%
%  Find which variant of r_shape has a '1' in the (1,1) position.
%
    for page = 1 : C
      ind = variantsR(1,1,page);
      if ( ind == 1 )
        keep = page;
        break
      end
    end
    r_shape_pari = variantsR(:,:,keep);
  end
  
    
  for colour = 1 : C
    b(colour) = length ( find ( r_shape_pari == colour ) );
  end
%
%  Add elements from constraints.
%
  count = 0;
  for i = C + 1 : m
    count = count + 1;
    b(i) = d(count);
  end
%
%  Construct the coefficient matrix a. First initialize the array that stores 
%  all n variants of the pariominoes. Each variant is packed into a
%  matrix with the same size as r_shape
%  
  for page = 1 : n
    variants(:,:,page) = zeros ( rowsR, colsR );
  end

  for count = 1 : p_num
    polyomino = p_shapes(:,:,count);  
    [ rows, cols ] = size ( polyomino );
    for j = 1 : cols
      for i = 1 : rows
        if ( polyomino(i,j) == 1 )
          for colour = 1 : C
            page = ( count - 1 ) * C + colour;
            variants(i,j,page) = mod ( i+j+page-1, C ) + 1;
          end
        end
      end
    end
  end
%
%  Count how many 1's, 2's, ..., C's there are in each variant and construct A.
%  
  for colour = 1 : C
    for page = 1 : n
      polyomino = variants(:,:,page);
      A(colour,page) = length ( find ( polyomino == colour ) );
    end
  end
%
%  Add elements due to constraints.
%
  for i = C + 1 : m
    low = ( i - C - 1) * C + 1;
    up = ( i - C ) * C;
    A(i,low:up) = 1;
  end
    
  return
end
            
 
