29 #ifndef __ParameterSpace_H
30 #define __ParameterSpace_H 1
34 #include <boost/lambda/lambda.hpp>
35 #include <boost/foreach.hpp>
39 #if defined(FEELPP_HAS_ANN_H)
45 #include <feel/feelcore/feel.hpp>
58 class ParameterSpace:
public boost::enable_shared_from_this<ParameterSpace<P> >
77 typedef boost::shared_ptr<parameterspace_type> parameterspace_ptrtype;
84 class Element :
public Eigen::Matrix<double,Dimension,1>
86 typedef Eigen::Matrix<double,Dimension,1> super;
89 typedef boost::shared_ptr<parameterspace_type> parameterspace_ptrtype;
136 template<
typename OtherDerived>
137 super&
operator=(
const Eigen::MatrixBase<OtherDerived>& other )
142 void setParameterSpace( parameterspace_ptrtype
const& space )
158 if ( !M_space->check() )
160 LOG(INFO) <<
"No need to check element since parameter space is no valid (yet)\n";
166 mpi::all_reduce( Environment::worldComm().localComm(), *
this, sum,
168 int proc_number = Environment::worldComm().globalRank();
169 if( (this->array()-sum.array()/Environment::numberOfProcessors()).abs().maxCoeff() > 1e-7 )
171 std::cout <<
"Parameter not identical on all processors: "<<
"current parameter on proc "<<proc_number<<
" : [";
172 for(
int i=0; i<this->size(); i++) std::cout <<std::setprecision(15)<< this->operator()(i) <<
", ";
173 std::cout<<std::setprecision(15)<<this->operator()(this->size()-1)<<
" ]";
174 std::cout <<std::setprecision(15)<<
" and test" << (this->array()-sum.array()/Environment::numberOfProcessors()).abs().maxCoeff() <<
"\n";
176 Environment::worldComm().barrier();
177 CHECK( (this->array()-sum.array()/Environment::numberOfProcessors()).abs().maxCoeff() < 1e-7 )
178 <<
"Parameter not identical on all processors(" << Environment::worldComm().masterRank() <<
"/" << Environment::numberOfProcessors() <<
")\n"
179 <<
"max error: " << (this->array()-sum.array()/Environment::numberOfProcessors()).abs().maxCoeff() <<
"\n"
180 <<
"error: " << std::setprecision(15) << (this->array()-sum.array()/Environment::numberOfProcessors()).abs() <<
"\n"
181 <<
"current parameter " << std::setprecision(15) << *
this <<
"\n"
182 <<
"sum parameter " << std::setprecision(15) << sum <<
"\n"
183 <<
"space min : " << M_space->min() <<
"\n"
184 <<
"space max : " << M_space->max() <<
"\n";
187 friend class boost::serialization::access;
188 template<
class Archive>
189 void save( Archive & ar,
const unsigned int version )
const
191 ar & boost::serialization::base_object<super>( *this );
195 template<
class Archive>
196 void load( Archive & ar,
const unsigned int version )
198 ar & boost::serialization::base_object<super>( *this );
201 BOOST_SERIALIZATION_SPLIT_MEMBER()
204 parameterspace_ptrtype M_space;
208 typedef boost::shared_ptr<
Element> element_ptrtype;
210 element_ptrtype elementPtr() { element_ptrtype e(
new element_type( this->shared_from_this() ) ); *e=element();
return e; }
213 return (
min()-
max()).array().abs().maxCoeff() > 1e-10;
221 typedef std::vector<Element> super;
225 typedef boost::shared_ptr<sampling_type> sampling_ptrtype;
228 typedef boost::shared_ptr<parameterspace_type> parameterspace_ptrtype;
231 typedef boost::shared_ptr<element_type> element_ptrtype;
233 #if defined(FEELPP_HAS_ANN_H)
234 typedef ANNkd_tree kdtree_type;
235 typedef boost::shared_ptr<kdtree_type> kdtree_ptrtype;
239 Sampling( parameterspace_ptrtype space,
int N = 1, sampling_ptrtype supersampling = sampling_ptrtype() )
243 M_supersampling( supersampling )
244 #if defined( FEELPP_HAS_ANN_H )
254 void setElements( std::vector< element_type > V )
256 CHECK( M_space ) <<
"Invalid(null pointer) parameter space for parameter generation\n";
261 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
263 BOOST_FOREACH(
auto mu, V )
264 super::push_back( mu );
267 boost::mpi::broadcast( Environment::worldComm() , *
this , Environment::worldComm().masterRank() );
275 void randomize(
int N )
277 CHECK( M_space ) <<
"Invalid(null pointer) parameter space for parameter generation\n";
285 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
287 for (
int i = 0; i < N; ++i )
294 boost::mpi::broadcast( Environment::worldComm() , *
this , Environment::worldComm().masterRank() );
302 void logEquidistribute(
int N )
307 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
310 for (
int i = 0; i < N; ++i )
312 double factor = double( i )/( N-1 );
317 boost::mpi::broadcast( Environment::worldComm() , *
this , Environment::worldComm().masterRank() );
324 void logEquidistributeProduct( std::vector<int> N )
329 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
331 int number_of_directions = N.size();
334 std::vector< std::vector< double > > components;
335 components.resize( number_of_directions );
337 for(
int direction=0; direction<number_of_directions; direction++)
340 components[direction]= coeff_vector ;
343 generateElementsProduct( components );
346 boost::mpi::broadcast( Environment::worldComm() , *
this , Environment::worldComm().masterRank() );
355 void mixEquiLogEquidistributeProduct( std::vector<int> Nlogequi , std::vector<int> Nequi )
360 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
363 int number_of_directions_equi = Nequi.size();
364 int number_of_directions_logequi = Nlogequi.size();
365 FEELPP_ASSERT( number_of_directions_equi == number_of_directions_logequi )( number_of_directions_equi )( number_of_directions_logequi )
366 .error(
"incompatible number of directions, you don't manipulate the same parameter space for log-equidistributed and equidistributed sampling" );
369 std::vector< std::vector< double > > components_equi;
370 components_equi.resize( number_of_directions_equi );
371 std::vector< std::vector< double > > components_logequi;
372 components_logequi.resize( number_of_directions_logequi );
374 for(
int direction=0; direction<number_of_directions_equi; direction++)
377 components_equi[direction]= coeff_vector_equi ;
379 components_logequi[direction]= coeff_vector_logequi ;
381 std::vector< std::vector< std::vector<double> > > vector_components(2);
382 vector_components[0]=components_logequi;
383 vector_components[1]=components_equi;
384 generateElementsProduct( vector_components );
387 boost::mpi::broadcast( Environment::worldComm() , *
this , Environment::worldComm().masterRank() );
396 void generateElementsProduct( std::vector< std::vector< double > > components )
398 std::vector< std::vector< std::vector< double > > > vector_components(1);
399 vector_components[0]=components;
400 generateElementsProduct( vector_components );
403 void generateElementsProduct( std::vector< std::vector< std::vector< double > > > vector_components )
405 int number_of_comp = vector_components.size();
406 for(
int comp=0; comp<number_of_comp; comp++)
408 int number_of_directions=vector_components[comp].size();
409 element_type mu( M_space );
412 std::vector<int> idx( number_of_directions );
413 for(
int direction=0; direction<number_of_directions; direction++)
416 int last_direction=number_of_directions-1;
419 auto min=M_space->min();
420 auto max=M_space->max();
424 bool already_exist =
true;
426 while( already_exist && !stop )
429 for(
int direction=0; direction<number_of_directions; direction++)
431 int index = idx[direction];
432 double value = vector_components[comp][direction][index];
433 mu(direction) = value ;
436 if( mu ==
min && comp>0 )
438 if( mu ==
max && comp>0 )
442 for(
int direction=0; direction<number_of_directions; direction++)
444 if( idx[direction] < vector_components[comp][direction].size()-1 )
452 if( direction == last_direction )
461 super::push_back( mu );
476 void writeOnFile( std::string file_name =
"list_of_parameters_taken" )
478 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
481 file.open( file_name,std::ios::out );
482 element_type mu( M_space );
483 int size = mu.size();
485 BOOST_FOREACH( mu, *
this )
487 file<<
" mu_"<<number<<
"= [ ";
488 for(
int i=0; i<size-1; i++)
489 file << mu[i]<<
" , ";
490 file<< mu[size-1] <<
" ] \n" ;
505 double readFromFile( std::string file_name=
"list_of_parameters_taken" )
507 std::ifstream file ( file_name );
512 while( ! file.eof() )
514 element_type mu( M_space );
524 super::push_back( mu );
539 std::vector<int> closestSamplingFromFile( std::string file_name=
"list_of_parameters_taken")
541 std::vector<int> index_vector;
542 std::ifstream file( file_name );
547 while( ! file.eof() )
549 element_type mu( M_space );
560 std::vector<int> local_index_vector;
561 auto neighbors = M_supersampling->searchNearestNeighbors( mu, 1 , local_index_vector );
562 auto closest_mu = neighbors->at(0);
563 int index = local_index_vector[0];
564 this->push_back( closest_mu , index );
567 index_vector.push_back( index );
577 void equidistribute(
int N )
582 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
585 for (
int i = 0; i < N; ++i )
587 double factor = double( i )/( N-1 );
592 boost::mpi::broadcast( Environment::worldComm() , *
this , Environment::worldComm().masterRank() );
599 void equidistributeProduct( std::vector<int> N )
604 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
606 int number_of_directions = N.size();
609 std::vector< std::vector< double > > components;
610 components.resize( number_of_directions );
612 for(
int direction=0; direction<number_of_directions; direction++)
615 components[direction] = coeff_vector ;
618 generateElementsProduct( components );
621 boost::mpi::broadcast( Environment::worldComm() , *
this , Environment::worldComm().masterRank() );
628 boost::tuple<element_type,size_type>
631 element_type mumin( M_space );
632 mumin = M_space->max();
634 element_type mu( M_space );
637 BOOST_FOREACH( mu, *
this )
639 if ( mu.norm() < mumin.norm() )
648 return boost::make_tuple( mumin, index );
654 boost::tuple<element_type,size_type>
657 element_type mumax( M_space );
658 mumax = M_space->min();
660 element_type mu( M_space );
663 BOOST_FOREACH( mu, *
this )
665 if ( mu.norm() > mumax.norm() )
674 return boost::make_tuple( mumax, index );
687 sampling_ptrtype
const& superSampling()
const
689 return M_supersampling;
695 void setSuperSampling( sampling_ptrtype
const& super )
697 M_supersampling = super;
706 sampling_ptrtype searchNearestNeighbors( element_type
const& mu,
size_type M , std::vector<int>& index_vector );
711 sampling_ptrtype complement()
const;
718 if ( M_supersampling ) M_superindices.push_back( index );
720 super::push_back( mu );
730 return M_superindices[ index ];
735 friend class boost::serialization::access;
736 template<
class Archive>
737 void save( Archive & ar,
const unsigned int version )
const
739 ar & boost::serialization::base_object<super>( *this );
741 ar & M_supersampling;
745 template<
class Archive>
746 void load( Archive & ar,
const unsigned int version )
748 ar & boost::serialization::base_object<super>( *this );
750 ar & M_supersampling;
753 BOOST_SERIALIZATION_SPLIT_MEMBER()
755 parameterspace_ptrtype M_space;
757 sampling_ptrtype M_supersampling;
759 #if defined( FEELPP_HAS_ANN_H )
760 kdtree_ptrtype M_kdtree;
764 typedef Sampling sampling_type;
765 typedef boost::shared_ptr<sampling_type> sampling_ptrtype;
767 sampling_ptrtype sampling() {
return sampling_ptrtype(
new sampling_type( this->shared_from_this() ) ); }
789 ParameterSpace( element_type
const& emin, element_type
const& emax )
794 M_min.setParameterSpace( this->shared_from_this() );
795 M_max.setParameterSpace( this->shared_from_this() );
804 static parameterspace_ptrtype
New()
840 element_type
const&
min()
const
848 element_type
const&
max()
const
858 return ( ( M_min.array().log() + M_max.array().log() )/2. ).log();
866 return ( ( M_min + M_max )/2. );
900 static element_type
logRandom( parameterspace_ptrtype space,
bool broadcast )
904 google::FlushLogFiles(google::GLOG_INFO);
907 element_type mu( space );
908 if ( !space->check() )
return mu;
909 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
913 mu = logRandom1( space );
917 boost::mpi::broadcast( Environment::worldComm() , mu , Environment::worldComm().masterRank() );
918 Environment::worldComm().barrier();
925 return logRandom1( space );
928 static element_type logRandom1( parameterspace_ptrtype space )
930 element_type mur( space );
931 mur.array() = element_type::Random().array().abs();
933 google::FlushLogFiles(google::GLOG_INFO);
934 element_type mu( space );
935 mu.array() = ( space->min().array().log()+mur.array()*( space->max().array().log()-space->min().array().log() ) ).array().exp();
943 static element_type
random( parameterspace_ptrtype space,
bool broadcast =
true )
945 std::srand( static_cast<unsigned>( std::time( 0 ) ) );
946 element_type mur( space );
949 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
951 mur.array() = element_type::Random().array().abs();
953 boost::mpi::broadcast( Environment::worldComm() , mur , Environment::worldComm().masterRank() );
957 mur.array() = element_type::Random().array().abs();
961 element_type mu( space );
962 mu.array() = space->min()+mur.array()*( space->max()-space->min() );
976 std::vector<double> result(N);
978 auto min_element = space->min();
979 auto max_element = space->max();
980 int space_dimension=space->dimension();
981 FEELPP_ASSERT( space_dimension > direction )( space_dimension )( direction ).error(
"bad dimension of vector containing number of samples on each direction" );
982 FEELPP_ASSERT( direction >= 0 )( direction ).error(
"the direction must be positive number" );
984 double min = min_element(direction);
985 double max = max_element(direction);
987 for(
int i=0; i<N; i++)
989 double factor = (double)(i)/(N-1);
990 result[i] = math::exp(math::log(min)+factor*( math::log(max)-math::log(min) ) );
1003 element_type mu( space );
1004 mu.array() = ( space->min().array().log()+factor*( space->max().array().log()-space->min().array().log() ) ).exp();
1017 std::vector<double> result(N);
1019 auto min_element = space->min();
1020 auto max_element = space->max();
1022 int mu_size = min_element.size();
1023 if( (mu_size < direction) || (direction < 0) )
1024 throw std::logic_error(
"[ParameterSpace::equidistributeInDirection] ERROR : bad dimension of vector containing number of samples on each direction" );
1026 double min = min_element(direction);
1027 double max = max_element(direction);
1029 for(
int i=0; i<N; i++)
1031 double factor = (double)(i)/(N-1);
1032 result[i] = min+factor*( max-
min );
1045 element_type mu( space );
1046 mu = space->min()+factor*( space->max()-space->min() );
1058 friend class boost::serialization::access;
1059 template<
class Archive>
1060 void save( Archive & ar,
const unsigned int version )
const
1066 template<
class Archive>
1067 void load( Archive & ar,
const unsigned int version )
1072 BOOST_SERIALIZATION_SPLIT_MEMBER()
1091 sampling_ptrtype complement;
1093 if ( M_supersampling )
1099 complement = sampling_ptrtype(
new sampling_type( M_space, 1, M_supersampling ) );
1100 complement->clear();
1102 if ( !M_superindices.empty() )
1104 for (
size_type i = 0, j = 0; i < M_supersampling->size(); ++i )
1106 if ( std::find( M_superindices.begin(), M_superindices.end(), i ) == M_superindices.end() )
1109 if ( j < M_supersampling->size()-this->size() )
1112 complement->push_back( M_supersampling->at( i ),i );
1125 boost::shared_ptr<typename ParameterSpace<P>::Sampling>
1128 std::vector<int>& index_vector)
1131 sampling_ptrtype neighbors(
new sampling_type( M_space, M ) );
1132 #if defined(FEELPP_HAS_ANN_H)
1134 if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1140 ANNpointArray data_pts;
1141 data_pts = annAllocPts( this->size(), M_space->Dimension );
1143 for (
size_type i = 0; i < this->size(); ++i )
1145 std::copy( this->at( i ).data(), this->at( i ).data()+M_space->Dimension, data_pts[i] );
1146 FEELPP_ASSERT( data_pts[i] != 0 )( i ) .error(
"invalid pointer" );
1150 M_kdtree = kdtree_ptrtype(
new kdtree_type( data_pts, this->size(), M_space->Dimension ) );
1155 if ( this->size() < M )
1158 std::vector<int> nnIdx( M );
1159 std::vector<double> dists( M );
1161 M_kdtree->annkSearch(
1162 const_cast<double*>( mu.data() ),
1170 if ( M_supersampling )
1172 neighbors->setSuperSampling( M_supersampling );
1174 if ( !M_superindices.empty() )
1175 neighbors->M_superindices.resize( M );
1182 neighbors->at( i ) = this->at( nnIdx[i] );
1183 index_vector.push_back( nnIdx[i] );
1186 if ( M_supersampling && !M_superindices.empty() )
1187 neighbors->M_superindices[i] = M_superindices[ nnIdx[i] ];
1190 annDeallocPts( data_pts );
1193 boost::mpi::broadcast( Environment::worldComm() , neighbors , Environment::worldComm().masterRank() );
1194 boost::mpi::broadcast( Environment::worldComm() , index_vector , Environment::worldComm().masterRank() );
size_type indexInSuperSampling(size_type index) const
given a local index, returns the index in the super sampling
Definition: parameterspace.hpp:728
element_type logMiddle() const
the log-middle point of the parameter space
Definition: parameterspace.hpp:856
Element(Element const &e)
Definition: parameterspace.hpp:103
static std::vector< double > equidistributeInDirection(parameterspace_ptrtype space, int direction, int N)
Returns a vector representing the values of equidistributed element of the parameter space in the giv...
Definition: parameterspace.hpp:1015
static parameterspace_ptrtype New()
Definition: parameterspace.hpp:804
static const uint16_type Dimension
dimension of the parameter space
Definition: parameterspace.hpp:68
element_type const & max() const
Definition: parameterspace.hpp:848
Element(parameterspace_ptrtype space)
Definition: parameterspace.hpp:111
void setMin(element_type const &min)
Definition: parameterspace.hpp:878
ParameterSpace()
default constructor
Definition: parameterspace.hpp:774
void push_back(element_type const &mu, size_type index)
add new parameter mu in sampling and store index in super sampling
Definition: parameterspace.hpp:716
static element_type logRandom(parameterspace_ptrtype space, bool broadcast)
Returns a log random element of the parameter space.
Definition: parameterspace.hpp:900
Parameter space sampling class.
Definition: parameterspace.hpp:219
element_type const & min() const
Definition: parameterspace.hpp:840
int dimension() const
Definition: parameterspace.hpp:832
static element_type logEquidistributed(double factor, parameterspace_ptrtype space)
Returns a log equidistributed element of the parameter space.
Definition: parameterspace.hpp:1001
static element_type equidistributed(double factor, parameterspace_ptrtype space)
Returns a equidistributed element of the parameter space.
Definition: parameterspace.hpp:1043
static element_type random(parameterspace_ptrtype space, bool broadcast=true)
Returns a log random element of the parameter space.
Definition: parameterspace.hpp:943
element_type middle() const
the middle point of the parameter space
Definition: parameterspace.hpp:864
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319
Parameter space class.
Definition: parameterspace.hpp:58
Element()
Definition: parameterspace.hpp:94
Elements & operator=(Elements const &e)
Definition: elements.hpp:335
ParameterSpace & operator=(ParameterSpace const &o)
copy operator
Definition: parameterspace.hpp:815
~Element()
Definition: parameterspace.hpp:119
~ParameterSpace()
destructor
Definition: parameterspace.hpp:798
ParameterSpace(ParameterSpace const &o)
copy constructor
Definition: parameterspace.hpp:783
Element & operator=(Element const &e)
Definition: parameterspace.hpp:125
void setMax(element_type const &max)
Definition: parameterspace.hpp:886
parameterspace_ptrtype parameterSpace() const
Retuns the parameter space.
Definition: parameterspace.hpp:151
element of a parameter space
Definition: parameterspace.hpp:84
static std::vector< double > logEquidistributeInDirection(parameterspace_ptrtype space, int direction, int N)
Returns a vector representing the values of log equidistributed element of the parameter space in the...
Definition: parameterspace.hpp:974