solver_mv_monomial_box_rep.hpp

Go to the documentation of this file.

/********************************************************************
 *   Author(s): A. Mantzaflaris, GALAAD, INRIA
 ********************************************************************/

#ifndef _realroot_SOLVE_MV_CONFRAC_BOXREP_H
#define _realroot_SOLVE_MV_CONFRAC_BOXREP_H
/********************************************************************/

# include <realroot/ring.hpp>
# include <realroot/polynomial_tensor.hpp>



# include <realroot/polynomial.hpp>
# include <realroot/homography.hpp>
# include <realroot/sign_variation.hpp>
# include <realroot/Interval.hpp>
# include <stack>

# include <realroot/solver_ucf.hpp>
# include <realroot/solver_uv_bspline.hpp>

// # include <realroot/solver_mv_monomial_tests.hpp>

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

# define INFO 1
//# define DIVISION
# define TODOUBLE as<double>


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

template<class POL>
class box_rep
{
    typedef typename POL::scalar_t C;
    typedef Seq<POL *> system_t;
    typedef std::stack<box_rep<POL> *>  STACK;
/*----------------------------------------------*/
  
  unsigned            dim;
  homography_mv<C>    hg;        // Homography
  //homography_mv<int>    hg;        // Homography


public:

  //int td;
  Seq<POL>            S;         // System
  
  void inline update_data()
    {

      //C t;
      
      // scale coefficients
#ifdef DIVISION
       for( unsigned i=0; i<S.size(); i++ ){
          t= *std::max_element( (S[i]).begin(), (S[i]).end() );
          S[i]= S[i]/t;
        }
#endif

      // Precondition
      Seq<double> c= this->point( this->middle() ) ;
      //double *fc, *jc, *ijc;
      Seq<POL> S1= S;

      // DPOL * J= jacobian(S1);

      // jc= eval_poly_matrix( c, J, dim);   // Jacobian evaluated on c
      // ijc= new double[dim*dim];
      // linear::LUinverse(ijc, jc, dim);

      // for( unsigned i=0; i<dim; i++ ){
      //   S[i]=0;
      //   for( unsigned j=0; j<dim; j++ )
      //     S[i] +=  ijc[i*dim+j] * S1[j] ;
      // }
    }
  
  box_rep() {};
  
  box_rep(Seq<POL>& sys, homography_mv<C>& h) 
    {
      hg = h ;
      dim= h.size();
      S= sys;
      
      update_data();
    };
  
    //point of subdivision in i-direction
  box_rep( box_rep<POL> b, int i) 
    {
          dim= b.nbvar();
          hg = b.homography();

          hg.shift_hom(b.middle(i) ,i);
          hg.colapse();

          S= b.system();          
    };

    
    box_rep( Seq<POL>& sys, unsigned d)
    {
        dim=d;
        hg  = homography_mv<C>(dim) ;
            //hg  = homography_mv<int>(dim) ;
        S= sys;    
    }

    //Destroy
    ~box_rep() {};    
    
    void restrict( Seq<Interval<C> > & dom0 )
    {
          //Seq<Interval<int> > dom;
          //let::assign(dom,dom0);

        unsigned i,j;
        
        for (i=0; i!=dim; i++)
        {
        hg.shift_hom     ( dom0[i].m      , i );
        hg.contract_hom  ( dom0[i].width(), i );
        hg.reciprocal_hom( 1             , i );
        hg.shift_hom     ( 1             , i );
        hg.reciprocal_hom( 1             , i );//mirror again
        for (j=0; j!=S.size(); j++)
        {
            shift      ( S[j].rep(), dom0[i].m      , i);
            contraction( S[j].rep(), dom0[i].width(), i);
            reciprocal ( S[j].rep(),                 i);
            shift      ( S[j].rep(), C(1),           i);
            reciprocal ( S[j].rep(),                 i);//mirror again
        }
        }
        update_data();
}
    
    inline Seq<POL> system() { return S; }
    inline homography_mv<C> homography() { return this->hg; }
    inline POL      system(unsigned i) { return S[i]; }
    inline unsigned  nbvar() { return dim; }
    inline unsigned  nbpol() { return S.size(); }
    
  template<class FT>
  FT volume()
    {
      Seq<Interval<FT> > s;
      FT v(1);
      
      s= domain<FT>();
      
      for (unsigned i=0; i!=s.size(); i++ )
        v *= s[i].width();
      
      return v;
    }
  
  //template<class C>
  //inline int signof_det(DPOL **p)
  //  {
  //    return p[0];
  //  }

  template<class C>
  int** submatrix(C **matrix, int order, int i, int j)
    {
      
      C **subm;
      int p, q;     // Indexes for matrix
      int a = 0, b; // Indexes for subm
      subm = new int* [order - 1];
      
      for(p = 0; p < order; p++) {
        if(p==i) continue;  //Skip ith row
        subm[a] = new C[order - 1];
        
        b = 0;
        
        for(q = 0; q< order; q++) {
          if(q==j) continue;    //Skip jth column
          subm[a][b++] = matrix[p][q];
        }
        a++; //Increment row index
      }
      return subm;
    }

  // Compute matrix determinant, recursively
  template<class C>  
  C det(C **matrix, int order)
    {
      if(order == 1)
        return **matrix; // Matrix is of order one
      
      int i;
      C ev(0);
      for(i = 0; i < order; i++)
        ev += (i%2==0?1:-1)* matrix[i][0] * det(submatrix(matrix, order, i, 0), order - 1);
      return ev;
    }


  inline int signof(DPOL * p)
    {
      
      Interval<double> ev;
      //int dim= p->nbvar();
      
      tensor::eval( ev , p->rep() , 
      this->template domain<double>() , dim );
      
      if (ev< .1e-15) return (0);
      return (ev>0?1:-1);
    }
  
    
  //inline int tdeg() { return td ; }
    homography_mv<C> hom() { return hg; }
    POL system(const int i) { return S[i]; }
    
  template<class FT>
  Seq<Interval<FT> > domain() 
    { 
      FT l, r;
      Seq<Interval<FT> >  s ;
      
      for ( unsigned i=0; i!=dim; ++i )
      {
        //lim to 0
        if ( hg[i].b!=0 && hg[i].d!=0 )
          l= as<FT>(hg[i].b)/as<FT>(hg[i].d);
        else if ( hg[i].d==0 )
          l= 10000000;
        else if ( hg[i].b==0 )
          l= 0 ;
        else 
          l= as<FT>(hg[i].a)/as<FT>(hg[i].c) ;
        
        //lim to inf
        if ( hg[i].a!=0 && hg[i].c!=0 )
          r= as<FT>(hg[i].a)/as<FT>(hg[i].c);
        else if ( hg[i].c==0 )
          r=10000000;
        else if ( hg[i].a==0 )
          r= 0 ;
        else 
          r= as<FT>(hg[i].b)/as<FT>(hg[i].d);
        
        if ( l<=r ) s << Interval<FT>(l,r);
        else        s << Interval<FT>(r,l);
      }     
      return s; 
    }
  
  Seq<C> middle()
    {
      Seq<C> m;
      for ( unsigned i=0; i!=dim; ++i )
        m <<  as<C>( floor( as<double>(hg[i].d/hg[i].c) )) ; //floor
      return (m);
    }
  
    C middle(int i)
    {
      C t;
      t= as<C>( floor( as<double>(hg[i].d/hg[i].c) )); //floor
      //if ( t==C(0) ) t=t+1;
      return (t+1); 
    }
  
  
    template<class FT>
    Seq<FT> point(Seq<FT> t)
    {
      assert( t.size()==dim );
      Seq<FT> m;
      for ( unsigned i=0; i!=dim; ++i )
        m <<  ( as<FT>(hg[i].a)*(t[i]) + as<FT>(hg[i].b) ) / 
          ( as<FT>(hg[i].c)*(t[i]) + as<FT>(hg[i].d) ) ;
      return m;
    }

  Seq<C> subdiv_center(const unsigned& i)
    {
      Seq<C> tt;
      // Seq<double> t;
      tt.resize( dim );
      tt[i]= this->middle(i);
      
      return tt;
    }
  
  Seq<double> subdiv_point(const unsigned& i)
    {
      return  this->template point<double>(this->subdiv_center(i));
    }

  template<class FT> inline
  Seq<FT> eval(Seq<FT> t)
    {
      assert( t.size()==dim );
      Seq<FT> res;
      FT ev;
      for ( unsigned i=0; i!=S.size() ; ++i )
      {
        tensor::eval( ev , S[i].rep(), t, dim );
        res<< ev;
      }
      return res;
    }
    
    //template<class FT> 
  inline bool is_root(Seq<C> t)
    {
      assert( t.size()==dim );
      C ev;
      for ( unsigned i=0; i!=S.size() ; ++i )
      {
        tensor::eval( ev , S[i].rep() , t , dim );
        if ( ev != 0 ) 
          return false;
      }
      return true;
    }
    
    
  bool reduce_domain()
    {
      C lb;
      unsigned i;
      bool flag= false; // true iff reduction takes place
      
      Seq<C> track;
      
      //Compute lower bound and shift system
      for ( i=0; i<dim ;++i)  // for all vars
      {
        lb= l_bound(i);
        track<< lb;
        if ( lb!=0 ) 
        {
          this->shift_box(lb,i);
          update_data();
          flag= true;
        }
      }
      return flag;
    };
  
  C l_bound( const int v )
    {
      C l(0) ;
      unsigned i;
      POL m0, m1 ;
      
      for (i=0; i!=S.size(); ++i )  // for all polys
      {
        m1= maxs( & S[i] , v);
        m0= mins( & S[i] , v);
        //std::cout<< m1<< " , "<< m0 << std::endl;
        
        if ( m0( C(0) ) > 0 ) 
          l= solver<ring<C,MonomialTensor>, CFfirstFloor>::template solve<C>(m0);
//                l= as<C>( solver<ring<double,Bernstein>, Bspline>::first_root(m0) );
        else if ( m1( C(0) ) < 0 ) 
          l= solver<ring<C,MonomialTensor>, CFfirstFloor>::template solve<C>(m1);
//                l= as<C>( solver<ring<double,Bernstein>, Bspline>::first_root(m1) );
      }
      return l;
    };
  
    
    POL mins( POL * f, int v)
    {
    POL h(0,f->rep().szs()[v]-1,0) ;//var is always x0
    tensor::mins(h.rep(), f->rep(), v );
    return h;
    };
    
    
    POL maxs( POL * f, int v)
    {
    POL h(0,f->rep().szs()[v]-1,0) ;//var is always x0
    tensor::maxs(h.rep(), f->rep(), v );
    return h;
    };
    
    
    void subdivide( STACK & stck )
    {
    int i;
    Seq<int> ind(dim);
    box_rep * b;

    for (;;) 
    {
        // copy box
        b = new box_rep( *this ) ;
        
        // transform box
        for (i = 0; i < dim ; ++i)  // for all vars
        if ( ind[i] )
            b->shift_box(1,i);
        else
            b->reverse_and_shift_box(i);
        b->update_data();       

        // push box in stack
        stck.push ( b );

        // next 
        for (i = 0; i < dim ; ++i) 
        if (++ind[i] <2 ) break; else ind[i]=0 ;
        if (i == dim ) break;
    }
    };


  void safe_split(const int &v, C m=C(1) )
    {
      box_rep * b;

      b = new box_rep( *this ) ;

      // evaluate polys at x_v=m

      //check sign
    };



    void subdivide( const int &v, STACK & stck, C m=C(1) )
    {
        box_rep * b;
    
        // right box
        b = new box_rep( *this ) ;
        b->shift_box( m ,v);
        b->update_data();       
        stck.push ( b );

        // left box
        b = new box_rep( *this ) ;
            b->contract_box(m,v);
        b->reverse_and_shift_box(v);
        //b->reverse_box( v); // produces thin boxes
        b->update_data();       
        stck.push ( b );
    };
    
    void subdivide( const int &v, C & m, STACK & stck )
    {
        box_rep * b;
        box_rep tmp(*this);
    
        // left box
        b = new box_rep( tmp ) ;
        b->contract_box( m ,v);
        b->reverse_and_shift_box(v);
        b->reverse_box( v); //mirror back
        b->update_data();       
        stck.push ( b );

        // right box
        b = new box_rep( tmp ) ;
        b->shift_box( m, v);
        b->update_data();       
        stck.push ( b );
    };
    
    void shift_box(const C& t, const int & v = 0)
    {
        unsigned i;
        
        for (i = 0; i < S.size() ; ++i) //for all polys do x=x+1
        shift( S[i].rep() , t, v);
        
        //update homography
        hg.shift_hom(t,v);
    };
    
    void contract_box(const C & t, const int & v )
    {
    unsigned i;
    
    for (i = 0; i < S.size() ; ++i)
        contraction( S[i].rep() , t, v);
    
    //update homography
    hg.contract_hom(t,v);
    };
    
    // x_v = 1/x_v 
    void reverse_box(const int & v )
    {
      for (unsigned i = 0; i < S.size() ; ++i) //for all polys
        reciprocal( S[i].rep() ,   v);
    
    //update homography
    hg.reciprocal_hom(1,v);
    };
    
    // x_v = 1/(x_v +1)
    void reverse_and_shift_box(const int & v )
    {
        for (unsigned i = 0; i < S.size() ; ++i) //for all polys
        {
        reciprocal( S[i].rep() ,   v);
        shift     ( S[i].rep() ,C(1), v);
        }
        
        //update homography
        hg.reciprocal_hom(1,v);
        hg.shift_hom     (1,v);
    };


    bool inline miranda_test(const int i,const int j)
    {
        POL u,l;

        tensor::face(l, S[i], j, 0);
        tensor::face(u, S[i], j, 1);

        return  ( no_variation(l)    && 
              no_variation(u)    &&
                      (l[0]>0)!=(u[0]>0) );
    };

    // assumes a square system for now
    template<class FT>
    bool include1(DPOL * J) 
    {
        Interval<double> ev;
        unsigned i,j,c,r ;
        POL u,l;

        bool t[dim][dim];

        tensor::eval( ev , J->rep() , 
        this->template domain<double>() , dim );
        if ( ev.m*ev.M < 0 ) 
        return false;

        for (i=0; i!=dim;++i)
        for (j=0; j!=dim;++j)
            t[i][j]= miranda_test(i,j);

        c=0;
        for (i=0; i!=dim;++i)
        for (j=0; j!=dim;++j)
            if (t[i][j]==true) 
            {   
            c++; break;
            }
        if (c<dim) return false;

        c=0;
        for (i=0; i!=dim;++i)
        for (j=0; j!=dim;++j)
            if (t[j][i]==true) 
            {   
            c++; break;
            }
        if (c<dim) return false;

//      std::cout<<"INCLUDE. ev="<<ev<<", c="<<c <<std::endl;

        return true;
    };

   bool include2(DPOL * J) 
    {
          Interval<double> ev;
      
          tensor::eval( ev , J->rep() , 
                        this->domain<double>(), dim );

//            if ( this->width<double>() > 0.001 )
              if ( ev.m*ev.M < 0 ) 
            return false;

              //td= this->topological_degree_2d<double>();
        //if ( (td==-1 || td==1) )
            //{ std::cout<<"INCLUDE. ev="<<ev<<", td="<<td <<std::endl;
            //this->print();  }

        return ( 0 );
    }

   bool include3(DPOL * J) 
    {
          //evaluate df_i(B) of the box = M_i
          // check if -Ji(c)*f(c) + (I-Ji(c)*M)*B
          // is contaied in B.

          Interval<double> ev;
      
          tensor::eval( ev , J->rep() , 
                        this->domain<double>(), dim );

//            if ( this->width<double>() > 0.001 )
              if ( ev.m*ev.M < 0 ) 
            return false;

              //td= this->topological_degree_2d<double>();
        //if ( (td==-1 || td==1) )
            //{ std::cout<<"INCLUDE. ev="<<ev<<", td="<<td <<std::endl;
            //this->print();  }

        return ( 0 );
    }


    template<class FT>
    bool exclude1( Seq<POL *>& S0) 
    {
      Interval<FT> ev;
      Seq<Interval<FT> > dom;
      dom= domain<FT>();
      
      for (unsigned i=0; i!=nbpol(); ++i)
      {
        tensor::eval( ev , S0[i]->rep(), 
                      dom , dim );  
        if ( ev.m*ev.M > 0 )
        {   
          //   std::cout<<i<<" ev"<< 
          //    dom<<"\nf= "<< *(S0[1])<<
          //   "\n f([])= " <<ev<<std::endl;
          //  if (td!=0)
          // std::cout<<"!!!!!!----td: "<<td <<std::endl;
          return true;
        }   
    
      }
      return false;
    };

     template<class FT>
     inline FT width() 
        {
            unsigned i;
            FT m=this->width<FT>(i);
        
            return m;
        };
    
     template<class FT>
    FT width(unsigned & t) 
    {
        unsigned i;
        FT w;
        Seq<Interval<FT> >  s ;

        s= domain<FT>();
        w= s[0].width(); t= 0;

        for ( i=0; i!=s.size(); ++i )
        if ( s[i].width() > w) 
        { w=s[i].width() ; t=i; }

        return w;
    };
    
    POL lface(const int & i, const int & v)
    {
        POL t;

        tensor::face(t, S[i],  v , 0 );
        rename_var( t.rep() , 1-v, 0 ) ; //1-v works for 2D only

        return t;
    };

    POL rface(const int & i, const int & v)
    {
        POL tmp;
        tensor::face(tmp, S[i],  v , 1 );
        rename_var( tmp.rep() , 1-v, 0 ) ; //1-v works for 2D only

        return tmp;
    };

    void print() 
    {
        std::cout << "-------------Box---------------"    << "\n" ;
        std::cout << system() << "\n";
        for (unsigned i=0; i!=dim;++i)
        std::cout<< "("<<hg[i].a <<"x + " << hg[i].b<<")/("<<hg[i].c<<"x+ "<<hg[i].d << ")"<<std::endl; ;
            //std::cout<<"td="<<tdeg()<<std::endl;
//  std::cout << this->template domain<QQ>()<<"\n" ;
    std::cout << this->template domain<double>()<<"\n" ;
    std::cout << "-------------------------------"   << "\n" ;
    };
    
  
};// box_rep
  
} //namespace realroot
} //namespace mmx 

# undef DPOL
#endif