diff --git a/CMakeLists.txt b/CMakeLists.txt index 715fc0cf..da961550 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -25,9 +25,13 @@ INCLUDE(FindPkgConfig) PKG_CHECK_MODULES(EO eo REQUIRED) PKG_CHECK_MODULES(MO mo REQUIRED) +FIND_PACKAGE(Boost 1.33.0) + INCLUDE_DIRECTORIES( ${EO_INCLUDE_DIRS} ${MO_INCLUDE_DIRS} + ${Boost_INCLUDE_DIRS} + /Dev/ometah-0.3/common ) ###################################################################################### diff --git a/application/cma_sa/CMakeLists.txt b/application/cma_sa/CMakeLists.txt index d68a49c0..c67ab0dd 100644 --- a/application/cma_sa/CMakeLists.txt +++ b/application/cma_sa/CMakeLists.txt @@ -1,7 +1,12 @@ PROJECT(cma_sa) +FIND_PACKAGE(Boost 1.33.0) + INCLUDE_DIRECTORIES(${CMAKE_CURRENT_SOURCE_DIR}) +INCLUDE_DIRECTORIES(${Boost_INCLUDE_DIRS}) +LINK_DIRECTORIES(${Boost_LIBRARY_DIRS}) + SET(RESOURCES cma_sa.param plot.py @@ -20,4 +25,4 @@ FILE(GLOB SOURCES *.cpp) SET(EXECUTABLE_OUTPUT_PATH ${DO_BINARY_DIR}) ADD_EXECUTABLE(${PROJECT_NAME} ${SOURCES}) -TARGET_LINK_LIBRARIES(${PROJECT_NAME} do ${EO_LIBRARIES} ${MO_LIBRARIES}) +TARGET_LINK_LIBRARIES(${PROJECT_NAME} do ${EO_LIBRARIES} ${MO_LIBRARIES} ${Boost_LIBRARIES}) diff --git a/application/cma_sa/main.cpp b/application/cma_sa/main.cpp index f0864a68..872636e7 100644 --- a/application/cma_sa/main.cpp +++ b/application/cma_sa/main.cpp @@ -52,7 +52,8 @@ int main(int ac, char** av) eoEvalFunc< EOT >* plainEval = new Sphere< EOT >(); state.storeFunctor(plainEval); - eoEvalFuncCounter< EOT > eval(*plainEval); + unsigned long max_eval = parser.getORcreateParam((unsigned long)0, "maxEval", "Maximum number of evaluations (0 = none)", 'E', "Stopping criterion").value(); // E + eoEvalFuncCounter< EOT > eval(*plainEval, max_eval); eoRndGenerator< double >* gen = new eoUniformGenerator< double >(-5, 5); //eoRndGenerator< double >* gen = new eoNormalGenerator< double >(0, 1); @@ -122,31 +123,20 @@ int main(int ac, char** av) eoCheckPoint< EOT >& checkpoint = do_make_checkpoint(parser, state, eval, monitoring_continue); - // appends some missing code to checkpoint + // eoPopStat< EOT >* popStat = new eoPopStat; + // state.storeFunctor(popStat); - // eoValueParam& plotPopParam = parser.createParam(false, "plotPop", "Plot sorted pop. every gen.", 0, "Graphical Output"); + // checkpoint.add(*popStat); - // if (plotPopParam.value()) // we do want plot dump - // { - // eoScalarFitnessStat* fitStat = new eoScalarFitnessStat; - // state.storeFunctor(fitStat); + // eoGnuplot1DMonitor* gnuplot = new eoGnuplot1DMonitor("gnuplot.txt"); + // state.storeFunctor(gnuplot); - // checkpoint.add(*fitStat); + // gnuplot->add(eval); + // gnuplot->add(*popStat); - // eoFileSnapshot* snapshot = new eoFileSnapshot("ResPop"); - // state.storeFunctor(snapshot); + //gnuplot->gnuplotCommand("set yrange [0:500]"); - // snapshot->add(*fitStat); - - // checkpoint.add(*snapshot); - // } - - // -------------------------- - - eoPopStat< EOT >* popStat = new eoPopStat; - state.storeFunctor(popStat); - - checkpoint.add(*popStat); + // checkpoint.add(*gnuplot); // eoMonitor* fileSnapshot = new doFileSnapshot< std::vector< std::string > >("ResPop"); // state.storeFunctor(fileSnapshot); @@ -206,6 +196,10 @@ int main(int ac, char** av) { do_run(*algo, pop); } + catch (eoReachedThresholdException& e) + { + eo::log << eo::warnings << e.what() << std::endl; + } catch (std::exception& e) { eo::log << eo::errors << "exception: " << e.what() << std::endl; diff --git a/application/cma_sa/plot.py b/application/cma_sa/plot.py old mode 100755 new mode 100644 diff --git a/copying b/copying new file mode 100644 index 00000000..90cbe47b --- /dev/null +++ b/copying @@ -0,0 +1 @@ +Private license diff --git a/matrix.hpp b/matrix.hpp new file mode 100644 index 00000000..89b77d85 --- /dev/null +++ b/matrix.hpp @@ -0,0 +1,560 @@ +/*************************************************************************** + * $Id: matrix.hpp,v 1.11 2006/05/13 10:05:53 nojhan Exp $ + * Copyright : Free Software Foundation + * Author : Johann Dréo + ****************************************************************************/ + +/* + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU Lesser General Public License as published by + * the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + */ + +#ifndef MATRIX +#define MATRIX + +#include +#include + +#include "Exception_oMetah.hpp" + +using namespace std; + +namespace ometah { + +//! Test if a vector is comprised in bounds +template +bool isInBounds( T aVector, T mins, T maxs) +{ + unsigned int i; + for(i=0; i maxs[i] ){ + return false; + } + } + return true; +} + + +//! Force a vector to be in bounds +template +T forceBounds( T aVector, T mins, T maxs) +{ + T CastedVector=aVector; + + unsigned int i; + for(i=0; i maxs[i] ){ + CastedVector[i]=maxs[i]; + } + } + return CastedVector; +} + +//! Create a 2D matrix filled with values +/* + if we want a vector > : + T stand for double + V stand for vector > +*/ +template +U matrixFilled( unsigned int dimL, unsigned int dimC, T fillValue ) +{ + unsigned int i; + + // make the vector possible at this step + typename U::value_type vec(dimC, fillValue); + + U mat; + for(i=0; i +vector > matrixFilled( unsigned int dimL, unsigned int dimC, T fillValue ) +{ + unsigned int i; + + // make the vector possible at this step + vector< T > vec(dimC, fillValue); + + vector > mat; + for(i=0; i +T multiply( T matA, T matB) +{ + + T newMat; + + unsigned int Al=matA.size(); + unsigned int Ac=matA[0].size(); + unsigned int Bl=matB.size(); + unsigned int Bc=matB[0].size(); + + newMat=matrixFilled( Al,Bc,0.0); + + if(Ac!=Bl) { + throw Exception_Size_Match("Cannot multiply matrices, sizes does not match", EXCEPTION_INFOS ); + } + + for( unsigned int i=0; i +U multiply(U aVector, T aNb) +{ + U res; + + res.reserve( aVector.size() ); + + unsigned int i; + for(i=0; i +T cholesky( T A) +{ + + // FIXME : vérifier que A est symétrique définie positive + + T B; + unsigned int Al=A.size(); + unsigned int Ac=A[0].size(); + B = matrixFilled(Al, Ac, 0.0); + + unsigned int i,j,k; + + // first column + i=0; + + // diagonal + j=0; + B[0][0]=sqrt(A[0][0]); + + // end of the column + for(j=1;j +T transpose( T &mat) +{ + unsigned int iSize=mat.size(); + unsigned int jSize=mat[0].size(); + + if ( iSize == 0 || jSize == 0 ) { + ostringstream msg; + msg << "ErrorSize: matrix not defined " + << "(iSize:" << iSize << ", jSize:" << jSize << ")"; + throw Exception_Size( msg.str(), EXCEPTION_INFOS ); + } + + typename T::value_type aVector; + T newMat; + + unsigned int i, j; + + for (j=0; j +vector mean( vector > mat) +{ + vector moyDim; + moyDim.reserve(mat.size()); + + unsigned int i,a; + a=mat.size(); + + for(i=0;i +T mean( vector aVector, unsigned int begin=0, unsigned int during=0) +{ + if (during==0) { + during = aVector.size() - begin; // if no end : take all + } + + T aSum, aMean; + + aSum = sum(aVector, begin, during); // Sum + aMean = aSum / (during - begin); // Mean + + return aMean; +} + +//! Calculate a variance-covariance matrix from a list of vector +/*! + For a population of p points on n dimensions : + if onRow==true, the matrix should have p rows and n columns. + if onRow==false, the matrix should have n rows and p columns. +*/ +template +U varianceCovariance( U pop, bool onRow = true) +{ +/* + // vector of means + typename U::value_type vecMeanCentered; + if(onRow) { + vecMeanCentered = mean( transpose(pop) ); // p rows and n columns => means of p + } else { + vecMeanCentered = mean( pop ); // n rows and p columns => means of n + } + + // centered population + // same size as the initial matrix + U popMeanCentered = matrixFilled(pop.size(),pop[0].size(), 0.0); + + // centering + // rows + for(unsigned int i=0;i covariance of p + } else { + popVar = multiply( popMeanCentered, popMeanCenteredT ); // if n rows and p columns => covariance of n + } + + // multiplication by 1/n : + for(unsigned int i=0;i +T sum(vector aVector, unsigned int begin=0, unsigned int during=0) +{ + if ( begin > aVector.size() || during > aVector.size() ) { + ostringstream msg; + msg << "ErrorSize: parameters are out of vector bounds " + << "(begin:" << begin << ", during:" << during + << ", size:" << aVector.size() << ")"; + throw Exception_Size_Index( msg.str(), EXCEPTION_INFOS ); + } + + if (during==0) { + during = aVector.size() - begin; + } + + T aSum=0; + + for (unsigned int j=begin; j +T stdev(vector aVector, unsigned int begin=0, unsigned int during=0) +{ + if ( begin > aVector.size() || during > aVector.size() ) { + ostringstream msg; + msg << "ErrorSize: parameters are out of vector bounds " + << "(begin:" << begin << ", during:" << during + << ", size:" << aVector.size() << ")"; + throw Exception_Size_Index( msg.str(), EXCEPTION_INFOS ); + } + + if (during==0) { + during = aVector.size() - begin; + } + + vector deviation; + T aMean, aDev, aStd; + + aMean = mean(aVector, begin, during); // mean + + for (unsigned int j=begin; j +typename T::value_type min(T aVector, unsigned int begin=0, unsigned int during=0) +{ + if ( begin > aVector.size() || during > aVector.size() ) { + ostringstream msg; + msg << "ErrorSize: parameters are out of vector bounds " + << "(begin:" << begin << ", during:" << during + << ", size:" << aVector.size() << ")"; + throw Exception_Size_Index( msg.str(), EXCEPTION_INFOS ); + } + + if (during==0) { + during = aVector.size() - begin; + } + + typename T::value_type aMin = aVector[begin]; + + for (unsigned int i=begin+1; i +vector mins(vector > aMatrix) +{ + vector mins; + + for( unsigned int i=0; i < aMatrix.size(); i++ ) { + mins.push_back( min(aMatrix[i]) ); + } + + return mins; +} + +//! Find the maximums values of a matrix, for each row +template +vector maxs(vector > aMatrix) +{ + vector maxs; + + for( unsigned int i=0; i < aMatrix.size(); i++ ) { + maxs.push_back( max(aMatrix[i]) ); + } + + return maxs; +} + +//! Find the maximum value of a vector +template +typename T::value_type max(T aVector, unsigned int begin=0, unsigned int during=0) +{ + if ( begin > aVector.size() || during > aVector.size() ) { + ostringstream msg; + msg << "ErrorSize: parameters are out of vector bounds " + << "(begin:" << begin << ", during:" << during + << ", size:" << aVector.size() << ")"; + throw Exception_Size_Index( msg.str(), EXCEPTION_INFOS ); + } + + if (during==0) { + during = aVector.size() - begin; + } + + typename T::value_type aMax = aVector[begin]; + + for (unsigned int i=begin+1; i aMax ) { + aMax = aVector[i]; + } + } + + return aMax; +} + +//! Substraction of two vectors, terms by terms +template +T substraction(T from, T that) +{ + T res; + + res.reserve(from.size()); + + for(unsigned int i=0; i +T addition(T from, T that) +{ + T res; + + res.reserve( from.size() ); + + for(unsigned int i=0; i +T absolute(T aVector) +{ + for(unsigned int i=0; i +vector gravityCenter( vector > points, vector weights ) +{ + + // if we have only one weight, we use it for all items + if ( weights.size() == 1 ) { + for ( unsigned int i=1; i < points.size(); i++ ) { + weights.push_back( weights[0] ); + } + } + + // if sizes does not match : error + if ( points.size() != weights.size() ) { + ostringstream msg; + msg << "ErrorSize: " + << "points size (" << points.size() << ")" + << " does not match weights size (" << weights.size() << ")"; + throw Exception_Size_Match( msg.str(), EXCEPTION_INFOS ); + } + + T weightsSum = sum(weights); + + vector > pointsT = transpose( points ); + + vector gravity; + + for ( unsigned int i=0; i < pointsT.size(); i++ ) { // dimensions + T g = 0; + for ( unsigned int j=0; j < pointsT[i].size(); j++ ) { // points + g += ( pointsT[i][j] * weights[j] ) / weightsSum; + } + gravity.push_back( g ); + } + + return gravity; +} + +} // ometah + +#endif // MATRIX diff --git a/readme b/readme new file mode 100644 index 00000000..64629f58 --- /dev/null +++ b/readme @@ -0,0 +1,57 @@ +This package contains the source code for BOPO problems. + +# Step 1 - Configuration +------------------------ +Rename the "install.cmake-dist" file as "install.cmake" and edit it, inserting the FULL PATH +to your ParadisEO distribution. +On Windows write your path with double antislash (ex: C:\\Users\\...) + + +# Step 2 - Build process +------------------------ +ParadisEO is assumed to be compiled. To download ParadisEO, please visit http://paradiseo.gforge.inria.fr/. +Go to the BOPO/build/ directory and lunch cmake: +(Unix) > cmake .. +(Windows) > cmake .. -G"Visual Studio 9 2008" + +Note for windows users: if you don't use VisualStudio 9, enter the name of your generator instead of "VisualStudio 9 2008". + + +# Step 3 - Compilation +---------------------- +In the bopo/build/ directory: +(Unix) > make +(Windows) Open the VisualStudio solution and compile it, compile also the target install. +You can refer to this tutorial if you don't know how to compile a solution: http://paradiseo.gforge.inria.fr/index.php?n=Paradiseo.VisualCTutorial + + +# Step 4 - Execution +--------------------- +A toy example is given to test the components. You can run these tests as following. +To define problem-related components for your own problem, please refer to the tutorials available on the website : http://paradiseo.gforge.inria.fr/. +In the bopo/build/ directory: +(Unix) > ctest +Windows users, please refer to this tutorial: http://paradiseo.gforge.inria.fr/index.php?n=Paradiseo.VisualCTutorial + +In the directory "application", there are several directory such as p_eoco which instantiate NSGAII on BOPO problems. + +(Unix) After compilation you can run the script "bopo/run.sh" and see results in "NSGAII.out". Parameters can be modified in the script. + +(Windows) Add argument "NSGAII.param" and execute the corresponding algorithms. +Windows users, please refer to this tutorial: http://paradiseo.gforge.inria.fr/index.php?n=Paradiseo.VisualCTutorial + + +# Documentation +--------------- +The API-documentation is available in doc/html/index.html + + +# Things to keep in mind when using BOPO +---------------------------------------- +* By default, the EO random generator's seed is initialized by the number of seconds since the epoch (with time(0)). It is available in the status file dumped at each execution. Please, keep in mind that if you start two run at the same second without modifying the seed, you will get exactly the same results. + +* Execution times are measured with the boost:timer, that measure wallclock time. Additionaly, it could not measure times larger than approximatively 596.5 hours (or even less). See http://www.boost.org/doc/libs/1_33_1/libs/timer/timer.htm + +* The q-quantile computation use averaging at discontinuities (in R, it correspond to the R-2 method, in SAS, SAS-5). For more explanations, see http://en.wikipedia.org/wiki/Quantile#Estimating_the_quantiles_of_a_population and http://stat.ethz.ch/R-manual/R-devel/library/stats/html/quantile.html + +* You can send a SIGUSR1 to a process to get some information (written down in the log file) on the current state of the search. diff --git a/src/doCMASA.h b/src/doCMASA.h index 1e28b783..7454f37c 100644 --- a/src/doCMASA.h +++ b/src/doCMASA.h @@ -86,8 +86,51 @@ public: _continuator(continuator), _cooling_schedule(cooling_schedule), _initial_temperature(initial_temperature), - _replacor(replacor) - {} + _replacor(replacor), + + _pop_results_destination("ResPop"), + + // directory where populations state are going to be stored. + _ofs_params("ResParams.txt"), + _ofs_params_var("ResVars.txt"), + + _bounds_results_destination("ResBounds") + { + + //------------------------------------------------------------- + // Temporary solution to store populations state at each + // iteration for plotting. + //------------------------------------------------------------- + + { + std::stringstream ss; + ss << "rm -rf " << _pop_results_destination; + ::system(ss.str().c_str()); + } + + ::mkdir(_pop_results_destination.c_str(), 0755); // create a first time the + + //------------------------------------------------------------- + + + //------------------------------------------------------------- + // Temporary solution to store bounds values for each distribution. + //------------------------------------------------------------- + + { + std::stringstream ss; + ss << "rm -rf " << _bounds_results_destination; + ::system(ss.str().c_str()); + } + +<<<<<<< HEAD + ::mkdir(bounds_results_destination.c_str(), 0755); // create once directory +======= + ::mkdir(_bounds_results_destination.c_str(), 0755); // create once directory + + //------------------------------------------------------------- + + } //! function that launches the CMASA algorithm. /*! @@ -114,44 +157,7 @@ public: //------------------------------------------------------------- int number_of_iterations = 0; - - //------------------------------------------------------------- - - - //------------------------------------------------------------- - // Temporary solution to store populations state at each - // iteration for plotting. - //------------------------------------------------------------- - - std::string pop_results_destination("ResPop"); - - { - std::stringstream ss; - ss << "rm -rf " << pop_results_destination; - ::system(ss.str().c_str()); - } - - ::mkdir(pop_results_destination.c_str(), 0755); // create a first time the - // directory where populations state are going to be stored. - std::ofstream ofs_params("ResParams.txt"); - std::ofstream ofs_params_var("ResVars.txt"); - - //------------------------------------------------------------- - - - //------------------------------------------------------------- - // Temporary solution to store bounds values for each distribution. - //------------------------------------------------------------- - - std::string bounds_results_destination("ResBounds"); - - { - std::stringstream ss; - ss << "rm -rf " << bounds_results_destination; - ::system(ss.str().c_str()); - } - - ::mkdir(bounds_results_destination.c_str(), 0755); // create once directory +>>>>>>> 36ec42d36204631eb4c25ae7b31a8728903697f8 //------------------------------------------------------------- @@ -166,10 +172,10 @@ public: for (int i = 0; i < size; ++i) { - hv.update( distrib.variance()[i] ); + //hv.update( distrib.varcovar()[i] ); } - ofs_params_var << hv << std::endl; + _ofs_params_var << hv << std::endl; } do @@ -260,7 +266,7 @@ public: { std::stringstream ss; - ss << pop_results_destination << "/" << number_of_iterations; + ss << _pop_results_destination << "/" << number_of_iterations; std::ofstream ofs(ss.str().c_str()); ofs << current_pop; } @@ -279,12 +285,16 @@ public: assert(size > 0); std::stringstream ss; - ss << bounds_results_destination << "/" << number_of_iterations; + ss << _bounds_results_destination << "/" << number_of_iterations; std::ofstream ofs(ss.str().c_str()); ofs << size << " "; std::copy(distrib.mean().begin(), distrib.mean().end(), std::ostream_iterator< double >(ofs, " ")); +<<<<<<< HEAD std::copy(distrib.variance().begin(), distrib.variance().end(), std::ostream_iterator< double >(ofs, " ")); +======= + //std::copy(distrib.varcovar().begin(), distrib.varcovar().end(), std::ostream_iterator< double >(ofs, " ")); +>>>>>>> 36ec42d36204631eb4c25ae7b31a8728903697f8 ofs << std::endl; } @@ -334,10 +344,10 @@ public: for (int i = 0; i < size; ++i) { - hv.update( distrib.variance()[i] ); + //hv.update( distrib.varcovar()[i] ); } - ofs_params_var << hv << std::endl; + _ofs_params_var << hv << std::endl; } //------------------------------------------------------------- @@ -383,6 +393,27 @@ private: //! A EOT replacor eoReplacement < EOT > & _replacor; + + + //------------------------------------------------------------- + // Temporary solution to store populations state at each + // iteration for plotting. + //------------------------------------------------------------- + + std::string _pop_results_destination; + std::ofstream _ofs_params; + std::ofstream _ofs_params_var; + + //------------------------------------------------------------- + + //------------------------------------------------------------- + // Temporary solution to store bounds values for each distribution. + //------------------------------------------------------------- + + std::string _bounds_results_destination; + + //------------------------------------------------------------- + }; #endif // !_doCMASA_h diff --git a/src/doEstimatorNormal.h b/src/doEstimatorNormal.h index 3bd2058f..010008a8 100644 --- a/src/doEstimatorNormal.h +++ b/src/doEstimatorNormal.h @@ -9,6 +9,8 @@ template < typename EOT > class doEstimatorNormal : public doEstimator< doNormal< EOT > > { public: + typedef typename EOT::AtomType AtomType; + doNormal< EOT > operator()(eoPop& pop) { unsigned int popsize = pop.size(); @@ -17,23 +19,12 @@ public: unsigned int dimsize = pop[0].size(); assert(dimsize > 0); - doCovMatrix cov(dimsize); + //doCovMatrix cov(dimsize); + doUblasCovMatrix< EOT > cov; - for (unsigned int i = 0; i < popsize; ++i) - { - cov.update(pop[i]); - } + cov.update(pop); - EOT mean(dimsize); - EOT covariance(dimsize); - - for (unsigned int d = 0; d < dimsize; ++d) - { - mean[d] = cov.get_mean(d); - covariance[d] = cov.get_var(d); - } - - return doNormal< EOT >(mean, covariance); + return doNormal< EOT >( cov.get_mean(), cov.get_varcovar() ); } }; diff --git a/src/doNormal.h b/src/doNormal.h index 6de887fc..f857da75 100644 --- a/src/doNormal.h +++ b/src/doNormal.h @@ -8,8 +8,10 @@ template < typename EOT > class doNormal : public doDistrib< EOT >, public doNormalParams< EOT > { public: - doNormal(EOT mean, EOT variance) - : doNormalParams< EOT >(mean, variance) + typedef typename EOT::AtomType AtomType; + + doNormal( const ublas::vector< AtomType >& mean, const ublas::symmetric_matrix< AtomType, ublas::lower >& varcovar ) + : doNormalParams< EOT >( mean, varcovar ) {} }; diff --git a/src/doNormalCenter.h b/src/doNormalCenter.h index 06113060..ff379ba2 100644 --- a/src/doNormalCenter.h +++ b/src/doNormalCenter.h @@ -12,7 +12,9 @@ public: void operator() ( doNormal< EOT >& distrib, EOT& mass ) { - distrib.mean() = mass; // vive les references!!! + ublas::vector< AtomType > mean( distrib.size() ); + std::copy( mass.begin(), mass.end(), mean.begin() ); + distrib.mean() = mean; } }; diff --git a/src/doNormalParams.h b/src/doNormalParams.h index 87a3a65a..557c5443 100644 --- a/src/doNormalParams.h +++ b/src/doNormalParams.h @@ -1,10 +1,19 @@ #ifndef _doNormalParams_h #define _doNormalParams_h +<<<<<<< HEAD +======= +#include +#include + +namespace ublas = boost::numeric::ublas; + +>>>>>>> 36ec42d36204631eb4c25ae7b31a8728903697f8 template < typename EOT > class doNormalParams { public: +<<<<<<< HEAD doNormalParams(EOT mean, EOT variance) : _mean(mean), _variance(variance) { @@ -18,12 +27,40 @@ public: unsigned int size() { assert(_mean.size() == _variance.size()); +======= + typedef typename EOT::AtomType AtomType; + + doNormalParams + ( + const ublas::vector< AtomType >& mean, + const ublas::symmetric_matrix< AtomType, ublas::lower >& varcovar + ) + : _mean(mean), _varcovar(varcovar) + { + assert(_mean.size() > 0); + assert(_mean.size() == _varcovar.size1()); + assert(_mean.size() == _varcovar.size2()); + } + + ublas::vector< AtomType >& mean(){return _mean;} + ublas::symmetric_matrix< AtomType, ublas::lower >& varcovar(){return _varcovar;} + + unsigned int size() + { + assert(_mean.size() == _varcovar.size1()); + assert(_mean.size() == _varcovar.size2()); +>>>>>>> 36ec42d36204631eb4c25ae7b31a8728903697f8 return _mean.size(); } private: +<<<<<<< HEAD EOT _mean; EOT _variance; +======= + ublas::vector< AtomType > _mean; + ublas::symmetric_matrix< AtomType, ublas::lower > _varcovar; +>>>>>>> 36ec42d36204631eb4c25ae7b31a8728903697f8 }; #endif // !_doNormalParams_h diff --git a/src/doSamplerNormal.h b/src/doSamplerNormal.h index 7ba4c2d5..413f7ada 100644 --- a/src/doSamplerNormal.h +++ b/src/doSamplerNormal.h @@ -6,6 +6,7 @@ #include "doSampler.h" #include "doNormal.h" #include "doBounder.h" +#include "doStats.h" /** * doSamplerNormal @@ -29,30 +30,61 @@ public: assert(size > 0); //------------------------------------------------------------- - // Point we want to sample to get higher a set of points - // (coordinates in n dimension) - // x = {x1, x2, ..., xn} + // Cholesky factorisation gererating matrix L from covariance + // matrix V. + // We must use cholesky.get_L() to get the resulting matrix. + // + // L = cholesky decomposition of varcovar //------------------------------------------------------------- - EOT solution; + Cholesky< EOT > cholesky; + cholesky.update( distrib.varcovar() ); + ublas::symmetric_matrix< AtomType, ublas::lower > L = cholesky.get_L(); //------------------------------------------------------------- //------------------------------------------------------------- - // Sampling all dimensions + // T = vector of size elements drawn in N(0,1) rng.normal(1.0) //------------------------------------------------------------- - for (unsigned int i = 0; i < size; ++i) + ublas::vector< AtomType > T( size ); + + for ( unsigned int i = 0; i < size; ++i ) { - solution.push_back( - rng.normal(distrib.mean()[i], - distrib.variance()[i]) - ); + T( i ) = rng.normal( 1.0 ); } //------------------------------------------------------------- + + //------------------------------------------------------------- + // LT = prod( L, trans(T) ) ? + // LT = prod( L, T ) + //------------------------------------------------------------- + + //ublas::symmetric_matrix< AtomType, ublas::lower > LT = ublas::prod( L, ublas::trans( T ) ); + ublas::vector< AtomType > LT = ublas::prod( L, T ); + + //------------------------------------------------------------- + + + //------------------------------------------------------------- + // solution = means + trans( LT ) ? + // solution = means + LT + //------------------------------------------------------------- + + ublas::vector< AtomType > mean = distrib.mean(); + + ublas::vector< AtomType > ublas_solution = mean + LT; + //ublas::vector< AtomType > ublas_solution = mean + ublas::trans( LT ); + + EOT solution( size ); + + std::copy( ublas_solution.begin(), ublas_solution.end(), solution.begin() ); + + //------------------------------------------------------------- + return solution; } }; diff --git a/src/doSamplerUniform.h b/src/doSamplerUniform.h index f49b6924..374d7a8e 100644 --- a/src/doSamplerUniform.h +++ b/src/doSamplerUniform.h @@ -20,8 +20,10 @@ public: unsigned int size = distrib.size(); assert(size > 0); + //------------------------------------------------------------- - // Point we want to sample to get higher a population + // Point we want to sample to get higher a set of points + // (coordinates in n dimension) // x = {x1, x2, ..., xn} //------------------------------------------------------------- diff --git a/src/doStats.h b/src/doStats.h index 94756fa3..50cc9c95 100644 --- a/src/doStats.h +++ b/src/doStats.h @@ -18,7 +18,14 @@ #ifndef _doStats_h #define _doStats_h +#include +#include +#include + #include +#include + +namespace ublas = boost::numeric::ublas; class doStats : public eoPrintable { @@ -110,9 +117,159 @@ protected: double _hv; }; +<<<<<<< HEAD class doCholesky : public doStats { +======= +template < typename EOT > +class doUblasCovMatrix : public doStats +{ +public: + typedef typename EOT::AtomType AtomType; + + virtual void update( const eoPop< EOT >& pop ) + { + unsigned int p_size = pop.size(); // population size + + assert(p_size > 0); + + unsigned int s_size = pop[0].size(); // solution size + + assert(s_size > 0); + + ublas::matrix< AtomType > sample( p_size, s_size ); + + for (unsigned int i = 0; i < p_size; ++i) + { + for (unsigned int j = 0; j < s_size; ++j) + { + sample(i, j) = pop[i][j]; + } + } + + _varcovar.resize(s_size, s_size); + + + //------------------------------------------------------------- + // variance-covariance matrix are symmetric (and semi-definite + // positive), thus a triangular storage is sufficient + //------------------------------------------------------------- + + ublas::symmetric_matrix< AtomType, ublas::lower > var(s_size, s_size); + + //------------------------------------------------------------- + + + //------------------------------------------------------------- + // variance-covariance matrix computation : A * transpose(A) + //------------------------------------------------------------- + + var = ublas::prod( sample, ublas::trans( sample ) ); + + //------------------------------------------------------------- + + + for (unsigned int i = 0; i < s_size; ++i) + { + // triangular LOWER matrix, thus j is not going further than i + for (unsigned int j = 0; j <= i; ++j) + { + // we want a reducted covariance matrix + _varcovar(i, j) = var(i, j) / p_size; + } + } + + _mean.resize(s_size); + + // unit vector + ublas::scalar_vector< AtomType > u( p_size, 1 ); + + // sum over columns + _mean = ublas::prod( ublas::trans( sample ), u ); + + // division by n + _mean /= p_size; + } + + const ublas::symmetric_matrix< AtomType, ublas::lower >& get_varcovar() const {return _varcovar;} + + const ublas::vector< AtomType >& get_mean() const {return _mean;} + +private: + ublas::symmetric_matrix< AtomType, ublas::lower > _varcovar; + ublas::vector< AtomType > _mean; +}; + +template < typename EOT > +class Cholesky : public doStats +{ +public: + typedef typename EOT::AtomType AtomType; + + virtual void update( const ublas::symmetric_matrix< AtomType, ublas::lower >& V) + { + unsigned int Vl = V.size1(); + + assert(Vl > 0); + + unsigned int Vc = V.size2(); + + assert(Vc > 0); + + _L.resize(Vl, Vc); + + unsigned int i,j,k; + + // first column + i=0; + + // diagonal + j=0; + _L(0, 0) = sqrt( V(0, 0) ); + + // end of the column + for ( j = 1; j < Vc; ++j ) + { + _L(j, 0) = V(0, j) / _L(0, 0); + } + + // end of the matrix + for ( i = 1; i < Vl; ++i ) // each column + { + + // diagonal + double sum = 0.0; + + for ( k = 0; k < i; ++k) + { + sum += _L(i, k) * _L(i, k); + } + + assert( ( V(i, i) - sum ) > 0 ); + + _L(i, i) = sqrt( V(i, i) - sum ); + + for ( j = i + 1; j < Vl; ++j ) // rows + { + // one element + sum = 0.0; + + for ( k = 0; k < i; ++k ) + { + sum += _L(j, k) * _L(i, k); + } + + _L(j, i) = (V(j, i) - sum) / _L(i, i); + } + } + } + + const ublas::symmetric_matrix< AtomType, ublas::lower >& get_L() const {return _L;} + +private: + ublas::symmetric_matrix< AtomType, ublas::lower > _L; +>>>>>>> 36ec42d36204631eb4c25ae7b31a8728903697f8 }; #endif // !_doStats_h diff --git a/src/todo b/src/todo new file mode 100644 index 00000000..13785cc8 --- /dev/null +++ b/src/todo @@ -0,0 +1,2 @@ +* deplacer les ecritures pour gnuplot dans des classes type eoContinue (eoMonitor) +* integrer ACP