MPolynomialReader.h

#pragma once

#if defined(MPolynomialReader_RECURSES)
#error Recursive header files inclusion detected in MPolynomialReader.h
#else // defined(MPolynomialReader_RECURSES)

#define MPolynomialReader_RECURSES

#if !defined MPolynomialReader_h

#define MPolynomialReader_h

// Inclusions
#include <iostream>
#include "DGtal/base/Common.h"
#include "DGtal/math/MPolynomial.h"
#include <boost/spirit/include/qi.hpp>
#include <boost/spirit/include/phoenix_core.hpp>
#include <boost/spirit/include/phoenix_operator.hpp>
#include <boost/spirit/include/phoenix_fusion.hpp>
#include <boost/spirit/include/phoenix_stl.hpp>
#include <boost/fusion/include/adapt_struct.hpp>
#include <boost/variant/recursive_variant.hpp>

namespace DGtal
{
  namespace detail {

    struct power_node {
      int k;
      int e;
    };

    struct monomial_node {
      double coef;
      std::vector<power_node> powers;
    };


    struct top_node;

    typedef 
    boost::variant< boost::recursive_wrapper<top_node>, monomial_node  >
    expr_node;

    struct top_node {
      std::vector<char> ops; 
      std::vector<expr_node> expressions; 
      int exp; 
    };

  }

}

BOOST_FUSION_ADAPT_STRUCT(
                          DGtal::detail::power_node,
                          (int, k)
                          (int, e)
)


BOOST_FUSION_ADAPT_STRUCT(
                          DGtal::detail::monomial_node,
                          (double, coef)
                          (std::vector<DGtal::detail::power_node>, powers)
)

BOOST_FUSION_ADAPT_STRUCT(
    DGtal::detail::top_node,
    (std::vector<char>, ops)
    (std::vector<DGtal::detail::expr_node>, expressions)
    (int, exp)
)

namespace DGtal
{
  namespace qi = boost::spirit::qi;
  namespace ascii = boost::spirit::ascii;
  namespace phoenix = boost::phoenix;

  template <typename Iterator>
  struct MPolynomialGrammar 
    : qi::grammar<Iterator, detail::top_node(), ascii::space_type>
  {
    MPolynomialGrammar()
      : MPolynomialGrammar::base_type(top)
    {
      using qi::eps;
      using qi::lit;
      using qi::int_;
      using qi::_val;
      using qi::_1;
      using qi::double_;
      using phoenix::at_c;
      using phoenix::push_back;

      top = // An expression is an additive or subtractive expression
        mulexpr [push_back(at_c<1>(_val), _1)] 
        >> *( ( lit('+') [ push_back(at_c<0>(_val), '+') ] 
                >> mulexpr [push_back(at_c<1>(_val), _1)]  )
              | ( ( lit('-') [ push_back(at_c<0>(_val), '-') ] 
                    >> mulexpr [push_back(at_c<1>(_val), _1)]  ) ) );
      
      mulexpr = // Each mul expression may be a product of sub-expressions
        subexpr [push_back(at_c<1>(_val), _1)] 
        >> *( ( lit('*') [ push_back(at_c<0>(_val), '*') ]
                >> mulexpr [push_back(at_c<1>(_val), _1)] ) );

      subexpr =  // a sub-expression is a monomial or some ( ) or some ( )^k
        monomial 
        | expexpr;

      expexpr = // ( expr ) or ( expr )^k
        lit('(') [ push_back(at_c<0>(_val), '^') ] 
        >> top [ push_back(at_c<1>(_val), _1) ]
        >> lit(')') 
        >> ( ( lit('^') >> int_ [ at_c<2>(_val) = _1 ] )
             | eps [ at_c<2>(_val) = 1 ]  ) ;
        
      monomial =   // coef.power(s)*, or power(s)+
        ( double_ [at_c<0>(_val) = _1] 
          >> *( genvariable [push_back(at_c<1>(_val), _1)] ) )
        | +( genvariable [push_back(at_c<1>(_val), _1)] 
             >> eps [at_c<0>(_val) = 1] );

      genvariable = // may be some X_k^m or x^m, y^m ...
        litvariable | variable;
      variable = // X_0 X_1 ... X_0^4 X_1^2 X_3^0 ...
        lit('X') >> lit('_') 
                 >> int_ [at_c<0>(_val) = _1]
                 >> ( ( lit('^') >> int_ [at_c<1>(_val) = _1] ) // X_k^e
                      | eps [at_c<1>(_val) = 1] // X_k
                      );
      litvariable = // x y z t x^4 y^5 z^2 ...
        ( lit('x') [at_c<0>(_val) = 0]
          | lit('y') [at_c<0>(_val) = 1] 
          | lit('z') [at_c<0>(_val) = 2] 
          | lit('t') [at_c<0>(_val) = 3] )
        >> ( ( lit('^') >> int_ [at_c<1>(_val) = _1] ) // x^3 z^4
             | eps [at_c<1>(_val) = 1] // x y z
             );

    }


    qi::rule<Iterator, detail::top_node(), ascii::space_type> top;
    qi::rule<Iterator, detail::top_node(), ascii::space_type> mulexpr;
    qi::rule<Iterator, detail::expr_node(), ascii::space_type> subexpr;
    qi::rule<Iterator, detail::top_node(), ascii::space_type> expexpr;
    qi::rule<Iterator, detail::monomial_node(), ascii::space_type> monomial;
    qi::rule<Iterator, detail::power_node(), ascii::space_type> variable;
    qi::rule<Iterator, detail::power_node(), ascii::space_type> genvariable;
    qi::rule<Iterator, detail::power_node(), ascii::space_type> litvariable;

  };
}


namespace DGtal
{
  // template class MPolynomialReader
  template <int n, typename TRing, 
            typename TAlloc = std::allocator<TRing>,
            typename TIterator = std::string::const_iterator>
  class MPolynomialReader
  {
  public:
    typedef TRing Ring;
    typedef TIterator Iterator;
    typedef TAlloc Alloc;
    typedef MPolynomial<n, Ring, Alloc > Polynomial;
    typedef MPolynomialGrammar<Iterator> Grammar;
    
    Grammar gpolynomial;
    
    MPolynomialReader() {}

    Iterator read( Polynomial & p, Iterator begin, Iterator end )
    {
      using qi::phrase_parse;
      using ascii::space;
      detail::top_node m;
      bool r = phrase_parse( begin, end, gpolynomial, space, m );
      if (r) p = make( m );
      return r ? begin : end;
    }



    // ----------------------- Interface --------------------------------------
  public:

    void selfDisplay ( std::ostream & out ) const;
    
    bool isValid() const;

     

    // ------------------------------ internals ------------------------------
  private:

    Polynomial make( const detail::power_node & pnode )
    {
      return Xe_k<n, Ring>( pnode.k, pnode.e );
    }

    Polynomial make( const detail::monomial_node & mnode )
    {
      Polynomial m;
      if ( mnode.powers.size() != 0 )
        {
          m = make( mnode.powers[ 0 ] );
          for ( unsigned int i = 1; i < mnode.powers.size(); ++i )
            m *= make( mnode.powers[ i ] );
        }
      else
        m = 1;
      return ( (Ring) mnode.coef ) * m;
    }

    struct ExprNodeMaker : boost::static_visitor<> {
      Polynomial myP;
      MPolynomialReader & myPR;
      ExprNodeMaker( MPolynomialReader & reader )
        : myPR( reader )
      {}
      void operator()( const detail::monomial_node & mnode)
      {
        myP = myPR.make( mnode );
      }
      void operator()( const detail::top_node & topnode)
      {
        myP = myPR.make( topnode );
      }
    };

    Polynomial make( const detail::top_node & topnode )
    {
      ASSERT( ! topnode.expressions.empty() );
      Polynomial p;
      ExprNodeMaker emaker( *this );
      if ( topnode.ops.empty() )
        {
          // Node is identity. Nothing special to do. 
          boost::apply_visitor( emaker, topnode.expressions[ 0 ] );
          p = emaker.myP;
        }
      else if ( topnode.ops[ 0 ] == '^' )
        {
          // Node is some exponent ( ... )^k. ^0 is admissible.
          boost::apply_visitor( emaker, topnode.expressions[ 0 ] );
          p = (Ring) 1; 
          for ( unsigned int i = 1; i <= (unsigned int)topnode.exp; ++i )
            p *= emaker.myP;
        }
      else
        {
          // Node is expr1 (*|+|-) expr2 (*|+|-) expr3 ... 
          // NB: either ops are in {+,-} or in {*} only.
          boost::apply_visitor( emaker, topnode.expressions[ 0 ] );
          p = emaker.myP;
          for ( unsigned int i = 0; i < (unsigned int)topnode.ops.size(); ++i )
            {
              boost::apply_visitor( emaker, topnode.expressions[ i+1 ] );
              switch ( topnode.ops[ i ] ) {
              case '+': p += emaker.myP; break;
              case '-': p -= emaker.myP; break;
              case '*': p *= emaker.myP; break;
              default: std::cerr << "[UNKNOWN-node]" << topnode.ops[ i ] << std::endl;
              }
            }
        }
      return p;
    }
                 
    // ------------------------- Datas --------------------------------------
  private:
      
      
    // ------------------------- Hidden services ----------------------------
  protected:
    

  }; // end of class MPolynomialReader


  template <int n, typename TRing, typename TAlloc, typename TIterator>
  std::ostream&
  operator<< ( std::ostream & out,
               const MPolynomialReader<n, TRing, TAlloc, TIterator> & object );

  template < int n, typename TRing, class TAlloc >
  std::istream&
  operator>> ( std::istream & in,
               MPolynomial<n,TRing,TAlloc> & aMPolynomial );

} // namespace DGtal


// Includes inline functions.
#include "DGtal/io/readers/MPolynomialReader.ih"

//                                                                           //

#endif // !defined MPolynomialReader_h

#undef MPolynomialReader_RECURSES
#endif // else defined(MPolynomialReader_RECURSES)