function t_puzzle_lp ( grid_index )

%*****************************************************************************80
%
%% t_puzzle_lp() describes a T puzzle as a linear programming problem.
%
%  Discussion:
%
%    Matlab's sparse storage option is used for the matrix A.
%
%    We are given a grid, and a set of tiles.
%
%    We seek configurations (rotation+reflection+translation) of each
%    tile which remain in the grid.
%
%    Each such configuration is a possible component of a covering.
%
%    Note that tile #4 has reflective symmetry.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    13 June 2022
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer GRID_INDEX: indicates the grid to use.  1 <= GRID_INDEX <= 4
%    1 = t.
%    2 = arrow
%    3 = rhombus
%    4 = fatarrow.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 't_puzzle_lp():\n' );
  fprintf ( 1, '  Define the T puzzle and a set of tiles.\n' );
  fprintf ( 1, '  Set up the linear system A*x=b for the tiling problem.\n' );
  fprintf ( 1, '  Use sparse storage for the matrix A.\n' );
  fprintf ( 1, '  Write that information to an LP file.\n' );
%
%  Specify the number of tiles.
%
  tile_num = 4;
%
%  Define the region R.
%    grid_word: boundary word for R.
%    grid_p: (I,J) coordinates of the base point of triangle #1 of R.
%
  [ grid_word, grid_p ] = t_grid_word ( grid_index );
%
%  For each triangle and node in the grid:
%    grid_tri: the indices of the nodes.
%    grid_type: the triangle type, 1, 2, 3 or 4.
%    grid_xy: coordinates of the nodes.
%
  [ grid_tri, grid_type, grid_xy ] = object_elements ( grid_word, grid_p );

  grid_element_num = size ( grid_type, 1 );
%
%  Allocate arrays:
%    A for the EQN_NUM x VAR_NUM linear system.
%    b for the EQN_NUM right hand side.
%    EQN_NUM = number of elements to cover and tiles to use.
%    VAR_NUM is determined by generating all distinct tile configurations.
%
  eqn_num1 = grid_element_num;
  nz_num1 = 0;
  eqn_num2 = tile_num;
  nz_num2 = 0;
%
%  On first pass, determine var_num and nz_num.
%  On second pass, allocate A and assign values to A and b.
%
  for pass = 1 : 2

    if ( pass == 2 )
      fprintf ( 1, '\n' );
      fprintf ( 1, '  Tile configuration counts:\n' )
      fprintf ( 1, '\n' );
      for tile_index = 1 : tile_num
        fprintf ( 1, '  %d  %d\n', tile_index, config(tile_index) );
      end
    end

    if ( pass == 2 )
      A1 = sparse ( [], [], [], eqn_num1, var_num, nz_num1 );
      b1 = ones ( eqn_num1, 1 );
      A2 = sparse ( [], [], [], eqn_num2, var_num, nz_num2 );
      b2 = ones ( eqn_num2, 1 );
    end
%
%  Starting pass 2, we need var_num from pass 1 for the "sparse()" command.
%  Then we need to reset it to 0!
%
    config = zeros ( tile_num, 1 );
    var_num = 0;
%
%  Loop over tiles.
%
    for tile_index = 1 : tile_num
%
%  Consider all counterclockwise tile rotations that are multiples of 90 degrees:
%
      for r = 0 : 3
%
%  Consider a tile or its reflection.
%
        if ( tile_index == 4 )
          reflect_max = 0;
        else
          reflect_max = 1;
        end

        for reflect = 0 : reflect_max
%
%  Get a fresh copy of the word and base point.
%
          [ tile_word, tile_p ] = t_tile_word ( tile_index );
%
%  Perform R counterclockwise rotations of 90 degrees.
%
          [ tile_word, tile_p ] = t_rotate_word ( tile_word, tile_p, r );
%
%  Perform REFLECT reflections across X axis.
%
          if ( reflect == 1 )
            [ tile_word, tile_p ] = t_reflect_word ( tile_word, tile_p );
          end
%
%  Get the TYPE of the triangles that are the tile elements.
%    tile_tri: the indices of the nodes.
%    tile_type: the triangle type, 1, 2, 3 or 4.
%    tile_xy: coordinates of the nodes.
%
          [ tile_tri, tile_type, tile_xy ] = object_elements ( tile_word, ...
            tile_p );
%
%  Get tile (i,j) coordinates of the base node of tile triangle #1.
%
          tile_node1 = tile_tri(1,1);
%
%  Key is the type of element #1.
%  We can only try to place this tile configuration starting at 
%  region elements of the same type as KEY.
%
          key = tile_type(1);
%
%  Consider triangle J in the grid.
%
          for j = 1 : grid_element_num
%
%  Does it match the type of triangle #1 of S?
%
            if ( grid_type(j) == key )
%
%  Determine the translation that would map triangle #1 of S to triangle J of R.
%
              grid_node1 = grid_tri(j,1);
              tile_p = grid_xy(grid_node1,1:2) - tile_xy(tile_node1,1:2);
%
%  Does this configuration remain in the region R?
%
              [ inside, match ] = region_contains_tile ( grid_word, grid_p, ...
                tile_word, tile_p );
%
%  The configuration constitutes a new variable.
%  Update the equation of each triangle of R that is covered by the configuration.
%  Update the constraint equation for the tile.
%
              if ( inside )

                config(tile_index) = config(tile_index) + 1;
                var_num = var_num + 1;
 
                if ( pass == 2 )
                  for k = 1 : length ( match )
                    eqn_tri = match(k);
                    A1(eqn_tri,var_num) = 1;
                    nz_num1 = nz_num1 + 1;
                  end
                  A2(tile_index,var_num) = 1;
                  nz_num2 = nz_num2 + 1;
                end

              end % end configuration remains in region.

            end % end key match

          end  % end element search

        end  % end reflection check

      end  % end rotation check

    end  % end tile check

  end % pass 1 and 2
%    
%  Create the LP file describing this linear system.
%
  [ name1, name2 ] = t_grid_name ( grid_index );

  filename = sprintf ( '%s.lp', name2 );
  label = sprintf ( '\\ %s lp file', name1 );
  lp_write ( filename, label, A1, b1, A2, b2 );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  linear programming system computed.\n' );
  fprintf ( 1, '  matrix A1 has %d rows, %d columns, %d nonzeros\n', ...
    eqn_num1, var_num, nz_num1 );
  fprintf ( 1, '  matrix A2 has %d rows, %d columns, %d nonzeros\n', ...
    eqn_num2, var_num, nz_num2 );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  The lp file has been saved as "%s"\n', filename );
%
%  Terminate.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 't_puzzle_lp():\n' );
  fprintf ( 1, '  Normal end of execution.\n' );

  return
end     

