realroot: include/realroot/solver_mv_bernstein.hpp Source File

00001 /********************************************************************
00002  *   This file is part of the source code of the realroot library.
00003  *   Author(s): J.P. Pavone, GALAAD, INRIA
00004  *   $Id: realroot.h,v 1.1 2005/07/11 10:03:51 jppavone Exp $
00005  ********************************************************************/
00006 #ifndef realroot_SOLVE_mslvbst_H
00007 #define realroot_SOLVE_mslvbst_H
00008 /********************************************************************/
00053 #include <set>
00054 #include <vector>
00055 #include <math.h>
00056 #include <realroot/texp_kernelof.hpp>
00057 #include <realroot/solver.hpp>
00058 
00059 #include <realroot/tensor_bernstein.hpp>
00060 
00061 #define REP(x) (x).rep()
00062 
00063 #include <realroot/GMP.hpp>
00064 #include <realroot/polynomial.hpp>
00065 #include <realroot/system_method.hpp>
00066 #include <realroot/system_support.hpp>
00067 #include <realroot/Seq.hpp>
00068 
00069 
00070 //====================================================================
00071 namespace mmx {
00072 //====================================================================
00073 namespace realroot
00074 {
00075 #define SBD_RD       0
00076 #define SBD_RDL      1
00077 #define SBD_RDS      2
00078 #define SBD_RDRDL    3
00079 #define SBD_RDRDLRDS 4
00080 #define SBD_RDRDS    5
00081 #define SBD_RDLRDS   6
00082 #define SBD_RDSRDL   7
00083 #ifndef SBD_MTH
00084 #define SBD_MTH SBD_RDRDL
00085 #endif
00086   
00087   template< class system, int mth > /* selection de methodes */
00088   struct select_mth
00089   { typedef realroot::method<system,strgy::newton<system>,rdslv::parallel,sbdrl::parametric> result_t; };
00090   template< class system > 
00091   struct select_mth<system,SBD_RD>
00092   { typedef realroot::method<system,strgy::simple<system>,rdslv::parallel,sbdrl::parametric> result_t; };
00093   template< class system > 
00094   struct select_mth<system,SBD_RDL>
00095   { typedef realroot::method<system,strgy::newton<system>,rdslv::parallel,sbdrl::parametric> result_t; };
00096   template< class system>
00097   struct select_mth<system,SBD_RDRDL>
00098   { typedef realroot::method<system,strgy::newton< system, strgy::simple<system> >,rdslv::parallel,sbdrl::parametric> 
00099     result_t; };
00100   template< class system >
00101   struct select_mth<system,SBD_RDRDLRDS>
00102   { typedef realroot::method<system,strgy::cvmatrix< system, strgy::newton<system,strgy::simple<system> > >,rdslv::parallel,sbdrl::parametric> result_t; };
00103   template< class system > 
00104     struct select_mth<system,SBD_RDS>
00105   { typedef realroot::method<system,strgy::cvmatrix<system>,rdslv::parallel,sbdrl::parametric> result_t; };
00106 
00107   template< class system > 
00108     struct select_mth<system,SBD_RDRDS>
00109   { typedef realroot::method<system,strgy::cvmatrix<system, strgy::simple<system> >,rdslv::parallel,sbdrl::parametric> result_t; };
00110   
00111   template< class system > 
00112     struct select_mth<system,SBD_RDLRDS>
00113   { typedef realroot::method<system,strgy::cvmatrix<system, strgy::newton<system> >,rdslv::parallel,sbdrl::parametric> result_t; };
00114 
00115   template< class system > 
00116     struct select_mth<system,SBD_RDSRDL>
00117   { typedef realroot::method<system,strgy::cvmatrix<system, strgy::newton<system> >,rdslv::parallel,sbdrl::parametric> result_t; };
00118   
00119   template< class T, class eenv_t, class MPOL >
00120   void fill_data( T * data, eenv_t * env, const MPOL& p, binomials<T>& benv  )
00121   {
00122     std::fill(data,data+env->sz(),0); 
00123     mpolfill(data,p.rep(),*env);
00124     convertm2b (data,*env,benv);
00125   };
00126   template< class monom >
00127   void scan_monom_env ( std::set< int >& env, const monom& m )
00128   { for ( unsigned i = 0; i < (m.rep()).size(); i++ ) if ( (m.rep())[i] != 0 ) env.insert(i); };
00129 }
00130 
00131   template<int MTH> struct ProjRd {};
00132   
00133 //====================================================================
00147 template<class C, int MTH>
00148 struct solver< C, ProjRd<MTH> >
00149 {
00150   typedef C coeff_t;
00151   typedef Seq<std::vector<C> >  Solutions;
00152   typedef realroot::system<C>   system_t;
00153   
00154   static C eps;
00155   static system_t* sys;
00156 
00158   // solver_multivariate_bernstein( const C& precision = (C)1e-6 ): eps(precision), sys(0){};
00159 
00160   
00161   static void set_precision(int prec) { using std::pow; eps = pow(2.0,-prec); }
00162 
00163   template< class MPDI >  static
00164   void init_bernstein( MPDI  first, MPDI last)
00165   {
00166 
00167     int    ne = 0;
00168     MPDI eptr;
00169     for ( eptr = first; eptr != last; eptr ++ ) ne++;
00170     int nvrs[ne];
00171     const int *szs[ne];
00172     int e;
00173     //    int mxvr = 0;
00174     std::set<int> gvrs;
00175     for ( e = 0, eptr = first; eptr != last; eptr ++, e++  ) { 
00176       nvrs[e] = eptr->rep().nbvar(); 
00177       szs[e]  = eptr->rep().szs();
00178       const int * tmp = eptr->rep().vrs();
00179       for ( int v = 0; v < eptr->rep().nvr(); v ++ ) gvrs.insert(tmp[v]);
00180     };
00181     int mxv = *(gvrs.rend());
00182     int lnames[ mxv+1 ];
00183     int nv = 0;
00184     for ( std::set<int>::const_iterator it = gvrs.begin(); it != gvrs.end(); lnames[*it++] = nv++ ) ;
00185     int * vrdat = new int[ nv * ne ];
00186     int * vrptr = vrdat;
00187     const int *  vrs[ne];
00188     e = 0;
00189     for ( eptr = first; eptr != last; vrptr+= eptr->rep().nbvar(), e++, eptr ++ ) {
00190       vrs[e] = vrptr;
00191       const int * tmp = eptr->rep().vrs();
00192       for ( int v = 0; v < eptr->rep().nvr();  v ++ )
00193         vrptr[v] = lnames[ tmp[v] ];    
00194     };
00195     //    std::cout << ne << " " << nv << std::endl;
00196     assert(ne>=nv);
00197     sys = new realroot::system<C>( ne, nv, nvrs, (const int**)vrs, (const int**)szs, eps );
00198     e=0; for ( MPDI eptr = first; eptr != last; add_equation( REP(*eptr), e++  ), eptr++ ) ;
00199     delete[] vrdat;
00200   };
00201   
00202   static void add_equation( const tensor::bernstein<C>& data, int e )
00203   {
00204     sys->set_bernstein2( data.begin(), e );
00205   };
00206   
00207   static void add_equation( const tensor::monomials<C>& data, int e ) 
00208   {
00209     sys->set_monomial2( data.begin(), e );
00210   };
00211   
00212   static int run( std::vector<C> & R, int& nv )
00213   {
00214     R.resize(0);
00215     typename realroot::select_mth<system_t,MTH>::result_t __mth;
00216     //int e = 0;for ( MPDI eptr = first; eptr != last; add( eptr->rep(), e++  ), eptr++ );
00217     //std::vector<C> tmp;
00218     realroot::system_ctrl< std::vector< C > >  intf(R);
00219     int nvars = sys->nvr();
00220     __mth.launch(sys,intf,(C*)0);
00221 
00222     int nsols = R.size() / (2*nvars);
00223     int nsols2 = realroot::clean_result( &R[0], nvars, nsols,(C)1e-3 );
00224     R.resize(nsols2 * nvars) ;
00225     int ns = nsols2;
00226     nv = nvars;
00227     return ns;
00228   };
00229   
00230   static int run( std::vector<C> & sol, C* dom, int& nv )
00231   {
00232     sol.resize(0);
00233     typename realroot::select_mth<system_t,MTH>::result_t __mth;
00234     realroot::system_ctrl< std::vector< C > >  intf(sol);
00235     int nvars = sys->nvr();
00236     //    C* dmns = &dom[0];
00237     // if ( (int)dom.size() >= 2*nvars )
00238 //       {
00239 //      dmns = new C[ 2*nvars ];
00240 //      for ( int i = 0; i < 2*nvars; i ++ ) let::assign(dmns[i],dom[i]);
00241 //       };
00242     __mth.launch(sys,intf,dom);
00243 
00244     int nsols = sol.size() / (2*nvars);
00245     int nsols2 = realroot::clean_result( &sol[0], nvars, nsols,(C)1e-3 );
00246     sol.resize(nsols2 * nvars) ;
00247     int ns = nsols2;
00248     nv = nvars;
00249     return ns;
00250   };
00251 
00252   //--------------------------------------------------------------------
00261   template< class output, class MPC, class Domain >
00262   static void solve_monomial( output& res, 
00263                               const MPC& l, 
00264                               const Domain  & dom,
00265                               bool verbose = false )
00266   {
00267     typedef C T;
00268     typedef typename MPC::value_type          mpoly_t;
00269     typedef MPC                               plist_t;
00270     typedef typename plist_t::const_iterator  pptr_t;
00271     typedef typename mpoly_t::const_iterator  mptr_t;
00272     //    typedef typename mpoly_t::monom_t         monom_t;
00273     //    typedef solver_multivariate_bernstein<C, MTH>           mth_t;
00274     typedef int                               sz_t;
00275     sz_t  nvars = 0, neqs, * nvrs, ** vrs, ** szs; 
00276     int e;
00277     pptr_t cpl; 
00278     // recuperation de l'environement global
00279     for ( neqs = 0, cpl = l.begin(); cpl != l.end(); cpl++, neqs++  )
00280       for ( mptr_t mn = cpl->begin(); mn != cpl->end(); mn++ )
00281         nvars = std::max(nvars,mn->nvars()+1);
00282     // recuperation de l'environement de chaque equation
00283     std::set<int> * seenvs = new std::set< int > [ neqs ];
00284     for ( e = 0, cpl = l.begin(); e < neqs; e++, cpl++  )
00285       for ( mptr_t mn = cpl->begin(); mn != cpl->end(); mn++ ) 
00286         realroot::scan_monom_env(seenvs[e],*mn);
00287     
00288     //  assert(sizeof(sz_t)==sizeof(sz_t*));
00289     sz_t  s = 0;
00290     for ( e = 0; e < neqs; s += seenvs[e].size(), e++ ) ;
00291     nvrs = new sz_t [neqs+2*s];
00292     vrs  = new sz_t*[2*neqs];
00293     szs  = vrs + neqs;
00294     sz_t * base = nvrs + neqs;
00295     std::set<int>::iterator seenvit;  
00296     for ( e = 0, cpl = l.begin(); e < neqs; base += 2*nvrs[e++], cpl++ )
00297       {
00298         nvrs[e] = seenvs[e].size();
00299         vrs [e] = base;
00300         szs [e] = base + nvrs[e];
00301         int k = 0;
00302         for ( seenvit = seenvs[e].begin(); k < nvrs[e]; seenvit++, k++  ) 
00303         {
00304           vrs[e][k]  = *seenvit;
00305           szs[e][k]  = degree( *cpl, *seenvit ) + 1;
00306         };
00307       };
00308     
00309     std::vector<C> result;
00310     C * dmns = 0;
00311     if ( (sz_t)dom.size() >= 2*nvars )
00312       {
00313         dmns = new C[ 2*nvars ];
00314         for ( int i = 0; i < 2*nvars; i ++ ) let::assign(dmns[i],dom[i]);
00315       };
00316     
00317     //  chrono tm;
00318     std::vector<T> tmp;
00319     
00320     {  /* remplissage et appel */
00321       using namespace realroot;
00322       binomials<T> bznv;
00323       typedef  realroot::system<T> system_t;
00324       system_t sys( neqs, nvars, nvrs, (const int**)vrs, (const int**)szs, (T)eps);
00325       for ( e = 0, cpl = l.begin(); e < neqs; e ++, cpl++  )
00326         realroot::fill_data( sys.data(e), sys.env(e), *cpl, bznv );
00327       realroot::system_ctrl< std::vector< T > >  intf(tmp);
00328       //    tm.start();
00329       typename realroot::select_mth<system_t,MTH>::result_t __mth;
00330       __mth.launch(&sys,intf,dmns);
00331       //tm.stop();
00332       if ( verbose )
00333         {
00334           std::cout << " iterations  = " << __mth.m_niter << std::endl;
00335           std::cout << "subdivisions = "<< __mth.m_nsbd  << std::endl;
00336         };
00337     };
00338     
00339     if ( !verbose && dmns ) delete[] dmns;
00340     delete[] seenvs;
00341     delete[] nvrs;
00342     delete[] vrs;
00343     if ( verbose )
00344       {
00345         std::cout << "<domain>";
00346         if  ( dmns )
00347           {
00348             std::cout << std::endl;
00349             for ( int i = 0; i < nvars; i ++ )
00350               std::cout << "\t" 
00351                         << "x" << i//monom_t::var_of_index(i) 
00352                         << " = [ " << dmns[2*i] 
00353                         << "," 
00354                         << dmns[2*i+1] 
00355                         << " ]";
00356           }
00357         else std::cout << "[0,1]^" << nvars;
00358         std::cout << "</domain>\n";
00359         //      std::cout << "\treductions   = " << bsys.nrest() << std::endl;
00360         //      std::cout << "\titerations   = " << bsys.niter() << std::endl;
00361         //      std::cout << "<solver_time>"; tm.write(std::cout); 
00362         //      std::cout <<"</solver_time>"<< std::endl;       
00363         if ( dmns ) delete[] dmns;
00364       };
00365     
00366     unsigned nsols = tmp.size() / (2*nvars);
00367     T boxes[nsols*2*nvars];
00368     std::copy(&tmp[0], &tmp[0]+2*nvars*nsols,boxes);
00369     int nsols2 = realroot::clean_result( &tmp[0], nvars, nsols,(T)1e-3 );
00370     if ( verbose ) { 
00371       //    xml::print(std::cout,nsols,"nb_box")<< std::endl;
00372       for ( unsigned i = 0; i < nsols; i ++ )
00373         {
00374           std::cout << "\t[ ";
00375           for (int j = 0; j < nvars; j ++ )
00376             std::cout << (boxes[2*(i*nvars+j)]+boxes[2*(i*nvars+j)+1])/2.0 << " ";
00377           std::cout << "]";
00378           std::cout << std::endl;
00379         };
00380       //    xml::print(std::cout,nsols2,"nb_sol")<< std::endl; 
00381       //    if ( nsols ) vctops::print(&tmp[0],nvars*nsols );
00382     };
00383     
00384     
00385     tmp.resize( nsols2 * nvars );
00386     //  cout <<"-----"<<endl;
00387     std::vector<T> vsol(nvars);
00388     //    Seq<vectdse<T> > res;
00389     for(int i=0;i<nsols2;++i)
00390       {
00391         for(int j=0; j<nvars; j++)
00392           vsol[j]=tmp[i*nvars+j];
00393         res.push_back(vsol);
00394       }
00395     //    return res;
00396   };
00397 
00398   template< class output, class MPC>
00399   static void solve_bernstein(output& sol, const MPC& l)
00400   {
00401      init_bernstein(l.begin(), l.end());
00402      std::vector<C> res; 
00403      int nv; //ns;
00404      //ns=
00405        run(sol,nv);
00406      //     std::cout<<"system: nv="<<nv<< " ns="<< ns<<std::endl;
00407 
00408   }
00409 
00410 
00411   template< class output, class MPL, class DOM>
00412   static void solve(output& sol, const MPL& eqs, const DOM& dom)
00413   {
00414      init_bernstein(eqs.begin(), eqs.end());
00415      int nvars = sys->nvr();
00416      int ns, nv; 
00417      DOM dmns(dom);
00418      if((int)dom.size()>=2*nvars)
00419        ns=run(sol,&dmns[0],nv);
00420      else
00421        ns=run(sol,nv);
00422 
00423      //     std::cout<<"system: nv="<<nv<< " ns="<< ns<<std::endl;
00424 
00425   }
00426   //  ~solver_multivariate_bernstein() { delete sys; };
00427 
00428 
00438   template<class MPC, class Domain>
00439   Seq<std::vector<C> > 
00440   static solve( const MPC& l, 
00441                 const Domain& D,
00442                 bool verbose = false ) {
00443 
00444     typedef typename MPC::value_type  POL;
00445     typedef typename POL::value_type coeff;
00446     typedef QQ Rational;
00447     Seq<Rational> box(D.size());
00448     for(unsigned i=0;i<D.size();i++) let::assign(box[i],D[i]);
00449     std::cout<<"Box: "<<box<<std::endl;
00450   
00451     Seq<std::vector<C> > sol; 
00452     solver<C, ProjRd<MTH> >::solve_monomial(sol,l,D,verbose);
00453     return sol;
00454   }
00455 };
00456 
00457   template<class C, int M> C solver<C,ProjRd<M> >::eps = std::pow(2.0,-solver_default_precision);
00458   template<class C, int M> realroot::system<C>* solver<C,ProjRd<M> >::sys = 0;
00459 
00460 //====================================================================
00461 // template<class C,int MTH= SBD_RDRDL> struct MVBST 
00462 // {
00463 //   static void set_precision(unsigned prec) 
00464 //   {
00465 //     solver_multivariate_bernstein<C,MTH>::set_precision(prec);
00466 //   }
00467 // };
00468 //====================================================================
00469 // template< class MPC, class Domain ,  class C, int MTH>
00470 // Seq<std::vector<C> > solve( const MPC& l, MVBST<C,MTH> slv, const Domain & D,
00471 //                          bool verbose = false  )
00472 // {
00473 //   Seq<std::vector<C> > sol; 
00474 //   solver_multivariate_bernstein<C,MTH>::solve(sol,l,D,verbose);
00475 //   return sol;
00476 // }
00477 //====================================================================
00478 #if 0
00479 template<class C, class MPC, class Domain , int MTH>
00480 Seq<std::vector<C> > 
00481 solve( const MPC& l, 
00482        const ProjRd<MTH> &, 
00483        const Domain & D ,
00484        bool verbose = false ) {
00485 
00486   typedef typename MPC::value_type  POL;
00487   typedef typename POL::value_type coeff;
00488   typedef tensor::bernstein<QQ> bpol;
00489   typedef QQ Rational;
00490   typedef polynomial<C,   with<Bernstein> > Polynomial;
00491   typedef polynomial<Rational, with<Bernstein> > PolynomialQQ;
00492 
00493   typedef polynomial<Rational, with<Sparse,DegRevLex> > PolynomialSparse;
00494 
00495   Seq<Rational> box(D.size()); 
00496   for(unsigned i=0;i<D.size();i++) 
00497     let::assign(box[i],D[i]);
00498   std::cout<<"Box: "<<box<<std::endl;
00499 
00500   Polynomial   p;
00501   PolynomialQQ q;
00502   Seq<Polynomial> L;
00503   typedef typename MPC::const_iterator iterator;
00504   for(iterator i=l.begin();i!=l.end();i++) {
00505     let::assign(q.rep(),i->rep(),box);
00506     Rational mx = array::max_abs(q);
00507     //    q/=mx;
00508     //    std::cout<<q<<std::endl;
00509     let::assign(p.rep(),q.rep());
00510     //    std::cout<<p.rep()<<std::endl;
00511     L<<p;
00512   }
00513 
00514   Seq<std::vector<C> > sol; 
00515   solver<C, ProjRd<MTH> >::solve_bernstein(sol,L,D,verbose);
00516   return sol;
00517 }
00518 #endif
00519 //====================================================================
00530 template<class C, class MPC, class Domain , int MTH>
00531 Seq<std::vector<C> > 
00532 solve( const MPC& l, const ProjRd<MTH> &, const Domain & D,
00533        bool verbose = false ) {
00534   typedef typename MPC::value_type  POL;
00535   typedef typename POL::value_type coeff;
00536   typedef QQ Rational;
00537   Seq<Rational> box(D.size());
00538   for(unsigned i=0;i<D.size();i++) let::assign(box[i],D[i]);
00539   std::cout<<"Box: "<<box<<std::endl;
00540   
00541   Seq<std::vector<C> > sol; 
00542   solver<C, ProjRd<MTH> >::solve_monomial(sol,l,D,verbose);
00543   return sol;
00544 }
00545 //====================================================================
00546 } //namespace mmx
00547 /********************************************************************/
00548 #endif //realroot_SOLVE_mslvbst_H