/Users/mourrain/Devel/mmx/realroot/include/realroot/bernstein_eenv.hpp

Go to the documentation of this file.

/********************************************************************
 *   This file is part of the source code of the realroot library.
 *   Author(s): J.P. Pavone, GALAAD, INRIA
 *   $Id: eenv.hpp,v 1.1 2005/07/11 10:03:55 jppavone Exp $
 ********************************************************************/
#ifndef realroot_SOLVE_SBDSLV_BERNSTEIN_EENV_H
#define realroot_SOLVE_SBDSLV_BERNSTEIN_EENV_H
//--------------------------------------------------------------------
#include <cstddef>
#include <iostream>
#include <realroot/loops_vctops.hpp>
#include <realroot/loops_brnops.hpp>
#include <realroot/loops_upoldse.hpp>
#include <realroot/bernstein_bzenv.hpp>
//--------------------------------------------------------------------
namespace mmx {
//--------------------------------------------------------------------
#ifndef _valueof_ 
#define _valueof_(x) std::cout << #x" = " << (x) << std::endl;
#endif

namespace bernstein
{
  struct eenv_base
  {
    typedef int  sz_t;
    sz_t  m_nvr, * m_szs, * m_str;
  protected:
    inline void compute_strides()
    {
      m_str[m_nvr-1] = 1;    
      for ( int i = m_nvr-2; i >= -1; i-- ) m_str[i] = m_szs[i+1]*m_str[i+1];
    };
  public:
    eenv_base( sz_t nvr = 0 ) { m_nvr = nvr; };
    inline void   set_szs( sz_t * szs ) { m_szs = szs; compute_strides();  };    
    inline sz_t sz( sz_t v  ) const { return m_szs[v]; };
    inline sz_t nvars()     const { return m_nvr; };
    inline sz_t data_size() const { return m_str[-1]; };

    template<typename real_t>
    void hodograph( real_t * dst, real_t * src, int v  )
    {
      int i,j;
      int k;
      int kstr = (m_szs[v]-1)*m_str[v];
      for ( k = i = 0; i < data_size(); i += m_str[v-1], k += kstr )
        for ( j = 0; j < m_str[v]; j++ )
          brnops::hodograph(dst+k+j,src+i+j,m_szs[v],m_str[v]);
    };
    
    template<typename X, typename real_t>
    void copy( X * dst, real_t const * const src ) { std::copy(src, src + data_size(), dst ); };
    
    template<typename real_t>
    real_t eval( real_t * data,  const real_t *  prm, real_t * chunk =  0 )
    {
      sz_t i,j;
      real_t * tmp;
      if ( chunk ) tmp = chunk; 
      else tmp = new real_t[ data_size() ];
      std::copy(data,data+data_size(), tmp );
      for ( int v = nvars()-1; v >= 0 ; v -- )
        for ( i = 0; i <data_size(); i += m_str[v-1] )
          brnops::decasteljau(tmp+i,m_szs[v],prm[v],m_str[v]);
      real_t val = tmp[0];
      if ( !chunk ) delete[] tmp;
      return val;
    };
    
    template<typename real_t, typename X>
    void spmeval( X& res, real_t* src, X * prm, unsigned * supp, unsigned nsupp ) const
    {
      sz_t v, lg, d;
      unsigned i,add,o;
      lg = 0;
      for ( v = 0; v < m_nvr; v ++ ) if ( m_szs[v] > m_szs[lg] ) lg = v;
      X powers [m_nvr][m_szs[lg]];
      X acc;
      for ( v = 0; v < m_nvr; v ++ )
        for ( powers[v][0] = 1.0, d = 1; d < m_szs[v]; powers[v][d] = prm[v]*powers[v][d-1], d ++ );
      res = 0;
      for ( i = 0; i < nsupp; res += acc, i ++ )
        for ( add = supp[i], acc = src[add], v = m_nvr-1; add; add /= m_szs[v], v -- )
          acc *= powers[v][add%m_szs[v]];
    };
    
    
    template<typename real_t,typename X>
    void meval( X& res, real_t * data,  const X *  prm  ) const
    {
      sz_t i;
      X * tmp = new X[ data_size() ];
      std::copy(data,data+data_size(),tmp);
      for ( int v = nvars()-1; v >= 0 ; v -- )
        for ( i = 0; i < data_size(); i += m_str[v-1] )
          *(tmp+i) = univariate::horner(tmp+i,m_szs[v],prm[v],m_str[v]);
      res = tmp[0];
      delete[] tmp;
    };
    
    template<typename real_t>
    real_t monoms_eval( real_t * data,  const real_t *  prm  ) const
    {
      sz_t i;
      real_t * tmp = new real_t[ data_size() ];
      std::copy(data,data+data_size(), tmp );
      for ( int v = nvars()-1; v >= 0 ; v -- )
        for ( i = 0; i < data_size(); i += m_str[v-1] )
          *(tmp+i) = upoldse::horner(tmp+i,m_szs[v],prm[v],m_str[v]);
      real_t val = tmp[0];
      delete[] tmp;
      return val;
    };
    
     
    template<typename real_t>
    void mins( real_t * _mins_, real_t const * const data, sz_t v ) const 
    {
      sz_t p,k,i,j;
      for ( p = 0, k = 0; k < m_szs[v]; k++, p += m_str[v] )
        {
          _mins_[k] = data[p];
          for ( i = 0; i < data_size(); i += m_str[v-1] )
            for ( j = i; j < i + m_str[v]; j ++ )
              if ( data[j+p] < _mins_[k] ) _mins_[k] = data[j+p];
        };
    };
    
    template<typename real_t>
    void maxs( real_t * _maxs_, real_t const * const data, sz_t v ) const 
    {
      sz_t p,k,i,j;

      for ( p = 0, k = 0; k < m_szs[v]; k++, p += m_str[v] )
        {
          _maxs_[k] = data[p];
          for ( i = 0; i < data_size(); i += m_str[v-1] )
            for ( j = i; j < i + m_str[v]; j ++ )
              if ( data[j+p] > _maxs_[k] ) _maxs_[k] = data[j+p];
        };
    };
        
    template<typename real_t>
    inline void maxs( real_t * _maxs_, real_t const  * const data ) const 
    { for ( sz_t v = 0; v < m_nvr; v ++ ) { maxs(_maxs_,data,v); _maxs_ += m_szs[v];}; };

    template<typename real_t>
    inline void mins( real_t * _mins_, real_t const  * const data ) const 
    { for ( sz_t v = 0; v < m_nvr; v ++ ) { mins(_mins_,data,v); _mins_ += m_szs[v];}; };

    template<typename real_t> 
    void split( real_t * left, real_t * right, int v, const real_t& t )
    {
      int i,j,k;
      for ( i = 0; i < data_size(); i += m_str[v-1] )
        for ( j = i; j < i + m_str[v]; j++ )
          split(left+j,right+j,m_szs[v],m_str[v],t);
    };
  
    template<typename real_t>
    void split2( real_t * left, real_t * right, int v )
    {
      int i,j;
      for ( i = 0; i < data_size(); i += m_str[v-1] )
        for ( j = i; j < i + m_str[v]; j++ )
          split2(left+j,right+j,m_szs[v],m_str[v]);
    };
    
    template<typename real_t>
    void lrestrict( real_t * data, int v, const real_t& mn )
    {
      sz_t i,j;
      for ( i = 0; i < data_size(); i += m_str[v-1] )
        for ( j = i; j < i + m_str[v]; j++ )
          brnops::lrestrict(data+j,m_szs[v],mn,m_str[v]);
    };
    
    template<typename real_t>
    void rrestrict( real_t * data, int v, const real_t& mx )
    {
      int i,j;
      for ( i = 0; i < data_size(); i += m_str[v-1] )
        for ( j = i; j < i + m_str[v]; j++ )
          brnops::rrestrict(data+j,m_szs[v],mx,m_str[v]);
    };
 
    template<typename real_t>
    bool sgnchg( real_t * data ) { return vctops::sgnchg(data,data_size()); };
  
    template<typename real_t>
    inline void scale( real_t * data ) { vctops::scale(data,data_size()); };
    

    /* conversion sur [0,1]^n */
    
    /* monomiale -> bernstein */
    template< typename real_t >     /* conversion partielle */
    void fromMonoms( real_t * data, sz_t v, bzenv<real_t> * env = bzenv<real_t>::_default_ )
    {
      sz_t i,j;
      for ( i = 0; i < data_size(); i += m_str[v-1] )
        for ( j = i; j < i + m_str[v]; j++ )
          env->fromMonoms(data+j,m_szs[v],m_str[v]);
    };
    template< typename real_t >    /* conversion complete */
    inline void fromMonoms( real_t * data, bzenv<real_t> * env = bzenv<real_t>::_default_ ) { for ( sz_t v = 0; v < m_nvr; fromMonoms(data,v,env), v++ ); };
    
    /* bernstein -> monomiale */
    template< typename real_t >   /* conversion partielle */
    void toMonoms( real_t * data, sz_t v, bzenv<real_t> * env = bzenv<real_t>::_default_ )
    {
      sz_t i,j;
      for ( i = 0; i < data_size(); i += m_str[v-1] )
        for ( j = i; j < i + m_str[v]; j++ )
          env->toMonoms(data+j,m_szs[v],m_str[v]);
    };
    template< typename real_t >  /* conversion complete */
    inline void toMonoms(  real_t * data, bzenv<real_t> * env = bzenv<real_t>::_default_ ) { for ( sz_t v = 0; v < m_nvr; toMonoms(data,v,env), v ++ ); };
    
//     template< typename real_t >
//     void mshift( real_t * data, real_t * prm )
//     {
//       sz_t i,j;
//       for ( i = 0; i < data_size(); i += m_str[v-1] )
//      for ( j = i; j < i + m_str[v]; j++ )
//        upoldse::shift(data+j,m_szs[v],m_str[v],prm[v]);
//     };

    template< typename real_t >
    real_t flatness( real_t * data, int v )
    {
      real_t m = 10;
      sz_t i,j;
      for ( i = 0; i < data_size(); i += m_str[v-1] )
        for ( j = i; j < i + m_str[v]; j++ )
          std::max(m,flatness(data+j,m_szs[v],m_str[v]) );
      return m;
    };
  };
  
  
  //  eenv_base + dictionnaire local/global pour les variables
  struct eenv : public eenv_base
  {
    typedef eenv_base::sz_t sz_t;
    struct add
    {
      
    };
    sz_t * m_vrs;
    sz_t * m_pszs;
    sz_t   m_mxvr;
    void _alloc_( sz_t nvr )
    {
      m_nvr  = nvr;  
      m_szs  = new sz_t[4*m_nvr+2]; 
      m_vrs  = m_szs + m_nvr;
      m_str  = m_vrs + m_nvr + 1;
      m_pszs = m_str + m_nvr + 1;
    };
    eenv(){};
    eenv( sz_t szu, sz_t szv )
    {
      _alloc_(2);
      m_szs[0] = szu; m_szs[1] = szv;
      m_vrs[0] = 0; m_vrs[1] = 1;
      m_mxvr = 1;
      compute_strides();
      m_pszs[-1] = 0;
      for ( sz_t v = 0; v < m_nvr; v ++ )
        m_pszs[v] = m_pszs[v-1] + m_szs[v];
    };
    
    eenv( sz_t nvr, sz_t const * szs, sz_t const * vrs )
    {
      _alloc_(nvr);
      std::copy(szs,szs+m_nvr,m_szs);
      std::copy(vrs,vrs+m_nvr,m_vrs);
      compute_strides();
      m_pszs[-1] = 0;
      m_mxvr = m_vrs[0];
      for (sz_t v = 0; v < m_nvr; v ++ ) if ( m_vrs[v] > m_mxvr ) m_mxvr = m_vrs[v];
      for ( sz_t v = 0; v < m_nvr; v ++ )
        m_pszs[v] = m_pszs[v-1] + m_szs[v];
    };
    
    void swap( eenv& e )
    {
      std::swap(m_nvr,e_m.hpp_nvr);
      std::swap(m_szs,e_m.hpp_szs);
      std::swap(m_vrs,e_m.hpp_vrs);
      std::swap(m_pszs,e_m.hpp_pszs);
      std::swap(m_mxvr,e_m.hpp_mxvr);
    };

    sz_t psz(){ return m_pszs[m_nvr-1]; };
    sz_t psft( sz_t v ) { return m_pszs[v-1];}
    
    ~eenv()
    {
      if ( nvars() ) delete[] m_szs;
      m_nvr = 0;
    };
    
    sz_t stride( sz_t v ) const { return m_str[v]; };
    sz_t var ( sz_t v ) const { return m_vrs[v]; };

    int indexofvar( sz_t gv ) const 
    {
      //      if ( gv < m_vrs[0] ||  gv > m_vrs[m_nvr-1] ) return -1;
      // remplacer par dichotomie ?
      for ( sz_t v = 0; v < m_nvr; v ++ )
        if ( m_vrs[v] == gv ) return v;
      return -1;
    };
    
    /* operations */
    /* environement de diff(e,v), v est local a e */ 
    void diff( eenv * e, sz_t v )
    {
      if ( nvars () ) this->~eenv();
      _alloc_( e->nvars() );
      std::copy( e->m_vrs, e->m_vrs+e->nvars(), m_vrs );
      m_mxvr = e->m_mxvr;
      for  ( sz_t i = 0; i < m_nvr; m_szs[i] = e->sz(i), i ++ );
      m_szs[v]--;
      compute_strides();
    };
    
    template<typename real_t>
    static void mdiff( eenv * res, eenv * a, real_t * dst, real_t const * const src, int v )
    {
      sz_t is = 0;
      for ( sz_t i = 0; i < res->data_size(); i += res->m_str[v-1], is += a->m_str[v-1] )
        for ( sz_t j = 0; j < a->m_str[v]; j ++ )
          upoldse::diff(dst+i+j,src+is+j,a->sz(v),a->m_str[v]);
    };
    
    template<typename real_t>
    static void monoms_derivative(eenv * res, eenv * a, real_t ** dst, real_t ** src, int v, int n = 1 )
    {
      res->diff(a,v);
      for ( int c = 0; c < n; c ++ )
        {
          dst[c] = new real_t[ res->data_size() ];
          mdiff( res, a, dst[c], src[c], v );
        };
    };
    
#include "bernstein_eenv_multiply.hpp"
#include "bernstein_eenv_simplify.hpp"
    
    template<typename real_t>
    bool print_monom( std::ostream& o, const real_t& c, sz_t * add, bool first )
    {
      if ( c ) 
        {
          if ( c < 0 )  o << "-";
          else 
            if ( ! first ) o << "+";
          o << std::abs(c);
          for ( sz_t v = 0; v < nvars(); v ++ )
            if ( add[var(v)] ) o << "*x" << var(v) << "^" << add[var(v)];
          return true;
        };
      return false;
    };
    
    template<typename real_t>
    void monoms_print( std::ostream& o, real_t * src )
    {
      unsigned i;
      unsigned mcount = 0;
      int c;
      sz_t add[m_mxvr+1];
      std::fill(add,add+m_mxvr+1,0);
      bool first = true;
      for ( i = 0; i < data_size(); i ++ )
        {
          bool b = print_monom(o,src[i],add,first);
          if ( first && b ) first = false;
          mcount += b;
          c = nvars()-1;
          add[var(c)]++;
          while ( c >0 && (add[var(c)] == m_szs[c])) 
            { 
              add[var(c)] = 0;  
              c--; 
              add[var(c)]++; 
            };
        }
      if ( mcount == 0 ) o << "0";
    };
        
    
    
    
    template<typename real_t>
    static void elevation( eenv * out, eenv * in )
    {
      
    };
    template<typename real_t>
    static void elevation(  eenv * out, eenv * in, real_t * dst, real_t * src, bzenv< real_t > * bznv = bzenv<real_t>::_default_  )
    {
      sz_t io = 0;
      sz_t v  = 0;
      for ( sz_t i = 0; i < in->data_size(); i += in->m_str[v-1], io += out->m_str[v-1] )
        for ( sz_t j = 0; j < in->m_str[v]; j ++ )
          bznv->elevate( dst + io + j, src + i + j , in->m_szs[v], in->m_str[v], out->m_str[v], out->m_szs[v]-in->m_szs[v] );
    };
        
    #include "bernstein_eenv_vrestrict.hpp"
      
  };
  
  
  inline std::ostream& operator<<( std::ostream& out, const eenv& env )
  {
    out << "eenv:";
    out << "\n\tszs = ";
    vctops::print(env_m.hpp_szs,env_m.hpp_nvr,1,out);
    out << "\n\tvrs = ";
    vctops::print(env_m.hpp_vrs,env_m.hpp_nvr,1,out);
    out << "\n";
    out << "\n\ttsz = " << env.data_size() << std::endl;
    return out;
  };
};
//--------------------------------------------------------------------
} //namespace mmx
/********************************************************************/
#endif //