system_dstack.h

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: dstack.h,v 1.1 2005/07/11 10:03:57 jppavone Exp $
 ********************************************************************/
#ifndef realroot_SOLVE_SBDSLV_SYSTEM_DOMAIN_HPP
#define realroot_SOLVE_SBDSLV_SYSTEM_DOMAIN_HPP
/********************************************************************/
#include <stdlib.h>
#include <realroot/system_epsilon.hpp>
#include <realroot/rounding_mode.hpp>
#include <realroot/system_nstack.h>
#include <realroot/Interval.hpp>
//--------------------------------------------------------------------
namespace mmx {
//--------------------------------------------------------------------
namespace realroot
{
  template< class creal_t >
  struct dstack : memory::nstack< typename Intervals<creal_t>::rdw_t  >
  {
    typedef   unsigned sz_t;
    typedef   typename Intervals<creal_t>::rdw_t interval_t;
    method_base * m_error;
    void receiver( method_base * mb ) { m_error = mb; };
    void error( const char * msg ) { m_error->error( msg ); };
    
    /* domaine précedent */

    inline interval_t * previous() { return this->top()+this->nvr(); };
    /* domaine courant   */
    inline interval_t * current () { return this->top(); };
    inline const interval_t * current() const { return this->top(); };
    sz_t    m_nvr;   /* nombre de variables   */
    sz_t    nvr() const { return m_nvr; };
    creal_t m_peps  ;/* précision à atteindre */
    const creal_t& peps() const { return m_peps; };
    creal_t m_spleps;/* seuil d'isolation      */
    creal_t m_rdprm; /* parametre de réduction */
    const creal_t& rdscr_prm() const { return m_rdprm; };
    sz_t    m_large; /* variable correspondant au plus grand intervalle */
    sz_t    m_small; /* variable correspondant au plus petit intervalle */
    
    inline sz_t    largest() { return m_large; };
    inline sz_t    smallest() { return m_small; };
    inline static creal_t size( const interval_t& i ) { return i.upper()-i.lower(); };
    inline creal_t prc()  {  return size(this->top()[ largest() ]); };

    /* initialisation de la pile avec le domaine de départ donnés par inits (2*nvr() valeurs) */
    void init( creal_t * inits )
    {
      this->popall(), this->push();
      if ( inits )
    for ( sz_t v = 0; v < this->nvr(); v ++ )
      new (this->top()+v)  interval_t(inits[2*v],inits[2*v+1]);
      else  
    std::fill( this->top(), this->top() + nvr(), interval_t((creal_t)0.0,(creal_t)1.0) );
    };
    
    void reset() 
    {
      const interval_t * curr(current());
      m_small = m_large = nvr()-1;
      creal_t mn(curr[nvr()-1].size());
      creal_t mx(curr[nvr()-1].size());
      for ( sz_t v = 0; v < nvr(); v += 2  )
    {
      creal_t a = curr[v].size();
      creal_t b = curr[v+1].size();
      if ( a > b  ) 
        {
          if ( b < mn ) { m_small = v + 1; mn = b; };
          if ( mx < a ) { m_large = v;     mx = a; };
        }
      else
        {
          if ( a < mn ) { m_small = v;     mn = a; };
          if ( mx < b ) { m_large = v + 1; mx = b; };
        };
    };
    };
    
    /* fixe la précision que l'on souhaite atteindre et la précision de l'isolation */
    void set_precision( const creal_t& peps, const creal_t& spleps = 0 )
    {
      m_peps   = peps;
      m_spleps = std::max(std::max(peps,spleps),(creal_t)1e-6);
    };
    
    dstack( sz_t _nvr_, const creal_t& peps, const creal_t& spleps ): 
      memory::nstack< interval_t >(_nvr_*2), 
      m_nvr(_nvr_)
    { 
      assert(_nvr_);
      m_nvr = _nvr_;
      //      m_rdprm = 0.9;//
      m_rdprm = (nvr()-(creal_t)0.5)/(nvr());
      /*      m_rdprm = 0.99; */
      //      std::cout << m_rdprm << std::endl;
      set_precision(peps,spleps);};

    /* garder en mémoire le domaine */ 
    inline void setmem() { std::copy(this->top(),this->top()+nvr(),this->top()+nvr()); };
    inline void mpush() {
      interval_t * prv = previous(); 
      this->push();
      std::copy(prv,prv+nvr(),this->top()+nvr());
    };
    inline void cpush() { 
      interval_t * prv = current(); 
      this->push();
      std::copy(prv,prv+2*nvr(),this->top());
    };
    
    /* échantillonage aléatoire du domaine  */
    void rndsupport( creal_t * dst, int nsmp = 1 )
    {
      interval_t * dm = this->top();
      creal_t * ptr;
      for ( unsigned v = 0; v < nvr(); v ++ )
    {
      creal_t a = dm[v].upper();
      creal_t b = dm[v].lower()-a;
      for ( ptr = dst + v; ptr < dst + nsmp * nvr(); ptr += nvr() )
        *ptr = a + ((creal_t)rand())/((creal_t)RAND_MAX)*b;
    };
    };
    

    static inline 
    creal_t _linterp_dw( const interval_t& intvl, const creal_t& t )
    {
      using namespace numerics::rdw;
      return dwadd(dwmul(upsub(creal_t(1),t),intvl.lower()),
           dwmul(t,intvl.upper())) ;
    };

    static inline 
    creal_t _linterp_up( const interval_t& intvl, const creal_t& t )
    { 
      using namespace numerics::rdw;
      return upadd(upmul(dwsub(creal_t(1),t),intvl.lower()),
           upmul(t,intvl.upper()));
    };
        
    /* mise a jour domaine local (src) vers global (dst) */
    void update( interval_t * dst, interval_t * src, interval_t * locals)
      //            ,bool debug = false )
      {
    numerics::rounding< creal_t > rnd( numerics::rnd_dw() );
    for ( sz_t v = 0; v < nvr(); v ++ )
      new(dst+v) interval_t( std::max(_linterp_dw(src[v],locals[v].lower()),
                     src[v].lower()),
                 std::min(_linterp_up(src[v],locals[v].upper()),
                    src[v].upper()) );
      };
    
    creal_t rdscr( interval_t const * const , interval_t const * const locals )
    {
      sz_t v;
      creal_t d0;
      for ( v = 0, d0=0.0 ; v < nvr(); d0 += size(locals[v]), v ++ ) ;
      d0/=nvr();
      //      std::cout << d0 << "/" << rdscr_prm() << std::endl;
      return d0;
    };
    
    
    /*  empilage des réductions */
    bool reduce( std::vector<interval_t> * ri, unsigned ) {
    unsigned nvrs = nvr();
    unsigned nsdm = 1, rs[nvrs];
    /* le nombre de sous-domaines (nsdm) présents dans la réduction est     */
    /* égal au produit du nombre d'intervalles trouvés pour chaque variable */
    for ( sz_t v = 0; v < nvrs; nsdm *= rs[v], v++ )
      {
        if (!ri[v].size()) ri[v].push_back(interval_t(0.0,1.0));
        rs[v] = ri[v].size();
      };
    /* allocation du nombre d'intervalles nécessaires représenter les sous-domaines */
    /* sdm[i] pointe sur sys->nvr() (nvrs) intervalles, pour i = 0..nsdm-1        */
    unsigned s,a,v; /* sous-domaine, addressage, variable */
    interval_t  *_tmp= new interval_t[ nsdm * nvrs ];
    interval_t** sdm = new interval_t*[ nsdm ];
    for ( sdm[0] = _tmp, s = 1; s < nsdm; sdm[s] = sdm[s-1] + nvrs, s ++) ;
    
    /* traduction entre l'addresse ou numéro (a) du sous-domaine et sa définition */
    /* sdm[a]    = [ r[0][x[0]], r[1][x[1]], r[2][x[2]], ..., r[nvrs][x[3]] ]     */
    /*         a = x[0]+rs[0]*(x[1]+rs[1]*( ...                                   */
    for ( s = 0; s < nsdm; s ++  )
      for ( a = s, v = 0; v < nvrs; sdm[s][v] = ri[v][a%rs[v]], a/=rs[v], v++ ) ;
    /* reduction simple insuffisante */
    if ( nsdm == 1 && rdscr(current(),sdm[0]) > rdscr_prm() ) return false;
    setmem();      /* mémorisation               */
    for ( s = 0; s < nsdm-1; s ++ )
      {
        update( this->top(), previous(), sdm[s] );
        mpush();
      };
    update( this->top(), previous(), sdm[s] );
    return true;
      };

    bool subdivide( sz_t v = (sz_t)(-1), creal_t m = (creal_t)0 )
    {
      creal_t prec(prc());
      if ( prec < m_spleps ) return false;
      /* split(-1,m) = split(largest(),m) */
      v = ((v==((sz_t)(-1)))?largest():v);
      /* si la subdivision est trop proche des bornes de l'intervalle */
      if ( m*(1.0-m) < 0.01 ) m = 0.5;     
      /* taille de l'intervalle a subdiviser trop faible */
      creal_t d(size(current()[v]));
      if ( d < m_spleps ) 
    { 
      v = largest(); 
      d = prec; 
    };
      /* calcul de la position en absolu */
      creal_t t(current()[v].lower() + m * d);
      /* mémorisation du domaine courant */
      interval_t * sv = current();
      cpush ();
      new (sv+v) interval_t(sv[v].lower(),t);
      new (current()+v) interval_t(t,current()[v].upper());
      
      return true;
    
    };
  
  };
}
//--------------------------------------------------------------------
} //namespace mmx
/********************************************************************/
#endif //