include/realroot/system_method.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: method.hpp,v 1.1 2005/07/11 10:03:55 jppavone Exp $
 ********************************************************************/
#ifndef realroot_SOLVE_SBDSLV_METHOD_HPP
#define realroot_SOLVE_SBDSLV_METHOD_HPP
// JW - the following file is not under svn control ...
//#include "../realroot2/system.H"
#include <realroot/system_system.h>
#include <realroot/strgy_simple.h>
#include <realroot/strgy_newton.h>
#include <realroot/strgy_cvmatrix.h>
#include <realroot/strgy_rdslv_parallel.hpp>
#include <realroot/strgy_sbdrl_parametric.hpp>
//--------------------------------------------------------------------
namespace mmx {
//--------------------------------------------------------------------
//#include "debug_stream.hpp"

namespace realroot
{
  
  /* erreur detectées par le system provoquées par un les parametres */
  const int E_CTRL   = 0;
  const int E_STRGY  = 1;
  const int E_RDSLV  = 2;
  const int E_SBDRL  = 3;
  const int E_INIT   = 4;
  
  /* ordre recu par le systeme depuis l'un des parametres */
  const int C_ACCEPT = 5;
  
  /* reponses possibles de la méthode de réduction */
  
  const int R_REJECT = 6;

  const int R_ISOK   = 7;
  const int R_WEAK   = 8;
  const int R_FAIL   = 9;
  const int R_ERROR  = 10;
  const int D_REJECT = 11;
  
  struct method_debug /* retour d'erreurs */
  { 
      virtual     ~method_debug(void)       {}
      virtual void output(const char * msg) {
          std::cerr << msg << std::endl;
          exit(1);
      }
  } ;

  /* méthode paramétree par strate'gie, solveur de réduction et règle de subdivision */
  template< class system, 
            class _strgy_ = strgy::newton<system> ,
            template<class> class _rdslv_ = rdslv::parallel, 
            template<class> class _sbdrl_ = sbdrl::parametric >
  struct method : method_base 
  {
    typedef typename system::creal_t    creal_t;
    typedef typename system::interval_t interval_t;
    typedef typename system::sz_t       sz_t;
    method_debug * m_dbg;
    method_debug   m_dbgdefault; 
    /* paramètres de la methode */
    _strgy_         m_strgy;
    _rdslv_<system> m_rdslv;
    _sbdrl_<system> m_sbdrl;
    
    int m_state;
    
    void accept()
    {
      /* le system a recu l'ordre d'une des methodes d'accepter le domaine */
      m_state = C_ACCEPT;
    };
    
    void error( const char * sysmsg )
    {
      char msg[ 200 ];
      switch( m_state )
        {
        case E_CTRL:
          sprintf(msg,"domain control:\n\t%s\n",sysmsg);
          break;
        case E_STRGY:
          sprintf(msg,"strategy:\n\t%s\n",sysmsg);
          break;
        case E_RDSLV:
          sprintf(msg,"reduction(solveur):\n\t%s\n",sysmsg);
          break;
        case E_SBDRL:
          sprintf(msg,"subdivision:\n\t%s\n",sysmsg);
          break;
        case E_INIT:
          sprintf(msg,"initialisation:\n\t%s\n",sysmsg);
          break;
        case R_ERROR:
          sprintf(msg,"reduction:\n\t%s\n",sysmsg);
          break;
        case R_FAIL:
          sprintf(msg,"%s","reduction: projection stack is empty !\n");
          break;
        case D_REJECT:
          sprintf(msg,"%s","rejection: pb in pop()\n");
        };
      m_dbg->output( msg );
    };    
    
    void check_pstack( system * sys )
    {
      for ( sz_t v = 0; v < sys->nvr(); v ++ )
        {
          //      std::cout << "projection : " << v << std::endl;
        };
    };
    /* la méthode de réduction consiste à appliquer le solveur de réduction sur   */
    /* toutes les projections disponibles, et à mettre à jour la pile de domaines */
    /* paramètres: le systeme utilisé et le taux de réduction esperé              */
    /* retourne:                                                                  */
    /*          R_ISOK   en cas de succès                                         */
    /*          R_WEAK   en cas de réduction trop faible (i.e. < ratio)           */
    /*          R_FAIL   si la pile de projections est vide.                      */
    /*          R_REJECT si le solveur de réduction rejette le domaine            */
    
    int reduction( system * sys )
    {
      int nvrTemp =  sys->nvr();

      // std::vector<interval_t> ri[nvrTemp];

      std::vector<interval_t> *ri = new std::vector<interval_t>[nvrTemp];

      m_state = R_FAIL;
      for ( sz_t v = 0; v < sys->nvr(); v ++ )
        {
          if ( sys->projections(v).nbp()  ==  0 ) continue;
          m_state = E_RDSLV;
          creal_t prc = sys->current()[v].size();
          creal_t eps = std::min(std::max(sys->peps()/(10*prc),(creal_t)1e-4),(creal_t)1e-2);
          if ( ! m_rdslv.process(ri[v],sys->projections(v),eps) ) return R_REJECT;
        };      
      
      if ( m_state == R_FAIL ) { 
        delete ri;
        return R_FAIL;
      };
      /* on propose les réductions au système, c'est à lui de disposer */
      int r = sys->reduce(ri,sys->nvr())?R_ISOK:R_WEAK;

      delete ri;

      return r;
    };
    
    void//int 
    subdivision( system * sys )
    {
      m_state = E_SBDRL;
      m_sbdrl.process( sys );
    };
    
    int m_niter;
    int m_nsbd ; 
    int m_seq;

    /* allumage de l'usine à gaz */
    template<class Prm, class Bounds> 
    void launch( system * sys, Prm& prm, 
                 Bounds * inits,//creal_t  * inits,
                 method_debug * dbg = 0  )
    {
      m_seq   = 0;
      m_niter = 0;
      m_nsbd  = 0;
      /* si la stream de debug n'est pas fourni on utilise celle donnée par défaut */
      m_dbg = (dbg)?dbg:&m_dbgdefault;
      unsigned answ;
      sys->receiver(this);       /* la methode ecoute les messages du systeme */
      m_state = E_INIT;
      sys->init( inits );  /* initialisation du domaine de recherche    */
      /* tant que la pile de domaines n'est pas vide */
      while ( sys->current() ) /* tant que la pile de domaines n'est pas vide */
        {
          // 12
          //      if ( m_niter == 12 ) std::cout << "POPO\n";
          m_niter++;
          m_seq ++;
          sys->reset();
          /* on passe le domaine à la classe de contrôle */
          m_state = E_CTRL;
          //      vctops::print(sys->current(),sys->nvr());
          answ = prm.check(sys->current(),sys->nvr());
          /* si elle retourne false on rejette le domaine */
          m_state = D_REJECT; if ( !answ ) { sys->pop(); continue; };
          /* on applique la strategie qui consiste à remplir la pile de projection */
          m_state = E_STRGY, answ = m_strgy.process(sys);
          /* pour une raison inconnue la strategie veut que l'on accepte le domaine*/
          if ( m_state == C_ACCEPT ) { 
            //      std::cout << " ACCEPT = STRGY\n";
            prm.output(sys->current(),sys->nvr()); sys->pop();  continue; };
          /* la strategie peut rejeter le domaine                                  */
          m_state = D_REJECT; if ( !answ ) { sys->pop(); continue; };
          /* phase de réduction du domaine par l'utilisation de la pile de projection          */

          if ( sys->thickness()  ) { 
            //      std::cout << "ACCEPT = THICK \n";
            prm.output(sys->current(),sys->nvr()); sys->pop();  continue; };

          switch(reduction(sys))
            {
            case R_REJECT:  /* le domaine est à rejeter */
              m_state = R_REJECT;
              sys->pop(); 
              continue;
            case R_ISOK:      /* le domaine est réduit    */
              m_state = R_ISOK;
              break;
            case R_WEAK:   /* la réduction n'est pas acceptée */
              break;
            case R_FAIL:   /* échec de la réduction    */
              m_state = R_FAIL;
              error("");
              break;
            };
          
          if ( m_state == R_ISOK ) { sys->dreset(); };
          
          if ( sys->prc() < sys->peps() )
            {
              //              std::cout << " ACCEPT = PREC\n";
              prm.output(sys->current(),sys->nvr());
              sys->pop();
            }
          else
            {
              if ( m_state != R_ISOK )
                {
                  m_seq = 0;
                  m_nsbd ++;
                  m_state = E_SBDRL;
                  if ( !m_sbdrl.process(sys) ) { 
                    //          std::cout << " ACCEPT = SBD\n";
                    prm.output(sys->current(),sys->nvr()); sys->pop(); 
                  };
                };
            };
        };

    };
  };
};
//--------------------------------------------------------------------
} //namespace mmx
/********************************************************************/
#endif //

---