# ifndef CPPAD_LOCAL_COLOR_SYMMETRIC_HPP
# define CPPAD_LOCAL_COLOR_SYMMETRIC_HPP
# include <cppad/configure.hpp>
# include <cppad/local/cppad_colpack.hpp>

// SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
// SPDX-FileCopyrightText: Bradley M. Bell <bradbell@seanet.com>
// SPDX-FileContributor: 2003-22 Bradley M. Bell
// ----------------------------------------------------------------------------

namespace CppAD { namespace local { // BEGIN_CPPAD_LOCAL_NAMESPACE
/*!
\file color_symmetric.hpp
Coloring algorithm for a symmetric sparse matrix.
*/
// --------------------------------------------------------------------------
/*!
CppAD algorithm for determining which rows of a symmetric sparse matrix can be
computed together.

\tparam SizeVector
is a simple vector class with elements of type size_t.

\tparam SetVector
is a vector_of_sets class.

\param pattern [in]
Is a representation of the sparsity pattern for the matrix.

\param row [in/out]
is a vector specifying which row indices to compute.

\param col [in/out]
is a vector, with the same size as row,
that specifies which column indices to compute.
\n
\n
Input:
For each  valid index k, the index pair
<code>(row[k], col[k])</code> must be present in the sparsity pattern.
It may be that some entries in the sparsity pattern do not need to be computed;
i.e, do not appear in the set of
<code>(row[k], col[k])</code> entries.
\n
\n
Output:
On output, some of row and column indices may have been swapped
\code
   std::swap( row[k], col[k] )
\endcode
So the the the color for row[k] can be used to compute entry
(row[k], col[k]).

\param color [out]
is a vector with size m.
The input value of its elements does not matter.
Upon return, it is a coloring for the rows of the sparse matrix.
Note that if color[i] == m, then there is no index k for which
row[k] == i (for the return value of row).
\n
\n
Fix any (i, j) in the sparsity pattern.
Suppose that there is a row index i1 with
i1 != i, color[i1] == color[i] and (i1, j) is in the sparsity pattern.
If follows that for all j1 with
j1 != j and color[j1] == color[j],
(j1, i ) is not in the sparsity pattern.
\n
\n
This routine tries to minimize, with respect to the choice of colors,
the maximum, with respect to k, of <code>color[ row[k] ]</code>.
*/
template <class SetVector>
void color_symmetric_cppad(
   const SetVector&        pattern   ,
   CppAD::vector<size_t>&  row       ,
   CppAD::vector<size_t>&  col       ,
   CppAD::vector<size_t>&  color     )
{
   size_t K = row.size();
   size_t m = pattern.n_set();
   CPPAD_ASSERT_UNKNOWN( m == pattern.end() );
   CPPAD_ASSERT_UNKNOWN( color.size() == m );
   CPPAD_ASSERT_UNKNOWN( col.size()   == K );

   // row, column pairs that appear in ( row[k], col[k] )
   CppAD::vector< std::set<size_t> > pair_needed(m);
   std::set<size_t>::iterator itr1, itr2;
   for(size_t k1 = 0;  k1 < K; k1++)
   {  CPPAD_ASSERT_UNKNOWN( pattern.is_element(row[k1], col[k1]) );
      pair_needed[ row[k1] ].insert( col[k1] );
      pair_needed[ col[k1] ].insert( row[k1] );
   }

   // order the rows decending by number of pairs needed
   CppAD::vector<size_t> key(m), order2row(m);
   for(size_t i1 = 0; i1 < m; i1++)
   {  CPPAD_ASSERT_UNKNOWN( pair_needed[i1].size() <= m );
      key[i1] = m - pair_needed[i1].size();
   }
   CppAD::index_sort(key, order2row);

   // mapping from order index to row index
   CppAD::vector<size_t> row2order(m);
   for(size_t o1 = 0; o1 < m; o1++)
      row2order[ order2row[o1] ] = o1;

   // initial coloring
   color.resize(m);
   size_t c1 = 0;
   for(size_t o1 = 0; o1 < m; o1++)
   {  size_t i1 = order2row[o1];
      if( pair_needed[i1].empty() )
         color[i1] = m;
      else
         color[i1] = c1++;
   }

   // which colors are forbidden for this row
   CppAD::vector<bool> forbidden(m);

   // must start with row zero so that we remove results computed for it
   for(size_t o1 = 0; o1 < m; o1++) // for each row that appears (in order)
   if( color[ order2row[o1] ] < m )
   {  size_t i1 = order2row[o1];
      c1 = color[i1];

      // initial all colors as ok for this row
      // (value of forbidden for c > c1 does not matter)
      for(size_t c2 = 0; c2 <= c1; c2++)
         forbidden[c2] = false;

      // -----------------------------------------------------
      // Forbid grouping with rows that would destroy results that are
      // needed for this row.
      itr1 = pair_needed[i1].begin();
      while( itr1 != pair_needed[i1].end() )
      {  // entry (i1, j1) is needed for this row
         size_t j1 = *itr1;

         // Forbid rows i2 != i1 that have non-zero sparsity at (i2, j1).
         // Note that this is the same as non-zero sparsity at (j1, i2)
         typename SetVector::const_iterator pattern_itr(pattern, j1);
         size_t i2 = *pattern_itr;
         while( i2 != pattern.end() )
         {  size_t c2 = color[i2];
            if( c2 < c1 )
               forbidden[c2] = true;
            i2 = *(++pattern_itr);
         }
         itr1++;
      }
      // -----------------------------------------------------
      // Forbid grouping with rows that this row would destroy results for
      for(size_t o2 = 0; o2 < o1; o2++)
      {  size_t i2 = order2row[o2];
         size_t c2 = color[i2];
         itr2 = pair_needed[i2].begin();
         while( itr2 != pair_needed[i2].end() )
         {  size_t j2 = *itr2;
            // row i2 needs pair (i2, j2).
            // Forbid grouping with i1 if (i1, j2) has non-zero sparsity
            if( pattern.is_element(i1, j2) )
               forbidden[c2] = true;
            itr2++;
         }
      }

      // pick the color with smallest index
      size_t c2 = 0;
      while( forbidden[c2] )
      {  c2++;
         CPPAD_ASSERT_UNKNOWN( c2 <= c1 );
      }
      color[i1] = c2;

      // no longer need results that are computed by this row
      itr1 = pair_needed[i1].begin();
      while( itr1 != pair_needed[i1].end() )
      {  size_t j1 = *itr1;
         if( row2order[j1] > o1 )
         {  itr2 = pair_needed[j1].find(i1);
            if( itr2 != pair_needed[j1].end() )
            {  pair_needed[j1].erase(itr2);
               if( pair_needed[j1].empty() )
                  color[j1] = m;
            }
         }
         itr1++;
      }
   }

   // determine which sparsity entries need to be reflected
   for(size_t k1 = 0; k1 < row.size(); k1++)
   {  size_t i1   = row[k1];
      size_t j1   = col[k1];
      itr1 = pair_needed[i1].find(j1);
      if( itr1 == pair_needed[i1].end() )
      {  row[k1] = j1;
         col[k1] = i1;
# ifndef NDEBUG
         itr1 = pair_needed[j1].find(i1);
         CPPAD_ASSERT_UNKNOWN( itr1 != pair_needed[j1].end() );
# endif
      }
   }
   return;
}

// --------------------------------------------------------------------------
/*!
Colpack algorithm for determining which rows of a symmetric sparse matrix
can be computed together.

\copydetails CppAD::local::color_symmetric_cppad
*/
template <class SetVector>
void color_symmetric_colpack(
   const SetVector&        pattern   ,
   CppAD::vector<size_t>&  row       ,
   CppAD::vector<size_t>&  col       ,
   CppAD::vector<size_t>&  color     )
{
# if ! CPPAD_HAS_COLPACK
   CPPAD_ASSERT_UNKNOWN(false);
   return;
# else
   size_t i, j, k;
   size_t m = pattern.n_set();
   CPPAD_ASSERT_UNKNOWN( m == pattern.end() );
   CPPAD_ASSERT_UNKNOWN( row.size() == col.size() );

   // Determine number of non-zero entries in each row
   CppAD::vector<size_t> n_nonzero(m);
   size_t n_nonzero_total = 0;
   for(i = 0; i < m; i++)
   {  n_nonzero[i] = 0;
      typename SetVector::const_iterator pattern_itr(pattern, i);
      j = *pattern_itr;
      while( j != pattern.end() )
      {  n_nonzero[i]++;
         j = *(++pattern_itr);
      }
      n_nonzero_total += n_nonzero[i];
   }

   // Allocate memory and fill in Adolc sparsity pattern
   CppAD::vector<unsigned int*> adolc_pattern(m);
   CppAD::vector<unsigned int>  adolc_memory(m + n_nonzero_total);
   size_t i_memory = 0;
   for(i = 0; i < m; i++)
   {  adolc_pattern[i]    = adolc_memory.data() + i_memory;
      CPPAD_ASSERT_KNOWN(
         std::numeric_limits<unsigned int>::max() >= n_nonzero[i],
         "Matrix is too large for colpack"
      );
      adolc_pattern[i][0] = static_cast<unsigned int>( n_nonzero[i] );
      typename SetVector::const_iterator pattern_itr(pattern, i);
      j = *pattern_itr;
      k = 1;
      while(j != pattern.end() )
      {
         CPPAD_ASSERT_KNOWN(
            std::numeric_limits<unsigned int>::max() >= j,
            "Matrix is too large for colpack"
         );
         adolc_pattern[i][k++] = static_cast<unsigned int>( j );
         j = *(++pattern_itr);
      }
      CPPAD_ASSERT_UNKNOWN( k == 1 + n_nonzero[i] );
      i_memory += k;
   }
   CPPAD_ASSERT_UNKNOWN( i_memory == m + n_nonzero_total );

   // Must use an external routine for this part of the calculation because
   // ColPack/ColPackHeaders.h has as 'using namespace std' at global level.
   cppad_colpack_symmetric(color, m, adolc_pattern);

   // determine which sparsity entries need to be reflected
   for(size_t k1 = 0; k1 < row.size(); k1++)
   {  size_t i1 = row[k1];
      size_t j1 = col[k1];
      bool reflect = false;
      for(size_t i2 = 0; i2 < m; i2++)
      if( (i1 != i2) & (color[i1]==color[i2]) )
      {  for(size_t k2 = 1; k2 <= adolc_pattern[i2][0]; k2++)
         {  size_t j2 = adolc_pattern[i2][k2];
            reflect |= (j1 == j2);
         }
      }
      if( reflect )
      {  row[k1] = j1;
         col[k1] = i1;
      }
   }
   return;
# endif // CPPAD_HAS_COLPACK
}

} } // END_CPPAD_LOCAL_NAMESPACE

# endif
