00001
00002
00003
00004
00005
00007 #ifndef SYNAPS_SOLVE_SBDSLV_H
00008 #define SYNAPS_SOLVE_SBDSLV_H
00009
00010 #include <set>
00011 #include<vector>
00012 #include <synaps/init.h>
00013 #include <synaps/base/xml.h>
00014 #include <synaps/mpol.h>
00015 #include <synaps/linalg/VectDse.h>
00016 #include <synaps/linalg/MatrDse.h>
00017 #include <synaps/msolve/sbdslv/method.hpp>
00018 #include <synaps/msolve/sbdslv/support.hpp>
00019 #include <synaps/base/chrono.h>
00020 #include <synaps/mpol/MPolDse.h>
00021
00022 __BEGIN_NAMESPACE_SYNAPS
00023
00024 namespace sbdslv
00025 {
00026 #define SBDSLV_RD 0
00027 #define SBDSLV_RDL 1
00028 #define SBDSLV_RDS 2
00029 #define SBDSLV_RDRDL 3
00030 #define SBDSLV_RDRDLRDS 4
00031 #define SBDSLV_RDRDS 5
00032 #define SBDSLV_RDLRDS 6
00033 #define SBDSLV_RDSRDL 7
00034 #ifndef SBDSLV_MTH
00035 #define SBDSLV_MTH SBDSLV_RDRDL
00036 #endif
00037
00038 template< class system, int mth >
00039 struct select_mth
00040 { typedef sbdslv::method<system,strgy::newton<system>,rdslv::parallel,sbdrl::parametric> result_t; };
00041 template< class system >
00042 struct select_mth<system,SBDSLV_RD>
00043 { typedef sbdslv::method<system,strgy::simple<system>,rdslv::parallel,sbdrl::parametric> result_t; };
00044 template< class system >
00045 struct select_mth<system,SBDSLV_RDL>
00046 { typedef sbdslv::method<system,strgy::newton<system>,rdslv::parallel,sbdrl::parametric> result_t; };
00047 template< class system>
00048 struct select_mth<system,SBDSLV_RDRDL>
00049 { typedef sbdslv::method<system,strgy::newton< system, strgy::simple<system> >,rdslv::parallel,sbdrl::parametric>
00050 result_t; };
00051 template< class system >
00052 struct select_mth<system,SBDSLV_RDRDLRDS>
00053 { typedef sbdslv::method<system,strgy::cvmatrix< system, strgy::newton<system,strgy::simple<system> > >,rdslv::parallel,sbdrl::parametric> result_t; };
00054 template< class system >
00055 struct select_mth<system,SBDSLV_RDS>
00056 { typedef sbdslv::method<system,strgy::cvmatrix<system>,rdslv::parallel,sbdrl::parametric> result_t; };
00057
00058 template< class system >
00059 struct select_mth<system,SBDSLV_RDRDS>
00060 { typedef sbdslv::method<system,strgy::cvmatrix<system, strgy::simple<system> >,rdslv::parallel,sbdrl::parametric> result_t; };
00061
00062 template< class system >
00063 struct select_mth<system,SBDSLV_RDLRDS>
00064 { typedef sbdslv::method<system,strgy::cvmatrix<system, strgy::newton<system> >,rdslv::parallel,sbdrl::parametric> result_t; };
00065
00066 template< class system >
00067 struct select_mth<system,SBDSLV_RDSRDL>
00068 { typedef sbdslv::method<system,strgy::cvmatrix<system, strgy::newton<system> >,rdslv::parallel,sbdrl::parametric> result_t; };
00069
00070 };
00071
00083 template< class T, int mth = SBDSLV_RDRDL>
00084 struct SBDSLV
00085 {
00086 typedef T coeff_t;
00087 T prec;
00089 SBDSLV( const T& precision = (T)1e-6 ): prec(precision){};
00090 };
00091
00092 namespace SolveSBDSLV
00093 {
00094 template< class T, class eenv_t, class MPOL >
00095 void fill_data( T * data, eenv_t * env, const MPOL& p, binomials<T>& benv )
00096 {
00097 typename MPOL::const_iterator mn; int i,a;
00098 std::fill(data,data+env->sz(),0);
00099 mpolfill(data,p,*env);
00100 convertm2b (data,*env,benv);
00101 };
00102 template< class monom >
00103 void scan_monom_env ( std::set< int >& env, const monom& m )
00104 { for ( unsigned i = 0; i < (m.rep()).size(); i++ ) if ( (m.rep())[i] != 0 ) env.insert(i); };
00105 }
00106
00107
00108
00115 template< class MPC, class C, int MTH, class Domains >
00116 Seq<VectDse<C> > solve( const MPC& l, const SBDSLV< C, MTH >& mth,
00117 const Domains & dom,
00118 bool verbose = false )
00119 {
00120 typedef C T;
00121 typedef typename MPC::value_type MPoly;
00122 typedef MPoly mpoly_t;
00123 typedef MPC plist_t;
00124 typedef typename plist_t::const_iterator pptr_t;
00125 typedef typename MPoly::const_iterator mptr_t;
00126 typedef typename MPoly::monom_t monom_t;
00127 typedef SBDSLV< C, MTH > mth_t;
00128 typedef int sz_t;
00129 sz_t nvars, neqs, * nvrs, ** vrs, ** szs;
00130 pptr_t cpl;
00131 int e;
00132 nvars = 0;
00133
00134 for ( neqs = 0, cpl = l.begin(); cpl != l.end(); cpl++, neqs++ )
00135 for ( mptr_t mn = cpl->begin(); mn != cpl->end(); mn++ )
00136 nvars = std::max(nvars,mn->nvars()+1);
00137
00138 std::set<int> * seenvs = new std::set< int > [ neqs ];
00139 for ( e = 0, cpl = l.begin(); e < neqs; e++, cpl++ )
00140 for ( mptr_t mn = cpl->begin(); mn != cpl->end(); mn++ )
00141 SolveSBDSLV::scan_monom_env(seenvs[e],*mn);
00142
00143
00144
00145 sz_t s = 0;
00146 for ( e = 0; e < neqs; s += seenvs[e].size(), e++ );
00147 nvrs = new sz_t [neqs+2*s];
00148 vrs = new sz_t*[2*neqs];
00149 szs = vrs + neqs;
00150 sz_t * base = nvrs + neqs;
00151 std::set<int>::iterator seenvit;
00152 for ( e = 0, cpl = l.begin(); e < neqs; base += 2*nvrs[e++], cpl++ )
00153 {
00154 nvrs[e] = seenvs[e].size();
00155 vrs [e] = base;
00156 szs [e] = base + nvrs[e];
00157 int k = 0;
00158 for ( seenvit = seenvs[e].begin(); k < nvrs[e]; seenvit++, k++ )
00159 {
00160 vrs[e][k] = *seenvit;
00161 szs[e][k] = MPOLDST::degree( *cpl, *seenvit ) + 1;
00162 };
00163 };
00164
00165 std::vector<C> result;
00166 C * dmns = 0;
00167 if ( dom.size() >= 2*nvars )
00168 {
00169 using namespace let;
00170 dmns = new C[ 2*nvars ];
00171 for ( int i = 0; i < 2*nvars; i ++ ) assign(dmns[i],dom[i]);
00172 };
00173
00174 chrono tm;
00175 std::vector<T> tmp;
00176
00177 {
00178 using namespace sbdslv;
00179 binomials<T> bznv;
00180 typedef sbdslv::system<T> system_t;
00181 system_t sys( neqs, nvars, nvrs, (const int**)vrs, (const int**)szs, (T)(mth.prec) );
00182 for ( e = 0, cpl = l.begin(); e < neqs; e ++, cpl++ )
00183 SolveSBDSLV::fill_data( sys.data(e), sys.env(e), *cpl, bznv );
00184 sbdslv::system_ctrl< std::vector< T > > intf(tmp);
00185 tm.start();
00186 typename sbdslv::select_mth<system_t,MTH>::result_t __mth;
00187 __mth.launch(&sys,intf,dmns);
00188 tm.stop();
00189 if ( verbose )
00190 {
00191 std::cout << " iterations = " << __mth.m_niter << std::endl;
00192 std::cout << "subdivisions = "<< __mth.m_nsbd << std::endl;
00193 };
00194 };
00195
00196 if ( !verbose && dmns ) delete[] dmns;
00197 delete[] seenvs;
00198 delete[] nvrs;
00199 delete[] vrs;
00200 if ( verbose )
00201 {
00202 std::cout << "<domain>";
00203 if ( dmns )
00204 {
00205 std::cout << std::endl;
00206 for ( int i = 0; i < nvars; i ++ )
00207 std::cout << "\t"
00208 << "x" << i
00209 << " = [ " << dmns[2*i]
00210 << ","
00211 << dmns[2*i+1]
00212 << " ]";
00213 }
00214 else std::cout << "[0,1]^" << nvars;
00215 std::cout << "</domain>\n";
00216
00217
00218 std::cout << "<solver_time>"; tm.write(std::cout);
00219 std::cout <<"</solver_time>"<< std::endl;
00220 if ( dmns ) delete[] dmns;
00221 };
00222
00223 int nsols = tmp.size() / (2*nvars);
00224 T boxes[nsols*2*nvars];
00225 std::copy(&tmp[0], &tmp[0]+2*nvars*nsols,boxes);
00226 int nsols2 = sbdslv::clean_result( &tmp[0], nvars, nsols,(T)1e-3 );
00227 if ( verbose ) {
00228 xml::print(std::cout,nsols,"nb_box")<< std::endl;
00229 for ( unsigned i = 0; i < nsols; i ++ )
00230 {
00231 std::cout << "\t[ ";
00232 for ( unsigned j = 0; j < nvars; j ++ )
00233 std::cout << (boxes[2*(i*nvars+j)]+boxes[2*(i*nvars+j)+1])/2.0 << " ";
00234 std::cout << "]";
00235 std::cout << std::endl;
00236 };
00237 xml::print(std::cout,nsols2,"nb_sol")<< std::endl;
00238
00239 };
00240
00241
00242 tmp.resize( nsols2 * nvars );
00243
00244 VectDse<T> vsol(nvars, AsSize());
00245 Seq<VectDse<T> > res;
00246 for(unsigned i=0;i<nsols2;++i)
00247 {
00248 for(unsigned j=0; j<nvars; j++)
00249 vsol[j]=tmp[i*nvars+j];
00250 res.push_back(vsol);
00251 }
00252 return res;
00253 };
00254
00255
00256 template< class MPC, class C, int MTH>
00257 Seq<VectDse<C> > solve( const MPC& l, const SBDSLV< C, MTH >& mth,
00258 bool verbose = false )
00259 {
00260 return solve(l,mth,VectDse<C>(),verbose);
00261 }
00262
00263
00264 template<class C, int METHOD>
00265 struct SbdSolver
00266 {
00267 typedef sbdslv::system<C> system_t;
00268 system_t * sys;
00269
00270
00271
00272 template< class MPDI >
00273 SbdSolver( MPDI first, MPDI last, const C& eps )
00274 {
00275
00276 int ne = 0;
00277 MPDI eptr;
00278 for ( eptr = first; eptr != last; eptr ++ ) ne++;
00279 int nvrs[ne];
00280 const int *szs[ne];
00281 int e;
00282 int mxvr = 0;
00283 std::set<int> gvrs;
00284 for ( e = 0, eptr = first; eptr != last; eptr ++, e++ ) {
00285 nvrs[e] = eptr->nvr();
00286 szs[e] = eptr->szs();
00287 const int * tmp = eptr->vrs();
00288 for ( int v = 0; v < eptr->nvr(); v ++ ) gvrs.insert(tmp[v]);
00289 };
00290 int mxv = *(gvrs.rend());
00291 int lnames[ mxv+1 ];
00292 int nv = 0;
00293 for ( std::set<int>::const_iterator it = gvrs.begin(); it != gvrs.end(); lnames[*it++] = nv++ );
00294 int * vrdat = new int[ nv * ne ];
00295 int * vrptr = vrdat;
00296 const int * vrs[ne];
00297 e = 0;
00298 for ( eptr = first; eptr != last; vrptr+= eptr->nvr(), e++, eptr ++ ) {
00299 vrs[e] = vrptr;
00300 const int * tmp = eptr->vrs();
00301 for ( int v = 0; v < eptr->nvr(); v ++ )
00302 vrptr[v] = lnames[ tmp[v] ];
00303 };
00304
00305 assert(ne>=nv);
00306 sys = new sbdslv::system<C>( ne, nv, nvrs, (const int**)vrs, (const int**)szs, eps );
00307 delete[] vrdat;
00308 };
00309
00310 void add( const mpoldse::bernstein<C>& data, int e )
00311 {
00312 sys->set_bernstein2( data.begin(), e );
00313 };
00314
00315 void add( const mpoldse::monomial<C>& data, int e )
00316 {
00317 sys->set_monomial2( data.begin(), e );
00318 };
00319
00320 template<class MPDI>
00321 void solve( std::vector<C> & R, int& ns, int& nv, MPDI first, MPDI last )
00322 {
00323 R.resize(0);
00324 typedef SBDSLV< C, METHOD > MTH;
00325 typename sbdslv::select_mth<system_t,METHOD>::result_t __mth;
00326 int e = 0;
00327 for ( MPDI eptr = first; eptr != last; add( eptr->rep(), e++ ), eptr++ );
00328
00329 sbdslv::system_ctrl< std::vector< C > > intf(R);
00330 int nvars = sys->nvr();
00331 __mth.launch(sys,intf,0);
00332
00333 int nsols = R.size() / (2*nvars);
00334 int nsols2 = sbdslv::clean_result( &R[0], nvars, nsols,(C)1e-3 );
00335 R.resize(nsols2) ;
00336 ns = nsols2;
00337 nv = nvars;
00338 };
00339
00340 ~SbdSolver() { delete sys; };
00341 };
00342
00343
00344
00345 __END_NAMESPACE_SYNAPS
00346
00347 #endif //SYNAPS_SOLVE_SBDSLV_H