solver_mv_monomial.hpp

Go to the documentation of this file.

/********************************************************************
 *   Author(s): A. Mantzaflaris, GALAAD, INRIA
 ********************************************************************/
#ifndef _realroot_SOLVE_MV_CONFRAC_H
#define _realroot_SOLVE_MV_CONFRAC_H
/********************************************************************/
# include <realroot/Interval.hpp>
#include <realroot/texp_rationalof.hpp>
#include <ctime>

//# include <realroot/IEEE754.hpp>
//# include <realroot/polynom.hpp>
//# include <realroot/ring_tens.hpp>

# include <realroot/univariate_bounds.hpp>
# include <realroot/linear_doolittle.hpp>

# include <realroot/homography.hpp>
# include <realroot/solver_mv_monomial_box_rep.hpp>
# include <realroot/solver_mv_monomial_tests.hpp>

# include <stack>
//# include <realroot/solver.hpp>

//====================================================================
# define TMPL  template<class POL>
# define TMPLF template<class FT, class POL>
# define SOLVER solver<C,MCFisolate >

# define BOX   box_rep<POL>
# define N_ITER 50000


# define ALLBOXES 0
# define VERBOSE  0

# define DPOL polynomial<double,with<MonomialTensor> >


# define C_INCLUDE include1(b,J) // Miranda test.  (degens?)
//# define C_INCLUDE include3(b,J,S0) // Rump's Test

//# define C_EXCLUDE exclude1(b,S0) // Interval arithmetic
# define C_EXCLUDE exclude3(b)  //Sign Inspection
//# define C_EXCLUDE exclude1(b,S0) || exclude3(b) //both


//====================================================================
namespace mmx {
namespace realroot
{

//====================================================================
// ---------------------- Solver class  ------------------------------
//====================================================================

TMPLF
class solver_mv_monomial
{
  typedef typename POL::scalar_t C;
  typedef std::stack< BOX * >  STACK;
  
  double    eps;// precision
  //DPOL   *  jac;// jacobian matrix with double coefficients
  DPOL      * J;// jacobian matrix with double coefficients
  Seq<POL *> S0;// initial polynomials 
public:
  
  solver_mv_monomial(double e=1e-7)
    {
      J=NULL;
      eps=e;
    }; 
  
  Seq<BOX> solve_system(BOX & sys)
    {
      unsigned c=0,cand=0, i, it=0, subs=0, ver=0;
      BOX * b = NULL;
      BOX * r = NULL;
      Seq<BOX> sol;
      Seq<BOX> psol;
      bool red, out;
      STACK boxes;

      FT v(0), bx;
      bx= sys.template volume <FT>();
      if (VERBOSE) sys.print();
      
      boxes.push( new BOX(sys) );
      
      while ( !boxes.empty() )
      {
        it++;
        b = boxes.top() ;
        boxes.pop();

        /*reduce the domain */
        out= false;
        red= true;
        while ( !out && red )
        {
          if ( C_EXCLUDE )
          { 
            out=true;
            if (VERBOSE) { 
              //std::cout<<"- Excluded:"<<std::endl;
              //b->print();
            }               
            if (ALLBOXES) //FOR AXEL
            {
              v+= b->template volume<FT>();
              r = new BOX( *b ) ;
              sol << (*r); 
              free(r);
            }
            c++;
          }
          else
          {
            red= b->reduce_domain(); 
            if (VERBOSE && red) { 
              //std::cout<<"- Reduced:"<<std::endl;
              //b->print();
            }            
            //red= false;
          }
        }

        if ( out ) 
        {
          free(b);
          continue;
        }       
        /*check for root */
        if ( C_INCLUDE )
        { 
          if (VERBOSE) { 
            std::cout<<"- Solution found:"<<std::endl;
            b->print();
          }
          
          ver++ ;
          sol << (*b);
          free(b);
          continue; 
        }
    
        /*Subdivide*/
        if ( it > N_ITER )
        {
          cand++;
          std::cout<<"MAX iters"<<" ("<<N_ITER<<") "
                   <<"reached!" << std::endl;
          b->print();
          sol << (*b);
          free(b);
        }
        else
        {
          if ( b->template width<double>(i) > eps )
          { 
            subs++;
          
            if (check_root( b->subdiv_point(i),eps) )
            {
              psol <<BOX( *b, i);
              //free(b);
              //continue; 
              //sol[sol.size()-1].print();
              ver++;
            }

            b->subdivide (i,boxes, b->middle(i) );
            //b->subdivide (i,boxes);
            //b->subdivide (i,boxes,   1 );
            //b->subdivide (boxes);
            free(b); 
          }
          else 
          {
            if ( C_EXCLUDE  ){
              if (ALLBOXES) sol << (*b); //FOR AXEL
            }else
            {   
              cand++;
              sol << (*b);
              if (VERBOSE) { 
                std::cout<<"- EPS reached:"<<std::endl;
                //b->print();
              }
            }
            free(b);
          }
        }       
      }/*while*/
      
      unsigned j=0;
      //if (0)
      if ( !ALLBOXES && S0.size()==2)
      while (j<sol.size())
      {
        if ( exclude2( &sol[j],J) )
        {
          sol.erase(j);
          cand--;
        }
        else 
        { //std::cout  <<"td="<<sol[j].td<<std::endl;
          ++j; 
        }
      }

      if (VERBOSE) {
      std::cout<< "Iterations= "<< it           <<std::endl;    
      std::cout<< "Excluded  = "<< c            <<std::endl;    
      std::cout<< "Verified  = "<< ver          <<std::endl;    
      std::cout<< "Subdivs   = "<< subs         <<std::endl;
      if (ALLBOXES)
        std::cout<< "Reduced volume= "<< 
          as<double>(  100*(bx-v)/bx )<< "\%"     <<std::endl;  
      std::cout<< "#stack=     "<< cand         <<std::endl;    
      }
      sol<< psol;

      return sol;
    }
  
  Seq< std::vector<Interval<FT> > > isolate(Seq<POL>& p, unsigned &d)
    {
      homography_mv<C> h(d);                
      BOX sys= BOX(p,h);
      
      /*Global data*/
      free(J);
      for (unsigned i=0;i<S0.size();++i) 
        free(S0[i]);
      S0=Seq<POL *>();
      for (unsigned i=0;i<sys.nbpol();++i) 
        S0 << new POL( sys.system(i) );     
      J= jacobian(S0);
      /*end Global data*/
      
      Seq< BOX> s(solve_system(sys) );
      
      Seq< std::vector<Interval<FT> > > a;
      
      for (unsigned i=0; i<s.size(); ++i)
        a << s[i].template domain<FT>().rep();
               
//      free(J);
//      for (unsigned i=0;i<S0.size();++i) 
//        free(S0[i]);  
      return a;
    }
  
  Seq< std::vector<Interval<FT> > > isolate(Seq<POL>& p, Seq<Interval<C> > & dom)
    {
          BOX * sys= new BOX(p, dom.size() );
          
          /*Global data*/
          free(J);
          for (unsigned i=0;i<S0.size();++i) 
            free(S0[i]);
          S0=Seq<POL *>();
          for (unsigned i=0;i<sys->nbpol();++i) 
            S0 << new POL( sys->system(i) );        
          J= jacobian(S0);
          /*end Global data*/

          sys->restrict(dom);
          Seq< BOX> s(solve_system( *sys) );
      
          Seq< std::vector<Interval<FT> > > a;

          for (unsigned i=0; i<s.size(); ++i)
            a << s[i].template domain<FT>().rep();
          
//          free(J);
//          for (unsigned i=0;i<S0.size();++i) 
//            free(S0[i]);
          return a;
    }
  
  
  Seq< std::vector<FT> > approximate(Seq<POL>& p, unsigned &d)
    {
      homography_mv<C> h(d);                
      BOX sys= BOX(p,h);
      
      /*Global data*/
      free(J);
      for (unsigned i=0;i<S0.size();++i) 
        free(S0[i]);
      S0=Seq<POL *>();
      for (unsigned i=0;i<sys->nbpol();++i) 
        S0 << new POL( sys->system(i) );        
      J= jacobian(S0);
      /*end Global data*/
      
      Seq<FT> t;
      Seq< BOX> s(solve_system(sys) );
      
      Seq< std::vector<FT> > a;
      
      // ... inf ?
      
//      free(J);
//      for (unsigned i=0;i<S0.size();++i) 
//        free(S0[i]);
      return a;
    }
  
  Seq< std::vector<FT> > approximate(Seq<POL>& p, Seq<Interval<C> > & dom)
    {
      BOX * sys= new BOX(p, dom.size() );
      
       /*Global data*/
       // free(J); //crash
      for (unsigned i=0;i<S0.size();++i) 
        free(S0[i]);
      S0=Seq<POL *>();
      for (unsigned i=0;i<sys->nbpol();++i) 
        S0 << new POL( sys->system(i) );        



      

      J= jacobian(S0);
      /*end Global data*/
      
      sys->restrict(dom);
      unsigned v; //, dim(dom.size() );
      Seq<FT> d;
      BOX * l, * b;
      Seq<C> t;
      //FT ev;
      Seq<BOX> s(solve_system(*sys) );
      
      Seq< std::vector<FT> > a;
      
      for (unsigned i=0;i<s.size(); ++i )
      {
        b= new BOX(s[i]);
    
        while ( b->template width<FT>(v) > eps )
        {
          t= b->middle();
          if ( b->is_root(t) ) 
          {
            d= b->template point<FT>(t);
            a << d.rep();
          }
          
          l= new BOX( *b ) ;
          l->shift_box( t[v] , v );
          
          if ( C_INCLUDE ) 
          {   
            free(b);
            b=l; 
            continue;
          }
          else 
          {   
            free(l);
            b->contract_box(t[i],v);
            b->reverse_and_shift_box(v);
            b->reverse_box(v);
          }
        }
        d= b->template domain<FT>();
        a << d.rep();
        free(b);
      }
      return a;
    }
  

  bool check_root(const Seq<double> t, const double eps) 
    {
      DPOL p(0);
      double ev;
      for (unsigned j=0; j!=S0.size(); ++j)
      {
        
        let::assign(p, *S0[j]);
        tensor::eval( ev , p.rep() , 
                      t , t.size() );

        //std::cout<<"check: "<< ev<<std::endl;
        if ( abs(ev) > 1e-10 )
          return false;
      } 
      std::cout<<"Root on split: "<< t  <<std::endl;
      return true;

    };

  
};// solver_mv_monomial
  
} //namespace realroot
  
//====================================================================
// --------------------------- INTERFACE -----------------------------
//====================================================================

    struct MCFisolate{};
    struct MCFapproximate{};

//--------------------------------------------------------------------

    template<class C>
    struct solver<C, MCFisolate > {

    typedef C   base_t;

    template<class FT, class POL> 
    static Seq< std::vector<Interval<FT> > >
    solve( Seq<POL>& p, Seq< Interval<base_t> > & dom);

    template<class FT, class POL>
    static Seq<std::vector<Interval<FT> > >
    solve( Seq<POL>& p);
};

template<class C> 
template<class FT, class POL> 
Seq< std::vector<Interval<FT> > >
SOLVER::solve( Seq<POL>& p, Seq<Interval<base_t> > & dom)
{

    realroot::solver_mv_monomial<FT,POL> slv(1e-3);
    return  slv.isolate(p, dom);
}

template<class C> 
template<class FT,class POL> 
Seq<std::vector<Interval<FT> > >
SOLVER::solve(Seq<POL>& p) {
    
    unsigned d(p[0].nbvar() );
    realroot::solver_mv_monomial<FT,POL>  slv(1e-3);
    return slv.isolate(p,d);
}


//--------------------------------------------------------------------
template<class C>
struct solver<C, MCFapproximate > {

    typedef C   base_t;
    template<class FT, class POL> static Seq< std::vector<FT> >
    solve( Seq<POL>& p, Seq<Interval<base_t> > & dom);

    template<class FT, class POL> static Seq< std::vector<FT> >
    solve( Seq<POL>& p);
};

template<class Ring> 
template<class FT, class POL> 
Seq< std::vector<FT> >
solver<Ring,MCFapproximate >::solve( Seq<POL>& p, Seq<Interval<base_t> > & dom)
{
    realroot::solver_mv_monomial<FT,POL> slv(1e-3);
    return slv.approximate(p, dom );
}

template<class C> 
template<class FT, class POL> 
Seq< std::vector<FT> >
solver<C,MCFapproximate >::solve( Seq<POL>& p)
{
    
    unsigned d( p[0].nbvar() );
    realroot::solver_mv_monomial<FT,POL> slv(1e-3);
    
    return slv.approximate(p,d);
}


} //namespace mmx 
//====================================================================
# undef TMPL
# undef SOLVER
# undef DPOL
//====================================================================
#endif