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 = 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 ) { 
    return R_FAIL;
      };
      /* on propose les réductions au système, c'est à lui de disposer */
      return  sys->reduce(ri,sys->nvr())?R_ISOK:R_WEAK;
    };
    
    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 //