|
realroot_doc 0.1.1
|
|
realroot_doc 0.1.1
|
00001 /******************************************************************** 00002 * This file is part of the source code of the realroot library. 00003 * Author(s): J.P. Pavone, GALAAD, INRIA 00004 * $Id: method.hpp,v 1.1 2005/07/11 10:03:55 jppavone Exp $ 00005 ********************************************************************/ 00006 #ifndef realroot_SOLVE_SBDSLV_METHOD_HPP 00007 #define realroot_SOLVE_SBDSLV_METHOD_HPP 00008 // JW - the following file is not under svn control ... 00009 //#include "../realroot2/system.H" 00010 #include <realroot/system_system.h> 00011 #include <realroot/strgy_simple.h> 00012 #include <realroot/strgy_newton.h> 00013 #include <realroot/strgy_cvmatrix.h> 00014 #include <realroot/strgy_rdslv_parallel.hpp> 00015 #include <realroot/strgy_sbdrl_parametric.hpp> 00016 //-------------------------------------------------------------------- 00017 namespace mmx { 00018 //-------------------------------------------------------------------- 00019 //#include "debug_stream.hpp" 00020 00021 namespace realroot 00022 { 00023 00024 /* erreur detectées par le system provoquées par un les parametres */ 00025 const int E_CTRL = 0; 00026 const int E_STRGY = 1; 00027 const int E_RDSLV = 2; 00028 const int E_SBDRL = 3; 00029 const int E_INIT = 4; 00030 00031 /* ordre recu par le systeme depuis l'un des parametres */ 00032 const int C_ACCEPT = 5; 00033 00034 /* reponses possibles de la méthode de réduction */ 00035 00036 const int R_REJECT = 6; 00037 00038 const int R_ISOK = 7; 00039 const int R_WEAK = 8; 00040 const int R_FAIL = 9; 00041 const int R_ERROR = 10; 00042 const int D_REJECT = 11; 00043 00044 struct method_debug /* retour d'erreurs */ 00045 { 00046 virtual ~method_debug(void) {} 00047 virtual void output(const char * msg) { 00048 std::cerr << msg << std::endl; 00049 exit(1); 00050 } 00051 } ; 00052 00053 /* méthode paramétree par strate'gie, solveur de réduction et règle de subdivision */ 00054 template< class system, 00055 class _strgy_ = strgy::newton<system> , 00056 template<class> class _rdslv_ = rdslv::parallel, 00057 template<class> class _sbdrl_ = sbdrl::parametric > 00058 struct method : method_base 00059 { 00060 typedef typename system::creal_t creal_t; 00061 typedef typename system::interval_t interval_t; 00062 typedef typename system::sz_t sz_t; 00063 method_debug * m_dbg; 00064 method_debug m_dbgdefault; 00065 /* paramètres de la methode */ 00066 _strgy_ m_strgy; 00067 _rdslv_<system> m_rdslv; 00068 _sbdrl_<system> m_sbdrl; 00069 00070 int m_state; 00071 00072 void accept() 00073 { 00074 /* le system a recu l'ordre d'une des methodes d'accepter le domaine */ 00075 m_state = C_ACCEPT; 00076 }; 00077 00078 void error( const char * sysmsg ) 00079 { 00080 char msg[ 200 ]; 00081 switch( m_state ) 00082 { 00083 case E_CTRL: 00084 sprintf(msg,"domain control:\n\t%s\n",sysmsg); 00085 break; 00086 case E_STRGY: 00087 sprintf(msg,"strategy:\n\t%s\n",sysmsg); 00088 break; 00089 case E_RDSLV: 00090 sprintf(msg,"reduction(solveur):\n\t%s\n",sysmsg); 00091 break; 00092 case E_SBDRL: 00093 sprintf(msg,"subdivision:\n\t%s\n",sysmsg); 00094 break; 00095 case E_INIT: 00096 sprintf(msg,"initialisation:\n\t%s\n",sysmsg); 00097 break; 00098 case R_ERROR: 00099 sprintf(msg,"reduction:\n\t%s\n",sysmsg); 00100 break; 00101 case R_FAIL: 00102 sprintf(msg,"%s","reduction: projection stack is empty !\n"); 00103 break; 00104 case D_REJECT: 00105 sprintf(msg,"%s","rejection: pb in pop()\n"); 00106 }; 00107 m_dbg->output( msg ); 00108 }; 00109 00110 void check_pstack( system * sys ) 00111 { 00112 for ( sz_t v = 0; v < sys->nvr(); v ++ ) 00113 { 00114 // std::cout << "projection : " << v << std::endl; 00115 }; 00116 }; 00117 /* la méthode de réduction consiste à appliquer le solveur de réduction sur */ 00118 /* toutes les projections disponibles, et à mettre à jour la pile de domaines */ 00119 /* paramètres: le systeme utilisé et le taux de réduction esperé */ 00120 /* retourne: */ 00121 /* R_ISOK en cas de succès */ 00122 /* R_WEAK en cas de réduction trop faible (i.e. < ratio) */ 00123 /* R_FAIL si la pile de projections est vide. */ 00124 /* R_REJECT si le solveur de réduction rejette le domaine */ 00125 00126 int reduction( system * sys ) 00127 { 00128 std::vector<interval_t> ri[ sys->nvr() ]; 00129 m_state = R_FAIL; 00130 for ( sz_t v = 0; v < sys->nvr(); v ++ ) 00131 { 00132 if ( sys->projections(v).nbp() == 0 ) continue; 00133 m_state = E_RDSLV; 00134 creal_t prc = sys->current()[v].size(); 00135 creal_t eps = std::min(std::max(sys->peps()/(10*prc),(creal_t)1e-4),(creal_t)1e-2); 00136 if ( ! m_rdslv.process(ri[v],sys->projections(v),eps) ) return R_REJECT; 00137 }; 00138 00139 if ( m_state == R_FAIL ) { 00140 return R_FAIL; 00141 }; 00142 /* on propose les réductions au système, c'est à lui de disposer */ 00143 return sys->reduce(ri,sys->nvr())?R_ISOK:R_WEAK; 00144 }; 00145 00146 void//int 00147 subdivision( system * sys ) 00148 { 00149 m_state = E_SBDRL; 00150 m_sbdrl.process( sys ); 00151 }; 00152 00153 int m_niter; 00154 int m_nsbd ; 00155 int m_seq; 00156 00157 /* allumage de l'usine à gaz */ 00158 template<class Prm, class Bounds> 00159 void launch( system * sys, Prm& prm, 00160 Bounds * inits,//creal_t * inits, 00161 method_debug * dbg = 0 ) 00162 { 00163 m_seq = 0; 00164 m_niter = 0; 00165 m_nsbd = 0; 00166 /* si la stream de debug n'est pas fourni on utilise celle donnée par défaut */ 00167 m_dbg = (dbg)?dbg:&m_dbgdefault; 00168 unsigned answ; 00169 sys->receiver(this); /* la methode ecoute les messages du systeme */ 00170 m_state = E_INIT; 00171 sys->init( inits ); /* initialisation du domaine de recherche */ 00172 /* tant que la pile de domaines n'est pas vide */ 00173 while ( sys->current() ) /* tant que la pile de domaines n'est pas vide */ 00174 { 00175 // 12 00176 // if ( m_niter == 12 ) std::cout << "POPO\n"; 00177 m_niter++; 00178 m_seq ++; 00179 sys->reset(); 00180 /* on passe le domaine à la classe de contrôle */ 00181 m_state = E_CTRL; 00182 // vctops::print(sys->current(),sys->nvr()); 00183 answ = prm.check(sys->current(),sys->nvr()); 00184 /* si elle retourne false on rejette le domaine */ 00185 m_state = D_REJECT; if ( !answ ) { sys->pop(); continue; }; 00186 /* on applique la strategie qui consiste à remplir la pile de projection */ 00187 m_state = E_STRGY, answ = m_strgy.process(sys); 00188 /* pour une raison inconnue la strategie veut que l'on accepte le domaine*/ 00189 if ( m_state == C_ACCEPT ) { 00190 // std::cout << " ACCEPT = STRGY\n"; 00191 prm.output(sys->current(),sys->nvr()); sys->pop(); continue; }; 00192 /* la strategie peut rejeter le domaine */ 00193 m_state = D_REJECT; if ( !answ ) { sys->pop(); continue; }; 00194 /* phase de réduction du domaine par l'utilisation de la pile de projection */ 00195 00196 if ( sys->thickness() ) { 00197 // std::cout << "ACCEPT = THICK \n"; 00198 prm.output(sys->current(),sys->nvr()); sys->pop(); continue; }; 00199 00200 switch(reduction(sys)) 00201 { 00202 case R_REJECT: /* le domaine est à rejeter */ 00203 m_state = R_REJECT; 00204 sys->pop(); 00205 continue; 00206 case R_ISOK: /* le domaine est réduit */ 00207 m_state = R_ISOK; 00208 break; 00209 case R_WEAK: /* la réduction n'est pas acceptée */ 00210 break; 00211 case R_FAIL: /* échec de la réduction */ 00212 m_state = R_FAIL; 00213 error(""); 00214 break; 00215 }; 00216 00217 if ( m_state == R_ISOK ) { sys->dreset(); }; 00218 00219 if ( sys->prc() < sys->peps() ) 00220 { 00221 // std::cout << " ACCEPT = PREC\n"; 00222 prm.output(sys->current(),sys->nvr()); 00223 sys->pop(); 00224 } 00225 else 00226 { 00227 if ( m_state != R_ISOK ) 00228 { 00229 m_seq = 0; 00230 m_nsbd ++; 00231 m_state = E_SBDRL; 00232 if ( !m_sbdrl.process(sys) ) { 00233 // std::cout << " ACCEPT = SBD\n"; 00234 prm.output(sys->current(),sys->nvr()); sys->pop(); 00235 }; 00236 }; 00237 }; 00238 }; 00239 00240 }; 00241 }; 00242 }; 00243 //-------------------------------------------------------------------- 00244 } //namespace mmx 00245 /********************************************************************/ 00246 #endif //
