Logo  0.95.0-final
Finite Element Embedded Library and Language in C++
Feel++ Feel++ on Github Feel++ community
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
crbscm.hpp
1 /* -*- mode: c++ -*-
2 
3  This file is part of the Feel library
4 
5  Author(s): Christophe Prud'homme <christophe.prudhomme@feelpp.org>
6  Date: 2009-08-07
7 
8  Copyright (C) 2009 Université Joseph Fourier (Grenoble I)
9 
10  This library is free software; you can redistribute it and/or
11  modify it under the terms of the GNU Lesser General Public
12  License as published by the Free Software Foundation; either
13  version 2.1 of the License, or (at your option) any later version.
14 
15  This library is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public
21  License along with this library; if not, write to the Free Software
22  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
23 */
29 #ifndef __CRBSCM_H
30 #define __CRBSCM_H 1
31 
32 #include <boost/archive/text_oarchive.hpp>
33 #include <boost/archive/text_iarchive.hpp>
34 #include <boost/timer.hpp>
35 #include <boost/tuple/tuple.hpp>
36 #include <boost/format.hpp>
37 #include <boost/assign/list_of.hpp>
38 #include <boost/foreach.hpp>
39 #include <feel/feelcore/feel.hpp>
40 #include <feel/feelcore/parameter.hpp>
43 #include <feel/feelcrb/crbdb.hpp>
45 #if defined(FEELPP_HAS_GLPK_H)
46 #include <glpk.h>
47 #endif /* FEELPP_HAS_GLPK_H */
48 
49 
50 namespace Feel
51 {
68 template<typename TruthModelType>
69 class CRBSCM : public CRBDB
70 {
71  typedef CRBDB super;
72 public:
73 
74 
78 
79 
81 
85 
86  typedef TruthModelType truth_model_type;
87  typedef truth_model_type model_type;
88  typedef boost::shared_ptr<truth_model_type> truth_model_ptrtype;
89 
90  typedef double value_type;
91  typedef boost::tuple<double,double> bounds_type;
92 
94  typedef boost::shared_ptr<parameterspace_type> parameterspace_ptrtype;
96  typedef typename parameterspace_type::element_ptrtype parameter_ptrtype;
98  typedef typename parameterspace_type::sampling_ptrtype sampling_ptrtype;
99 
100  typedef boost::tuple<double, parameter_type, size_type, double, double> relative_error_type;
101 
102  typedef typename model_type::backend_type backend_type;
103  typedef typename model_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
104  typedef typename model_type::vector_ptrtype vector_ptrtype;
105  typedef typename model_type::beta_vector_type beta_vector_type;
106 
107 
108  typedef Eigen::VectorXd y_type;
109  typedef Eigen::VectorXd vector_type;
110  typedef std::vector< std::vector<y_type> > y_set_type;
111  typedef std::vector< std::vector<boost::tuple<double,double> > > y_bounds_type;
113 
117 
120  :
121  super(),
122  M_is_initialized( false ),
123  M_model(),
124  M_tolerance( 1e-2 ),
125  M_iter_max( 3 ),
126  M_Mplus( 10 ),
127  M_Malpha( 4 ),
128  M_Dmu( new parameterspace_type ),
129  M_Xi( M_Dmu ),
130  M_C( M_Dmu, 1, M_Xi ),
131  M_C_complement( M_Dmu, 1, M_Xi ),
132  M_scm_for_mass_matrix( false ),
133  M_mu_ref( M_Dmu->element() ),
134  M_use_scm( true )
135  {
136  }
137 
139  CRBSCM( std::string const& name,
140  po::variables_map const& vm )
141  :
142  super( "scm",
143  ( boost::format( "%1%" ) % name ).str(),
144  ( boost::format( "%1%" ) % name ).str(),
145  vm ),
146  M_is_initialized( false ),
147  M_model(),
148  M_tolerance( vm["crb.scm.tol"].template as<double>() ),
149  M_iter_max( vm["crb.scm.iter-max"].template as<int>() ),
150  M_Mplus( vm["crb.scm.Mplus"].template as<int>() ),
151  M_Malpha( vm["crb.scm.Malpha"].template as<int>() ),
152  M_Dmu( new parameterspace_type ),
153  M_Xi( new sampling_type( M_Dmu ) ),
154  M_C( new sampling_type( M_Dmu, 1, M_Xi ) ),
155  M_C_complement( new sampling_type( M_Dmu, 1, M_Xi ) ),
156  M_vm( vm ),
157  M_scm_for_mass_matrix( false ),
158  M_mu_ref( M_Dmu->element() ),
159  M_use_scm( vm["crb.scm.use-scm"].template as<bool>() )
160  {
161  if ( this->loadDB() )
162  std::cout << "Database " << this->lookForDB() << " available and loaded\n";
163  }
164 
165 
167  CRBSCM( std::string const& name,
168  po::variables_map const& vm ,
169  truth_model_ptrtype const & model )
170  :
171  super( "scm",
172  ( boost::format( "%1%" ) % name ).str(),
173  ( boost::format( "%1%" ) % name ).str(),
174  vm ),
175  M_is_initialized( false ),
176  M_model(),
177  M_tolerance( vm["crb.scm.tol"].template as<double>() ),
178  M_iter_max( vm["crb.scm.iter-max"].template as<int>() ),
179  M_Mplus( vm["crb.scm.Mplus"].template as<int>() ),
180  M_Malpha( vm["crb.scm.Malpha"].template as<int>() ),
181  M_Dmu( new parameterspace_type ),
182  M_Xi( new sampling_type( M_Dmu ) ),
183  M_C( new sampling_type( M_Dmu, 1, M_Xi ) ),
184  M_C_complement( new sampling_type( M_Dmu, 1, M_Xi ) ),
185  M_vm( vm ),
186  M_scm_for_mass_matrix( false ),
187  M_mu_ref( M_Dmu->element() ),
188  M_use_scm( vm["crb.scm.use-scm"].template as<bool>() )
189  {
190  this->setTruthModel( model );
191  if ( this->loadDB() )
192  std::cout << "Database " << this->lookForDB() << " available and loaded\n";
193  }
194 
196  CRBSCM( std::string const& name,
197  po::variables_map const& vm ,
198  truth_model_ptrtype const & model ,
199  bool scm_for_mass_matrix )
200  :
201  super( "scm",
202  ( boost::format( "%1%" ) % name ).str(),
203  ( boost::format( "%1%" ) % name ).str(),
204  vm ),
205  M_is_initialized( false ),
206  M_model(),
207  M_tolerance( vm["crb.scm.tol"].template as<double>() ),
208  M_iter_max( vm["crb.scm.iter-max"].template as<int>() ),
209  M_Mplus( vm["crb.scm.Mplus"].template as<int>() ),
210  M_Malpha( vm["crb.scm.Malpha"].template as<int>() ),
211  M_Dmu( new parameterspace_type ),
212  M_Xi( new sampling_type( M_Dmu ) ),
213  M_C( new sampling_type( M_Dmu, 1, M_Xi ) ),
214  M_C_complement( new sampling_type( M_Dmu, 1, M_Xi ) ),
215  M_vm( vm ),
216  M_scm_for_mass_matrix( scm_for_mass_matrix ),
217  M_mu_ref( M_Dmu->element() ),
218  M_use_scm( vm["crb.scm.use-scm"].template as<bool>() )
219  {
220  this->setTruthModel( model );
221  if ( this->loadDB() )
222  std::cout << "Database " << this->lookForDB() << " available and loaded\n";
223  }
224 
226  CRBSCM( CRBSCM const & o )
227  :
228  super( o ),
229  M_is_initialized( o.M_is_initialized ),
230  M_tolerance( o.M_tolerance ),
231  M_iter_max( o.M_iter_max ),
232  M_Mplus( o.M_Mplus ),
233  M_Malpha( o.M_Malpha ),
234  M_Dmu( o.M_Dmu ),
235  M_Xi( o.M_Xi ),
236  M_C( o.M_C ),
237  M_C_complement( o.M_C_complement ),
238  M_vm( o.M_vm ),
239  M_scm_for_mass_matrix( o.M_scm_for_mass_matrix ),
240  M_mu_ref( o.M_mu_ref ),
241  M_use_scm( o.M_use_scm )
242  {
243  }
244 
247  {}
248 
250 
254 
256  CRBSCM& operator=( CRBSCM const & o )
257  {
258  if ( this != &o )
259  {
260  }
261 
262  return *this;
263  }
265 
269 
271  size_type KMax() const
272  {
273  return M_C->size();
274  }
275 
277  int maxIter() const
278  {
279  return M_iter_max;
280  }
281 
283  parameterspace_ptrtype Dmu() const
284  {
285  return M_Dmu;
286  }
287 
289  int Mplus() const
290  {
291  return M_Mplus;
292  }
293 
295  int Malpha()
296  {
297  return M_Malpha;
298  }
299 
301 
305 
307  void setTolerance( double tolerance )
308  {
309  M_tolerance = tolerance;
310  }
311 
313  void setMplus( int Mplus )
314  {
315  M_Mplus = Mplus;
316  }
317 
319  void setMalpha( int Malpha )
320  {
321  M_Malpha = Malpha;
322  }
323 
325  void setTruthModel( truth_model_ptrtype const& model )
326  {
327  M_model = model;
328  M_Dmu = M_model->parameterSpace();
329 
330  if ( ! loadDB() )
331  {
332  M_Xi = sampling_ptrtype( new sampling_type( M_Dmu ) );
333  M_C = sampling_ptrtype( new sampling_type( M_Dmu ) );
334  M_C_complement = sampling_ptrtype( new sampling_type( M_Dmu ) );
335  }
336  }
337 
339  void setMaxIter( int K )
340  {
341  M_iter_max = K;
342  }
343 
345  void setScmForMassMatrix( bool b )
346  {
347  M_scm_for_mass_matrix = b;
348  }
349 
350  int mMax( int q ) const
351  {
352  if ( M_scm_for_mass_matrix )
353  return M_model->mMaxM( q );
354  else
355  return M_model->mMaxA( q );
356  }
357 
358  int nb_decomposition_terms_q ( void ) const
359  {
360  if ( M_scm_for_mass_matrix )
361  return M_model->Qm();
362  else
363  return M_model->Qa();
364  }
365 
366  //total number of decomposition terms
367  int nb_decomposition_terms_qm ( void ) const
368  {
369  int count=0;
370  for(int q=0;q<nb_decomposition_terms_q();q++)
371  {
372  for(int m=0;m<mMax(q);m++)
373  count++;
374  }
375  return count;
376  }
377 
379 
383 
387  void generate();
388 
393  void checkC( size_type K ) const;
394 
395  boost::tuple<value_type,value_type> ex( parameter_type const& mu ) const;
396 
406  boost::tuple<value_type,value_type> lb( parameter_type const& mu, size_type K = invalid_size_type_value, int indexmu = -1 ) const;
407  boost::tuple<value_type,value_type> lbSCM( parameter_type const& mu, size_type K = invalid_size_type_value, int indexmu = -1 ) const;
408  boost::tuple<value_type,value_type> lbNoSCM( parameter_type const& mu, size_type K = invalid_size_type_value, int indexmu = -1 ) const;
418  boost::tuple<value_type,value_type> lb( parameter_ptrtype const& mu, size_type K = invalid_size_type_value, int indexmu = -1 ) const
419  {
420  return lb( *mu, K, indexmu );
421  }
422 
423 
433  boost::tuple<value_type,value_type> ub( parameter_type const& mu, size_type K = invalid_size_type_value ) const;
434 
444  boost::tuple<value_type,value_type> ub( parameter_ptrtype const& mu, size_type K = invalid_size_type_value ) const
445  {
446  return ub( *mu, K );
447  }
448 
449 
453  std::vector<boost::tuple<double,double,double> > offline();
454  std::vector<boost::tuple<double,double,double> > offlineSCM();
455  std::vector<boost::tuple<double,double,double> > offlineNoSCM();
456 
460  bounds_type
461  online() const
462  {
463  return bounds_type();
464  }
465 
470  relative_error_type maxRelativeError( size_type K ) const;
471 
475  void computeYBounds();
476 
481  std::vector<double> run( parameter_type const& mu, int K );
482 
492  void run( const double * X, unsigned long N, double * Y, unsigned long P );
493 
497  void saveDB();
498 
502  bool loadDB();
503 
504  bool rebuildDB();
505 
506  bool doScmForMassMatrix();
508 
509 
510  sampling_ptrtype c() const
511  {
512  return M_C;
513  }
514 
515 protected:
516 
517 private:
518 private:
519  friend class boost::serialization::access;
520 
521  template<class Archive>
522  void save( Archive & ar, const unsigned int version ) const;
523 
524  template<class Archive>
525  void load( Archive & ar, const unsigned int version ) ;
526 
527  BOOST_SERIALIZATION_SPLIT_MEMBER()
528 
529 private:
530 
531  bool M_is_initialized;
532 
533 
534  truth_model_ptrtype M_model;
535 
536  double M_tolerance;
537  int M_iter_max;
538  size_type M_Mplus;
539  size_type M_Malpha;
540 
541  // parameter space
542  parameterspace_ptrtype M_Dmu;
543 
544  // sampling of parameter space
545  sampling_ptrtype M_Xi;
546 
547  // sampling of parameter space (subset of M_Xi later on)
548  sampling_ptrtype M_C, M_C_complement;
549 
551  std::map<size_type,value_type> M_C_eigenvalues;
552  std::vector<double> M_eig;
553  mutable std::map<size_type,std::map<size_type,value_type> > M_C_alpha_lb;
554 
556  y_set_type M_Y_ub;
557 
559  y_bounds_type M_y_bounds;
560  std::vector< std::vector<double> > M_y_bounds_0;
561  std::vector< std::vector<double> > M_y_bounds_1;
562 
563  po::variables_map M_vm;
564 
565  bool M_scm_for_mass_matrix;
566  bool M_print_matrix;
567 
568  parameter_type M_mu_ref;
569  bool M_use_scm;
570 };
571 
572 po::options_description crbSCMOptions( std::string const& prefix = "" );
573 
574 
575 
576 template<typename TruthModelType>
577 std::vector<boost::tuple<double,double,double> >
578 CRBSCM<TruthModelType>::offline()
579 {
580  if( M_use_scm )
581  return this->offlineSCM();
582  else
583  return this->offlineNoSCM();
584 }
585 
586 template<typename TruthModelType>
587 std::vector<boost::tuple<double,double,double> >
589 {
590  sparse_matrix_ptrtype inner_prod,sym,Matrix;
591  M_mu_ref = M_model->refParameter();
592  M_model->computeAffineDecomposition();
593  if ( M_scm_for_mass_matrix )
594  {
595  inner_prod = M_model->innerProductForMassMatrix();
596  boost::tie( Matrix, boost::tuples::ignore, boost::tuples::ignore ) = M_model->update( M_mu_ref );
597  }
598  else
599  {
600  inner_prod = M_model->innerProduct();
601  boost::tie( boost::tuples::ignore, Matrix, boost::tuples::ignore ) = M_model->update( M_mu_ref );
602  }
603  sym = M_model->newMatrix();sym->close();
604  Matrix->symmetricPart( sym );
605  // solve for eigenvalue problem at mu_ref
606  SolverEigen<double>::eigenmodes_type modes;
607 
608  modes=
609  eigs( _matrixA=sym,
610  _matrixB=inner_prod,
611  _solver=( EigenSolverType )M_vm["crb.scm.solvereigen-solver-type"].template as<int>(),
612  _spectrum=SMALLEST_REAL,
613  //_spectrum=LARGEST_MAGNITUDE,
614  _transform=SINVERT,
615  _ncv=M_vm["crb.scm.solvereigen-ncv"].template as<int>(),
616  _nev=M_vm["crb.scm.solvereigen-nev"].template as<int>(),
617  _tolerance=M_vm["crb.scm.solvereigen-tol"].template as<double>(),
618  _maxit=M_vm["crb.scm.solvereigen-maxiter"].template as<int>()
619  );
620  double eigen_value = modes.begin()->second.template get<0>();
621  //std::cout<<"-------------------------------------------"<<std::endl;
622  //std::cout<<"eigenvalue ( min ) for mu_ref : "<<eigen_value<<std::endl;
623  //std::cout<<"-------------------------------------------"<<std::endl;
624 
625  //store the eigen value in M_C_eigenvalues
626  M_C_eigenvalues[0] = modes.begin()->second.template get<0>();
627 
628  saveDB();
629 
630  //only to have the same signature that offlineSCM
631  std::vector<boost::tuple<double,double,double> > ckconv;
632  return ckconv;
633 }
634 
635 template<typename TruthModelType>
636 std::vector<boost::tuple<double,double,double> >
637 CRBSCM<TruthModelType>::offlineSCM()
638 {
639  std::ofstream os_y( "y.m" );
640  std::ofstream os_C( "C.m" );
641 
642 
643  std::vector<boost::tuple<double,double,double> > ckconv;
644  // do the affine decomposition
645  M_model->computeAffineDecomposition();
646 
647  // random sampling
648  M_Xi->randomize( M_vm["crb.scm.sampling-size"].template as<int>() );
649  //M_Xi->logEquidistribute( M_vm["crb.scm.sampling-size"].template as<int>() );
650  M_C->setSuperSampling( M_Xi );
651  parameter_type mu( M_Dmu );
652 
653 
654  M_print_matrix = M_vm["crb.scm.print-matrix"].template as<bool>() ;
655 
656 #if 0
657 
658  //Test coercicivty
659  //std::cout << "Testing coercivity at sample points\n" ;
660  for ( int i=0; i<M_Xi->size(); i++ )
661  {
662  double testalpha = ex( M_Xi->at( i ) );
663  //std::cout << "mu[" << i+1 << "]\n" << M_Xi->at(i) << "\nalpha = " << testalpha << "\n" ;
664  }
665 
666 #endif
667  double relative_error = 1e30;
668 
669 
670  // compute the bounds for y in R^Q
671  this->computeYBounds();
672 
673  // empty sets
674  M_C->clear();
675  M_Y_ub.clear();
676  M_C_alpha_lb.clear();
677 
678  size_type index;
679 
680  bool use_predefined_C = option(_name="crb.scm.use-predefined-C").template as<bool>();
681  int N_log_equi = this->vm()["crb.scm.use-logEquidistributed-C"].template as<int>() ;
682  int N_equi = this->vm()["crb.scm.use-equidistributed-C"].template as<int>() ;
683  std::vector<int> index_vector;
684 
685  if( N_log_equi > 0 || N_equi > 0 )
686  use_predefined_C = true;
687 
688  if( use_predefined_C )
689  {
690  std::string file_name = ( boost::format("SamplingC") ).str();
691  std::ifstream file ( file_name );
692  if( ! file )
693  throw std::logic_error( "[CRBSCM::offline] ERROR the file SamplingC doesn't exist so it's impossible to known which parameters you want to use to build the database" );
694  else
695  {
696  M_C->clear();
697  index_vector = M_C->closestSamplingFromFile(file_name);
698  int sampling_size = index_vector.size();
699  M_iter_max = sampling_size;
700  }
701  mu = M_C->at( 0 ); // first element
702  index = index_vector[0];
703 
704  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
705  std::cout<<"[CRBSCM::offline] read sampling C ( sampling size : "<<M_iter_max<<" )"<<std::endl;
706  }
707  else
708  {
709  // start with M_C = { arg min mu, mu \in Xi }
710  boost::tie( mu, index ) = M_Xi->min();
711  M_C->push_back( mu, index );
712  }
713 
714  M_C_complement = M_C->complement();
715  //std::cout << " -- start with mu = " << mu << "\n";
716  //std::cout << " -- C size : " << M_C->size() << "\n";
717  //std::cout << " -- C complement size : " << M_C_complement->size() << "\n";
718 
719  sparse_matrix_ptrtype B,symmMatrix,Matrix;
720  std::vector<vector_ptrtype> F;
721 
722  // dimension of Y_UB and U_LB
723  int K = 1;
724 
725 
726  std::vector< std::vector<sparse_matrix_ptrtype> > Matrixq;
727 
728  if ( M_scm_for_mass_matrix )
729  boost::tie( Matrixq, boost::tuples::ignore, boost::tuples::ignore ) = M_model->computeAffineDecomposition();
730  else
731  boost::tie( boost::tuples::ignore, Matrixq, boost::tuples::ignore ) = M_model->computeAffineDecomposition();
732 
733  int Qmax = nb_decomposition_terms_q();
734 
735  while ( relative_error > M_tolerance && K <= M_iter_max )
736  {
737  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
738  {
739  std::cout << "============================================================\n";
740  std::cout << "K=" << K << "\n";
741  }
742 
743  os_C << M_C->at( K-1 ) << "\n";
744 
745  // resize y_ub
746  M_Y_ub.resize( K );
747  M_Y_ub[K-1].resize( Qmax );
748  for(int q=0; q<Qmax; q++)
749  M_Y_ub[K-1][q].resize( mMax(q) );
750 
751  M_model->solve( mu );
752 
753 
754  // for a given parameter \p mu assemble the left and right hand side
755  if ( M_scm_for_mass_matrix )
756  {
757  B = M_model->innerProductForMassMatrix();
758  boost::tie( Matrix, boost::tuples::ignore, F ) = M_model->update( mu );
759  }
760  else
761  {
762  B = M_model->innerProduct();
763  boost::tie( boost::tuples::ignore, Matrix, F ) = M_model->update( mu );
764  }
765 
766  std::string mu_str;
767  for(int i=0;i<mu.size();i++)
768  mu_str= mu_str + (boost::format("_%1%") %mu[i]).str() ;
769 
770  std::string file_name;
771 
772  symmMatrix = M_model->newMatrix();symmMatrix->close();
773  Matrix->symmetricPart( symmMatrix );
774 
775  if( M_print_matrix && ( Environment::worldComm().globalSize() == 1 ) )
776  {
777  if( M_scm_for_mass_matrix)
778  {
779  file_name = "offline_Matrix_M" + mu_str +".m";
780  Matrix->printMatlab( file_name );
781  file_name = "offline_symmMatrix_M" + mu_str +".m";
782  symmMatrix->printMatlab( file_name );
783  }
784  else
785  {
786  file_name = "offline_Matrix_A" + mu_str +".m";
787  Matrix->printMatlab( file_name );
788  file_name = "offline_symmMatrix_A" + mu_str +".m";
789  symmMatrix->printMatlab( file_name );
790  }
791  file_name = "offline_B" + mu_str +".m";
792  B->printMatlab( file_name );
793  }
794 
795  //
796  // Build Y_UB
797  //
798  vector_ptrtype eigenvector;
799  // solve for eigenvalue problem at \p mu
800  SolverEigen<double>::eigenmodes_type modes;
801 
802  modes=
803  eigs( _matrixA=symmMatrix,
804  _matrixB=B,
805  _solver=( EigenSolverType )M_vm["crb.scm.solvereigen-solver-type"].template as<int>(),
806  _spectrum=SMALLEST_REAL,
807  //_spectrum=LARGEST_MAGNITUDE,
808  _transform=SINVERT,
809  _ncv=M_vm["crb.scm.solvereigen-ncv"].template as<int>(),
810  _nev=M_vm["crb.scm.solvereigen-nev"].template as<int>(),
811  _tolerance=M_vm["crb.scm.solvereigen-tol"].template as<double>(),
812  _maxit=M_vm["crb.scm.solvereigen-maxiter"].template as<int>()
813  );
814 
815  if ( modes.empty() )
816  {
817  LOG(INFO) << "eigs failed to converge\n";
818  return ckconv;
819  }
820 
821  LOG( INFO ) << "[fe eig] mu=" << std::setprecision( 4 ) << mu ;
822  LOG( INFO ) << "[fe eig] eigmin : " << std::setprecision( 16 ) << modes.begin()->second.template get<0>() ;
823  LOG( INFO ) << "[fe eig] ndof:" << M_model->functionSpace()->nDof() ;
824 
825  // extract the eigenvector associated with the smallest eigenvalue
826  eigenvector = modes.begin()->second.template get<2>();
827 
828  // store the eigenvalue associated with \p mu
829  M_C_eigenvalues[index] = modes.begin()->second.template get<0>();
830  typedef std::pair<size_type,value_type> key_t;
831 
832  //BOOST_FOREACH( key_t eig, M_C_eigenvalues )
833  //std::cout << "[fe eig] stored/map eig=" << eig.second <<" ( " << eig.first << " ) " << "\n";
834 
835  M_eig.push_back( M_C_eigenvalues[index] );
836  //BOOST_FOREACH( value_type eig, M_eig )
837  //std::cout << "[fe eig] stored/vec eig=" << eig << "\n";
838 
839 
840  /*
841  * now apply eigenvector to the Aq to compute
842  * y( eigenvector ) = ( a_q( eigenvector, eigenvector )/ ||eigenvector||^2
843  */
844  for ( size_type q = 0; q < nb_decomposition_terms_q() ; ++q )
845  {
846  for ( size_type m_eim = 0; m_eim < mMax(q) ; ++m_eim )
847  {
848  value_type aqmw = Matrixq[q][m_eim]->energy( eigenvector, eigenvector );
849  value_type bw = B->energy( eigenvector, eigenvector );
850  LOG( INFO ) << "[scm_offline] q=" << q << " aqmw = " << aqmw << ", bw = " << bw ;
851  M_Y_ub[K-1][ q ](m_eim) = aqmw/bw;
852  }
853 
854  //checkC( K );
855  //std::cout << "[scm_offline] M_Y_ub[" << K-1 << "]=" << M_Y_ub[K-1] << "\n";
856 
857  // save Y
858  os_y << M_Y_ub[K-1][q] << "\n";
859  }//q
860 
861 
862  if( M_scm_for_mass_matrix )
863  LOG( INFO ) <<"scm is done for mass matrix";
864  else
865  LOG( INFO )<<"scm is done for a( . , . ; mu )";
866 
867  double minerr, meanerr;
868  if( use_predefined_C )
869  {
870  relative_error = M_tolerance+10;
871  minerr = relative_error;
872  meanerr = relative_error;
873  if( K < M_iter_max )
874  {
875  mu = M_C->at( K );
876  index = index_vector[ K ];
877  }
878  }
879  else
880  boost::tie( relative_error, mu, index, minerr, meanerr ) = maxRelativeError( K );
881 #if 0
882  std::cout << " -- max relative error = " << relative_error
883  << " at mu = " << mu
884  << " at index = " << index << "\n";
885 #endif
886  //ofs << K << " " << std::setprecision(16) << relative_error << std::endl;
887  ckconv.push_back( boost::make_tuple( relative_error, minerr, meanerr ) );
888 
889  // could be that the max relative error is smaller than the tolerance if
890  // the coercivity constant is independant of the parameter set
891  if ( relative_error > M_tolerance && K < M_iter_max && ! use_predefined_C )
892  {
893  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
894  std::cout << " -- inserting mu - index : "<<index<<" - in C (" << M_C->size() << ")\n";
895  M_C->push_back( mu, index );
896 
897  //for ( size_type _i =0; _i < M_C->size(); ++_i )
898  //std::cout << " -- mu [" << _i << "]=" << M_C->at( _i ) << std::endl;
899 
900  M_C_complement = M_C->complement();
901 
902  //for ( size_type _i =0; _i < M_C_complement->size(); ++_i )
903  //std::cout << " -- mu complement [" << _i << "]=" << M_C_complement->at( _i ) << std::endl;
904  }
905 
906  ++K;
907  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
908  std::cout << "============================================================\n";
909  }
910 
911  //before call saveDB we have to split the vector of tuple M_y_bounds
912  M_y_bounds_0.resize( Qmax );
913  M_y_bounds_1.resize( Qmax );
914  for ( int q=0; q<Qmax; q++ )
915  {
916  for(int m=0; m<mMax(q); m++)
917  {
918  M_y_bounds_0[q].push_back( M_y_bounds[q][m].template get<0>() );
919  M_y_bounds_1[q].push_back( M_y_bounds[q][m].template get<1>() );
920  }
921  }
922 
923  saveDB();
924  return ckconv;
925 }
926 template<typename TruthModelType>
927 void
929 {
930  y_type err( M_Xi->size() );
931 
932  // check that for each mu in C_K the lb and ub are very close
933  for ( int k = 0; k < K; ++k )
934  {
935  std::cout << "[checkC] k= " << k << " K=" << K << " index in superSampling: " << M_C->indexInSuperSampling( k ) <<"\n";
936  parameter_type const& mu = M_C->at( k );
937  size_type index = M_C->indexInSuperSampling( k );
938  double _lb = lb( mu, K, index ).template get<0>();
939  double _lblb = lb( mu, std::max( K-1,size_type( 0 ) ), index ).template get<0>();
940 
941  if ( _lblb - _lb > 1e-10 )
942  {
943  LOG(INFO) << "[lberror] the lower bound is decreasing\n"
944  << "[lberror] _lblb = " << _lblb << "\n"
945  << "[lberror] _lb = " << _lb << "\n";
946  }
947 
948  double _ub = ub( mu, K ).template get<0>();
949  double _ex = M_C_eigenvalues.find( index )->second;
950  std::cout << "_lb = " << _lb << " _lblb=" << _lblb << " ub = " << _ub << " _ex = " << _ex << "\n";
951 
952  err( k ) = 1. - _lb/_ub;
953 
954  if ( std::abs( err( k ) ) > 1e-6 )
955  {
956  std::ostringstream ostr;
957  ostr << "error is not small as it should be " << k << " " << err( k ) << "\n";
958  ostr << "[checkC] relative error : " << mu << "\n"
959  << "[checkC] lb=" << std::setprecision( 16 ) << _lb << "\n"
960  << "[checkC] ub=" << std::setprecision( 16 ) << _ub << "\n"
961  << "[checkC] relerr(" << k << ")=" << std::setprecision( 16 ) << err( k ) << "\n"
962  << "[checkC] exact(" << k << ")=" << std::setprecision( 16 ) << _ex << "\n";
963  throw std::logic_error( ostr.str() );
964  }
965 
966 #if 1
967  std::cout << "[checkC] relative error : " << mu << "\n"
968  << "[checkC] lb=" << std::setprecision( 16 ) << _lb << "\n"
969  << "[checkC] ub=" << std::setprecision( 16 ) << _ub << "\n"
970  << "[checkC] relerr(" << k << ")=" << std::setprecision( 16 ) << err( k ) << "\n"
971  << "[checkC] exerr(" << k << ")=" << std::setprecision( 16 ) << _ex << "\n";
972 #endif
973  }
974 }
975 
976 
977 template<typename TruthModelType>
978 boost::tuple<typename CRBSCM<TruthModelType>::value_type,
979  typename CRBSCM<TruthModelType>::value_type>
980  CRBSCM<TruthModelType>::ex( parameter_type const& mu ) const
981 {
982  boost::timer ti;
983  sparse_matrix_ptrtype Matrix,symmMatrix;
984  sparse_matrix_ptrtype M ;
985  std::vector<vector_ptrtype> F;
986 
987  if ( M_scm_for_mass_matrix )
988  {
989  M = M_model->innerProductForMassMatrix();
990  boost::tie( Matrix, boost::tuples::ignore, F ) = M_model->update( mu );
991  }
992  else
993  {
994  M = M_model->innerProduct();
995  boost::tie( boost::tuples::ignore, Matrix, F ) = M_model->update( mu );
996  }
997 
998  symmMatrix = M_model->newMatrix();
999  symmMatrix->close();
1000  Matrix->symmetricPart( symmMatrix );
1001 
1002  SolverEigen<double>::eigenmodes_type modesmin=
1003  eigs( _matrixA=Matrix,
1004  _matrixB=M,
1005  _solver=( EigenSolverType )M_vm["crb.scm.solvereigen-solver-type"].template as<int>(),
1006  //_spectrum=LARGEST_MAGNITUDE,
1007  //_spectrum=LARGEST_REAL,
1008  _spectrum=SMALLEST_REAL,
1009  //_spectrum=SMALLEST_MAGNITUDE,
1010  _transform=SINVERT,
1011  _ncv=M_vm["crb.scm.solvereigen-ncv"].template as<int>(),
1012  _nev=M_vm["crb.scm.solvereigen-nev"].template as<int>(),
1013  _tolerance=M_vm["crb.scm.solvereigen-tol"].template as<double>(),
1014  _maxit=M_vm["crb.scm.solvereigen-maxiter"].template as<int>()
1015  );
1016 
1017  if ( modesmin.empty() )
1018  {
1019  std::cout << "no eigenmode converged: increase --solvereigen-ncv\n";
1020  return 0.;
1021  }
1022 
1023  double eigmin = modesmin.begin()->second.template get<0>();
1024  //std::cout << std::setprecision(4) << mu << " eigmin = "
1025  // << std::setprecision(16) << eigmin << "\n";
1026 
1027  return boost::make_tuple( eigmin, ti.elapsed() );
1028 #if 0
1029  SolverEigen<double>::eigenmodes_type modesmax=
1030  eigs( _matrixA=Matrix,
1031  _matrixB=M,
1032  _solver=( EigenSolverType )M_vm["crb.scm.solvereigen-solver-type"].as<int>(),
1033  _spectrum=LARGEST_MAGNITUDE,
1034  _ncv=M_vm["crb.scm.solvereigen-ncv"].as<int>(),
1035  _nev=M_vm["crb.scm.solvereigen-nev"].as<int>(),
1036  _tolerance=M_vm["crb.scm.solvereigen-tol"].as<double>(),
1037  _maxit=M_vm["crb.scm.solvereigen-maxiter"].as<int>()
1038  );
1039 
1040  if ( modesmax.empty() )
1041  {
1042  std::cout << "no eigenmode converged: increase --solvereigen-ncv\n";
1043  return;
1044  }
1045 
1046  double eigmax = modesmax.rbegin()->second.template get<0>();
1047 
1048  //std::cout << "[crbmodel] " << std::setprecision(4) << mu << " "
1049  // << std::setprecision(16) << eigmin << " " << eigmax
1050  // << " " << Xh->nDof() << "\n";
1051 
1052 #endif
1053 
1054 }
1055 
1056 
1057 
1058 template<typename TruthModelType>
1059 boost::tuple<typename CRBSCM<TruthModelType>::value_type, double>
1060 CRBSCM<TruthModelType>::lb( parameter_type const& mu ,size_type K ,int indexmu ) const
1061 {
1062  if( M_use_scm )
1063  return this->lbSCM( mu, K , indexmu );
1064  else
1065  return this->lbNoSCM( mu, K , indexmu );
1066 }
1067 
1068 template<typename TruthModelType>
1069 boost::tuple<typename CRBSCM<TruthModelType>::value_type, double>
1070 CRBSCM<TruthModelType>::lbNoSCM( parameter_type const& mu ,size_type K ,int indexmu ) const
1071 {
1072 
1073  boost::mpi::timer ti;
1074 
1075  vector_type vec_min_coeff;
1076  auto all_beta = M_model->computeBetaQm( mu );
1077  auto all_beta_ref = M_model->computeBetaQm( M_mu_ref );
1078  beta_vector_type beta_mu;
1079  beta_vector_type beta_mu_ref;
1080  if ( M_scm_for_mass_matrix )
1081  {
1082  beta_mu = all_beta.template get<0>();
1083  beta_mu_ref = all_beta_ref.template get<0>();
1084  }
1085  else
1086  {
1087  beta_mu = all_beta.template get<1>();
1088  beta_mu_ref = all_beta_ref.template get<1>();
1089  }
1090 
1091  int Q = this->nb_decomposition_terms_q();
1092  vec_min_coeff.resize(Q);
1093  for(int q=0; q<Q; q++)
1094  {
1095  //for each q, search the min coeff
1096  vector_type vec_local_min_coeff;
1097  int M = beta_mu[q].size();
1098  vec_local_min_coeff.resize( M );
1099  for(int m=0; m<M; m++)
1100  {
1101  double __beta = beta_mu[q][m]/beta_mu_ref[q][m];
1102  vec_local_min_coeff(m) = __beta ;
1103  }
1104  double local_min = vec_local_min_coeff.minCoeff();
1105  vec_min_coeff(q) = local_min ;
1106  }
1107  double min = vec_min_coeff.minCoeff();
1108  double lower_bound = min * M_C_eigenvalues.find(0)->second;
1109 
1110  return boost::make_tuple( lower_bound, ti.elapsed() );
1111 }
1112 
1113 template<typename TruthModelType>
1114 boost::tuple<typename CRBSCM<TruthModelType>::value_type, double>
1115 CRBSCM<TruthModelType>::lbSCM( parameter_type const& mu ,size_type K ,int indexmu ) const
1116 {
1117 
1118 #if defined(FEELPP_HAS_GLPK_H)
1119  if ( K == invalid_size_type_value ) K = this->KMax();
1120 
1121  if ( K > this->KMax() ) K = this->KMax();
1122 
1123  boost::mpi::timer ti;
1124 
1125  // value if K==0
1126  if ( K <= 0 ) return 0.0;
1127 
1128  if ( indexmu >= 0 && ( M_C_alpha_lb[indexmu].find( K ) != M_C_alpha_lb[indexmu].end() ) )
1129  return M_C_alpha_lb[indexmu][K] ;
1130 
1131  //if ( K == std::max(size_type(0),M_C->size()-(size_type)M_vm["crb.scm.level"].template as<int>() ) ) return 0.0;
1132  //int level = M_vm["crb.scm.level"].template as<int>();
1133 
1134  //std::cout << "[CRBSCM::lb] Alphalb size " << M_C_alpha_lb.size() << "\n";
1135 
1136  beta_vector_type beta_qm;
1137  glp_prob *lp;
1138 
1139 
1140  lp = glp_create_prob();
1141  glp_set_prob_name( lp, ( boost::format( "scm_%1%" ) % K ).str().c_str() );
1142  glp_set_obj_dir( lp, GLP_MIN );
1143 
1144  int Malpha = std::min( M_Malpha,std::min( K,M_Xi->size() ) );
1145 
1146  int Mplus = std::min( M_Mplus,M_Xi->size()-K );
1147 
1148 
1149  if ( M_vm["crb.scm.strategy"].template as<int>()==2 )
1150  Mplus = std::min( M_Mplus,std::min( K, M_Xi->size()-K ) );
1151 
1152  // we have exactely Qa*(M+ + Malpha) entries in the matrix
1153 
1154 
1155  int nnz = nb_decomposition_terms_qm()*( Malpha+Mplus );
1156 
1157  int ia[1+10000], ja[1+10000];
1158  double ar[1+10000];
1159  int nnz_index = 1;
1160 
1161  // set the auxiliary variables: we have first Malpha of them from C_K and
1162  // Mplus from its complement
1163  glp_add_rows( lp,Malpha+Mplus );
1164  // search the the M_Malpha closest points in C_K, M_Malpha must be < K
1165  // index_vector will contain index of neighbors
1166  std::vector<int> index_vector;
1167  sampling_ptrtype C_neighbors = M_C->searchNearestNeighbors( mu, Malpha , index_vector);
1168 
1169  //std::cout << "[CRBSCM::lb] C_neighbors size = " << C_neighbors->size() << "\n";
1170 
1171  //std::cout << "[CRBSCM::lb] add rows associated with C_K\n";
1172  // first the constraints associated with C_K: make sure that Malpha is not
1173  // greater by taking the min of M_Malpha and K
1174 
1175  for ( int m=0; m < Malpha; ++m )
1176  {
1177 
1178  //std::cout << "[CRBSCM::lb] add row " << m << " from C_K\n";
1179  parameter_type mup = C_neighbors->at( m );
1180 
1181  // update the theta_q associated with mup
1182  if ( M_scm_for_mass_matrix )
1183  {
1184  boost::tie( beta_qm , boost::tuples::ignore , boost::tuples::ignore ) = M_model->computeBetaQm( mup );
1185  }
1186  else
1187  {
1188  boost::tie( boost::tuples::ignore, beta_qm, boost::tuples::ignore ) = M_model->computeBetaQm( mup );
1189  }
1190 
1191  //std::cout << "[CRBSCM::lb] thetaq = " << theta_q << "\n";
1192 
1193  //std::cout << "[CRBSCM::lb] row name : " << (boost::format( "c_%1%_%2%" ) % K % m).str() << "\n";
1194  glp_set_row_name( lp, m+1, ( boost::format( "c_%1%_%2%" ) % K % m ).str().c_str() );
1195 
1196  //std::cout << "[CRBSCM::lb] row bounds\n";
1197  //std::cout << "[CRBSCM::lb] index in super sampling: " << C_neighbors->indexInSuperSampling( m ) << "\n";
1198  //std::cout << "[CRBSCM::lb] eig value: " << M_C_eigenvalues.find( C_neighbors->indexInSuperSampling( m ) )->second << "\n";
1199  glp_set_row_bnds( lp,
1200  m+1,
1201  GLP_LO,
1202  M_C_eigenvalues.find( C_neighbors->indexInSuperSampling( m ) )->second,
1203  0.0 );
1204 
1205 
1206  //std::cout << "[CRBSCM::lb] constraints matrix\n";
1207  int count=1;
1208  for ( int q = 0; q < nb_decomposition_terms_q(); ++q )
1209  {
1210  for ( int m_eim = 0; m_eim < mMax( q ); ++m_eim, ++nnz_index, ++count )
1211  {
1212  //std::cout << "[CRBSCM::lb] constraints matrix q,m = " << q << ","<< m <<"\n";
1213  ia[nnz_index]=m+1;
1214  ja[nnz_index]=count;
1215  ar[nnz_index]=beta_qm[q][m_eim];
1216  }
1217  }
1218 
1219  }
1220  //std::cout << "[CRBSCM::lb] add rows associated with C_K done. nnz=" << nnz_index << "\n";
1221 
1222  // search the the Mplus closest points in Xi\C_K
1223  sampling_ptrtype Xi_C_neighbors = M_C_complement->searchNearestNeighbors( mu, Mplus , index_vector);
1224 
1225  //std::cout << "[CRBSCM::lb] C_complement size = " << Xi_C_neighbors->size() << "\n";
1226 
1227  //std::cout << "[CRBSCM::lb] add rows associated with Xi\\C_K\n";
1228  //std::cout << "[CRBSCM::lb] Mplus =" << Mplus << " , neighbors= " << Xi_C_neighbors->size() << "\n";
1229  // second the monotonicity constraints associated with the complement of C_K
1230  for ( int m=0; m < Mplus; ++m )
1231  //for( int m=0;m < Xi_C_neighbors->size(); ++m )
1232  {
1233 
1234  //std::cout << "[CRBSCM::lb] add row " << m << " from Xi\\C_K\n";
1235  parameter_type mup = Xi_C_neighbors->at( m );
1236 
1237  double _lb;
1238 
1239 #if 1
1240  //std::cout << "[CRBSCM::lb] Entering find alphamap\n" << Xi_C_neighbors->indexInSuperSampling( m ) << std::endl ;
1241  //std::cout << "[CRBSCM::lb] Size of alphalb " << M_C_alpha_lb.find( Xi_C_neighbors->indexInSuperSampling( m ) )->second.size() << "\n" ;
1242 
1243  if ( M_C_alpha_lb[Xi_C_neighbors->indexInSuperSampling( m ) ].find( K-1 ) != M_C_alpha_lb[Xi_C_neighbors->indexInSuperSampling( m ) ].end() )
1244  {
1245  //std::cout << "[CRBSCM::lb] lb déjŕ calculée\n" ;
1246  }
1247  else
1248  {
1249  //std::cout << "[CRBSCM::lb] Nouvelle lb\n" ;
1250  M_C_alpha_lb[ Xi_C_neighbors->indexInSuperSampling( m )][K-1] = lb( mup, K-1, Xi_C_neighbors->indexInSuperSampling( m ) ).template get<0>();
1251  }
1252 
1253  _lb = M_C_alpha_lb[ Xi_C_neighbors->indexInSuperSampling( m )][K-1];
1254 #endif
1255 
1256  // update the theta_q associated with mup
1257  if ( M_scm_for_mass_matrix )
1258  {
1259  boost::tie( beta_qm, boost::tuples::ignore, boost::tuples::ignore ) = M_model->computeBetaQm( mup );
1260  }
1261 
1262  else
1263  {
1264  boost::tie( boost::tuples::ignore, beta_qm, boost::tuples::ignore ) = M_model->computeBetaQm( mup );
1265  }
1266 
1267  glp_set_row_name( lp, Malpha+m+1, ( boost::format( "xi_c_%1%_%2%" ) % K % m ).str().c_str() );
1268 
1269  switch ( M_vm["crb.scm.strategy"].template as<int>() )
1270  {
1271  // Patera
1272  case 0:
1273  {
1274  //std::cout << "Patera\n" ;
1275  glp_set_row_bnds( lp, Malpha+m+1, GLP_LO, 0.0, 0.0 );
1276  }
1277  break;
1278  // Maday
1279 
1280  case 1:
1281 
1282  // Prud'homme
1283  case 2:
1284  {
1285  //std::cout << "Prud'homme\n" ;
1286  // std::cout << "Lb = " << _lb << std::endl ;
1287  glp_set_row_bnds( lp, Malpha+m+1, GLP_LO, _lb, 0.0 );
1288  }
1289  break;
1290  }
1291 
1292  int count=1;
1293  int nb_already_done = nb_decomposition_terms_qm();
1294  for ( int q = 0; q < nb_decomposition_terms_q(); ++q )
1295  {
1296  for ( int m_eim = 0; m_eim < this->mMax(q); ++m_eim, ++nnz_index,++count )
1297  {
1298  //std::cout << "[CRBSCM::lb] constraints matrix q = " << q <<" , "<<m<< "\n";
1299  ia[nnz_index]=Malpha+m+1;
1300  ja[nnz_index]=count;
1301  ar[nnz_index]=beta_qm[q][m_eim];
1302  }
1303  }
1304  }
1305 
1306  //std::cout << "[CRBSCM::lb] add rows associated with C_K complement done. nnz=" << nnz_index << "\n";
1307 
1308  // set the structural variables, we have M_model->Qa() of them
1309  if ( M_scm_for_mass_matrix )
1310  {
1311  boost::tie( beta_qm,boost::tuples::ignore, boost::tuples::ignore ) = M_model->computeBetaQm( mu );
1312  }
1313  else
1314  {
1315  boost::tie( boost::tuples::ignore, beta_qm, boost::tuples::ignore ) = M_model->computeBetaQm( mu );
1316  }
1317 
1318  //nb_columns
1319 #if 0
1320  int nb_col = 0;
1321  for(int q=0; q<nb_decomposition_terms_q(); q++)
1322  {
1323  for(int m_eim=0; m_eim<this->mMax(q); m_eim++)
1324  nb_col++;
1325  }
1326 #endif
1327  // glp_add_cols( lp, nb_col );
1328  glp_add_cols( lp, nb_decomposition_terms_qm() );
1329  int count=0;
1330  for ( int q = 0; q < nb_decomposition_terms_q(); ++q )
1331  {
1332  for(int m_eim=0; m_eim<this->mMax(q); ++m_eim,++count)
1333  {
1334  glp_set_col_name( lp, count+1, ( boost::format( "y_%1%-%2%" ) % q %m_eim ).str().c_str() );
1335  glp_set_col_bnds( lp, count+1, GLP_DB,
1336  M_y_bounds[q][m_eim].template get<0>(),
1337  M_y_bounds[q][m_eim].template get<1>() );
1338  glp_set_obj_coef( lp, count+1, beta_qm[q][m_eim] );
1339  }
1340  }
1341 
1342  // now we need to build the matrix
1343  glp_load_matrix( lp, nnz, ia, ja, ar );
1344  glp_smcp parm;
1345  glp_init_smcp( &parm );
1346  parm.msg_lev = GLP_MSG_ERR;
1347 
1348  // use the simplex method and solve the LP
1349  glp_simplex( lp, &parm );
1350 
1351  // retrieve the minimum
1352  double Jobj = glp_get_obj_val( lp );
1353  //std::cout << "Jobj = " << Jobj << "\n";
1354 #if 0
1355 
1356  for ( size_type q = 0; q < M_model->Qa(); ++q )
1357  {
1358  double y = glp_get_col_prim( lp,q+1 );
1359  //std::cout << "y" << q << " = " << y << "\n";
1360  }
1361 
1362 #endif
1363  glp_delete_prob( lp );
1364 
1365  if ( indexmu >= 0 )
1366  {
1367 
1368  M_C_alpha_lb[ indexmu ][K] = Jobj;
1369 
1370  }
1371 #else // FEELPP_HAS_GLPK_H
1372 
1373  double Jobj = 0;
1374  boost::timer ti;
1375 #endif /* FEELPP_HAS_GLPK_H */
1376 
1377  return boost::make_tuple( Jobj, ti.elapsed() );
1378 
1379 }
1380 
1381 
1382 template<typename TruthModelType>
1383 boost::tuple<typename CRBSCM<TruthModelType>::value_type, double>
1385 {
1386  if ( K == invalid_size_type_value ) K = this->KMax();
1387 
1388  if ( K > this->KMax() ) K = this->KMax();
1389 
1390  boost::timer ti;
1391  beta_vector_type beta_qm;
1392 
1393  if ( M_scm_for_mass_matrix )
1394  {
1395  boost::tie( beta_qm ,boost::tuples::ignore, boost::tuples::ignore ) = M_model->computeBetaQm( mu );
1396  }
1397 
1398  else
1399  {
1400  boost::tie( boost::tuples::ignore, beta_qm, boost::tuples::ignore ) = M_model->computeBetaQm( mu );
1401  }
1402 
1403  //std::cout << "[CRBSCM<TruthModelType>::ub] theta_q = " << theta_q << "\n";
1404  y_type y( K );
1405 
1406  //for ( size_type k = 0; k < K; ++k )
1407  //{
1408  // y( k ) = ( beta_qm.array()*M_Y_ub[k].array() ).sum();
1409  //}
1410  for ( size_type k = 0; k < K; ++k )
1411  {
1412  y( k ) = 0;
1413  for(int q=0; q<nb_decomposition_terms_q(); q++)
1414  {
1415  //y( k ) = ( beta_qm[q].array()*M_Y_ub[k][q].array() ).sum();
1416  for(int m=0; m<mMax(q); m++)
1417  y( k ) += beta_qm[q][m]*M_Y_ub[k][q][m];
1418  }
1419  }
1420 
1421  //std::cout << "[CRBSCM<TruthModelType>::ub] y = " << y << "\n";
1422  return boost::make_tuple( y.minCoeff(), ti.elapsed() );
1423 
1424 }
1425 
1426 template<typename TruthModelType>
1427 typename CRBSCM<TruthModelType>::relative_error_type
1429 {
1430  //std::cout << "==================================================\n";
1431 #if 0
1432  y_type err( M_Xi->size() );
1433 
1434  for ( size_type k = 0; k < M_Xi->size(); ++k )
1435  {
1436  parameter_type const& mu = M_Xi->at( k );
1437 #else
1438  y_type err( M_C_complement->size() );
1439 
1440  for ( size_type k = 0; k < M_C_complement->size(); ++k )
1441  {
1442  parameter_type const& mu = M_C_complement->at( k );
1443 #endif
1444 
1445  //std::cout << "[maxRelativeError] Calcul de lb pour mu[" << k << "]\n" ;
1446  double _lb = lb( mu, K, k ).template get<0>();
1447  //std::cout << "[maxRelativeError] Calcul de lblb pour mu[" << k << "]\n" ;
1448  double _lblb = lb( mu, std::max( K-1,size_type( 0 ) ), k ).template get<0>();
1449 
1450  if ( _lblb - _lb > 1e-10 )
1451  {
1452  //LOG(INFO) << "[lberror] the lower bound is decreasing\n"
1453  // << "[lberror] _lblb = " << _lblb << "\n"
1454  // << "[lberror] _lb = " << _lb << "\n";
1455  }
1456 
1457  double _ub = ub( mu, K ).template get<0>();
1458 
1459  err( k ) = 1. - _lb/_ub;
1460 #if 0
1461  //std::cout << "[maxRelativeError] k = " << k << "\n";
1462  //std::cout << "[maxRelativeError] parameter : " << mu << "\n";
1463  //std::cout << "[maxRelativeError] lb = " << std::setprecision(16) << _lb << "\n";
1464  //std::cout << "[maxRelativeError] ub = " << std::setprecision(16) << _ub << "\n";
1465  //std::cout << "[maxRelativeError] rel err(" << k << ")=" << std::setprecision(16) << err(k) << "\n";
1466 #endif
1467  }
1468 
1469 
1470  Eigen::MatrixXf::Index index;
1471  double maxerr = err.array().abs().maxCoeff( &index );
1472  //std::cout << "[maxRelativeError] K=" << K << " max Error = " << maxerr << " at index = " << index << "\n";
1473  Eigen::MatrixXf::Index indexmin;
1474  double minerr = err.array().abs().minCoeff( &indexmin );
1475  //std::cout << "[maxRelativeError] K=" << K << " min Error = " << minerr << " at index = " << indexmin << "\n";
1476  double meanerr = err.array().abs().sum()/err.size();
1477  //std::cout << "[maxRelativeError] K=" << K << " mean Error = " << meanerr << "\n";
1478  //std::cout << "==================================================\n";
1479  //return boost::make_tuple( maxerr, M_Xi->at( index ), index, minerr, meanerr );
1480 
1481 #if 0
1482  parameter_type const& mu = M_C_complement->at( index );
1483  double _lb = lb( mu, K, index ).template get<0>();
1484  double _ub = ub( mu, K ).template get<0>();
1485  std::cout<<"_lb = "<<_lb<<" and _ub = "<<_ub<<std::endl;
1486 #endif
1487 
1488  return boost::make_tuple( maxerr, M_C_complement->at( index ), M_C_complement->indexInSuperSampling( index ), minerr, meanerr );
1489 }
1490 
1491 template<typename TruthModelType>
1492 void
1494 {
1495  //std::cout << "************************************************************\n";
1496  LOG(INFO) << "[CRBSCM<TruthModelType>::computeYBounds()] start...\n";
1497  int Qmax = nb_decomposition_terms_q();
1498  M_y_bounds.resize(Qmax);
1499  sparse_matrix_ptrtype Matrix, symmMatrix=M_model->newMatrix(), B;
1500 
1501 
1502  if ( M_scm_for_mass_matrix )
1503  B=M_model->innerProductForMassMatrix();
1504  else
1505  B=M_model->innerProduct();
1506 
1507  B->close();
1508 
1509 
1510  // solve 2 * Q_a eigenproblems
1511  for ( int q = 0; q < nb_decomposition_terms_q() ; ++q )
1512  {
1513 
1514  for ( int m = 0; m < mMax(q) ; ++m )
1515  {
1516  //std::cout << "================================================================================\n";
1517  //std::cout << "[ComputeYBounds] = q = " << q << " / " << M_model->Qa() << "\n";
1518  LOG(INFO) << "[CRBSCM<TruthModelType>::computeYBounds()] q = " << q << "/" << M_model->Qa() << "\n";
1519  // for a given parameter \p mu assemble the left and right hand side
1520  std::ostringstream os;
1521 
1522  if ( M_scm_for_mass_matrix )
1523  Matrix = M_model->Mqm( q , m );
1524  else
1525  Matrix = M_model->Aqm( q , m );
1526 
1527  Matrix->close();
1528  Matrix->symmetricPart( symmMatrix );
1529  if( Environment::worldComm().globalSize() == 1 )
1530  {
1531  os << "yb_Matrix" << q << " - "<< m << ".m";
1532  Matrix->printMatlab( os.str() );
1533  os.str( "" );
1534  os << "yb_symmMatrix" << q << " - "<< m <<".m";
1535  symmMatrix->printMatlab( os.str() );
1536  }
1537  if ( symmMatrix->l1Norm()==0.0 )
1538  {
1539  std::cout << "matrix is null\n" ;
1540  M_y_bounds[q].push_back( boost::make_tuple( 0.0, 1e-10 ) );
1541  }
1542 
1543  else
1544  {
1545 
1546 #if 0
1547  // solve for eigenvalue problem at \p mu
1548  boost::tie( nconv, eigenvalue_lb, boost::tuples::ignore, boost::tuples::ignore ) =
1549  eigs( _matrixA=symmA,
1550  _matrixB=B,
1551  _spectrum=SMALLEST_MAGNITUDE,
1552  _ncv=15 );
1553 
1554  //std::cout << " -- lower bounds q=" << q
1555  // << " nconv=" << nconv
1556  // << " eigenvalue_lb = " << eigenvalue_lb << "\n";
1557  double eigmin=eigenvalue_lb;
1558 
1559  boost::tie( nconv, eigenvalue_ub, boost::tuples::ignore, boost::tuples::ignore ) =
1560  eigs( _matrixA=symmA,
1561  _matrixB=B,
1562  _spectrum=LARGEST_MAGNITUDE );
1563 
1564  //std::cout << " -- upper bounds q=" << q
1565  // << " nconv=" << nconv
1566  // << " eigenvalue_ub = " << eigenvalue_ub << "\n";
1567  double eigmax=eigenvalue_ub;
1568 #else
1569 
1570  SolverEigen<double>::eigenmodes_type modes;
1571 #if 1
1572  // solve for eigenvalue problem at \p mu
1573  modes =
1574  eigs( _matrixA=symmMatrix,
1575  _matrixB=B,
1576  //_problem=(EigenProblemType)PGNHEP,
1577  _problem=( EigenProblemType )GHEP,
1578  _solver=( EigenSolverType )M_vm["crb.scm.solvereigen-solver-type"].template as<int>(),
1579  //_spectrum=SMALLEST_REAL,
1580  _spectrum=SMALLEST_MAGNITUDE,
1581  _ncv=M_vm["crb.scm.solvereigen-ncv"].template as<int>(),
1582  _nev=M_vm["crb.scm.solvereigen-nev"].template as<int>(),
1583  _tolerance=M_vm["crb.scm.solvereigen-tol"].template as<double>(),
1584  _maxit=M_vm["crb.scm.solvereigen-maxiter"].template as<int>()
1585  );
1586 
1587 #endif
1588 
1589  if ( modes.empty() )
1590  {
1591  LOG(INFO) << "[Computeybounds] eigmin did not converge for q=" << q << " (set to 0)\n";
1592  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1593  std::cout << "[Computeybounds] eigmin did not converge for q=" << q << " (set to 0)"<<std::endl;
1594  }
1595 
1596  double eigmin = modes.empty()?0:modes.begin()->second.template get<0>();
1597 #if 1
1598  modes=
1599  eigs( _matrixA=symmMatrix,
1600  _matrixB=B,
1601  //_problem=(EigenProblemType)PGNHEP,
1602  _problem=( EigenProblemType )GHEP,
1603  _solver=( EigenSolverType )M_vm["crb.scm.solvereigen-solver-type"].template as<int>(),
1604  _spectrum=LARGEST_REAL,
1605  //_spectrum=LARGEST_MAGNITUDE,
1606  _ncv=M_vm["crb.scm.solvereigen-ncv"].template as<int>(),
1607  //_ncv=20,
1608  _nev=M_vm["crb.scm.solvereigen-nev"].template as<int>(),
1609  _tolerance=M_vm["crb.scm.solvereigen-tol"].template as<double>(),
1610  //_tolerance=1e-7,
1611  _maxit=M_vm["crb.scm.solvereigen-maxiter"].template as<int>()
1612  );
1613 
1614 #endif
1615 
1616  if ( modes.empty() )
1617  {
1618  LOG(INFO) << "[Computeybounds] eigmax did not converge for q=" << q << " (set to 0)\n";
1619  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1620  std::cout << "[Computeybounds] eigmax did not converge for q=" << q << " (set to 0)"<<std::endl;
1621  }
1622 
1623  double eigmax = modes.empty()?0:modes.rbegin()->second.template get<0>();
1624  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1625  std::cout<<"[computeYBounds] bounds for (q,m) = ("<<q<<","<<m<<") [ "<<eigmin<<" ; "<<eigmax<<"]"<<std::endl;
1626  LOG(INFO)<<"[computeYBounds] bounds for (q,m) = ("<<q<<","<<m<<") [ "<<eigmin<<" ; "<<eigmax<<"]\n";
1627  //std::cout << "[Computeybounds] q= " << q << " eigmin=" << std::setprecision(16) << eigmin << " eigmax=" << std::setprecision(16) << eigmax << "\n";
1628  //std::cout << std::setprecision(16) << eigmin << " " << eigmax << "\n";
1629 
1630 #endif
1631 
1632  if ( eigmin==0.0 && eigmax==0.0 ) throw std::logic_error( "eigs null\n" );
1633  M_y_bounds[q].push_back( boost::make_tuple( eigmin, eigmax ) );
1634  }//m
1635  }//q
1636  }
1637 
1638  LOG(INFO) << "[CRBSCM<TruthModelType>::computeYBounds()] stop.\n";
1639  //std::cout << "************************************************************\n";
1640 }
1641 
1642 template<typename TruthModelType>
1643 std::vector<double>
1645 {
1646  std::cout << "------------------------------------------------------------\n";
1647  double alpha_lb,alpha_lbti;
1648  boost::tie( alpha_lb, alpha_lbti ) = this->lb( mu, K );
1649  double alpha_ub,alpha_ubti;
1650  if( M_use_scm )
1651  boost::tie( alpha_ub, alpha_ubti ) = this->ub( mu, K );
1652  else
1653  {
1654  alpha_ub = alpha_lb;
1655  alpha_ubti = alpha_lbti;
1656  }
1657  double alpha_ex, alpha_exti;
1658  boost::tie( alpha_ex, alpha_exti ) = this->ex( mu );
1659  LOG( INFO ) << "alpha_lb=" << alpha_lb << " alpha_ub=" << alpha_ub << " alpha_ex=" << alpha_ex << "\n";
1660  LOG( INFO ) << ( alpha_ex-alpha_lb )/( alpha_ub-alpha_lb ) << "\n";
1661  LOG( INFO ) << K << " "
1662  << std::setprecision( 16 ) << alpha_lb << " "
1663  << std::setprecision( 3 ) << alpha_lbti << " "
1664  << std::setprecision( 16 ) << alpha_ub << " "
1665  << std::setprecision( 3 ) << alpha_ubti << " "
1666  << std::setprecision( 16 ) << alpha_ex << " "
1667  << std::setprecision( 16 ) << alpha_exti << " "
1668  << std::setprecision( 16 ) << ( alpha_ub-alpha_lb )/( alpha_ub ) << " "
1669  << std::setprecision( 16 ) << ( alpha_ex-alpha_lb )/( alpha_ex ) << " "
1670  << std::setprecision( 16 ) << ( alpha_ub-alpha_ex )/( alpha_ex ) << " "
1671  ;
1672  LOG( INFO ) << "------------------------------------------------------------\n";
1673  double rel_diff = (alpha_ex - alpha_lb)/alpha_ex;
1674  return boost::assign::list_of( alpha_lb )( alpha_lbti )( alpha_ub )( alpha_ubti )( alpha_ex )( alpha_exti )( rel_diff );
1675 }
1676 
1677 template<typename TruthModelType>
1678 void
1679 CRBSCM<TruthModelType>::run( const double * X, unsigned long N, double * Y, unsigned long P )
1680 {
1681  parameter_type mu( M_Dmu );
1682 
1683  // the last parameter is the max error
1684  for ( unsigned long p= 0; p < N-3; ++p )
1685  mu( p ) = X[p];
1686 
1687  for ( unsigned long i=0; i<N; i++ ) std::cout<<"X["<<i<<"] = "<<X[i]<<std::endl;
1688 
1689  double meshSize = X[N-3];
1690  //M_model->setMeshSize( meshSize );
1691 
1692  size_type K = this->KMax();
1693  double alpha_lb,lbti;
1694  boost::tie( alpha_lb, lbti ) = this->lb( mu, K );
1695  double alpha_ub,ubti;
1696  boost::tie( alpha_ub, ubti ) = this->ub( mu, K );
1697  double alpha_ex, alpha_exti;
1698  boost::tie( alpha_ex, alpha_exti ) = this->ex( mu );
1699  std::cout << "lb=" << alpha_lb << " ub=" << alpha_ub << " ex=" << alpha_ex << "\n";
1700  std::cout << ( alpha_ex-alpha_lb )/( alpha_ub-alpha_lb ) << "\n";
1701  Y[0] = alpha_lb;
1702  Y[1] = alpha_ub;
1703 
1704 }
1705 
1706 template<typename TruthModelType>
1707 template<class Archive>
1708 void
1709 CRBSCM<TruthModelType>::save( Archive & ar, const unsigned int version ) const
1710 {
1711  if( M_use_scm )
1712  {
1713  ar & M_use_scm;
1714  ar & M_Malpha;
1715  ar & M_Mplus;
1716  ar & M_C_alpha_lb;
1717  ar & M_C;
1718  ar & M_C_complement;
1719  ar & M_C_eigenvalues;
1720  ar & M_y_bounds_0;
1721  ar & M_y_bounds_1;
1722  ar & M_Y_ub;
1723  ar & M_Xi;
1724  }
1725  else
1726  {
1727  ar & M_use_scm;
1728  ar & M_mu_ref;
1729  ar & M_C_eigenvalues;
1730  }
1731 }
1732 
1733 template<typename TruthModelType>
1734 template<class Archive>
1735 void
1736 CRBSCM<TruthModelType>::load( Archive & ar, const unsigned int version )
1737 {
1738 
1739  ar & M_use_scm ;
1740  bool use_scm = option(_name="crb.scm.use-scm").template as<bool>() ;
1741  bool rebuild = option(_name="crb.scm.rebuild-database").template as<bool>() ;
1742  if( M_use_scm != use_scm && rebuild==false)
1743  {
1744  if( use_scm )
1745  throw std::logic_error( "[CRBSCM::load] ERROR the database created is not appropriate to use SCM. Use option crb.scm.rebuild-database=true");
1746  else
1747  throw std::logic_error( "[CRBSCM::load] ERROR the database was created to use SCM. Use option crb.scm.rebuild-database=true");
1748  }
1749 
1750  if( M_use_scm )
1751  {
1752  ar & M_Malpha;
1753  ar & M_Mplus;
1754  ar & M_C_alpha_lb;
1755  ar & M_C;
1756  ar & M_C_complement;
1757  ar & M_C_eigenvalues;
1758  ar & M_y_bounds_0;
1759  ar & M_y_bounds_1;
1760  ar & M_Y_ub;
1761  ar & M_Xi;
1762  int Qmax = this->nb_decomposition_terms_q();
1763  M_y_bounds.resize( Qmax );
1764  for ( int q=0; q<Qmax; q++ )
1765  {
1766  for(int m=0; m<mMax(q); m++)
1767  M_y_bounds[q].push_back( boost::make_tuple( M_y_bounds_0[q][m] , M_y_bounds_1[q][m] ) );
1768  }
1769  }
1770  else
1771  {
1772  ar & M_mu_ref;
1773  ar & M_C_eigenvalues;
1774  }
1775 
1776 }
1777 
1778 template<typename TruthModelType>
1779 bool
1780 CRBSCM<TruthModelType>::doScmForMassMatrix()
1781 {
1782  bool b = this->vm()["crb.scm.do-scm-for-mass-matrix"].template as<bool>();
1783  return b;
1784 }
1785 
1786 template<typename TruthModelType>
1787 bool
1788 CRBSCM<TruthModelType>::rebuildDB()
1789 {
1790  bool rebuild = this->vm()["crb.scm.rebuild-database"].template as<bool>();
1791  return rebuild;
1792 }
1793 
1794 template<typename TruthModelType>
1795 void
1797 {
1798  fs::ofstream ofs( this->dbLocalPath() / this->dbFilename() );
1799 
1800  if ( ofs )
1801  {
1802  boost::archive::text_oarchive oa( ofs );
1803  // write class instance to archive
1804  oa << *this;
1805  // archive and stream closed when destructors are called
1806  }
1807 }
1808 template<typename TruthModelType>
1809 bool
1811 {
1812  if ( this->rebuildDB() )
1813  return false;
1814 
1815  fs::path db = this->lookForDB();
1816 
1817  if ( db.empty() )
1818  return false;
1819 
1820  if ( !fs::exists( db ) )
1821  return false;
1822 
1823  std::cout << "Loading " << db << "...\n";
1824  fs::ifstream ifs( db );
1825  if ( ifs )
1826  {
1827  boost::archive::text_iarchive ia( ifs );
1828  // write class instance to archive
1829  ia >> *this;
1830  std::cout << "Loading " << db << " done...\n";
1831  this->setIsLoaded( true );
1832  // archive and stream closed when destructors are called
1833  return true;
1834  }
1835 
1836  return false;
1837 }
1838 
1839 } // Feel
1840 
1841 namespace boost
1842 {
1843 namespace serialization
1844 {
1845 template< typename T>
1846 struct version< Feel::CRBSCM<T> >
1847 {
1848  // at the moment the version of the CRBSCM DB is 0. if any changes is done
1849  // to the format it is mandatory to increase the version number below
1850  // and use the new version number of identify the new entries in the DB
1851  typedef mpl::int_<0> type;
1852  typedef mpl::integral_c_tag tag;
1853  static const unsigned int value = version::type::value;
1854 };
1855 template<typename T> const unsigned int version<Feel::CRBSCM<T> >::value;
1856 }
1857 }
1858 #endif /* __CRBSCM_H */
1859 

Generated on Fri Oct 25 2013 14:24:07 for Feel++ by doxygen 1.8.4