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>
45 #if defined(FEELPP_HAS_GLPK_H)
68 template<
typename TruthModelType>
86 typedef TruthModelType truth_model_type;
87 typedef truth_model_type model_type;
88 typedef boost::shared_ptr<truth_model_type> truth_model_ptrtype;
90 typedef double value_type;
91 typedef boost::tuple<double,double> bounds_type;
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;
100 typedef boost::tuple<double, parameter_type, size_type, double, double> relative_error_type;
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;
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;
122 M_is_initialized( false ),
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() ),
140 po::variables_map
const&
vm )
143 ( boost::format(
"%1%" ) % name ).str(),
144 ( boost::format(
"%1%" ) % name ).str(),
146 M_is_initialized( false ),
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>() ),
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>() )
162 std::cout <<
"Database " << this->
lookForDB() <<
" available and loaded\n";
168 po::variables_map
const&
vm ,
169 truth_model_ptrtype
const & model )
172 ( boost::format(
"%1%" ) % name ).str(),
173 ( boost::format(
"%1%" ) % name ).str(),
175 M_is_initialized( false ),
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>() ),
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>() )
192 std::cout <<
"Database " << this->
lookForDB() <<
" available and loaded\n";
197 po::variables_map
const&
vm ,
198 truth_model_ptrtype
const & model ,
199 bool scm_for_mass_matrix )
202 ( boost::format(
"%1%" ) % name ).str(),
203 ( boost::format(
"%1%" ) % name ).str(),
205 M_is_initialized( false ),
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>() ),
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>() )
222 std::cout <<
"Database " << this->
lookForDB() <<
" available and loaded\n";
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 ),
237 M_C_complement( o.M_C_complement ),
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 )
283 parameterspace_ptrtype
Dmu()
const
309 M_tolerance = tolerance;
328 M_Dmu = M_model->parameterSpace();
334 M_C_complement = sampling_ptrtype(
new sampling_type( M_Dmu ) );
347 M_scm_for_mass_matrix = b;
350 int mMax(
int q )
const
352 if ( M_scm_for_mass_matrix )
353 return M_model->mMaxM( q );
355 return M_model->mMaxA( q );
358 int nb_decomposition_terms_q (
void )
const
360 if ( M_scm_for_mass_matrix )
361 return M_model->Qm();
363 return M_model->Qa();
367 int nb_decomposition_terms_qm (
void )
const
370 for(
int q=0;q<nb_decomposition_terms_q();q++)
372 for(
int m=0;m<mMax(q);m++)
395 boost::tuple<value_type,value_type> ex( parameter_type
const& mu )
const;
420 return lb( *mu, K, indexmu );
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();
463 return bounds_type();
481 std::vector<double>
run( parameter_type
const& mu,
int K );
492 void run(
const double * X,
unsigned long N,
double * Y,
unsigned long P );
506 bool doScmForMassMatrix();
510 sampling_ptrtype c()
const
519 friend class boost::serialization::access;
521 template<
class Archive>
522 void save( Archive & ar,
const unsigned int version )
const;
524 template<
class Archive>
525 void load( Archive & ar,
const unsigned int version ) ;
527 BOOST_SERIALIZATION_SPLIT_MEMBER()
531 bool M_is_initialized;
534 truth_model_ptrtype M_model;
542 parameterspace_ptrtype M_Dmu;
545 sampling_ptrtype M_Xi;
548 sampling_ptrtype M_C, M_C_complement;
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;
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;
563 po::variables_map M_vm;
565 bool M_scm_for_mass_matrix;
568 parameter_type M_mu_ref;
572 po::options_description crbSCMOptions( std::
string const& prefix = "" );
576 template<typename TruthModelType>
577 std::vector<boost::tuple<
double,
double,
double> >
581 return this->offlineSCM();
583 return this->offlineNoSCM();
586 template<
typename TruthModelType>
587 std::vector<boost::tuple<double,double,double> >
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 )
595 inner_prod = M_model->innerProductForMassMatrix();
596 boost::tie( Matrix, boost::tuples::ignore, boost::tuples::ignore ) = M_model->update( M_mu_ref );
600 inner_prod = M_model->innerProduct();
601 boost::tie( boost::tuples::ignore, Matrix, boost::tuples::ignore ) = M_model->update( M_mu_ref );
603 sym = M_model->newMatrix();sym->close();
604 Matrix->symmetricPart( sym );
606 SolverEigen<double>::eigenmodes_type modes;
611 _solver=(
EigenSolverType )M_vm[
"crb.scm.solvereigen-solver-type"].
template as<int>(),
612 _spectrum=SMALLEST_REAL,
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>()
620 double eigen_value = modes.begin()->second.template get<0>();
626 M_C_eigenvalues[0] = modes.begin()->second.template get<0>();
631 std::vector<boost::tuple<double,double,double> > ckconv;
635 template<
typename TruthModelType>
636 std::vector<boost::tuple<double,double,double> >
637 CRBSCM<TruthModelType>::offlineSCM()
639 std::ofstream os_y(
"y.m" );
640 std::ofstream os_C(
"C.m" );
643 std::vector<boost::tuple<double,double,double> > ckconv;
645 M_model->computeAffineDecomposition();
648 M_Xi->randomize( M_vm[
"crb.scm.sampling-size"].
template as<int>() );
650 M_C->setSuperSampling( M_Xi );
651 parameter_type mu( M_Dmu );
654 M_print_matrix = M_vm[
"crb.scm.print-matrix"].template as<bool>() ;
660 for (
int i=0; i<M_Xi->size(); i++ )
662 double testalpha = ex( M_Xi->at( i ) );
667 double relative_error = 1e30;
676 M_C_alpha_lb.clear();
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;
685 if( N_log_equi > 0 || N_equi > 0 )
686 use_predefined_C =
true;
688 if( use_predefined_C )
690 std::string file_name = ( boost::format(
"SamplingC") ).str();
691 std::ifstream file ( file_name );
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" );
697 index_vector = M_C->closestSamplingFromFile(file_name);
698 int sampling_size = index_vector.size();
699 M_iter_max = sampling_size;
702 index = index_vector[0];
704 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
705 std::cout<<
"[CRBSCM::offline] read sampling C ( sampling size : "<<M_iter_max<<
" )"<<std::endl;
710 boost::tie( mu, index ) = M_Xi->min();
711 M_C->push_back( mu, index );
714 M_C_complement = M_C->complement();
719 sparse_matrix_ptrtype B,symmMatrix,Matrix;
720 std::vector<vector_ptrtype> F;
726 std::vector< std::vector<sparse_matrix_ptrtype> > Matrixq;
728 if ( M_scm_for_mass_matrix )
729 boost::tie( Matrixq, boost::tuples::ignore, boost::tuples::ignore ) = M_model->computeAffineDecomposition();
731 boost::tie( boost::tuples::ignore, Matrixq, boost::tuples::ignore ) = M_model->computeAffineDecomposition();
733 int Qmax = nb_decomposition_terms_q();
735 while ( relative_error > M_tolerance && K <= M_iter_max )
737 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
739 std::cout <<
"============================================================\n";
740 std::cout <<
"K=" << K <<
"\n";
743 os_C << M_C->at( K-1 ) <<
"\n";
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) );
751 M_model->solve( mu );
755 if ( M_scm_for_mass_matrix )
757 B = M_model->innerProductForMassMatrix();
758 boost::tie( Matrix, boost::tuples::ignore, F ) = M_model->update( mu );
762 B = M_model->innerProduct();
763 boost::tie( boost::tuples::ignore, Matrix, F ) = M_model->update( mu );
767 for(
int i=0;i<mu.size();i++)
768 mu_str= mu_str + (boost::format(
"_%1%") %mu[i]).str() ;
770 std::string file_name;
772 symmMatrix = M_model->newMatrix();symmMatrix->close();
773 Matrix->symmetricPart( symmMatrix );
775 if( M_print_matrix && ( Environment::worldComm().globalSize() == 1 ) )
777 if( M_scm_for_mass_matrix)
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 );
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 );
791 file_name =
"offline_B" + mu_str +
".m";
792 B->printMatlab( file_name );
798 vector_ptrtype eigenvector;
800 SolverEigen<double>::eigenmodes_type modes;
803 eigs( _matrixA=symmMatrix,
805 _solver=(
EigenSolverType )M_vm[
"crb.scm.solvereigen-solver-type"].
template as<int>(),
806 _spectrum=SMALLEST_REAL,
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>()
817 LOG(INFO) <<
"eigs failed to converge\n";
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() ;
826 eigenvector = modes.begin()->second.template get<2>();
829 M_C_eigenvalues[index] = modes.begin()->second.template get<0>();
830 typedef std::pair<size_type,value_type> key_t;
835 M_eig.push_back( M_C_eigenvalues[index] );
844 for (
size_type q = 0; q < nb_decomposition_terms_q() ; ++q )
846 for (
size_type m_eim = 0; m_eim < mMax(q) ; ++m_eim )
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;
858 os_y << M_Y_ub[K-1][q] <<
"\n";
862 if( M_scm_for_mass_matrix )
863 LOG( INFO ) <<
"scm is done for mass matrix";
865 LOG( INFO )<<
"scm is done for a( . , . ; mu )";
867 double minerr, meanerr;
868 if( use_predefined_C )
870 relative_error = M_tolerance+10;
871 minerr = relative_error;
872 meanerr = relative_error;
876 index = index_vector[ K ];
880 boost::tie( relative_error, mu, index, minerr, meanerr ) =
maxRelativeError( K );
882 std::cout <<
" -- max relative error = " << relative_error
884 <<
" at index = " << index <<
"\n";
887 ckconv.push_back( boost::make_tuple( relative_error, minerr, meanerr ) );
891 if ( relative_error > M_tolerance && K < M_iter_max && ! use_predefined_C )
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 );
900 M_C_complement = M_C->complement();
907 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
908 std::cout <<
"============================================================\n";
912 M_y_bounds_0.resize( Qmax );
913 M_y_bounds_1.resize( Qmax );
914 for (
int q=0; q<Qmax; q++ )
916 for(
int m=0; m<mMax(q); m++)
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>() );
926 template<
typename TruthModelType>
930 y_type err( M_Xi->size() );
933 for (
int k = 0; k < K; ++k )
935 std::cout <<
"[checkC] k= " << k <<
" K=" << K <<
" index in superSampling: " << M_C->indexInSuperSampling( k ) <<
"\n";
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>();
941 if ( _lblb - _lb > 1e-10 )
943 LOG(INFO) <<
"[lberror] the lower bound is decreasing\n"
944 <<
"[lberror] _lblb = " << _lblb <<
"\n"
945 <<
"[lberror] _lb = " << _lb <<
"\n";
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";
952 err( k ) = 1. - _lb/_ub;
954 if ( std::abs( err( k ) ) > 1e-6 )
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() );
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";
977 template<
typename TruthModelType>
978 boost::tuple<typename CRBSCM<TruthModelType>::value_type,
979 typename CRBSCM<TruthModelType>::value_type>
983 sparse_matrix_ptrtype Matrix,symmMatrix;
984 sparse_matrix_ptrtype M ;
985 std::vector<vector_ptrtype> F;
987 if ( M_scm_for_mass_matrix )
989 M = M_model->innerProductForMassMatrix();
990 boost::tie( Matrix, boost::tuples::ignore, F ) = M_model->update( mu );
994 M = M_model->innerProduct();
995 boost::tie( boost::tuples::ignore, Matrix, F ) = M_model->update( mu );
998 symmMatrix = M_model->newMatrix();
1000 Matrix->symmetricPart( symmMatrix );
1002 SolverEigen<double>::eigenmodes_type modesmin=
1003 eigs( _matrixA=Matrix,
1005 _solver=(
EigenSolverType )M_vm[
"crb.scm.solvereigen-solver-type"].
template as<int>(),
1008 _spectrum=SMALLEST_REAL,
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>()
1017 if ( modesmin.empty() )
1019 std::cout <<
"no eigenmode converged: increase --solvereigen-ncv\n";
1023 double eigmin = modesmin.begin()->second.template get<0>();
1027 return boost::make_tuple( eigmin, ti.elapsed() );
1029 SolverEigen<double>::eigenmodes_type modesmax=
1030 eigs( _matrixA=Matrix,
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>()
1040 if ( modesmax.empty() )
1042 std::cout <<
"no eigenmode converged: increase --solvereigen-ncv\n";
1046 double eigmax = modesmax.rbegin()->second.template get<0>();
1058 template<
typename TruthModelType>
1059 boost::tuple<typename CRBSCM<TruthModelType>::value_type,
double>
1063 return this->lbSCM( mu, K , indexmu );
1065 return this->lbNoSCM( mu, K , indexmu );
1068 template<
typename TruthModelType>
1069 boost::tuple<typename CRBSCM<TruthModelType>::value_type,
double>
1073 boost::mpi::timer ti;
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 )
1082 beta_mu = all_beta.template get<0>();
1083 beta_mu_ref = all_beta_ref.template get<0>();
1087 beta_mu = all_beta.template get<1>();
1088 beta_mu_ref = all_beta_ref.template get<1>();
1091 int Q = this->nb_decomposition_terms_q();
1092 vec_min_coeff.resize(Q);
1093 for(
int q=0; q<Q; q++)
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++)
1101 double __beta = beta_mu[q][m]/beta_mu_ref[q][m];
1102 vec_local_min_coeff(m) = __beta ;
1104 double local_min = vec_local_min_coeff.minCoeff();
1105 vec_min_coeff(q) = local_min ;
1107 double min = vec_min_coeff.minCoeff();
1108 double lower_bound = min * M_C_eigenvalues.find(0)->second;
1110 return boost::make_tuple( lower_bound, ti.elapsed() );
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
1118 #if defined(FEELPP_HAS_GLPK_H)
1121 if ( K > this->
KMax() ) K = this->
KMax();
1123 boost::mpi::timer ti;
1126 if ( K <= 0 )
return 0.0;
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] ;
1136 beta_vector_type beta_qm;
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 );
1144 int Malpha = std::min( M_Malpha,std::min( K,M_Xi->size() ) );
1146 int Mplus = std::min( M_Mplus,M_Xi->size()-K );
1149 if ( M_vm[
"crb.scm.strategy"].
template as<int>()==2 )
1150 Mplus = std::min( M_Mplus,std::min( K, M_Xi->size()-K ) );
1155 int nnz = nb_decomposition_terms_qm()*( Malpha+
Mplus );
1157 int ia[1+10000], ja[1+10000];
1163 glp_add_rows( lp,Malpha+Mplus );
1166 std::vector<int> index_vector;
1167 sampling_ptrtype C_neighbors = M_C->searchNearestNeighbors( mu, Malpha , index_vector);
1175 for (
int m=0; m <
Malpha; ++m )
1179 parameter_type mup = C_neighbors->at( m );
1182 if ( M_scm_for_mass_matrix )
1184 boost::tie( beta_qm , boost::tuples::ignore , boost::tuples::ignore ) = M_model->computeBetaQm( mup );
1188 boost::tie( boost::tuples::ignore, beta_qm, boost::tuples::ignore ) = M_model->computeBetaQm( mup );
1194 glp_set_row_name( lp, m+1, ( boost::format(
"c_%1%_%2%" ) % K % m ).str().c_str() );
1199 glp_set_row_bnds( lp,
1202 M_C_eigenvalues.find( C_neighbors->indexInSuperSampling( m ) )->second,
1208 for (
int q = 0; q < nb_decomposition_terms_q(); ++q )
1210 for (
int m_eim = 0; m_eim < mMax( q ); ++m_eim, ++nnz_index, ++count )
1214 ja[nnz_index]=count;
1215 ar[nnz_index]=beta_qm[q][m_eim];
1223 sampling_ptrtype Xi_C_neighbors = M_C_complement->searchNearestNeighbors( mu, Mplus , index_vector);
1230 for (
int m=0; m <
Mplus; ++m )
1235 parameter_type mup = Xi_C_neighbors->at( m );
1243 if ( M_C_alpha_lb[Xi_C_neighbors->indexInSuperSampling( m ) ].find( K-1 ) != M_C_alpha_lb[Xi_C_neighbors->indexInSuperSampling( m ) ].end() )
1250 M_C_alpha_lb[ Xi_C_neighbors->indexInSuperSampling( m )][K-1] =
lb( mup, K-1, Xi_C_neighbors->indexInSuperSampling( m ) ).
template get<0>();
1253 _lb = M_C_alpha_lb[ Xi_C_neighbors->indexInSuperSampling( m )][K-1];
1257 if ( M_scm_for_mass_matrix )
1259 boost::tie( beta_qm, boost::tuples::ignore, boost::tuples::ignore ) = M_model->computeBetaQm( mup );
1264 boost::tie( boost::tuples::ignore, beta_qm, boost::tuples::ignore ) = M_model->computeBetaQm( mup );
1267 glp_set_row_name( lp, Malpha+m+1, ( boost::format(
"xi_c_%1%_%2%" ) % K % m ).str().c_str() );
1269 switch ( M_vm[
"crb.scm.strategy"].
template as<int>() )
1275 glp_set_row_bnds( lp, Malpha+m+1, GLP_LO, 0.0, 0.0 );
1287 glp_set_row_bnds( lp, Malpha+m+1, GLP_LO, _lb, 0.0 );
1293 int nb_already_done = nb_decomposition_terms_qm();
1294 for (
int q = 0; q < nb_decomposition_terms_q(); ++q )
1296 for (
int m_eim = 0; m_eim < this->mMax(q); ++m_eim, ++nnz_index,++count )
1299 ia[nnz_index]=Malpha+m+1;
1300 ja[nnz_index]=count;
1301 ar[nnz_index]=beta_qm[q][m_eim];
1309 if ( M_scm_for_mass_matrix )
1311 boost::tie( beta_qm,boost::tuples::ignore, boost::tuples::ignore ) = M_model->computeBetaQm( mu );
1315 boost::tie( boost::tuples::ignore, beta_qm, boost::tuples::ignore ) = M_model->computeBetaQm( mu );
1321 for(
int q=0; q<nb_decomposition_terms_q(); q++)
1323 for(
int m_eim=0; m_eim<this->mMax(q); m_eim++)
1328 glp_add_cols( lp, nb_decomposition_terms_qm() );
1330 for (
int q = 0; q < nb_decomposition_terms_q(); ++q )
1332 for(
int m_eim=0; m_eim<this->mMax(q); ++m_eim,++count)
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] );
1343 glp_load_matrix( lp, nnz, ia, ja, ar );
1345 glp_init_smcp( &parm );
1346 parm.msg_lev = GLP_MSG_ERR;
1349 glp_simplex( lp, &parm );
1352 double Jobj = glp_get_obj_val( lp );
1356 for (
size_type q = 0; q < M_model->Qa(); ++q )
1358 double y = glp_get_col_prim( lp,q+1 );
1363 glp_delete_prob( lp );
1368 M_C_alpha_lb[ indexmu ][K] = Jobj;
1371 #else // FEELPP_HAS_GLPK_H
1377 return boost::make_tuple( Jobj, ti.elapsed() );
1382 template<
typename TruthModelType>
1383 boost::tuple<typename CRBSCM<TruthModelType>::value_type,
double>
1388 if ( K > this->
KMax() ) K = this->
KMax();
1391 beta_vector_type beta_qm;
1393 if ( M_scm_for_mass_matrix )
1395 boost::tie( beta_qm ,boost::tuples::ignore, boost::tuples::ignore ) = M_model->computeBetaQm( mu );
1400 boost::tie( boost::tuples::ignore, beta_qm, boost::tuples::ignore ) = M_model->computeBetaQm( mu );
1413 for(
int q=0; q<nb_decomposition_terms_q(); q++)
1416 for(
int m=0; m<mMax(q); m++)
1417 y( k ) += beta_qm[q][m]*M_Y_ub[k][q][m];
1422 return boost::make_tuple( y.minCoeff(), ti.elapsed() );
1426 template<
typename TruthModelType>
1427 typename CRBSCM<TruthModelType>::relative_error_type
1432 y_type err( M_Xi->size() );
1434 for (
size_type k = 0; k < M_Xi->size(); ++k )
1438 y_type err( M_C_complement->size() );
1440 for (
size_type k = 0; k < M_C_complement->size(); ++k )
1446 double _lb =
lb( mu, K, k ).template get<0>();
1448 double _lblb =
lb( mu, std::max( K-1,
size_type( 0 ) ), k ).template get<0>();
1450 if ( _lblb - _lb > 1e-10 )
1457 double _ub =
ub( mu, K ).template get<0>();
1459 err( k ) = 1. - _lb/_ub;
1470 Eigen::MatrixXf::Index index;
1471 double maxerr = err.array().abs().maxCoeff( &index );
1473 Eigen::MatrixXf::Index indexmin;
1474 double minerr = err.array().abs().minCoeff( &indexmin );
1476 double meanerr = err.array().abs().sum()/err.size();
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;
1488 return boost::make_tuple( maxerr, M_C_complement->at( index ), M_C_complement->indexInSuperSampling( index ), minerr, meanerr );
1491 template<
typename TruthModelType>
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;
1502 if ( M_scm_for_mass_matrix )
1503 B=M_model->innerProductForMassMatrix();
1505 B=M_model->innerProduct();
1511 for (
int q = 0; q < nb_decomposition_terms_q() ; ++q )
1514 for (
int m = 0; m < mMax(q) ; ++m )
1518 LOG(INFO) <<
"[CRBSCM<TruthModelType>::computeYBounds()] q = " << q <<
"/" << M_model->Qa() <<
"\n";
1520 std::ostringstream os;
1522 if ( M_scm_for_mass_matrix )
1523 Matrix = M_model->Mqm( q , m );
1525 Matrix = M_model->Aqm( q , m );
1528 Matrix->symmetricPart( symmMatrix );
1529 if( Environment::worldComm().globalSize() == 1 )
1531 os <<
"yb_Matrix" << q <<
" - "<< m <<
".m";
1532 Matrix->printMatlab( os.str() );
1534 os <<
"yb_symmMatrix" << q <<
" - "<< m <<
".m";
1535 symmMatrix->printMatlab( os.str() );
1537 if ( symmMatrix->l1Norm()==0.0 )
1539 std::cout <<
"matrix is null\n" ;
1540 M_y_bounds[q].push_back( boost::make_tuple( 0.0, 1e-10 ) );
1548 boost::tie( nconv, eigenvalue_lb, boost::tuples::ignore, boost::tuples::ignore ) =
1549 eigs( _matrixA=symmA,
1551 _spectrum=SMALLEST_MAGNITUDE,
1557 double eigmin=eigenvalue_lb;
1559 boost::tie( nconv, eigenvalue_ub, boost::tuples::ignore, boost::tuples::ignore ) =
1560 eigs( _matrixA=symmA,
1562 _spectrum=LARGEST_MAGNITUDE );
1567 double eigmax=eigenvalue_ub;
1570 SolverEigen<double>::eigenmodes_type modes;
1574 eigs( _matrixA=symmMatrix,
1578 _solver=(
EigenSolverType )M_vm[
"crb.scm.solvereigen-solver-type"].
template as<int>(),
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>()
1589 if ( modes.empty() )
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;
1596 double eigmin = modes.empty()?0:modes.begin()->second.template get<0>();
1599 eigs( _matrixA=symmMatrix,
1603 _solver=(
EigenSolverType )M_vm[
"crb.scm.solvereigen-solver-type"].
template as<int>(),
1604 _spectrum=LARGEST_REAL,
1606 _ncv=M_vm[
"crb.scm.solvereigen-ncv"].
template as<int>(),
1608 _nev=M_vm[
"crb.scm.solvereigen-nev"].
template as<int>(),
1609 _tolerance=M_vm[
"crb.scm.solvereigen-tol"].
template as<double>(),
1611 _maxit=M_vm[
"crb.scm.solvereigen-maxiter"].
template as<int>()
1616 if ( modes.empty() )
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;
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";
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 ) );
1638 LOG(INFO) <<
"[CRBSCM<TruthModelType>::computeYBounds()] stop.\n";
1642 template<
typename TruthModelType>
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;
1651 boost::tie( alpha_ub, alpha_ubti ) = this->
ub( mu, K );
1654 alpha_ub = alpha_lb;
1655 alpha_ubti = alpha_lbti;
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 ) <<
" "
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 );
1677 template<
typename TruthModelType>
1684 for (
unsigned long p= 0; p < N-3; ++p )
1687 for (
unsigned long i=0; i<N; i++ ) std::cout<<
"X["<<i<<
"] = "<<X[i]<<std::endl;
1689 double meshSize = X[N-3];
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";
1706 template<
typename TruthModelType>
1707 template<
class Archive>
1718 ar & M_C_complement;
1719 ar & M_C_eigenvalues;
1729 ar & M_C_eigenvalues;
1733 template<
typename TruthModelType>
1734 template<
class Archive>
1736 CRBSCM<TruthModelType>::load( Archive & ar,
const unsigned int version )
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)
1745 throw std::logic_error(
"[CRBSCM::load] ERROR the database created is not appropriate to use SCM. Use option crb.scm.rebuild-database=true");
1747 throw std::logic_error(
"[CRBSCM::load] ERROR the database was created to use SCM. Use option crb.scm.rebuild-database=true");
1756 ar & M_C_complement;
1757 ar & M_C_eigenvalues;
1762 int Qmax = this->nb_decomposition_terms_q();
1763 M_y_bounds.resize( Qmax );
1764 for (
int q=0; q<Qmax; q++ )
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] ) );
1773 ar & M_C_eigenvalues;
1778 template<
typename TruthModelType>
1780 CRBSCM<TruthModelType>::doScmForMassMatrix()
1782 bool b = this->
vm()[
"crb.scm.do-scm-for-mass-matrix"].template as<bool>();
1786 template<
typename TruthModelType>
1788 CRBSCM<TruthModelType>::rebuildDB()
1790 bool rebuild = this->
vm()[
"crb.scm.rebuild-database"].template as<bool>();
1794 template<
typename TruthModelType>
1802 boost::archive::text_oarchive oa( ofs );
1808 template<
typename TruthModelType>
1812 if ( this->rebuildDB() )
1820 if ( !fs::exists( db ) )
1823 std::cout <<
"Loading " << db <<
"...\n";
1824 fs::ifstream ifs( db );
1827 boost::archive::text_iarchive ia( ifs );
1830 std::cout <<
"Loading " << db <<
" done...\n";
1831 this->setIsLoaded(
true );
1843 namespace serialization
1845 template<
typename T>
1846 struct version< Feel::
CRBSCM<T> >
1851 typedef mpl::int_<0> type;
1852 typedef mpl::integral_c_tag tag;
1853 static const unsigned int value = version::type::value;
1855 template<
typename T>
const unsigned int version<Feel::CRBSCM<T> >::value;