/Users/mourrain/Devel/mmx/realroot/include/realroot/solver_mv_bernstein.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: realroot.h,v 1.1 2005/07/11 10:03:51 jppavone Exp $
 ********************************************************************/
#ifndef realroot_SOLVE_mslvbst_H
#define realroot_SOLVE_mslvbst_H
/********************************************************************/
#include <set>
#include <vector>
#include <math.h>
#include <realroot/texp_kernelof.hpp>
#include <realroot/solver.hpp>

#include <realroot/tensor_bernstein.hpp>

#define REP(x) (x).rep()

#include <realroot/GMP.hpp>
#include <realroot/polynomial.hpp>
#include <realroot/system_method.hpp>
#include <realroot/system_support.hpp>
#include <realroot/Seq.hpp>


//====================================================================
namespace mmx {
//====================================================================
namespace realroot
{
#define SBD_RD       0
#define SBD_RDL      1
#define SBD_RDS      2
#define SBD_RDRDL    3
#define SBD_RDRDLRDS 4
#define SBD_RDRDS    5
#define SBD_RDLRDS   6
#define SBD_RDSRDL   7
#ifndef SBD_MTH
#define SBD_MTH SBD_RDRDL
#endif
  
  template< class system, int mth > /* selection de methodes */
  struct select_mth
  { typedef realroot::method<system,strgy::newton<system>,rdslv::parallel,sbdrl::parametric> result_t; };
  template< class system > 
  struct select_mth<system,SBD_RD>
  { typedef realroot::method<system,strgy::simple<system>,rdslv::parallel,sbdrl::parametric> result_t; };
  template< class system > 
  struct select_mth<system,SBD_RDL>
  { typedef realroot::method<system,strgy::newton<system>,rdslv::parallel,sbdrl::parametric> result_t; };
  template< class system>
  struct select_mth<system,SBD_RDRDL>
  { typedef realroot::method<system,strgy::newton< system, strgy::simple<system> >,rdslv::parallel,sbdrl::parametric> 
    result_t; };
  template< class system >
  struct select_mth<system,SBD_RDRDLRDS>
  { typedef realroot::method<system,strgy::cvmatrix< system, strgy::newton<system,strgy::simple<system> > >,rdslv::parallel,sbdrl::parametric> result_t; };
  template< class system > 
    struct select_mth<system,SBD_RDS>
  { typedef realroot::method<system,strgy::cvmatrix<system>,rdslv::parallel,sbdrl::parametric> result_t; };

  template< class system > 
    struct select_mth<system,SBD_RDRDS>
  { typedef realroot::method<system,strgy::cvmatrix<system, strgy::simple<system> >,rdslv::parallel,sbdrl::parametric> result_t; };
  
  template< class system > 
    struct select_mth<system,SBD_RDLRDS>
  { typedef realroot::method<system,strgy::cvmatrix<system, strgy::newton<system> >,rdslv::parallel,sbdrl::parametric> result_t; };

  template< class system > 
    struct select_mth<system,SBD_RDSRDL>
  { typedef realroot::method<system,strgy::cvmatrix<system, strgy::newton<system> >,rdslv::parallel,sbdrl::parametric> result_t; };
  
  template< class T, class eenv_t, class MPOL >
  void fill_data( T * data, eenv_t * env, const MPOL& p, binomials<T>& benv  )
  {
    std::fill(data,data+env->sz(),0); 
    mpolfill(data,p.rep(),*env);
    convertm2b (data,*env,benv);
  };
  template< class monom >
  void scan_monom_env ( std::set< int >& env, const monom& m )
  { for ( unsigned i = 0; i < (m.rep()).size(); i++ ) if ( (m.rep())[i] != 0 ) env.insert(i); };
}

  template<int MTH> struct ProjRd {};
  
//====================================================================
template<class C, int MTH>
struct solver< C, ProjRd<MTH> >
{
  typedef C coeff_t;
  typedef Seq<std::vector<C> >  Solutions;
  typedef realroot::system<C>   system_t;
  
  static C eps;
  static system_t* sys;

  // solver_multivariate_bernstein( const C& precision = (C)1e-6 ): eps(precision), sys(0){};

  
  static void set_precision(int prec) { using std::pow; eps = pow(2.0,-prec); }

  template< class MPDI >  static
  void init_bernstein( MPDI  first, MPDI last)
  {

    int    ne = 0;
    MPDI eptr;
    for ( eptr = first; eptr != last; eptr ++ ) ne++;
    int nvrs[ne];
    const int *szs[ne];
    int e;
    //    int mxvr = 0;
    std::set<int> gvrs;
    for ( e = 0, eptr = first; eptr != last; eptr ++, e++  ) { 
      nvrs[e] = eptr->rep().nbvar(); 
      szs[e]  = eptr->rep().szs();
      const int * tmp = eptr->rep().vrs();
      for ( int v = 0; v < eptr->rep().nvr(); v ++ ) gvrs.insert(tmp[v]);
    };
    int mxv = *(gvrs.rend());
    int lnames[ mxv+1 ];
    int nv = 0;
    for ( std::set<int>::const_iterator it = gvrs.begin(); it != gvrs.end(); lnames[*it++] = nv++ ) ;
    int * vrdat = new int[ nv * ne ];
    int * vrptr = vrdat;
    const int *  vrs[ne];
    e = 0;
    for ( eptr = first; eptr != last; vrptr+= eptr->rep().nbvar(), e++, eptr ++ ) {
      vrs[e] = vrptr;
      const int * tmp = eptr->rep().vrs();
      for ( int v = 0; v < eptr->rep().nvr();  v ++ )
        vrptr[v] = lnames[ tmp[v] ];    
    };
    //    std::cout << ne << " " << nv << std::endl;
    assert(ne>=nv);
    sys = new realroot::system<C>( ne, nv, nvrs, (const int**)vrs, (const int**)szs, eps );
    e=0; for ( MPDI eptr = first; eptr != last; add_equation( REP(*eptr), e++  ), eptr++ ) ;
    delete[] vrdat;
  };
  
  static void add_equation( const tensor::bernstein<C>& data, int e )
  {
    sys->set_bernstein2( data.begin(), e );
  };
  
  static void add_equation( const tensor::monomials<C>& data, int e ) 
  {
    sys->set_monomial2( data.begin(), e );
  };
  
  static int run( std::vector<C> & R, int& nv )
  {
    R.resize(0);
    typename realroot::select_mth<system_t,MTH>::result_t __mth;
    //int e = 0;for ( MPDI eptr = first; eptr != last; add( eptr->rep(), e++  ), eptr++ );
    //std::vector<C> tmp;
    realroot::system_ctrl< std::vector< C > >  intf(R);
    int nvars = sys->nvr();
    __mth.launch(sys,intf,(C*)0);

    int nsols = R.size() / (2*nvars);
    int nsols2 = realroot::clean_result( &R[0], nvars, nsols,(C)1e-3 );
    R.resize(nsols2 * nvars) ;
    int ns = nsols2;
    nv = nvars;
    return ns;
  };
  
  static int run( std::vector<C> & sol, C* dom, int& nv )
  {
    sol.resize(0);
    typename realroot::select_mth<system_t,MTH>::result_t __mth;
    realroot::system_ctrl< std::vector< C > >  intf(sol);
    int nvars = sys->nvr();
    //    C* dmns = &dom[0];
    // if ( (int)dom.size() >= 2*nvars )
//       {
//      dmns = new C[ 2*nvars ];
//      for ( int i = 0; i < 2*nvars; i ++ ) let::assign(dmns[i],dom[i]);
//       };
    __mth.launch(sys,intf,dom);

    int nsols = sol.size() / (2*nvars);
    int nsols2 = realroot::clean_result( &sol[0], nvars, nsols,(C)1e-3 );
    sol.resize(nsols2 * nvars) ;
    int ns = nsols2;
    nv = nvars;
    return ns;
  };

  //--------------------------------------------------------------------
  template< class output, class MPC, class Domains >
  static void solve_monomial( output& res, 
                              const MPC& l, const Domains & dom, 
                              bool verbose = false )
  {
    typedef C T;
    typedef typename MPC::value_type          mpoly_t;
    typedef MPC                               plist_t;
    typedef typename plist_t::const_iterator  pptr_t;
    typedef typename mpoly_t::const_iterator  mptr_t;
    //    typedef typename mpoly_t::monom_t         monom_t;
    //    typedef solver_multivariate_bernstein<C, MTH>           mth_t;
    typedef int                               sz_t;
    sz_t  nvars = 0, neqs, * nvrs, ** vrs, ** szs; 
    int e;
    pptr_t cpl; 
    // recuperation de l'environement global
    for ( neqs = 0, cpl = l.begin(); cpl != l.end(); cpl++, neqs++  )
      for ( mptr_t mn = cpl->begin(); mn != cpl->end(); mn++ )
        nvars = std::max(nvars,mn->nvars()+1);
    // recuperation de l'environement de chaque equation
    std::set<int> * seenvs = new std::set< int > [ neqs ];
    for ( e = 0, cpl = l.begin(); e < neqs; e++, cpl++  )
      for ( mptr_t mn = cpl->begin(); mn != cpl->end(); mn++ ) 
        realroot::scan_monom_env(seenvs[e],*mn);
    
    //  assert(sizeof(sz_t)==sizeof(sz_t*));
    sz_t  s = 0;
    for ( e = 0; e < neqs; s += seenvs[e].size(), e++ ) ;
    nvrs = new sz_t [neqs+2*s];
    vrs  = new sz_t*[2*neqs];
    szs  = vrs + neqs;
    sz_t * base = nvrs + neqs;
    std::set<int>::iterator seenvit;  
    for ( e = 0, cpl = l.begin(); e < neqs; base += 2*nvrs[e++], cpl++ )
      {
        nvrs[e] = seenvs[e].size();
        vrs [e] = base;
        szs [e] = base + nvrs[e];
        int k = 0;
        for ( seenvit = seenvs[e].begin(); k < nvrs[e]; seenvit++, k++  ) 
        {
          vrs[e][k]  = *seenvit;
          szs[e][k]  = degree( *cpl, *seenvit ) + 1;
        };
      };
    
    std::vector<C> result;
    C * dmns = 0;
    if ( (sz_t)dom.size() >= 2*nvars )
      {
        dmns = new C[ 2*nvars ];
        for ( int i = 0; i < 2*nvars; i ++ ) let::assign(dmns[i],dom[i]);
      };
    
    //  chrono tm;
    std::vector<T> tmp;
    
    {  /* remplissage et appel */
      using namespace realroot;
      binomials<T> bznv;
      typedef  realroot::system<T> system_t;
      system_t sys( neqs, nvars, nvrs, (const int**)vrs, (const int**)szs, (T)eps);
      for ( e = 0, cpl = l.begin(); e < neqs; e ++, cpl++  )
        realroot::fill_data( sys.data(e), sys.env(e), *cpl, bznv );
      realroot::system_ctrl< std::vector< T > >  intf(tmp);
      //    tm.start();
      typename realroot::select_mth<system_t,MTH>::result_t __mth;
      __mth.launch(&sys,intf,dmns);
      //tm.stop();
      if ( verbose )
        {
          std::cout << " iterations  = " << __mth.m_niter << std::endl;
          std::cout << "subdivisions = "<< __mth.m_nsbd  << std::endl;
        };
    };
    
    if ( !verbose && dmns ) delete[] dmns;
    delete[] seenvs;
    delete[] nvrs;
    delete[] vrs;
    if ( verbose )
      {
        std::cout << "<domain>";
        if  ( dmns )
          {
            std::cout << std::endl;
            for ( int i = 0; i < nvars; i ++ )
              std::cout << "\t" 
                        << "x" << i//monom_t::var_of_index(i) 
                        << " = [ " << dmns[2*i] 
                        << "," 
                        << dmns[2*i+1] 
                        << " ]";
          }
        else std::cout << "[0,1]^" << nvars;
        std::cout << "</domain>\n";
        //      std::cout << "\treductions   = " << bsys.nrest() << std::endl;
        //      std::cout << "\titerations   = " << bsys.niter() << std::endl;
        //      std::cout << "<solver_time>"; tm.write(std::cout); 
        //      std::cout <<"</solver_time>"<< std::endl;       
        if ( dmns ) delete[] dmns;
      };
    
    unsigned nsols = tmp.size() / (2*nvars);
    T boxes[nsols*2*nvars];
    std::copy(&tmp[0], &tmp[0]+2*nvars*nsols,boxes);
    int nsols2 = realroot::clean_result( &tmp[0], nvars, nsols,(T)1e-3 );
    if ( verbose ) { 
      //    xml::print(std::cout,nsols,"nb_box")<< std::endl;
      for ( unsigned i = 0; i < nsols; i ++ )
        {
          std::cout << "\t[ ";
          for (int j = 0; j < nvars; j ++ )
            std::cout << (boxes[2*(i*nvars+j)]+boxes[2*(i*nvars+j)+1])/2.0 << " ";
          std::cout << "]";
          std::cout << std::endl;
        };
      //    xml::print(std::cout,nsols2,"nb_sol")<< std::endl; 
      //    if ( nsols ) vctops::print(&tmp[0],nvars*nsols );
    };
    
    
    tmp.resize( nsols2 * nvars );
    //  cout <<"-----"<<endl;
    std::vector<T> vsol(nvars);
    //    Seq<vectdse<T> > res;
    for(int i=0;i<nsols2;++i)
      {
        for(int j=0; j<nvars; j++)
          vsol[j]=tmp[i*nvars+j];
        res.push_back(vsol);
      }
    //    return res;
  };

  template< class output, class MPC>
  static void solve_bernstein(output& sol, const MPC& l)
  {
     init_bernstein(l.begin(), l.end());
     std::vector<C> res; 
     int ns, nv; 
     ns=run(sol,nv);
     //     std::cout<<"system: nv="<<nv<< " ns="<< ns<<std::endl;

  }


  template< class output, class MPL, class DOM>
  static void solve(output& sol, const MPL& eqs, const DOM& dom)
  {
     init_bernstein(eqs.begin(), eqs.end());
     int nvars = sys->nvr();
     int ns, nv; 
     DOM dmns(dom);
     if((int)dom.size()>=2*nvars)
       ns=run(sol,&dmns[0],nv);
     else
       ns=run(sol,nv);

     //     std::cout<<"system: nv="<<nv<< " ns="<< ns<<std::endl;

  }
  //  ~solver_multivariate_bernstein() { delete sys; };

};

  template<class C, int M> C solver<C,ProjRd<M> >::eps = std::pow(2.0,-solver_default_precision);
  template<class C, int M> realroot::system<C>* solver<C,ProjRd<M> >::sys = 0;

//====================================================================
// template<class C,int MTH= SBD_RDRDL> struct MVBST 
// {
//   static void set_precision(unsigned prec) 
//   {
//     solver_multivariate_bernstein<C,MTH>::set_precision(prec);
//   }
// };
//====================================================================
// template< class MPC, class Domains,  class C, int MTH>
// Seq<std::vector<C> > solve( const MPC& l, MVBST<C,MTH> slv, const Domains& D, 
//                          bool verbose = false  )
// {
//   Seq<std::vector<C> > sol; 
//   solver_multivariate_bernstein<C,MTH>::solve(sol,l,D,verbose);
//   return sol;
// }
//====================================================================
#if 0
template<class C, class MPC, class Domains, int MTH>
Seq<std::vector<C> > 
solve( const MPC& l, 
       const ProjRd<MTH> &, const Domains& D , 
       bool verbose = false ) {

  typedef typename MPC::value_type  POL;
  typedef typename POL::value_type coeff;
  typedef tensor::bernstein<QQ> bpol;
  typedef QQ Rational;
  typedef polynomial<C,   with<Bernstein> > Polynomial;
  typedef polynomial<Rational, with<Bernstein> > PolynomialQQ;

  typedef polynomial<Rational, with<Sparse,DegRevLex> > PolynomialSparse;

  Seq<Rational> box(D.size()); 
  for(unsigned i=0;i<D.size();i++) 
    let::assign(box[i],D[i]);
  std::cout<<"Box: "<<box<<std::endl;

  Polynomial   p;
  PolynomialQQ q;
  Seq<Polynomial> L;
  typedef typename MPC::const_iterator iterator;
  for(iterator i=l.begin();i!=l.end();i++) {
    let::assign(q.rep(),i->rep(),box);
    Rational mx = array::max_abs(q);
    //    q/=mx;
    //    std::cout<<q<<std::endl;
    let::assign(p.rep(),q.rep());
    //    std::cout<<p.rep()<<std::endl;
    L<<p;
  }

  Seq<std::vector<C> > sol; 
  solver<C, ProjRd<MTH> >::solve_bernstein(sol,L,D,verbose);
  return sol;
}
#endif
//====================================================================
template<class C, class MPC, class Domains, int MTH>
Seq<std::vector<C> > 
solve( const MPC& l, const ProjRd<MTH> &, const Domains& D , 
                            bool verbose = false ) {
  typedef typename MPC::value_type  POL;
  typedef typename POL::value_type coeff;
  //  typedef typename kernelof<coeff>::ieee C;
    typedef QQ Rational;
        Seq<Rational> box(D.size());
        for(unsigned i=0;i<D.size();i++)
                let::assign(box[i],D[i]);
        std::cout<<"Box: "<<box<<std::endl;
        
        Seq<std::vector<C> > sol; 
        solver<C, ProjRd<MTH> >::solve_monomial(sol,l,D,verbose);
  return sol;
}
//====================================================================
} //namespace mmx
/********************************************************************/
#endif //realroot_SOLVE_mslvbst_H