MultivariatePolynomial.hpp

Go to the documentation of this file.

/*
 * Copyright(c):
 * Signal Image and Communications (SIC) Department 
 * http://www.sic.sp2mi.univ-poitiers.fr/ 
 *  - University of Poitiers, France http://www.univ-poitiers.fr 
 *  - XLIM  Institute UMR CNRS 7252  http://www.xlim.fr/ 
 * 
 * and 
 * 
 * D2 Fluid, Thermic and Combustion 
 *  - University of Poitiers, France http://www.univ-poitiers.fr 
 *  - PPRIME Institute -  UPR CNRS 3346  http://www.pprime.fr
 *  - ISAE-ENSMA http://www.ensma.fr
 *
 * Contributor(s):
 * The SLIP team,
 * Benoit Tremblais <tremblais_AT_sic.univ-poitiers.fr>,
 * Laurent David <laurent.david_AT_lea.univ-poitiers.fr>, 
 * Ludovic Chatellier <ludovic.chatellier_AT_univ-poitiers.fr>, 
 * Lionel Thomas <lionel.thomas_AT_univ-poitiers.fr>, 
 * Denis Arrivault <arrivault_AT_sic.univ-poitiers.fr>, 
 * Julien Dombre <julien.dombre_AT_univ-poitiers.fr>.
 *
 * Description:
 * The Simple Library of Image Processing (SLIP) is a new image processing 
 * library. It is written in the C++ language following as much as possible 
 * the ISO/ANSI C++ standard. It is consequently compatible with any system  
 * satisfying the ANSI C++ complience. It works on different Unix , Linux ,  
 * Mircrosoft Windows and Mac OS X plateforms. SLIP is a research library that  
 * was created by the Signal, Image and Communications (SIC) departement of  
 * the XLIM, UMR 7252 CNRS Institute in collaboration with the Fluids, Thermic  
 * and Combustion departement of the P', UPR 3346 CNRS Institute of the  
 * University of Poitiers.
 * 
 * The SLIP Library source code has been registered to the APP (French Agency 
 * for the Protection of Programs) by the University of Poitiers and CNRS, 
 * under  registration number IDDN.FR.001.300034.000.S.P.2010.000.21000.

 * http://www.sic.sp2mi.univ-poitiers.fr/slip/
 *
 * This software is governed by the CeCILL-C license under French law and
 * abiding by the rules of distribution of free software.  You can  use, 
 * modify and/ or redistribute the software under the terms of the CeCILL-C 
 * license as circulated by CEA, CNRS and INRIA at the following URL 
 * http://www.cecill.info.
 * As a counterpart to the access to the source code and  rights to copy,
 * modify and redistribute granted by the license, users are provided only
 * with a limited warranty  and the software's author,  the holder of the
 * economic rights,  and the successive licensors  have only  limited
 * liability.
 * 
 * In this respect, the user's attention is drawn to the risks associated
 * with loading,  using,  modifying and/or developing or reproducing the
 * software by the user in light of its specific status of free software,
 * that may mean  that it is complicated to manipulate,  and  that  also
 * therefore means  that it is reserved for developers  and  experienced
 * professionals having in-depth computer knowledge. Users are therefore
 * encouraged to load and test the software's suitability as regards their
 * requirements in conditions enabling the security of their systems and/or
 * data to be ensured and,  more generally, to use and operate it in the 
 * same conditions as regards security. 
 *
 * The fact that you are presently reading this means that you have had
 * knowledge of the CeCILL-C license and that you accept its terms.
 */



#ifndef SLIP_MULTIVARIATEPOLYNOMIAL_HPP
#define SLIP_MULTIVARIATEPOLYNOMIAL_HPP

#include <iostream>
#include <iterator>
#include <cassert>
#include <numeric>
#include <algorithm>
#include <map>
#include <cmath>
#include <complex>
#include <string>
#include "Block.hpp"
#include "Array2d.hpp"
#include "Polynomial.hpp"
#include "linear_least_squares.hpp"
#include "statistics.hpp"
#include "polynomial_algo.hpp"

#include <boost/serialization/access.hpp>
#include <boost/serialization/split_member.hpp>
#include <boost/serialization/complex.hpp>
#include <boost/serialization/version.hpp>
#include <boost/serialization/base_object.hpp>
#include <boost/serialization/map.hpp>

namespace slip
{
template <std::size_t DIM>
struct Monomial;

template <std::size_t DIM>
std::ostream& operator<<(std::ostream & out, 
                                                 const Monomial<DIM>& m);

template <std::size_t DIM>
bool operator==(const Monomial<DIM>& x, 
                                const Monomial<DIM>& y);

template <std::size_t DIM>
bool operator!=(const Monomial<DIM>& x, 
                                const Monomial<DIM>& y);

template <std::size_t DIM>
bool operator<(const Monomial<DIM>& x, 
                       const Monomial<DIM>& y);

template <std::size_t DIM>
bool operator>(const Monomial<DIM>& x, 
                       const Monomial<DIM>& y);

template <std::size_t DIM>
Monomial<DIM> operator*(const Monomial<DIM>& x, 
                                                const Monomial<DIM>& y);


template <std::size_t DIM>
struct Monomial
{
  slip::block<unsigned char,DIM> powers;

  std::string name() const {return "Monomial";}
  std::size_t dim() const {return DIM;}
  std::size_t degree() const 
  {
    return std::accumulate(powers.begin(),powers.end(),std::size_t(0));
  }

  bool is_a_constant() const
  {
    return (std::size_t(std::count(powers.begin(),powers.end(),(unsigned char)(0))) == DIM);
  }

 
  
  friend bool operator== <>(const Monomial<DIM>& x, 
                                                    const Monomial<DIM>& y);

  friend bool operator!= <>(const Monomial<DIM>& x, 
                                                    const Monomial<DIM>& y);


  friend bool operator< <>(const Monomial<DIM>& x, 
                                                   const Monomial<DIM>& y);

  friend bool operator> <>(const Monomial<DIM>& x, 
                                                   const Monomial<DIM>& y);

private:
    friend class boost::serialization::access;
    template<class Archive>
    void save(Archive & ar, const unsigned int version) const
    {
      ar & powers;
    }
  template<class Archive>
    void load(Archive & ar, const unsigned int version)
    {
      ar & powers;
    }
   BOOST_SERIALIZATION_SPLIT_MEMBER()
  

};

}//end slip::


namespace slip{

  template<std::size_t DIM>
  Monomial<DIM> operator*(const Monomial<DIM>& monomial1, 
                                                  const Monomial<DIM>& monomial2);
}
  /* @{ */
namespace slip
{
  template <std::size_t DIM>
  inline
  std::ostream& operator<<(std::ostream & out, 
                                                   const Monomial<DIM>& m)
  {
    if(m.is_a_constant())
      {
      }
    else
      {
                for(std::size_t i = 0; i < (m.powers).size(); ++i)
                  out<<"x"<<(i+1)<<"^"<<(int)m.powers[i]<<" ";
      }
    return out;
  }
/* @} */
  /* @{ */
  template <std::size_t DIM>
  inline
  bool operator==(const Monomial<DIM>& x, 
                                  const Monomial<DIM>& y)
  {
    return (x.powers == y.powers);
  }
  
  template <std::size_t DIM>
  inline
  bool operator!=(const Monomial<DIM>& x, 
                                  const Monomial<DIM>& y)
  {
    return !(x == y);
  }

/* @} */

  /* @{ */
  template <std::size_t DIM>
  inline
  bool operator<(const Monomial<DIM>& x, 
                                 const Monomial<DIM>& y)
  {
    return std::lexicographical_compare((x.powers).rbegin(), (x.powers).rend(),
                                                                                (y.powers).rbegin(), (y.powers).rend());
    //return x.powers < y.powers;
  }
  
  
  template <std::size_t DIM>
  inline
  bool operator>(const Monomial<DIM>& x, 
                                 const Monomial<DIM>& y)
  {
    return (y < x);
  }

/* @} */
 
  template<std::size_t DIM>
  inline
  Monomial<DIM> operator*(const Monomial<DIM>& monomial1,
                                                  const Monomial<DIM>& monomial2)
  {
    Monomial<DIM> tmp;
    std::transform((monomial1.powers).begin(),
                                   (monomial1.powers).end(),
                                   (monomial2.powers).begin(),
                                   (tmp.powers).begin(),
                                   std::plus<unsigned char>());
    return tmp;
  }

}//end slip::
namespace slip
{

template <typename T, std::size_t DIM>
class MultivariatePolynomial;

template <typename T, std::size_t DIM>
std::ostream& operator<<(std::ostream & out, 
                                                 const MultivariatePolynomial<T,DIM>& v);

template <typename T, std::size_t DIM>
class MultivariatePolynomial:public std::map<slip::Monomial<DIM>,T>
{
public :

  typedef std::map<slip::Monomial<DIM>,T> base;
  typedef typename std::map<slip::Monomial<DIM>,T>::key_type key_type;
  typedef typename std::map<slip::Monomial<DIM>,T>::mapped_type mapped_type;
  typedef typename std::map<slip::Monomial<DIM>,T>::value_type value_type;
  typedef typename std::map<slip::Monomial<DIM>,T>::key_compare key_compare;
  typedef typename std::map<slip::Monomial<DIM>,T>::value_compare value_compare;
  typedef MultivariatePolynomial<T,DIM> self;
  typedef const MultivariatePolynomial<T,DIM> const_self;

  typedef typename std::map<slip::Monomial<DIM>,T>::pointer pointer;
  typedef  const typename std::map<slip::Monomial<DIM>,T>::const_pointer const_pointer;
  typedef typename std::map<slip::Monomial<DIM>,T>::reference reference;
  typedef typename std::map<slip::Monomial<DIM>,T>::const_reference const_reference;

  typedef typename std::map<slip::Monomial<DIM>,T>::difference_type difference_type;
  typedef typename std::map<slip::Monomial<DIM>,T>::size_type size_type;

  typedef typename std::map<slip::Monomial<DIM>,T>::iterator iterator;
  typedef typename std::map<slip::Monomial<DIM>,T>::const_iterator const_iterator;

  typedef typename std::map<slip::Monomial<DIM>,T>::reverse_iterator reverse_iterator;
  typedef typename std::map<slip::Monomial<DIM>,T>::const_reverse_iterator const_reverse_iterator;


  
  MultivariatePolynomial():
    base()
  {
    slip::Monomial<DIM> constant;
    constant.powers.fill((unsigned char)0);
    this->insert(constant,T());
  }
  
  MultivariatePolynomial(const std::size_t degree):
    base()
  {

    this->create(degree);
  }

  template<typename InputIterator>
  MultivariatePolynomial(const std::size_t degree,
                                                 InputIterator first,
                                                 InputIterator last):
    base()
  {
    this->create(degree);
    assert(int(last-first)== int(this->total_nb_coeff()));
    
    for(iterator it = this->begin(); it!=this->end();++it)
      {
                (*it).second = *first++;
      }
  }

  template<class InputIterator>
  MultivariatePolynomial(InputIterator first, InputIterator last):
    base(first,last)
  {}

  

  MultivariatePolynomial(const self& rhs):
    base(rhs)
  {}

  ~MultivariatePolynomial()
  {}

  friend std::ostream& operator<< <>(std::ostream & out, 
                                                                     const self& v);
  void insert(const slip::Monomial<DIM>& monomial, const T& coefficient);

  self& operator+=(const T& val);
  self& operator-=(const T& val);
  self& operator*=(const T& val);
  self& operator/=(const T& val);

  self  operator-() const;


  self& operator+=(const self& rhs);
  self& operator-=(const self& rhs);
  self& operator*=(const self& rhs);

    

  
  self derivative(const std::size_t dim) const;

  void derivative(const std::size_t dim);
  

  self derivative(const std::size_t dim,
                                  const std::size_t order) const;

  void derivative(const std::size_t dim,
                                   const std::size_t order);


  self integral(const std::size_t dim,
                                 const T& K = T()) const;

  void integral(const std::size_t dim,
                                const T& K = T());

  self integral(const std::size_t dim,
                                 const std::size_t order,
                                 const T& K = T()) const;

  
  void integral(const std::size_t dim,
                                 const std::size_t order,
                                 const T& K = T());

  template <typename Vector1, typename Vector2>
  double fit(Vector1 X,
                     Vector2 Y,
                     const std::size_t order)
  {
                
    //input data
    //computes the number of coefficients estimate
    slip::Array<std::size_t> V(DIM+order);
    slip::iota(V.begin(),V.end(),1,1);
    std::size_t coeff_size  = (std::accumulate(V.begin()+order,V.end(),1,std::multiplies<int>()))/(std::accumulate(V.begin(),V.begin()+DIM,1,std::multiplies<int>()));

    slip::Array<T> Coeff(coeff_size);
    //output data
    double chi_square = 0.0;
    slip::EvalPowerNdBasis<typename Vector1::value_type,typename slip::Array<T>::iterator> power_basis;
    slip::svd_least_square_nd(X,
                                                      Y,
                                                      Coeff,
                                                      order,
                                                      power_basis);
    this->create(order);
    typename slip::MultivariatePolynomial<T,DIM>::iterator it = this->begin();
    typename slip::Array<T>::const_iterator ita = Coeff.begin();
    for(;it != this->end();++it,++ita)
      {
                (*it).second = *ita;
      }
    return chi_square;
  }
  

  template<typename InputIterator>
  typename std::iterator_traits<InputIterator>::value_type
  evaluate(InputIterator first,
                   InputIterator last) const
  {
    assert((last-first) == DIM);
    typedef  typename slip::MultivariatePolynomial<T,DIM>::const_iterator const_it;
 
    typedef  typename std::iterator_traits<InputIterator>::value_type Return;
    slip::Array2d<Return> values(DIM,this->degree()+1);
   
    //pre-computation of the variables' powers
    for(std::size_t i = 0; i < values.rows(); ++i)
      {
                slip::eval_power_basis(*first++,this->degree(),values.row_begin(i),values.row_end(i));
      }
    
 
      Return result = Return();
    for(const_it first_m = this->begin();first_m!=this->end();++first_m)
      {
                Return tmp = values[0][((*first_m).first).powers[0]];
                for(std::size_t i = 1; i < DIM; ++i)
                  {
                    tmp*=values[i][((*first_m).first).powers[i]];
                  }
                result+= ((*first_m).second) * tmp;
      }

     return result;
  }


  MultivariatePolynomial<T,DIM> partial_evaluate(const std::size_t number,
                                                                                                 const T& val);
 

  self compose(const std::vector<self>& VP);


  std::string name() const;


  size_type degree() const;

  size_type order() const;

  size_type total_nb_coeff() const;

  template<typename RandomAccessIterator>
  void all_coeff(RandomAccessIterator first,
                                 RandomAccessIterator last) const
  {
    typedef typename std::iterator_traits<RandomAccessIterator>::difference_type _Distance;
    assert((last - first) == _Distance(this->total_nb_coeff()));
    assert(DIM > 0);

    //array of all the powers combinations sort by lexicographical order
  
    slip::Array2d<std::size_t> powers_a(slip::power(this->order()+1,DIM),DIM);
    for(std::size_t d = 0; d < powers_a.dim2(); ++d)
      {
                std::size_t step = slip::power((this->order()+1),d);
                for(std::size_t ord = 0, ord2 = 0; ord <= (powers_a.dim1() - step); ord+=step,++ord2)
                  {
                    std::size_t tmp = ord2%(this->order()+1);
                    for(std::size_t k = 0; k < step; ++k)
                      {
                                powers_a[ord+k][d] = tmp;
                      }
                  }
      }

    std::map<slip::Monomial<DIM>,T> tmp_map;
    slip::Array2d<std::size_t> powers_cut(this->total_nb_coeff(),DIM);

    for(std::size_t ord = 0; ord < powers_a.dim1(); ++ord)
      {
                 Monomial<DIM> m;
                 std::copy(powers_a.row_begin(ord),powers_a.row_end(ord),(m.powers).begin());
                 if(m.degree() <= this->order())
                   {
                     tmp_map[m] = T();
                   }
      }
   
    
    typedef typename slip::MultivariatePolynomial<T,DIM>::const_iterator const_it;
    
    for(const_it it = this->begin();it!=this->end();++it)
      {
                Monomial<DIM> m = (*it).first;
                T coefficient = (*it).second;
                tmp_map[m] += coefficient;
      }

    for(const_it it = tmp_map.begin(); it!=tmp_map.end();++it)
      {
                *first++ = (*it).second;
      }
    
  }

  
private:
  void create(const std::size_t degree);
 friend class boost::serialization::access;
    template<class Archive>
    void save(Archive & ar, const unsigned int version) const
    {
      ar & boost::serialization::base_object<std::map<slip::Monomial<DIM>,T> >(*this);
    }
  template<class Archive>
    void load(Archive & ar, const unsigned int version)
    {
      ar & boost::serialization::base_object<std::map<slip::Monomial<DIM>,T> >(*this);
    }
   BOOST_SERIALIZATION_SPLIT_MEMBER()
};


template<typename T, std::size_t DIM>
MultivariatePolynomial<T,DIM> operator+(const MultivariatePolynomial<T,DIM>& P1, 
                                                                    const MultivariatePolynomial<T,DIM>& P2);

template<typename T, std::size_t DIM>
MultivariatePolynomial<T,DIM> operator+(const MultivariatePolynomial<T,DIM>& P, 
                                                                                const T& val);

template<typename T, std::size_t DIM>
MultivariatePolynomial<T,DIM> operator+(const T& val, 
                                                                                const MultivariatePolynomial<T,DIM>& P);


template<typename T, std::size_t DIM>
MultivariatePolynomial<T,DIM> operator-(const MultivariatePolynomial<T,DIM>& P1, 
                                                                                const MultivariatePolynomial<T,DIM>& P2);

  template<typename T, std::size_t DIM>
  MultivariatePolynomial<T,DIM> operator-(const MultivariatePolynomial<T,DIM>& P, 
                                                                                  const T& val);

template<typename T, std::size_t DIM>
MultivariatePolynomial<T,DIM> operator-(const T& val, 
                                                                                const MultivariatePolynomial<T,DIM>& P);

template<typename T, std::size_t DIM>
MultivariatePolynomial<T,DIM> operator*(const MultivariatePolynomial<T,DIM>& P1, 
                                                                                const MultivariatePolynomial<T,DIM>& P2);

template<typename T, std::size_t DIM>
MultivariatePolynomial<T,DIM> operator*(const MultivariatePolynomial<T,DIM>& P, 
                                                                    const T& val);

template<typename T, std::size_t DIM>
MultivariatePolynomial<T,DIM> operator*(const T& val, 
                                                                    const MultivariatePolynomial<T,DIM>& P);


template<typename T, std::size_t DIM>
MultivariatePolynomial<T,DIM> operator/(const MultivariatePolynomial<T,DIM>& P, 
                                                                    const T& val);
}//slip::



namespace slip
{

  template<typename T, std::size_t DIM>
  inline
  void 
  MultivariatePolynomial<T,DIM>::insert(const slip::Monomial<DIM>& monomial, 
                                                                                const T& coefficient)
  {
    if(monomial.is_a_constant()) 
      {
                (*this)[monomial] += coefficient;
      }    
    else
      {
                if(coefficient != T())
                  {
                    (*this)[monomial] += coefficient;
                  }
      }
    
  }

 
  template <typename T, std::size_t DIM>
  inline
  std::ostream& operator<<(std::ostream & out, 
                                                   const MultivariatePolynomial<T,DIM>& m)
  {
    typename slip::MultivariatePolynomial<T,DIM>::const_iterator it = m.begin();
    out.setf(std::ios_base::scientific);
    out.setf(std::ios_base::internal);
    
    
    for(;it!=m.end();++it)
      {
                Monomial<DIM> m = (*it).first;
                T coefficient = (*it).second;
                if(m.is_a_constant())
                  {
                    out<<coefficient<<" ";
                  }
                else
                  {
                    if(coefficient >= T(0))
                      {
                                out<<"+ "<<coefficient;
                      }
                    else
                      {
                                out<<"- "<<std::abs(coefficient);
                      }
                    out<<" "<<(*it).first<<" ";
                  }
      }
    out<<std::endl;
    return out;
  }


template <std::size_t DIM>
  inline
  std::ostream& operator<<(std::ostream & out, 
                                                   const MultivariatePolynomial<std::complex<double>,DIM>& m)
  {
    typename slip::MultivariatePolynomial<std::complex<double>,DIM>::const_iterator it = m.begin();
    out.setf(std::ios_base::scientific);
    out.setf(std::ios_base::internal);
    
    
    for(;it!=m.end();++it)
      {
                Monomial<DIM> m = (*it).first;
                std::complex<double> coefficient = (*it).second;
                if(m.is_a_constant())
                  {
                    out<<coefficient<<" ";
                  }
                else
                  {
                   
                    out<<"+ "<<coefficient;
                    out<<" "<<(*it).first<<" ";
                  }
      }
    out<<std::endl;
    return out;
  }


  template <std::size_t DIM>
  inline
  std::ostream& operator<<(std::ostream & out, 
                                                   const MultivariatePolynomial<std::complex<float>,DIM>& m)
  {
    typename slip::MultivariatePolynomial<std::complex<float>,DIM>::const_iterator it = m.begin();
    out.setf(std::ios_base::scientific);
    out.setf(std::ios_base::internal);
    
    
    for(;it!=m.end();++it)
      {
                Monomial<DIM> m = (*it).first;
                std::complex<float> coefficient = (*it).second;
                if(m.is_a_constant())
                  {
                    out<<coefficient<<" ";
                  }
                else
                  {
                   
                    out<<"+ "<<coefficient;
                    out<<" "<<(*it).first<<" ";
                  }
      }
    out<<std::endl;
    return out;
  }

 template <std::size_t DIM>
  inline
  std::ostream& operator<<(std::ostream & out, 
                                                   const MultivariatePolynomial<std::complex<long double>,DIM>& m)
  {
    typename slip::MultivariatePolynomial<std::complex<long double>,DIM>::const_iterator it = m.begin();
    out.setf(std::ios_base::scientific);
    out.setf(std::ios_base::internal);
    
    
    for(;it!=m.end();++it)
      {
                Monomial<DIM> m = (*it).first;
                std::complex<long double> coefficient = (*it).second;
                if(m.is_a_constant())
                  {
                    out<<coefficient<<" ";
                  }
                else
                  {
                   
                    out<<"+ "<<coefficient;
                    out<<" "<<(*it).first<<" ";
                  }
      }
    out<<std::endl;
    return out;
  }


 // template<typename T, std::size_t DIM>
 //  inline
 //  MultivariatePolynomial<T,DIM>& MultivariatePolynomial<T,DIM>::derivative(const std::size_t dim) 
 //  {
 //    assert(dim > 0);
 //    assert(dim <= DIM);

 //    if(this->degree()!=0)
 //      {
 //             slip::MultivariatePolynomial<T,DIM>* tmp = new slip::MultivariatePolynomial<T,DIM>(*this);
 //             this->clear();
 //             typename slip::MultivariatePolynomial<T,DIM>::iterator it = tmp->begin();
 //             for(;it!=tmp->end();++it)
 //               {
 //                 slip::Monomial<DIM> tmp_m = (*it).first;
 //                 T tmp_coeff = (*it).second * tmp_m.powers[dim-1];
 //                 int p =  (int)tmp_m.powers[dim-1] - 1;
 //                 if(p >= 0)
 //                   {
 //                             tmp_m.powers[dim-1]-=1;
 //                             this->insert(tmp_m,tmp_coeff);
 //                   }
 //               }
 //              tmp->clear();
 //              delete tmp;
 //      }
 //     else
 //      {
 //             slip::Monomial<DIM> constant;
 //             constant.powers.fill((unsigned char)0);
 //             this->insert(constant,T());
 //      }
 //      return *this;
 //  }

 // template<typename T, std::size_t DIM>
 //  inline
 //  MultivariatePolynomial<T,DIM>& MultivariatePolynomial<T,DIM>::derivative(const std::size_t dim,
 // const std::size_t order) 
 //  {
 //    for(std::size_t i = 0; i < order; ++i)
 //      {
 //             this->derivative(dim);
 //      }
 //    return *this;
 //  }


 template<typename T, std::size_t DIM>
 inline
 void MultivariatePolynomial<T,DIM>::derivative(const std::size_t dim) 
  {
    assert(dim > 0);
    assert(dim <= DIM);

    if(this->degree()!=0)
      {
                slip::MultivariatePolynomial<T,DIM>* tmp = new slip::MultivariatePolynomial<T,DIM>(*this);
                this->clear();
                typename slip::MultivariatePolynomial<T,DIM>::iterator it = tmp->begin();
                for(;it!=tmp->end();++it)
                  {
                    slip::Monomial<DIM> tmp_m = (*it).first;
                    T tmp_coeff = (*it).second * tmp_m.powers[dim-1];
                    int p =  (int)tmp_m.powers[dim-1] - 1;
                    if(p >= 0)
                      {
                                tmp_m.powers[dim-1]-=1;
                                this->insert(tmp_m,tmp_coeff);
                      }
                  }
                 tmp->clear();
                 delete tmp;
      }
     else
      {
                slip::Monomial<DIM> constant;
                constant.powers.fill((unsigned char)0);
                this->insert(constant,T());
      }
    // return *this;
  }

 template<typename T, std::size_t DIM>
  inline
  MultivariatePolynomial<T,DIM> MultivariatePolynomial<T,DIM>::derivative(const std::size_t dim) const 
  {
    assert(dim > 0);
    assert(dim <= DIM);

    slip::MultivariatePolynomial<T,DIM> tmp(*this);
    tmp.derivative(dim);
    return tmp;
  }


 template<typename T, std::size_t DIM>
 inline
 void MultivariatePolynomial<T,DIM>::derivative(const std::size_t dim,
                                                                                                const std::size_t order) 
  {
    for(std::size_t i = 0; i < order; ++i)
      {
                this->derivative(dim);
      }
  }

 template<typename T, std::size_t DIM>
  inline
  MultivariatePolynomial<T,DIM> MultivariatePolynomial<T,DIM>::derivative(const std::size_t dim,
 const std::size_t order) const
  {
     slip::MultivariatePolynomial<T,DIM> tmp(*this);
     tmp.derivative(dim,order);
     return tmp;
  }



  // template<typename T, std::size_t DIM>
 //  inline
 //  MultivariatePolynomial<T,DIM>& MultivariatePolynomial<T,DIM>::integral(const std::size_t dim,
 //                                                                                                                                              const T& K) 
 //  {
 //    assert(dim > 0);
 //    assert(dim <= DIM);
 //   slip::MultivariatePolynomial<T,DIM>* tmp = new slip::MultivariatePolynomial<T,DIM>(*this);
 //   this->clear();
 //    typename slip::MultivariatePolynomial<T,DIM>::iterator it = tmp->begin();
 //    for(;it!=tmp->end();++it)
 //      {
 //             slip::Monomial<DIM> tmp_m;
 //             tmp_m.powers = (it->first).powers;
 //             T tmp_coeff = it->second / T(tmp_m.powers[dim-1] + 1);
 //             tmp_m.powers[dim-1]+=1;
 //             if(tmp_coeff != 0)
 //               {
 //                 (*this)[tmp_m] = tmp_coeff;
 //               }
 //       }
  
 //    slip::Monomial<DIM> constant;
 //    constant.powers.fill((unsigned char)0);
 //    this->insert(constant,K);
 //    tmp->clear();
 //    delete tmp;
 //    return *this;
 //  }

 // template<typename T, std::size_t DIM>
 //  inline
 //  MultivariatePolynomial<T,DIM>& 
 // MultivariatePolynomial<T,DIM>::integral(const std::size_t dim,
 //                                                                              const std::size_t order,
 //                                                                              const T& K) 
 //  {
 //     for(std::size_t i = 0; i < order; ++i)
 //      {
 //             this->integral(dim,K);
 //      }
 //    return *this;
 //  }


 template<typename T, std::size_t DIM>
  inline
  void MultivariatePolynomial<T,DIM>::integral(const std::size_t dim,
                                                                                       const T& K)
  {
    assert(dim > 0);
    assert(dim <= DIM);
   slip::MultivariatePolynomial<T,DIM>* tmp = new slip::MultivariatePolynomial<T,DIM>(*this);
   this->clear();
    typename slip::MultivariatePolynomial<T,DIM>::iterator it = tmp->begin();
    for(;it!=tmp->end();++it)
      {
                slip::Monomial<DIM> tmp_m;
                tmp_m.powers = (it->first).powers;
                T tmp_coeff = it->second / T(tmp_m.powers[dim-1] + 1);
                tmp_m.powers[dim-1]+=1;
                if(tmp_coeff != 0)
                  {
                    (*this)[tmp_m] = tmp_coeff;
                  }
       }
  
    slip::Monomial<DIM> constant;
    constant.powers.fill((unsigned char)0);
    this->insert(constant,K);
    tmp->clear();
    delete tmp;
  }

  template<typename T, std::size_t DIM>
  inline
  MultivariatePolynomial<T,DIM> MultivariatePolynomial<T,DIM>::integral(const std::size_t dim,
                                                                                                                                                 const T& K) const
  {
    MultivariatePolynomial<T,DIM> tmp(*this);
    tmp.integral(dim,K);
    return tmp;
  }

 template<typename T, std::size_t DIM>
  inline
  void
 MultivariatePolynomial<T,DIM>::integral(const std::size_t dim,
                                                                                 const std::size_t order,
                                                                                 const T& K) 
  {
     for(std::size_t i = 0; i < order; ++i)
      {
                this->integral(dim,K);
      }
  }

  template<typename T, std::size_t DIM>
  inline
  MultivariatePolynomial<T,DIM> 
  MultivariatePolynomial<T,DIM>::integral(const std::size_t dim,
                                                                                  const std::size_t order,
                                                                                  const T& K) const
  {
   MultivariatePolynomial<T,DIM> tmp(*this);
   tmp.integral(dim,order,K);
   return tmp;
  }



  template<typename T, std::size_t DIM>
  inline
  MultivariatePolynomial<T,DIM> 
  MultivariatePolynomial<T,DIM>::partial_evaluate(const std::size_t number,
                                                                                                  const T& val)
  {
    assert(number <= DIM);
    MultivariatePolynomial<T,DIM> pol;   
    typename slip::MultivariatePolynomial<T,DIM>::const_iterator it = this->begin();
   
    for(;it!=this->end();++it)
      {
                Monomial<DIM> m((*it).first); 
                T coefficient = (*it).second * std::pow(val,((*it).first).powers[number-1]);
                m.powers[number-1] = 0;
                pol.insert(m,coefficient);
      }
    return pol;
    
  }
  
  template<typename T, std::size_t DIM>
  inline
  MultivariatePolynomial<T,DIM> 
  MultivariatePolynomial<T,DIM>::compose(const std::vector<self>& VP)
  {
    assert(VP.size() == DIM);
  //compute Polynomial powers
  const std::size_t P_degree = this->degree();
  std::vector<std::size_t> V(DIM+P_degree);
  slip::iota(V.begin(),V.end(),1,1);
   std::size_t result_size = (std::accumulate(V.begin()+P_degree,V.end(),1,std::multiplies<int>()))/(std::accumulate(V.begin(),V.begin()+DIM,1,std::multiplies<int>()));
  std::vector<slip::MultivariatePolynomial<T,DIM> > PolyPowers(result_size);
  slip::eval_multi_poly_power_nd_basis(VP,P_degree,PolyPowers.begin(),PolyPowers.end());
    
  std::vector<T> all_coeff(result_size);
  this->all_coeff(all_coeff.begin(),all_coeff.end());

  slip::MultivariatePolynomial<T,DIM> Result =
    std::inner_product(PolyPowers.begin(),PolyPowers.end(),
                                       all_coeff.begin(),slip::MultivariatePolynomial<T,DIM>());
  
  
 
  return Result;
  }

  template<typename T, std::size_t DIM>
  inline
  std::string 
  MultivariatePolynomial<T,DIM>::name() const
  {
    return "MultivariatePolynomial";
  }

  template<typename T, std::size_t DIM>
  inline
  typename MultivariatePolynomial<T,DIM>::size_type MultivariatePolynomial<T,DIM>::degree() const {
    if(this->size() != 0)
      {
                 typename slip::MultivariatePolynomial<T,DIM>::const_iterator it = this->end();
                 it--;
                 
                 return (it->first).degree();
      }
    else
      return 0;
  }

  template<typename T, std::size_t DIM>
  inline
  typename MultivariatePolynomial<T,DIM>::size_type MultivariatePolynomial<T,DIM>::order() const 
  {
    return this->degree();
  }

  template<typename T, std::size_t DIM>
  inline
  typename MultivariatePolynomial<T,DIM>::size_type MultivariatePolynomial<T,DIM>::total_nb_coeff() const 
  {
    std::size_t nb_coeff = 1;
    std::size_t fact = 1;
    std::size_t degree = this->degree();
    for(std::size_t i = 1; i <= DIM;++i)
      {
                nb_coeff *= degree + i;
                fact *= i;
      }
    nb_coeff/=fact;
    return nb_coeff;
  }

  

  template<typename T, std::size_t DIM>
  inline
  MultivariatePolynomial<T,DIM>& 
  MultivariatePolynomial<T,DIM>::operator+=(const T& val)
  {
    typename slip::MultivariatePolynomial<T,DIM>::iterator it = this->begin();

    if(((*it).first).is_a_constant())
      {
                (*it).second += val;
      }
    else
      {
                slip::Monomial<DIM> constant;
                constant.powers.fill((unsigned char)0);
                this->insert(constant,val);
      }
    return *this;
  }

  template<typename T, std::size_t DIM>
  inline
  MultivariatePolynomial<T,DIM>& 
  MultivariatePolynomial<T,DIM>::operator-=(const T& val)
  {
    typename slip::MultivariatePolynomial<T,DIM>::iterator it = this->begin();

    if(((*it).first).is_a_constant())
      {
                (*it).second -= val;
      }
    else
      {
                slip::Monomial<DIM> constant;
                constant.powers.fill((unsigned char)0);
                this->insert(constant,-val);
      }
    return *this;
  }

  template<typename T, std::size_t DIM>
  inline
  MultivariatePolynomial<T,DIM>& 
  MultivariatePolynomial<T,DIM>::operator*=(const T& val)
  {
    typename slip::MultivariatePolynomial<T,DIM>::iterator it = this->begin();
    typename slip::MultivariatePolynomial<T,DIM>::iterator last = this->end();
    for(;it!=last;++it)
      {
                (*it).second *= val;            
      }
    return *this;
  }

  template<typename T, std::size_t DIM>
  inline
  MultivariatePolynomial<T,DIM>& 
  MultivariatePolynomial<T,DIM>::operator/=(const T& val)
  {
    typename slip::MultivariatePolynomial<T,DIM>::iterator it = this->begin();
    typename slip::MultivariatePolynomial<T,DIM>::iterator last = this->end();
    for(;it!=last;++it)
      {
                (*it).second /= val;            
      }
    return *this;
  }

  template<typename T, std::size_t DIM>
  inline
  MultivariatePolynomial<T,DIM> 
  MultivariatePolynomial<T,DIM>::operator-() const
  {
    MultivariatePolynomial<T,DIM> tmp(*this);
    typename slip::MultivariatePolynomial<T,DIM>::const_iterator it = this->begin();
    typename slip::MultivariatePolynomial<T,DIM>::const_iterator last = this->end();
    typename slip::MultivariatePolynomial<T,DIM>::iterator it_tmp = tmp.begin();
   
    for(;it!=last;++it,++it_tmp)
      {
                (*it_tmp).second = -(*it).second;               
      }
    return tmp;
  }


  template<typename T, std::size_t DIM>
  inline
  MultivariatePolynomial<T,DIM>& MultivariatePolynomial<T,DIM>::operator+=(const MultivariatePolynomial<T,DIM>& rhs)
  {
    typename slip::MultivariatePolynomial<T,DIM>::const_iterator it_rhs = rhs.begin();
    typename slip::MultivariatePolynomial<T,DIM>::const_iterator last_rhs = rhs.end();
    
    for(;it_rhs!=last_rhs;++it_rhs)
      {
                typename slip::MultivariatePolynomial<T,DIM>::iterator it_tmp = 
                  this->find((*it_rhs).first);
                if(it_tmp!=this->end())
                  {
                   (*it_tmp).second += (*it_rhs).second;
                    if((*it_tmp).second == T(0))
                      {
                                this->erase(it_tmp);
                      }
                  }
                else
                  {
                    this->insert((*it_rhs).first,(*it_rhs).second);
                  }
      }
    return *this;
  }
  
 template<typename T, std::size_t DIM>
  inline
  MultivariatePolynomial<T,DIM>& MultivariatePolynomial<T,DIM>::operator-=(const MultivariatePolynomial<T,DIM>& rhs)
  {
    typename slip::MultivariatePolynomial<T,DIM>::const_iterator it_rhs = rhs.begin();
    typename slip::MultivariatePolynomial<T,DIM>::const_iterator last_rhs = rhs.end();
    
    for(;it_rhs!=last_rhs;++it_rhs)
      {
                typename slip::MultivariatePolynomial<T,DIM>::iterator it_tmp = 
                  this->find((*it_rhs).first);
                if(it_tmp!=this->end())
                  {
                   (*it_tmp).second -= (*it_rhs).second;
                    if((*it_tmp).second == T(0) && !(*it_tmp).first.is_a_constant())
                      {
                                this->erase(it_tmp);
                      }
                  }
                else
                  {
                    this->insert((*it_rhs).first,-(*it_rhs).second);
                  }
      }
  
    return *this;
  }
  


  template<typename T, std::size_t DIM>
  inline
  MultivariatePolynomial<T,DIM>& 
  MultivariatePolynomial<T,DIM>::operator*=(const MultivariatePolynomial<T,DIM>& rhs)
  {
    MultivariatePolynomial<T,DIM> tmp;

    for(typename MultivariatePolynomial<T,DIM>::const_iterator const_iter1 = 
      this->begin();const_iter1 != this->end(); ++const_iter1)
      {
                
                for(typename MultivariatePolynomial<T,DIM>::const_iterator const_iter2 = rhs.begin(); const_iter2 != rhs.end();++const_iter2)
                  {
                    tmp[(*const_iter1).first * (*const_iter2).first] +=
                        (*const_iter1).second * (*const_iter2).second;
                    if( tmp[(*const_iter1).first * (*const_iter2).first] == T())
                      {
                                tmp.erase((*const_iter1).first * (*const_iter2).first);
                      }
                  }
      }
    
    //this->clear();
    *this = tmp;
    
    return *this;
  }


  template<typename T, std::size_t DIM>
  inline
  void MultivariatePolynomial<T,DIM>::create(const std::size_t order)
  {
    //array of all the powers combinations sort by lexicographical order
     slip::Array2d<std::size_t> powers(slip::power(order+1,DIM),DIM);
    for(std::size_t d = 0; d < powers.dim2(); ++d)
      {
                std::size_t step = slip::power((order+1),d);
                for(std::size_t ord = 0, ord2 = 0; ord <= (powers.dim1() - step); ord+=step,++ord2)
                  {
                    std::size_t tmp = ord2%(order+1);
                    for(std::size_t k = 0; k < step; ++k)
                      {
                                powers[ord+k][d] = tmp;
                      }
                  }
      }

    for(std::size_t k = 0; k < powers.dim1(); ++k)
      {
                slip::Monomial<DIM> m;
                
                std::copy(powers.row_begin(k),powers.row_end(k),(m.powers).begin());
                if (m.degree() <= order)
                  {
                    this->insert(m,T(1));
                  }
      }
  }

}//slip::

  /* @{ */
namespace slip{
template<typename T,std::size_t DIM>
inline
MultivariatePolynomial<T,DIM> operator+(const MultivariatePolynomial<T,DIM>& P1, 
                                                                                const MultivariatePolynomial<T,DIM>& P2)
{
  MultivariatePolynomial<T,DIM> tmp(P1);
  tmp+=P2;
  return tmp;
}

template<typename T,std::size_t DIM>
inline
MultivariatePolynomial<T,DIM> operator+(const MultivariatePolynomial<T,DIM>& P1, 
                                                                                const T& val)
{
  MultivariatePolynomial<T,DIM> tmp(P1);
  tmp+=val;
  return tmp;
}

template<typename T,std::size_t DIM>
inline
MultivariatePolynomial<T,DIM> operator+(const T& val,
                                                                                const MultivariatePolynomial<T,DIM>& P1)
{
  return P1+val;
}


template<typename T,std::size_t DIM>
inline
MultivariatePolynomial<T,DIM> operator-(const MultivariatePolynomial<T,DIM>& P1, 
                                                                                const MultivariatePolynomial<T,DIM>& P2)
{
  MultivariatePolynomial<T,DIM> tmp(P1);
  tmp-=P2;
  return tmp;
}

template<typename T,std::size_t DIM>
inline
MultivariatePolynomial<T,DIM> operator-(const MultivariatePolynomial<T,DIM>& P1, 
                                                                                const T& val)
{
  MultivariatePolynomial<T,DIM> tmp(P1);
  tmp-=val;
  return tmp;
}

template<typename T,std::size_t DIM>
inline
MultivariatePolynomial<T,DIM> operator-(const T& val,
                                                                                const MultivariatePolynomial<T,DIM>& P1)
{
  return -(P1-val);
}


template<typename T,std::size_t DIM>
inline
MultivariatePolynomial<T,DIM> operator*(const MultivariatePolynomial<T,DIM>& P1, 
                                                                                const MultivariatePolynomial<T,DIM>& P2)
{
  MultivariatePolynomial<T,DIM> tmp(P1);
  tmp*=P2;
  return tmp;
}

template<typename T,std::size_t DIM>
inline
MultivariatePolynomial<T,DIM> operator*(const MultivariatePolynomial<T,DIM>& P1, 
                                                                                const T& val)
{
  MultivariatePolynomial<T,DIM> tmp(P1);
  tmp*=val;
  return tmp;
}

template<typename T,std::size_t DIM>
inline
MultivariatePolynomial<T,DIM> operator*(const T& val,
                                                                                const MultivariatePolynomial<T,DIM>& P1)
{
  return P1*val;
}


template<typename T,std::size_t DIM>
inline
MultivariatePolynomial<T,DIM> operator/(const MultivariatePolynomial<T,DIM>& P1, 
                                                                                const T& val)
{
  MultivariatePolynomial<T,DIM> tmp(P1);
  tmp/=val;
  return tmp;
}
/* @} */
}//slip::
#endif //SLIP_MULTIVARIATEPOLYNOMIAL_HPP