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
parameterspace.hpp
Go to the documentation of this file.
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-11
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 __ParameterSpace_H
30 #define __ParameterSpace_H 1
31 
32 #include <vector>
33 #include <algorithm>
34 #include <boost/lambda/lambda.hpp>
35 #include <boost/foreach.hpp>
36 
37 #include <Eigen/Core>
38 
39 #if defined(FEELPP_HAS_ANN_H)
40 #include <ANN/ANN.h>
41 #endif /* FEELPP_HAS_ANN_H */
42 
43 
44 
45 #include <feel/feelcore/feel.hpp>
46 #include <feel/feelmesh/kdtree.hpp>
47 
48 namespace Feel
49 {
57 template<int P>
58 class ParameterSpace: public boost::enable_shared_from_this<ParameterSpace<P> >
59 {
60 public:
61 
62 
66 
68  static const uint16_type Dimension = P;
69 
71 
75 
77  typedef boost::shared_ptr<parameterspace_type> parameterspace_ptrtype;
78 
80 
84  class Element : public Eigen::Matrix<double,Dimension,1>
85  {
86  typedef Eigen::Matrix<double,Dimension,1> super;
87  public:
89  typedef boost::shared_ptr<parameterspace_type> parameterspace_ptrtype;
90 
95  :
96  super(),
97  M_space()
98  {}
99 
103  Element( Element const& e )
104  :
105  super( e ),
106  M_space( e.M_space )
107  {}
111  Element( parameterspace_ptrtype space )
112  :
113  super(),
114  M_space( space )
115  {}
120  {}
121 
125  Element& operator=( Element const& e )
126  {
127  if ( this != &e )
128  {
129  super::operator=( e );
130  M_space = e.M_space;
131  }
132 
133  return *this;
134  }
135 
136  template<typename OtherDerived>
137  super& operator=( const Eigen::MatrixBase<OtherDerived>& other )
138  {
139  return super::operator=( other );
140  }
141 
142  void setParameterSpace( parameterspace_ptrtype const& space )
143  {
144  M_space = space;
145  }
146 
147 
151  parameterspace_ptrtype parameterSpace() const
152  {
153  return M_space;
154  }
155 
156  void check() const
157  {
158  if ( !M_space->check() )
159  {
160  LOG(INFO) << "No need to check element since parameter space is no valid (yet)\n";
161  return;
162  }
163  Element sum;
164  sum.setZero();
165  // verify that the element is the same on all processors
166  mpi::all_reduce( Environment::worldComm().localComm(), *this, sum,
167  []( Element const& m1, Element const& m2 ) { return m1+m2; } );
168  int proc_number = Environment::worldComm().globalRank();
169  if( (this->array()-sum.array()/Environment::numberOfProcessors()).abs().maxCoeff() > 1e-7 )
170  {
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";
175  }
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";
185  }
186  private:
187  friend class boost::serialization::access;
188  template<class Archive>
189  void save( Archive & ar, const unsigned int version ) const
190  {
191  ar & boost::serialization::base_object<super>( *this );
192  ar & M_space;
193  }
194 
195  template<class Archive>
196  void load( Archive & ar, const unsigned int version )
197  {
198  ar & boost::serialization::base_object<super>( *this );
199  ar & M_space;
200  }
201  BOOST_SERIALIZATION_SPLIT_MEMBER()
202 
203  private:
204  parameterspace_ptrtype M_space;
205  };
206 
207  typedef Element element_type;
208  typedef boost::shared_ptr<Element> element_ptrtype;
209  element_type element() { return parameterspace_type::logRandom( this->shared_from_this(), true ); }
210  element_ptrtype elementPtr() { element_ptrtype e( new element_type( this->shared_from_this() ) ); *e=element(); return e; }
211  bool check() const
212  {
213  return (min()-max()).array().abs().maxCoeff() > 1e-10;
214  }
219  class Sampling : public std::vector<Element>
220  {
221  typedef std::vector<Element> super;
222  public:
223 
224  typedef Sampling sampling_type;
225  typedef boost::shared_ptr<sampling_type> sampling_ptrtype;
226 
228  typedef boost::shared_ptr<parameterspace_type> parameterspace_ptrtype;
229 
230  typedef typename parameterspace_type::Element element_type;
231  typedef boost::shared_ptr<element_type> element_ptrtype;
232 
233 #if defined(FEELPP_HAS_ANN_H)
234  typedef ANNkd_tree kdtree_type;
235  typedef boost::shared_ptr<kdtree_type> kdtree_ptrtype;
236 #endif /* FEELPP_HAS_ANN_H */
237 
238 
239  Sampling( parameterspace_ptrtype space, int N = 1, sampling_ptrtype supersampling = sampling_ptrtype() )
240  :
241  super( N ),
242  M_space( space ),
243  M_supersampling( supersampling )
244 #if defined( FEELPP_HAS_ANN_H )
245  ,M_kdtree()
246 #endif
247  {}
248 
249 
254  void setElements( std::vector< element_type > V )
255  {
256  CHECK( M_space ) << "Invalid(null pointer) parameter space for parameter generation\n";
257 
258  // first empty the set
259  this->clear();
260 
261  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
262  {
263  BOOST_FOREACH( auto mu, V )
264  super::push_back( mu );
265  }
266 
267  boost::mpi::broadcast( Environment::worldComm() , *this , Environment::worldComm().masterRank() );
268 
269  }
270 
275  void randomize( int N )
276  {
277  CHECK( M_space ) << "Invalid(null pointer) parameter space for parameter generation\n";
278 
279  // first empty the set
280  this->clear();
281  //std::srand(static_cast<unsigned>(std::time(0)));
282 
283  // fill with log Random elements from the parameter space
284  //only with one proc and then broadcast
285  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
286  {
287  for ( int i = 0; i < N; ++i )
288  {
289  super::push_back( parameterspace_type::logRandom( M_space, false ) );
290  }
291 
292  }
293 
294  boost::mpi::broadcast( Environment::worldComm() , *this , Environment::worldComm().masterRank() );
295 
296  }
297 
302  void logEquidistribute( int N )
303  {
304  // first empty the set
305  this->clear();
306 
307  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
308  {
309  // fill with log Random elements from the parameter space
310  for ( int i = 0; i < N; ++i )
311  {
312  double factor = double( i )/( N-1 );
313  super::push_back( parameterspace_type::logEquidistributed( factor,
314  M_space ) );
315  }
316  }
317  boost::mpi::broadcast( Environment::worldComm() , *this , Environment::worldComm().masterRank() );
318  }
319 
324  void logEquidistributeProduct( std::vector<int> N )
325  {
326  // first empty the set
327  this->clear();
328 
329  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
330  {
331  int number_of_directions = N.size();
332 
333  //contains values of parameters on each direction
334  std::vector< std::vector< double > > components;
335  components.resize( number_of_directions );
336 
337  for(int direction=0; direction<number_of_directions; direction++)
338  {
339  std::vector<double> coeff_vector = parameterspace_type::logEquidistributeInDirection( M_space , direction , N[direction] );
340  components[direction]= coeff_vector ;
341  }
342 
343  generateElementsProduct( components );
344 
345  }//end of master proc
346  boost::mpi::broadcast( Environment::worldComm() , *this , Environment::worldComm().masterRank() );
347  }
348 
349 
355  void mixEquiLogEquidistributeProduct( std::vector<int> Nlogequi , std::vector<int> Nequi )
356  {
357  // first empty the set
358  this->clear();
359 
360  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
361  {
362 
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" );
367 
368  //contains values of parameters on each direction
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 );
373 
374  for(int direction=0; direction<number_of_directions_equi; direction++)
375  {
376  std::vector<double> coeff_vector_equi = parameterspace_type::equidistributeInDirection( M_space , direction , Nequi[direction] );
377  components_equi[direction]= coeff_vector_equi ;
378  std::vector<double> coeff_vector_logequi = parameterspace_type::logEquidistributeInDirection( M_space , direction , Nlogequi[direction] );
379  components_logequi[direction]= coeff_vector_logequi ;
380  }
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 );
385 
386  }//end of master proc
387  boost::mpi::broadcast( Environment::worldComm() , *this , Environment::worldComm().masterRank() );
388  }
389 
390 
391  /*
392  * \param a vector of vector containing values of parameters on each direction
393  * example if components has 2 vectors : [1,2,3] and [100,200] then it creates the samping :
394  * (1,100) - (2,100) - (3,100) - (1,200) - (2,200) - (3,200)
395  */
396  void generateElementsProduct( std::vector< std::vector< double > > components )
397  {
398  std::vector< std::vector< std::vector< double > > > vector_components(1);
399  vector_components[0]=components;
400  generateElementsProduct( vector_components );
401  }
402 
403  void generateElementsProduct( std::vector< std::vector< std::vector< double > > > vector_components )
404  {
405  int number_of_comp = vector_components.size();
406  for(int comp=0; comp<number_of_comp; comp++)
407  {
408  int number_of_directions=vector_components[comp].size();
409  element_type mu( M_space );
410 
411  //initialization
412  std::vector<int> idx( number_of_directions );
413  for(int direction=0; direction<number_of_directions; direction++)
414  idx[direction]=0;
415 
416  int last_direction=number_of_directions-1;
417 
418  bool stop=false;
419  auto min=M_space->min();
420  auto max=M_space->max();
421 
422  do
423  {
424  bool already_exist = true;
425  //construction of the parameter
426  while( already_exist && !stop )
427  {
428  already_exist=false;
429  for(int direction=0; direction<number_of_directions; direction++)
430  {
431  int index = idx[direction];
432  double value = vector_components[comp][direction][index];
433  mu(direction) = value ;
434  }
435 
436  if( mu == min && comp>0 )
437  already_exist=true;
438  if( mu == max && comp>0 )
439  already_exist=true;
440 
441  //update values or stop
442  for(int direction=0; direction<number_of_directions; direction++)
443  {
444  if( idx[direction] < vector_components[comp][direction].size()-1 )
445  {
446  idx[direction]++;
447  break;
448  }
449  else
450  {
451  idx[direction]=0;
452  if( direction == last_direction )
453  stop = true;
454  }
455  }//end loop over directions
456 
457  }//while
458 
459 
460  //add the new parameter
461  super::push_back( mu );
462 
463  } while( ! stop );
464 
465  }//end loop on vector_component
466  }
467 
468 
476  void writeOnFile( std::string file_name = "list_of_parameters_taken" )
477  {
478  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
479  {
480  std::ofstream file;
481  file.open( file_name,std::ios::out );
482  element_type mu( M_space );
483  int size = mu.size();
484  int number = 0;
485  BOOST_FOREACH( mu, *this )
486  {
487  file<<" mu_"<<number<<"= [ ";
488  for(int i=0; i<size-1; i++)
489  file << mu[i]<<" , ";
490  file<< mu[size-1] << " ] \n" ;
491  number++;
492  }
493  file.close();
494  }
495  }
496 
505  double readFromFile( std::string file_name= "list_of_parameters_taken" )
506  {
507  std::ifstream file ( file_name );
508  double mui;
509  std::string str;
510  int number=0;
511  file>>str;
512  while( ! file.eof() )
513  {
514  element_type mu( M_space );
515  file>>str;
516  int i=0;
517  while ( str!="]" )
518  {
519  file >> mui;
520  mu[i] = mui;
521  file >> str;
522  i++;
523  }
524  super::push_back( mu );
525  number++;
526  file>>str;
527  }
528  file.close();
529  return number;
530  }
531 
532 
539  std::vector<int> closestSamplingFromFile( std::string file_name="list_of_parameters_taken")
540  {
541  std::vector<int> index_vector;
542  std::ifstream file( file_name );
543  double mui;
544  std::string str;
545  int number=0;
546  file>>str;
547  while( ! file.eof() )
548  {
549  element_type mu( M_space );
550  file>>str;
551  int i=0;
552  while ( str!="]" )
553  {
554  file >> mui;
555  mu[i] = mui;
556  file >> str;
557  i++;
558  }
559  //search the closest neighbor of mu in the supersampling
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 );
565  number++;
566  file>>str;
567  index_vector.push_back( index );
568  }
569  file.close();
570  return index_vector;
571  }
572 
577  void equidistribute( int N )
578  {
579  // first empty the set
580  this->clear();
581 
582  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
583  {
584  // fill with log Random elements from the parameter space
585  for ( int i = 0; i < N; ++i )
586  {
587  double factor = double( i )/( N-1 );
588  super::push_back( parameterspace_type::equidistributed( factor,
589  M_space ) );
590  }
591  }
592  boost::mpi::broadcast( Environment::worldComm() , *this , Environment::worldComm().masterRank() );
593  }
594 
599  void equidistributeProduct( std::vector<int> N )
600  {
601  // first empty the set
602  this->clear();
603 
604  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
605  {
606  int number_of_directions = N.size();
607 
608  //contains values of parameters on each direction
609  std::vector< std::vector< double > > components;
610  components.resize( number_of_directions );
611 
612  for(int direction=0; direction<number_of_directions; direction++)
613  {
614  std::vector<double> coeff_vector = parameterspace_type::equidistributeInDirection( M_space , direction , N[direction] );
615  components[direction] = coeff_vector ;
616  }
617 
618  generateElementsProduct( components );
619 
620  }//end of master proc
621  boost::mpi::broadcast( Environment::worldComm() , *this , Environment::worldComm().masterRank() );
622  }
623 
624 
628  boost::tuple<element_type,size_type>
629  min() const
630  {
631  element_type mumin( M_space );
632  mumin = M_space->max();
633 
634  element_type mu( M_space );
635  int index = 0;
636  int i = 0;
637  BOOST_FOREACH( mu, *this )
638  {
639  if ( mu.norm() < mumin.norm() )
640  {
641  mumin = mu;
642  index = i;
643  }
644 
645  ++i;
646  }
647  mumin.check();
648  return boost::make_tuple( mumin, index );
649  }
650 
654  boost::tuple<element_type,size_type>
655  max() const
656  {
657  element_type mumax( M_space );
658  mumax = M_space->min();
659 
660  element_type mu( M_space );
661  int index = 0;
662  int i = 0;
663  BOOST_FOREACH( mu, *this )
664  {
665  if ( mu.norm() > mumax.norm() )
666  {
667  mumax = mu;
668  index = i;
669  }
670 
671  ++i;
672  }
673  mumax.check();
674  return boost::make_tuple( mumax, index );
675  }
679  parameterspace_ptrtype parameterSpace() const
680  {
681  return M_space;
682  }
683 
687  sampling_ptrtype const& superSampling() const
688  {
689  return M_supersampling;
690  }
691 
695  void setSuperSampling( sampling_ptrtype const& super )
696  {
697  M_supersampling = super;
698  }
699 
706  sampling_ptrtype searchNearestNeighbors( element_type const& mu, size_type M , std::vector<int>& index_vector );
707 
711  sampling_ptrtype complement() const;
712 
716  void push_back( element_type const& mu, size_type index )
717  {
718  if ( M_supersampling ) M_superindices.push_back( index );
719 
720  super::push_back( mu );
721  }
722 
729  {
730  return M_superindices[ index ];
731  }
732  private:
733  Sampling() {}
734  private:
735  friend class boost::serialization::access;
736  template<class Archive>
737  void save( Archive & ar, const unsigned int version ) const
738  {
739  ar & boost::serialization::base_object<super>( *this );
740  ar & M_space;
741  ar & M_supersampling;
742  ar & M_superindices;
743  }
744 
745  template<class Archive>
746  void load( Archive & ar, const unsigned int version )
747  {
748  ar & boost::serialization::base_object<super>( *this );
749  ar & M_space;
750  ar & M_supersampling;
751  ar & M_superindices;
752  }
753  BOOST_SERIALIZATION_SPLIT_MEMBER()
754  private:
755  parameterspace_ptrtype M_space;
756 
757  sampling_ptrtype M_supersampling;
758  std::vector<size_type> M_superindices;
759 #if defined( FEELPP_HAS_ANN_H )
760  kdtree_ptrtype M_kdtree;
761 #endif
762  };
763 
764  typedef Sampling sampling_type;
765  typedef boost::shared_ptr<sampling_type> sampling_ptrtype;
766 
767  sampling_ptrtype sampling() { return sampling_ptrtype( new sampling_type( this->shared_from_this() ) ); }
768 
772 
775  :
776  M_min(),
777  M_max()
778  {
779  M_min.setZero();
780  M_max.setZero();
781  }
784  :
785  M_min( o.M_min ),
786  M_max( o.M_max )
787  {}
788 
789  ParameterSpace( element_type const& emin, element_type const& emax )
790  :
791  M_min( emin ),
792  M_max( emax )
793  {
794  M_min.setParameterSpace( this->shared_from_this() );
795  M_max.setParameterSpace( this->shared_from_this() );
796  }
799  {}
800 
804  static parameterspace_ptrtype New()
805  {
806  return parameterspace_ptrtype(new parameterspace_type);
807  }
809 
813 
816  {
817  if ( this != &o )
818  {
819  M_min = o.M_min;
820  M_max = o.M_max;
821  }
822 
823  return *this;
824  }
826 
830 
832  int dimension() const
833  {
834  return Dimension;
835  }
836 
840  element_type const& min() const
841  {
842  return M_min;
843  }
844 
848  element_type const& max() const
849  {
850  return M_max;
851  }
852 
856  element_type logMiddle() const
857  {
858  return ( ( M_min.array().log() + M_max.array().log() )/2. ).log();
859  }
860 
864  element_type middle() const
865  {
866  return ( ( M_min + M_max )/2. );
867  }
868 
870 
874 
878  void setMin( element_type const& min )
879  {
880  M_min = min;
881  }
882 
886  void setMax( element_type const& max )
887  {
888  M_max = max;
889  }
890 
892 
896 
900  static element_type logRandom( parameterspace_ptrtype space, bool broadcast )
901  {
902  //LOG(INFO) << "call logRandom...\n";
903  //LOG(INFO) << "call logRandom broadcast: " << broadcast << "...\n";
904  google::FlushLogFiles(google::GLOG_INFO);
905  if ( broadcast )
906  {
907  element_type mu( space );
908  if ( !space->check() ) return mu;
909  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
910  {
911  //LOG(INFO) << "generate random mu...\n";
912  //google::FlushLogFiles(google::GLOG_INFO);
913  mu = logRandom1( space );
914  }
915  //LOG(INFO) << "broadcast...\n";
916  //google::FlushLogFiles(google::GLOG_INFO);
917  boost::mpi::broadcast( Environment::worldComm() , mu , Environment::worldComm().masterRank() );
918  Environment::worldComm().barrier();
919  //LOG(INFO) << "check...\n";
920  mu.check();
921  return mu;
922  }
923  else
924  {
925  return logRandom1( space );
926  }
927  }
928  static element_type logRandom1( parameterspace_ptrtype space )
929  {
930  element_type mur( space );
931  mur.array() = element_type::Random().array().abs();
932  //LOG(INFO) << "random1 generate random mur= " << mur << " \n";
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();
936  //LOG(INFO) << "random1 generate random mu= " << mu << " \n";
937  return mu;
938  }
939 
943  static element_type random( parameterspace_ptrtype space, bool broadcast = true )
944  {
945  std::srand( static_cast<unsigned>( std::time( 0 ) ) );
946  element_type mur( space );
947  if ( broadcast )
948  {
949  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
950  {
951  mur.array() = element_type::Random().array().abs();
952  }
953  boost::mpi::broadcast( Environment::worldComm() , mur , Environment::worldComm().masterRank() );
954 
955  }
956  else
957  mur.array() = element_type::Random().array().abs();
958  //std::cout << "mur= " << mur << "\n";
959  //mur.setRandom()/RAND_MAX;
960  //std::cout << "mur= " << mur << "\n";
961  element_type mu( space );
962  mu.array() = space->min()+mur.array()*( space->max()-space->min() );
963  if ( broadcast )
964  mu.check();
965  return mu;
966  }
967 
974  static std::vector<double> logEquidistributeInDirection( parameterspace_ptrtype space, int direction , int N)
975  {
976  std::vector<double> result(N);
977 
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" );
983 
984  double min = min_element(direction);
985  double max = max_element(direction);
986 
987  for(int i=0; i<N; i++)
988  {
989  double factor = (double)(i)/(N-1);
990  result[i] = math::exp(math::log(min)+factor*( math::log(max)-math::log(min) ) );
991  }
992 
993  return result;
994  }
995 
996 
1001  static element_type logEquidistributed( double factor, parameterspace_ptrtype space )
1002  {
1003  element_type mu( space );
1004  mu.array() = ( space->min().array().log()+factor*( space->max().array().log()-space->min().array().log() ) ).exp();
1005  return mu;
1006  }
1007 
1008 
1015  static std::vector<double> equidistributeInDirection( parameterspace_ptrtype space, int direction , int N)
1016  {
1017  std::vector<double> result(N);
1018 
1019  auto min_element = space->min();
1020  auto max_element = space->max();
1021 
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" );
1025 
1026  double min = min_element(direction);
1027  double max = max_element(direction);
1028 
1029  for(int i=0; i<N; i++)
1030  {
1031  double factor = (double)(i)/(N-1);
1032  result[i] = min+factor*( max-min );
1033  }
1034 
1035  return result;
1036  }
1037 
1038 
1043  static element_type equidistributed( double factor, parameterspace_ptrtype space )
1044  {
1045  element_type mu( space );
1046  mu = space->min()+factor*( space->max()-space->min() );
1047  return mu;
1048  }
1049 
1050 
1052 
1053 
1054 
1055 protected:
1056 
1057 private:
1058  friend class boost::serialization::access;
1059  template<class Archive>
1060  void save( Archive & ar, const unsigned int version ) const
1061  {
1062  ar & M_min;
1063  ar & M_max;
1064  }
1065 
1066  template<class Archive>
1067  void load( Archive & ar, const unsigned int version )
1068  {
1069  ar & M_min;
1070  ar & M_max;
1071  }
1072  BOOST_SERIALIZATION_SPLIT_MEMBER()
1073 
1074  private:
1075 
1077  element_type M_min;
1078 
1080  element_type M_max;
1081 };
1082 
1083 template<int P> const uint16_type ParameterSpace<P>::Dimension;
1084 
1085 template<int P>
1086 //typename ParameterSpace<P>::sampling_ptrtype
1087 boost::shared_ptr<typename ParameterSpace<P>::Sampling>
1088 ParameterSpace<P>::Sampling::complement() const
1089 {
1090  //std::cout << "[ParameterSpace::Sampling::complement] start\n";
1091  sampling_ptrtype complement;
1092 
1093  if ( M_supersampling )
1094  {
1095  // std::cout << "[ParameterSpace::Sampling::complement] super sampling available\n";
1096  // std::cout << "[ParameterSpace::Sampling::complement] this->size = " << this->size() << std::endl;
1097  // std::cout << "[ParameterSpace::Sampling::complement] supersampling->size = " << M_supersampling->size() << std::endl;
1098  // std::cout << "[ParameterSpace::Sampling::complement] complement->size = " << M_supersampling->size()-this->size() << std::endl;
1099  complement = sampling_ptrtype( new sampling_type( M_space, 1, M_supersampling ) );
1100  complement->clear();
1101 
1102  if ( !M_superindices.empty() )
1103  {
1104  for ( size_type i = 0, j = 0; i < M_supersampling->size(); ++i )
1105  {
1106  if ( std::find( M_superindices.begin(), M_superindices.end(), i ) == M_superindices.end() )
1107  {
1108  //std::cout << "inserting " << j << "/" << i << " in complement of size (" << complement->size() << " out of " << M_supersampling->size() << std::endl;
1109  if ( j < M_supersampling->size()-this->size() )
1110  {
1111  //complement->at( j ) = M_supersampling->at( i );
1112  complement->push_back( M_supersampling->at( i ),i );
1113  ++j;
1114  }
1115  }
1116  }
1117  }
1118 
1119  return complement;
1120  }
1121 
1122  return complement;
1123 }
1124 template<int P>
1125 boost::shared_ptr<typename ParameterSpace<P>::Sampling>
1127  size_type _M ,
1128  std::vector<int>& index_vector)
1129 {
1130  size_type M=_M;
1131  sampling_ptrtype neighbors( new sampling_type( M_space, M ) );
1132 #if defined(FEELPP_HAS_ANN_H)
1133 
1134  if( Environment::worldComm().globalRank() == Environment::worldComm().masterRank() )
1135  {
1136  //std::cout << "[ParameterSpace::Sampling::searchNearestNeighbors] start\n";
1137  //if ( !M_kdtree )
1138  //{
1139  //std::cout << "[ParameterSpace::Sampling::searchNearestNeighbors] building data points\n";
1140  ANNpointArray data_pts;
1141  data_pts = annAllocPts( this->size(), M_space->Dimension );
1142 
1143  for ( size_type i = 0; i < this->size(); ++i )
1144  {
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" );
1147  }
1148 
1149  //std::cout << "[ParameterSpace::Sampling::searchNearestNeighbors] building tree in R^" << M_space->Dimension << "\n";
1150  M_kdtree = kdtree_ptrtype( new kdtree_type( data_pts, this->size(), M_space->Dimension ) );
1151  //}
1152 
1153 
1154  // make sure that the sampling set is big enough
1155  if ( this->size() < M )
1156  M=this->size();
1157 
1158  std::vector<int> nnIdx( M );
1159  std::vector<double> dists( M );
1160  double eps = 0;
1161  M_kdtree->annkSearch( // search
1162  const_cast<double*>( mu.data() ), // query point
1163  M, // number of near neighbors
1164  nnIdx.data(), // nearest neighbors (returned)
1165  dists.data(), // distance (returned)
1166  eps ); // error bound
1167 
1168 
1169 
1170  if ( M_supersampling )
1171  {
1172  neighbors->setSuperSampling( M_supersampling );
1173 
1174  if ( !M_superindices.empty() )
1175  neighbors->M_superindices.resize( M );
1176  }
1177 
1178  //std::cout << "[parameterspace::sampling::searchNearestNeighbors] neighbor size = " << neighbors->size() <<std::endl;
1179  for ( size_type i = 0; i < M; ++i )
1180  {
1181  //std::cout << "[parameterspace::sampling::searchNearestNeighbors] neighbor: " <<i << " distance = " << dists[i] << "\n";
1182  neighbors->at( i ) = this->at( nnIdx[i] );
1183  index_vector.push_back( nnIdx[i] );
1184 
1185  //std::cout<<"[parameterspace::sampling::searchNearestNeighbors] M_superindices.size() : "<<M_superindices.size()<<std::endl;
1186  if ( M_supersampling && !M_superindices.empty() )
1187  neighbors->M_superindices[i] = M_superindices[ nnIdx[i] ];
1188  //std::cout << "[parameterspace::sampling::searchNearestNeighbors] " << neighbors->at( i ) << "\n";
1189  }
1190  annDeallocPts( data_pts );
1191  }//end of procMaster
1192 
1193  boost::mpi::broadcast( Environment::worldComm() , neighbors , Environment::worldComm().masterRank() );
1194  boost::mpi::broadcast( Environment::worldComm() , index_vector , Environment::worldComm().masterRank() );
1195 
1196 
1197 #endif /* FEELPP_HAS_ANN_H */
1198 
1199  return neighbors;
1200 }
1201 } // Feel
1202 #endif /* __ParameterSpace_H */

Generated on Sun Oct 20 2013 08:25:03 for Feel++ by doxygen 1.8.4