OrderedAlphabet.cpp

#include "DGtal/base/OrderedAlphabet.h"
#include "DGtal/arithmetic/ModuloComputer.h"

using namespace std;

// class OrderedAlphabet

// Standard services - public :

DGtal::OrderedAlphabet::~OrderedAlphabet()
{
  if ( myOrder != 0 )
    delete[ ] myOrder;
}

string
DGtal::OrderedAlphabet::orderedAlphabet() const
{
  char *tbl;
  tbl = (char *)malloc((myNb + 1)*sizeof(char));
  for ( unsigned int i = 0; i < myNb; ++i )
    {
      tbl[ myOrder[ i ] ] = myFirst + i;
    }
  tbl[ myNb ] = '\0';
  string s( tbl );
  free(tbl);
  return s;
}

void
DGtal::OrderedAlphabet::shiftLeft()
{
  unsigned int* p = myOrder;
  unsigned int* q = myOrder + myNb;
  while ( p != q )
    {
      unsigned int k = *p;
      *p = ( k == 0 ) ? myNb - 1 : k - 1; 
      ++p;
    }
}

void 
DGtal::OrderedAlphabet::shiftRight()
{
  unsigned int* p = myOrder;
  unsigned int* q = myOrder + myNb;
  while ( p != q )
    {
      unsigned int k = *p + 1;
      *p = ( k == myNb ) ? 0 : k; 
      ++p;
    }
}

void 
DGtal::OrderedAlphabet::reverse()
{
  unsigned int* p = myOrder;
  unsigned int* q = myOrder + myNb;
  while ( p != q )
    {
      *p = myNb - 1 - *p;
      ++p;
    }
}

void 
DGtal::OrderedAlphabet::reverseAround12()
{
  unsigned int* p = myOrder;
  unsigned int* q = myOrder + myNb;
  while ( p != q )
    {
      *p = ( myNb + 3 - *p ) % myNb;
      ++p;
    }
}




void 
DGtal::OrderedAlphabet::firstLyndonFactor
( size_t & len, size_t & nb,
  const std::string & w, 
  index_t s, index_t e ) const
{
  index_t i = s;
  index_t j = s+1;
  while ( ( j < e ) && ( lessOrEqual( w[ i ], w[ j ] ) ) )
    {
      if ( equal( w[ i ], w[ j ] ) )
  ++i;
      else 
  i = s;
      ++j;
    }
  len = (size_t) j - i;
  nb = ( (size_t) ( j - s ) ) / len;
}


void
DGtal::OrderedAlphabet::firstLyndonFactorMod
( size_t & len, size_t & nb,
  const std::string & w, 
  index_t s, index_t e ) const
{
  size_t modulo = (DGtal::OrderedAlphabet::size_t)w.size();
  ModuloComputer< Integer > mc( modulo );
  index_t i = s;
  index_t j = mc.next( s );
  while ( ( j != e ) && ( lessOrEqual( w[ i ], w[ j ] ) ) )
    {
      if ( equal( w[ i ], w[ j ] ) )
  mc.increment( i );
      else 
  i = s;
      mc.increment( j );
    }
  len = j >= i ? (size_t) ( j - i )
    : (size_t) ( j + modulo - i );
  nb = ( (size_t) ( ( j + modulo - s ) % modulo ) ) / len;
}


bool DGtal::OrderedAlphabet::duvalPP( size_t & len, size_t & nb, 
    const std::string & w, index_t s, index_t e) const
{
  ASSERT( ( order( w[ s ] ) == 1 ) || ( order( w[ s ] ) == 2 ) );
  index_t i = s;
  index_t j = s+1;
  index_t p = s;
  index_t q = s+1;
  while ( ( j < e ) && ( lessOrEqual( w[ i ], w[ j ] ) ) )
  {
    if ( equal( w[ i ], w[ j ] ) )
    {
      if ( j == q ) 
        q += (p-s+1);
      ++i;
    }
    else
    {
      if ( ( j != q ) || ( order ( w[ j ] ) != 2 ) )
      {
        len = j-s; nb = 0;
        return false;
      }
      index_t tmp = p; 
      p = q;
      q += q - tmp;
      i = s;
    }
    ++j;
  }
  len = (size_t) j - i;
  nb = ( (size_t) (j-s) ) / len;
  return true;
}

bool
DGtal::OrderedAlphabet::duvalPPtoDSS
( size_t & len, size_t & nb,
  unsigned int & n1,  unsigned int & n2, 
  unsigned int & lf1, unsigned int & lf2,
  const std::string & w, 
  index_t s, index_t e
  ) const
{
  ASSERT( ( order( w[ s ] ) == 1 ) || ( order( w[ s ] ) == 2 ) );
  index_t i = s;
  index_t j = s+1;
  index_t p = s;
  index_t q = s+1;
  unsigned int slope1 = (order( w[ i ] ) == 1) ? 1 : 0;
  unsigned int slope2 = (order( w[ i ] ) == 2) ? 1 : 0;
  lf1 = n1 = slope1;
  lf2 = n2 = slope2;
  nb = 1;
  //cerr << "input : " << w << endl;
  
  // Convex is not used so I comment it
  while ( ( j < e ) && ( lessOrEqual( w[ i ], w[ j ] ) ) ) {

    //cerr << "i=" << i << " j=" << j << " p=" << p << " q=" 
    //  << q << " nb=" << nb << " n1=" << n1 << " n2=" << n2 
    //  << " lf1=" << lf1 << " lf2=" << lf2 << endl;g

    //This 'if/else if' is added to compute the vector defined by
    //the Christoffel word, this is usefull in order to compute the
    //leaning points.
    if (order( w[ j ] ) == 1)
      ++slope1;
    else if (order( w[ j ] ) == 2)
      ++slope2;

    if ( equal( w[ i ], w[ j ] ) ) {
      if ( j == q ) {
        q += (p-s+1);

      //A repetition is read when j==s+(nb+1)*(p-s+1)-1
      } 
      if ( j == nb * (p-s+1) + p ) {
        ++nb;
      }
      ++i;
    } else {
      if ( ( j != q ) || ( order ( w[ j ] ) != 2 ) ) {
        //   convex = false;
        break;
      }
      index_t tmp = p; 
      p = q;
      q += q - tmp;
      i = s;

      // Extra information to compute the leaning points
      lf1 = lf1 + (nb-1)*n1;
      lf2 = lf2 + (nb-1)*n2;
      n1 = slope1;
      n2 = slope2;
      nb = 1;
    }
    ++j;
  }
  len = (size_t) j - i;
  //cerr << "Termine" << endl;
  //cerr << "i=" << i << " j=" << j << " p=" << p << " q=" 
  //  << q << " nb=" << nb << " n1=" << n1 << " n2=" << n2 
  //  << " lf1=" << lf1 << " lf2=" << lf2 << endl;

  if ( nb != (size_t) (j-s) / len)
    cout << "ASSERT(" << nb << " == (" << j << "-" << s << ") / "<<len << ")" << endl;
  ASSERT( nb == (size_t) (j-s) / len);
  //nb = (size_t) (j-s) / len;
  return true;
}


bool 
DGtal::OrderedAlphabet::duvalPPMod( size_t & len, size_t & nb,
               const std::string & w, 
               index_t s, index_t e ) const
{
  ASSERT( ( order( w[ s ] ) == 1 )
          || ( order( w[ s ] ) == 2 ) );
  size_t modulo = (DGtal::OrderedAlphabet::size_t)w.size();
  ModuloComputer< Integer > mc( modulo );
  ModuloComputer< Integer >::UnsignedInteger i = s;
  index_t j = mc.next( s );
  unsigned int p = 1;
  unsigned int q = 2;
  while ( ( j != e ) && ( lessOrEqual( w[ i ], w[ j ] ) ) )
    {
      // cerr << "i=" << i << " j=" << j << " p=" << p << " q=" << q << endl;
      if ( equal( w[ i ], w[ j ] ) )
  {
    if ( j == mc.cast( s + q - 1 ) )
      q += p;
    mc.increment( i );
  }
      else 
  {
    if ( ( j != mc.cast( s + q - 1 ) ) || ( order ( w[ j ] ) != 2 ) )
      {
        len = j; nb = 0;
        return false;
      }
    unsigned int tmp = p; 
    p = q;
    q += q - tmp;
    i = s;
  }
      mc.increment( j );
    }
  len = j >= i ? (size_t) ( j - i )
    : (size_t) ( j + modulo - i );
  nb = ( (size_t) ( ( j + modulo - s ) % modulo ) ) / len;
  return true;
}


// Interface - public :

void
DGtal::OrderedAlphabet::selfDisplay ( std::ostream & out ) const
{
    out << "[OrderedAlphabet]";
    out << " " << orderedAlphabet() << endl;
}

bool
DGtal::OrderedAlphabet::isValid() const
{
    return true;
}



// ----------------------- MLP services ------------------------------

DGtal::OrderedAlphabet::size_t
DGtal::OrderedAlphabet::nextEdge( size_t & nb_a1,
            size_t & nb_a2,
            std::string & w,
            index_t & s,
            bool & cvx )
{
  ModuloComputer< Integer > mc( (const unsigned int)w.size() );
  size_t l;
  size_t len;
  size_t nb;
  bool inC = duvalPPMod( len, nb, w, s, s );
  if ( ! inC ) 
    // case : change of convexity
    {
//       cerr << "Convexity change" <<  orderedAlphabet() ;
//       cerr << " (" << w[ len ] << "==" << letter( 2 ) << ")";
      // JOL : temporary change of letter w[ s ]
      char tmp = w[ s ];
      index_t tmp_s = s;
      w[ s ] = letter( 2 ); // a3
      this->reverseAround12();
      cvx = ! cvx;
//       cerr << " => " << orderedAlphabet() << endl;
      l = nextEdge( nb_a1, nb_a2, w, s, cvx );
      // JOL : former letter is put back in string.
      w[ tmp_s ] = tmp;
    }
  else if ( ( len == 1 ) && ( order( w[ s ] ) == 1 ) ) 
    // case u=a1 => Quadrant change
    {
//       cerr << "Quadrant change " << orderedAlphabet() ;
      this->shiftRight();
//       cerr << " => " << orderedAlphabet() << endl;
      s = mc.cast( s + nb );
      nb_a1 = 0;
      nb_a2 = nb - 1;
      l = nb;
    }
  else 
    { // standard (convex) case.
//       cerr << "standard (convex) case " << orderedAlphabet() << endl;
      l = len * nb;
      char a2 = letter( 2 ); 
      nb_a1 = len;
      nb_a2 = 0;
      index_t ss = s;
      s = mc.cast( s + l );
      while ( len != 0 )
  {
    if ( w[ ss ] == a2 ) ++nb_a2;
    mc.increment( ss );
    --len;
  }
      nb_a1 -= nb_a2;
      nb_a1 *= nb;
      nb_a2 *= nb;
    }
  return l;
}
            

// Internals - private :

//                                                                           //