diff --git a/application/cma_sa/main.cpp b/application/cma_sa/main.cpp index 872636e7e..a452ead86 100644 --- a/application/cma_sa/main.cpp +++ b/application/cma_sa/main.cpp @@ -14,9 +14,13 @@ #include "Rosenbrock.h" #include "Sphere.h" -typedef eoReal EOT; +typedef eoReal EOT; -int main(int ac, char** av) +//typedef doUniform< EOT > Distrib; +//typedef doNormalMono< EOT > Distrib; +typedef doNormalMulti< EOT > Distrib; + +int main(int ac, char** av) { eoParserLogger parser(ac, av); @@ -37,19 +41,24 @@ int main(int ac, char** av) eoSelect< EOT >* selector = new eoDetSelect< EOT >(0.1); state.storeFunctor(selector); - //doEstimator< doUniform< EOT > >* estimator = new doEstimatorUniform< EOT >(); - doEstimator< doNormal< EOT > >* estimator = new doEstimatorNormal< EOT >(); + doEstimator< Distrib >* estimator = + //new doEstimatorUniform< EOT >(); + //new doEstimatorNormalMono< EOT >(); + new doEstimatorNormalMulti< EOT >(); state.storeFunctor(estimator); eoSelectOne< EOT >* selectone = new eoDetTournamentSelect< EOT >(); state.storeFunctor(selectone); - //doModifierMass< doUniform< EOT > >* modifier = new doUniformCenter< EOT >(); - doModifierMass< doNormal< EOT > >* modifier = new doNormalCenter< EOT >(); + doModifierMass< Distrib >* modifier = + //new doUniformCenter< EOT >(); + //new doNormalMonoCenter< EOT >(); + new doNormalMultiCenter< EOT >(); state.storeFunctor(modifier); - //eoEvalFunc< EOT >* plainEval = new BopoRosenbrock< EOT, double, const EOT& >(); - eoEvalFunc< EOT >* plainEval = new Sphere< EOT >(); + eoEvalFunc< EOT >* plainEval = + //new BopoRosenbrock< EOT, typename EOT::AtomType, const EOT& >(); + new Sphere< EOT >(); state.storeFunctor(plainEval); unsigned long max_eval = parser.getORcreateParam((unsigned long)0, "maxEval", "Maximum number of evaluations (0 = none)", 'E', "Stopping criterion").value(); // E @@ -66,6 +75,12 @@ int main(int ac, char** av) state.storeFunctor(init); + doStats< Distrib >* stats = + //new doStatsUniform< EOT >(); + //new doStatsNormalMono< EOT >(); + new doStatsNormalMulti< EOT >(); + state.storeFunctor(stats); + //----------------------------------------------------------------------------- @@ -98,8 +113,10 @@ int main(int ac, char** av) *gen); state.storeFunctor(bounder); - //doSampler< doUniform< EOT > >* sampler = new doSamplerUniform< EOT >(); - doSampler< doNormal< EOT > >* sampler = new doSamplerNormal< EOT >( *bounder ); + doSampler< Distrib >* sampler = + //new doSamplerUniform< EOT >(); + //new doSamplerNormalMono< EOT >( *bounder ); + new doSamplerNormalMulti< EOT >( *bounder ); state.storeFunctor(sampler); @@ -165,8 +182,7 @@ int main(int ac, char** av) // CMASA algorithm configuration //----------------------------------------------------------------------------- - //doAlgo< doUniform< EOT > >* algo = new doCMASA< doUniform< EOT > > - doAlgo< doNormal< EOT > >* algo = new doCMASA< doNormal< EOT > > + doAlgo< Distrib >* algo = new doCMASA< Distrib > (*selector, *estimator, *selectone, *modifier, *sampler, checkpoint, eval, *continuator, *cooling_schedule, initial_temperature, *replacor); diff --git a/src/do b/src/do index 478a3f216..10248481a 100644 --- a/src/do +++ b/src/do @@ -6,24 +6,27 @@ #include "doDistrib.h" #include "doUniform.h" -#include "doNormal.h" +#include "doNormalMono.h" +#include "doNormalMulti.h" #include "doEstimator.h" #include "doEstimatorUniform.h" -#include "doEstimatorNormal.h" +#include "doEstimatorNormalMono.h" +#include "doEstimatorNormalMulti.h" #include "doModifier.h" #include "doModifierDispersion.h" #include "doModifierMass.h" #include "doUniformCenter.h" -#include "doNormalCenter.h" +#include "doNormalMonoCenter.h" +#include "doNormalMultiCenter.h" #include "doSampler.h" #include "doSamplerUniform.h" -#include "doSamplerNormal.h" +#include "doSamplerNormalMono.h" +#include "doSamplerNormalMulti.h" #include "doVectorBounds.h" -#include "doNormalParams.h" #include "doBounder.h" #include "doBounderNo.h" @@ -31,5 +34,8 @@ #include "doBounderRng.h" #include "doStats.h" +#include "doStatsUniform.h" +#include "doStatsNormalMono.h" +#include "doStatsNormalMulti.h" #endif // !_do_ diff --git a/src/do.cpp b/src/do.cpp new file mode 100644 index 000000000..b33d98fe5 --- /dev/null +++ b/src/do.cpp @@ -0,0 +1 @@ +#include "do" diff --git a/src/doCMASA.h b/src/doCMASA.h index 31d1411a7..ec7de9d8b 100644 --- a/src/doCMASA.h +++ b/src/doCMASA.h @@ -33,6 +33,7 @@ #include "doEstimator.h" #include "doModifierMass.h" #include "doSampler.h" +#include "doHyperVolume.h" #include "doStats.h" using namespace boost::numeric::ublas; @@ -41,9 +42,12 @@ template < typename D > class doCMASA : public doAlgo< D > { public: - //! Alias for the type + //! Alias for the type EOT typedef typename D::EOType EOT; + //! Alias for the atom type + typedef typename EOT::AtomType AtomType; + //! Alias for the fitness typedef typename EOT::Fitness Fitness; @@ -74,7 +78,8 @@ public: moSolContinue < EOT > & continuator, moCoolingSchedule & cooling_schedule, double initial_temperature, - eoReplacement< EOT > & replacor + eoReplacement< EOT > & replacor, + doStats< D > & stats ) : _selector(selector), _estimator(estimator), @@ -87,6 +92,7 @@ public: _cooling_schedule(cooling_schedule), _initial_temperature(initial_temperature), _replacor(replacor), + _stats(stats), _pop_results_destination("ResPop"), @@ -164,14 +170,14 @@ public: double size = distrib.size(); assert(size > 0); - doHyperVolume hv; + doHyperVolume< EOT > hv; for (int i = 0; i < size; ++i) { //hv.update( distrib.varcovar()[i] ); } - _ofs_params_var << hv << std::endl; + // _ofs_params_var << hv << std::endl; } do @@ -332,14 +338,14 @@ public: double size = distrib.size(); assert(size > 0); - doHyperVolume hv; + doHyperVolume< EOT > hv; for (int i = 0; i < size; ++i) { //hv.update( distrib.varcovar()[i] ); } - _ofs_params_var << hv << std::endl; + //_ofs_params_var << hv << std::endl; } //------------------------------------------------------------- @@ -386,6 +392,8 @@ private: //! A EOT replacor eoReplacement < EOT > & _replacor; + //! Stats to print distrib parameters out + doStats< D > & _stats; //------------------------------------------------------------- // Temporary solution to store populations state at each diff --git a/src/doCheckPoint.h b/src/doCheckPoint.h new file mode 100644 index 000000000..e569e08d6 --- /dev/null +++ b/src/doCheckPoint.h @@ -0,0 +1,78 @@ +#ifndef _doCheckPoint_h +#define _doCheckPoint_h + +#include "doContinue.h" +#include "doStat.h" + +//! eoCheckPoint< EOT > classe fitted to Distribution Object library + +template < typename D > +class doCheckPoint : public doContinue< D > +{ +public: + typedef typename D::EOType EOType; + + doCheckPoint(doContinue< D >& _cont) + { + _continuators.push_back( &_cont ); + } + + bool operator()(const D& distrib) + { + for ( unsigned int i = 0, size = _stats.size(); i < size; ++i ) + { + (*_stats[i])( distrib ); + } + + bool bContinue = true; + for ( unsigned int i = 0, size = _continuators.size(); i < size; ++i ) + { + if ( !(*_continuators[i]( distrib )) ) + { + bContinue = false; + } + } + + if ( !bContinue ) + { + for ( unsigned int i = 0, size = _stats.size(); i < size; ++i ) + { + _stats[i]->lastCall( distrib ); + } + } + + return bContinue; + } + + void add(doContinue< D >& cont) { _continuators.push_back( &cont ); } + void add(doStatBase< D >& stat) { _stats.push_back( &stat ); } + + virtual std::string className(void) const { return "doCheckPoint"; } + + std::string allClassNames() const + { + std::string s("\n" + className() + "\n"); + + s += "Stats\n"; + for ( unsigned int i = 0, size = _stats.size(); i < size; ++i ) + { + s += _stats[i]->className() + "\n"; + } + s += "\n"; + + s += "Continuators\n"; + for ( unsigned int i = 0, size = _continuators.size(); i < size; ++i ) + { + s += _continuators[i]->className() + "\n"; + } + s += "\n"; + + return s; + } + +private: + std::vector< doContinue< D >* > _continuators; + std::vector< doStatBase< D >* > _stats; +}; + +#endif // !_doCheckPoint_h diff --git a/src/doContinue.h b/src/doContinue.h new file mode 100644 index 000000000..f7a7448d5 --- /dev/null +++ b/src/doContinue.h @@ -0,0 +1,27 @@ +#ifndef _doContinue_h +#define _doContinue_h + +#include +#include +#include + +//! eoContinue< EOT > classe fitted to Distribution Object library + +template < typename D > +class doContinue : public eoUF< const D&, bool >, public eoPersistent +{ +public: + virtual std::string className(void) const { return "doContinue"; } + + void readFrom(std::istream&) + { + /* It should be implemented by subclasses ! */ + } + + void printOn(std::ostream&) const + { + /* It should be implemented by subclasses ! */ + } +}; + +#endif // !_doContinue_h diff --git a/src/doEstimatorNormal.h b/src/doEstimatorNormal.h deleted file mode 100644 index 010008a85..000000000 --- a/src/doEstimatorNormal.h +++ /dev/null @@ -1,31 +0,0 @@ -#ifndef _doEstimatorNormal_h -#define _doEstimatorNormal_h - -#include "doEstimator.h" -#include "doUniform.h" -#include "doStats.h" - -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(); - assert(popsize > 0); - - unsigned int dimsize = pop[0].size(); - assert(dimsize > 0); - - //doCovMatrix cov(dimsize); - doUblasCovMatrix< EOT > cov; - - cov.update(pop); - - return doNormal< EOT >( cov.get_mean(), cov.get_varcovar() ); - } -}; - -#endif // !_doEstimatorNormal_h diff --git a/src/doEstimatorNormalMono.h b/src/doEstimatorNormalMono.h new file mode 100644 index 000000000..c6e168b77 --- /dev/null +++ b/src/doEstimatorNormalMono.h @@ -0,0 +1,70 @@ +#ifndef _doEstimatorNormalMono_h +#define _doEstimatorNormalMono_h + +#include "doEstimator.h" +#include "doNormalMono.h" + +template < typename EOT > +class doEstimatorNormalMono : public doEstimator< doNormalMono< EOT > > +{ +public: + typedef typename EOT::AtomType AtomType; + + class Variance + { + public: + Variance() : _sumvar(0){} + + void update(AtomType v) + { + _n++; + + AtomType d = v - _mean; + + _mean += 1 / _n * d; + _sumvar += (_n - 1) / _n * d * d; + } + + AtomType get_mean() const {return _mean;} + AtomType get_var() const {return _sumvar / (_n - 1);} + AtomType get_std() const {return sqrt( get_var() );} + + private: + AtomType _n; + AtomType _mean; + AtomType _sumvar; + }; + +public: + doNormalMono< EOT > operator()(eoPop& pop) + { + unsigned int popsize = pop.size(); + assert(popsize > 0); + + unsigned int dimsize = pop[0].size(); + assert(dimsize > 0); + + std::vector< Variance > var( dimsize ); + + for (unsigned int i = 0; i < popsize; ++i) + { + for (unsigned int d = 0; d < dimsize; ++d) + { + var[d].update( pop[i][d] ); + } + } + + EOT mean( dimsize ); + EOT variance( dimsize ); + + for (unsigned int d = 0; d < dimsize; ++d) + { + mean[d] = var[d].get_mean(); + variance[d] = var[d].get_var(); + } + + return doNormalMono< EOT >( mean, variance ); + } +}; + +#endif // !_doEstimatorNormalMono_h diff --git a/src/doEstimatorNormalMulti.h b/src/doEstimatorNormalMulti.h new file mode 100644 index 000000000..404a71d2c --- /dev/null +++ b/src/doEstimatorNormalMulti.h @@ -0,0 +1,109 @@ +#ifndef _doEstimatorNormalMulti_h +#define _doEstimatorNormalMulti_h + +#include "doEstimator.h" +#include "doUniform.h" +#include "doStats.h" + +template < typename EOT > +class doEstimatorNormalMulti : public doEstimator< doNormalMulti< EOT > > +{ +public: + class CovMatrix + { + public: + typedef typename EOT::AtomType AtomType; + + 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; + }; + +public: + typedef typename EOT::AtomType AtomType; + + doNormalMulti< EOT > operator()(eoPop& pop) + { + unsigned int popsize = pop.size(); + assert(popsize > 0); + + unsigned int dimsize = pop[0].size(); + assert(dimsize > 0); + + CovMatrix cov; + + cov.update(pop); + + return doNormalMulti< EOT >( cov.get_mean(), cov.get_varcovar() ); + } +}; + +#endif // !_doEstimatorNormalMulti_h diff --git a/src/doHyperVolume.h b/src/doHyperVolume.h new file mode 100644 index 000000000..c47f457b5 --- /dev/null +++ b/src/doHyperVolume.h @@ -0,0 +1,25 @@ +#ifndef _doHyperVolume_h +#define _doHyperVolume_h + +template < typename EOT > +class doHyperVolume +{ +public: + typedef typename EOT::AtomType AtomType; + + doHyperVolume() : _hv(1) {} + + void update(AtomType v) + { + _hv *= ::sqrt( v ); + + assert( _hv <= std::numeric_limits< AtomType >::max() ); + } + + AtomType get_hypervolume() const { return _hv; } + +protected: + AtomType _hv; +}; + +#endif // !_doHyperVolume_h diff --git a/src/doNormalCenter.h b/src/doNormalCenter.h deleted file mode 100644 index ff379ba22..000000000 --- a/src/doNormalCenter.h +++ /dev/null @@ -1,21 +0,0 @@ -#ifndef _doNormalCenter_h -#define _doNormalCenter_h - -#include "doModifierMass.h" -#include "doNormal.h" - -template < typename EOT > -class doNormalCenter : public doModifierMass< doNormal< EOT > > -{ -public: - typedef typename EOT::AtomType AtomType; - - void operator() ( doNormal< EOT >& distrib, EOT& mass ) - { - ublas::vector< AtomType > mean( distrib.size() ); - std::copy( mass.begin(), mass.end(), mean.begin() ); - distrib.mean() = mean; - } -}; - -#endif // !_doNormalCenter_h diff --git a/src/doNormalMonoCenter.h b/src/doNormalMonoCenter.h new file mode 100644 index 000000000..48b448a42 --- /dev/null +++ b/src/doNormalMonoCenter.h @@ -0,0 +1,19 @@ +#ifndef _doNormalMonoCenter_h +#define _doNormalMonoCenter_h + +#include "doModifierMass.h" +#include "doNormalMono.h" + +template < typename EOT > +class doNormalMonoCenter : public doModifierMass< doNormalMono< EOT > > +{ +public: + typedef typename EOT::AtomType AtomType; + + void operator() ( doNormalMono< EOT >& distrib, EOT& mass ) + { + distrib.mean() = mass; + } +}; + +#endif // !_doNormalMonoCenter_h diff --git a/src/doNormalMultiCenter.h b/src/doNormalMultiCenter.h new file mode 100644 index 000000000..65e9ec7cc --- /dev/null +++ b/src/doNormalMultiCenter.h @@ -0,0 +1,21 @@ +#ifndef _doNormalMultiCenter_h +#define _doNormalMultiCenter_h + +#include "doModifierMass.h" +#include "doNormalMulti.h" + +template < typename EOT > +class doNormalMultiCenter : public doModifierMass< doNormalMulti< EOT > > +{ +public: + typedef typename EOT::AtomType AtomType; + + void operator() ( doNormalMulti< EOT >& distrib, EOT& mass ) + { + ublas::vector< AtomType > mean( distrib.size() ); + std::copy( mass.begin(), mass.end(), mean.begin() ); + distrib.mean() = mean; + } +}; + +#endif // !_doNormalMultiCenter_h diff --git a/src/doSamplerNormalMono.h b/src/doSamplerNormalMono.h new file mode 100644 index 000000000..fa45b2e8b --- /dev/null +++ b/src/doSamplerNormalMono.h @@ -0,0 +1,65 @@ +#ifndef _doSamplerNormalMono_h +#define _doSamplerNormalMono_h + +#include + +#include "doSampler.h" +#include "doNormalMono.h" +#include "doBounder.h" +#include "doStats.h" + +/** + * doSamplerNormalMono + * This class uses the NormalMono distribution parameters (bounds) to return + * a random position used for population sampling. + */ +template < typename EOT > +class doSamplerNormalMono : public doSampler< doNormalMono< EOT > > +{ +public: + typedef typename EOT::AtomType AtomType; + + doSamplerNormalMono( doBounder< EOT > & bounder ) + : doSampler< doNormalMono< EOT > >( bounder ) + {} + + EOT sample( doNormalMono< EOT >& distrib ) + { + unsigned int size = distrib.size(); + assert(size > 0); + + + //------------------------------------------------------------- + // Point we want to sample to get higher a set of points + // (coordinates in n dimension) + // x = {x1, x2, ..., xn} + //------------------------------------------------------------- + + EOT solution; + + //------------------------------------------------------------- + + + //------------------------------------------------------------- + // Sampling all dimensions + //------------------------------------------------------------- + + for (unsigned int i = 0; i < size; ++i) + { + AtomType mean = distrib.mean()[i]; + AtomType variance = distrib.variance()[i]; + AtomType random = rng.normal(mean, variance); + + assert(variance >= 0); + + solution.push_back(random); + } + + //------------------------------------------------------------- + + + return solution; + } +}; + +#endif // !_doSamplerNormalMono_h diff --git a/src/doSamplerNormal.h b/src/doSamplerNormalMulti.h similarity index 54% rename from src/doSamplerNormal.h rename to src/doSamplerNormalMulti.h index 413f7ada6..ae54945b2 100644 --- a/src/doSamplerNormal.h +++ b/src/doSamplerNormalMulti.h @@ -1,34 +1,106 @@ -#ifndef _doSamplerNormal_h -#define _doSamplerNormal_h +#ifndef _doSamplerNormalMulti_h +#define _doSamplerNormalMulti_h + +#include +#include +#include #include #include "doSampler.h" -#include "doNormal.h" +#include "doNormalMulti.h" #include "doBounder.h" #include "doStats.h" /** - * doSamplerNormal + * doSamplerNormalMulti * This class uses the Normal distribution parameters (bounds) to return * a random position used for population sampling. */ template < typename EOT > -class doSamplerNormal : public doSampler< doNormal< EOT > > +class doSamplerNormalMulti : public doSampler< doNormalMulti< EOT > > { public: typedef typename EOT::AtomType AtomType; - doSamplerNormal( doBounder< EOT > & bounder ) - : doSampler< doNormal< EOT > >( bounder ) + class Cholesky + { + public: + 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; + }; + + doSamplerNormalMulti( doBounder< EOT > & bounder ) + : doSampler< doNormalMulti< EOT > >( bounder ) {} - EOT sample( doNormal< EOT >& distrib ) + EOT sample( doNormalMulti< EOT >& distrib ) { unsigned int size = distrib.size(); assert(size > 0); + //------------------------------------------------------------- // Cholesky factorisation gererating matrix L from covariance // matrix V. @@ -37,7 +109,7 @@ public: // L = cholesky decomposition of varcovar //------------------------------------------------------------- - Cholesky< EOT > cholesky; + Cholesky cholesky; cholesky.update( distrib.varcovar() ); ublas::symmetric_matrix< AtomType, ublas::lower > L = cholesky.get_L(); @@ -59,25 +131,21 @@ public: //------------------------------------------------------------- - // 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 ); @@ -89,4 +157,4 @@ public: } }; -#endif // !_doSamplerNormal_h +#endif // !_doSamplerNormalMulti_h diff --git a/src/doStat.h b/src/doStat.h new file mode 100644 index 000000000..46b213ba0 --- /dev/null +++ b/src/doStat.h @@ -0,0 +1,33 @@ +#ifndef _doStat_h +#define _doStat_h + +#include + +template < typename D > +class doStatBase : public eoUF< const D&, void > +{ +public: + // virtual void operator()( const D& ) = 0 (provided by eoUF< A1, R >) + + virtual void lastCall( const D& ) {} + virtual std::string className() const { return "doStatBase"; } +}; + +template < typename D > +class doStat : doStatBase< D > +{ +public: + typedef typename D::EOType EOType; + typedef typename EOType::AtomType AtomType; + +public: + doStat(){} + + + D& distrib() { return _distrib; } + +private: + D& _distrib; +}; + +#endif // !_doStat_h diff --git a/src/doStatNormalMono.h b/src/doStatNormalMono.h new file mode 100644 index 000000000..bc0a5b606 --- /dev/null +++ b/src/doStatNormalMono.h @@ -0,0 +1,21 @@ +#ifndef _doStatNormalMono_h +#define _doStatNormalMono_h + +#include "doStat.h" +#include "doNormalMono.h" + +template < typename EOT > +class doStatNormalMono : public doStat< doNormalMono< EOT > > +{ +public: + doStatNormalMono( doNormalMono< EOT >& distrib ) + : doStat< doNormalMono< EOT > >( distrib ) + {} + + virtual void printOn(std::ostream& os) const + { + os << this->distrib().mean() << " " << this->distrib().variance(); + } +}; + +#endif // !_doStatNormalMono_h diff --git a/src/doStatNormalMulti.h b/src/doStatNormalMulti.h new file mode 100644 index 000000000..7d57ee5a7 --- /dev/null +++ b/src/doStatNormalMulti.h @@ -0,0 +1,24 @@ +#ifndef _doStatNormalMulti_h +#define _doStatNormalMulti_h + +#include + +#include "doStat.h" +#include "doNormalMulti.h" +#include "doDistrib.h" + +template < typename EOT > +class doStatNormalMulti : public doStat< doNormalMulti< EOT > > +{ +public: + doStatNormalMulti( doNormalMulti< EOT >& distrib ) + : doStat< doNormalMulti< EOT > >( distrib ) + {} + + virtual void printOn(std::ostream& os) const + { + os << this->distrib().mean() << this->distrib().varcovar(); + } +}; + +#endif // !_doStatNormalMulti_h diff --git a/src/doStatUniform.h b/src/doStatUniform.h new file mode 100644 index 000000000..abeeb11d7 --- /dev/null +++ b/src/doStatUniform.h @@ -0,0 +1,21 @@ +#ifndef _doStatUniform_h +#define _doStatUniform_h + +#include "doStat.h" +#include "doUniform.h" + +template < typename EOT > +class doStatUniform : public doStat< doUniform< EOT > > +{ +public: + doStatUniform( doUniform< EOT >& distrib ) + : doStat< doUniform< EOT > >( distrib ) + {} + + virtual void printOn(std::ostream& os) const + { + os << this->distrib().min() << " " << this->distrib().max(); + } +}; + +#endif // !_doStatUniform_h diff --git a/src/doStats.cpp b/src/doStats.cpp deleted file mode 100644 index 80c9869fd..000000000 --- a/src/doStats.cpp +++ /dev/null @@ -1,192 +0,0 @@ -/* - * Copyright (C) 2005 Maarten Keijzer - * - * This program is free software; you can redistribute it and/or modify - * it under the terms of version 2 of the GNU General Public License as - * published by the Free Software Foundation. - * - * 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 General Public License for more details. - * - * You should have received a copy of the GNU 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 - */ - -#include -#include - -#include -#include - -#include "doStats.h" - -doStats::doStats() - : _n(0) -{} - -void doStats::printOn(std::ostream& _os) const -{ - _os << "Not implemented yet! "; -} - -doMean::doMean() - : _mean(0) -{} - -void doMean::update(double v) -{ - _n++; - - double d = v - _mean; - - _mean += 1 / _n * d; -} - -double doMean::get_mean() const -{ - return _mean; -} - -void doMean::printOn(std::ostream& _os) const -{ - _os << get_mean(); -} - -doVar::doVar() - : _sumvar(0) -{} - -void doVar::update(double v) -{ - _n++; - - double d = v - _mean; - - _mean += 1 / _n * d; - _sumvar += (_n - 1) / _n * d * d; -} - -double doVar::get_var() const -{ - return _sumvar / (_n - 1); -} - -double doVar::get_std() const -{ - return ::sqrt( get_var() ); -} - -void doVar::printOn(std::ostream& _os) const -{ - _os << get_var(); -} - -doCov::doCov() - : _meana(0), _meanb(0), _sumcov(0) -{} - -void doCov::update(double a, double b) -{ - ++_n; - - double da = a - _meana; - double db = b - _meanb; - - _meana += 1 / _n * da; - _meanb += 1 / _n * db; - - _sumcov += (_n - 1) / _n * da * db; -} - -double doCov::get_meana() const -{ - return _meana; -} - -double doCov::get_meanb() const -{ - return _meanb; -} - -double doCov::get_cov() const -{ - return _sumcov / (_n - 1); -} - -void doCov::printOn(std::ostream& _os) const -{ - _os << get_cov(); -} - -doCovMatrix::doCovMatrix(unsigned dim) - : _mean(dim), _sumcov(dim, std::vector< double >( dim )) -{} - -void doCovMatrix::update(const std::vector& v) -{ - assert(v.size() == _mean.size()); - - _n++; - - for (unsigned int i = 0; i < v.size(); ++i) - { - double d = v[i] - _mean[i]; - - _mean[i] += 1 / _n * d; - _sumcov[i][i] += (_n - 1) / _n * d * d; - - for (unsigned j = i; j < v.size(); ++j) - { - double e = v[j] - _mean[j]; // _mean[j] is not updated yet - - double upd = (_n - 1) / _n * d * e; - - _sumcov[i][j] += upd; - _sumcov[j][i] += upd; - } - } -} - -double doCovMatrix::get_mean(int i) const -{ - return _mean[i]; -} - -double doCovMatrix::get_var(int i) const -{ - return _sumcov[i][i] / (_n - 1); -} - -double doCovMatrix::get_std(int i) const -{ - return ::sqrt( get_var(i) ); -} - -double doCovMatrix::get_cov(int i, int j) const -{ - return _sumcov[i][j] / (_n - 1); -} - -doHyperVolume::doHyperVolume() - : _hv(1) -{} - -void doHyperVolume::update(double v) -{ - _hv *= ::sqrt(v); - - assert( _hv <= std::numeric_limits< double >::max() ); -} - -double doHyperVolume::get_hypervolume() const -{ - return _hv; -} - -void doHyperVolume::printOn(std::ostream& _os) const -{ - _os << get_hypervolume(); -} diff --git a/src/doStats.h b/src/doStats.h deleted file mode 100644 index 50cc9c95e..000000000 --- a/src/doStats.h +++ /dev/null @@ -1,275 +0,0 @@ -/* - * Copyright (C) 2005 Maarten Keijzer - * - * This program is free software; you can redistribute it and/or modify - * it under the terms of version 2 of the GNU General Public License as - * published by the Free Software Foundation. - * - * 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 General Public License for more details. - * - * You should have received a copy of the GNU 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 _doStats_h -#define _doStats_h - -#include -#include -#include - -#include -#include - -namespace ublas = boost::numeric::ublas; - -class doStats : public eoPrintable -{ -public: - doStats(); - - virtual void printOn(std::ostream&) const; - -protected: - double _n; -}; - -class doMean : public doStats -{ -public: - doMean(); - - virtual void update(double); - virtual void printOn(std::ostream&) const; - - double get_mean() const; - -protected: - double _mean; -}; - -class doVar : public doMean -{ -public: - doVar(); - - virtual void update(double); - virtual void printOn(std::ostream&) const; - - double get_var() const; - double get_std() const; - -protected: - double _sumvar; -}; - -/** Single covariance between two variates */ -class doCov : public doStats -{ -public: - doCov(); - - virtual void update(double, double); - virtual void printOn(std::ostream&) const; - - double get_meana() const; - double get_meanb() const; - double get_cov() const; - -protected: - double _meana; - double _meanb; - double _sumcov; -}; - -class doCovMatrix : public doStats -{ -public: - doCovMatrix(unsigned dim); - - virtual void update(const std::vector&); - - double get_mean(int) const; - double get_var(int) const; - double get_std(int) const; - double get_cov(int, int) const; - -protected: - std::vector< double > _mean; - std::vector< std::vector< double > > _sumcov; -}; - -class doHyperVolume : public doStats -{ -public: - doHyperVolume(); - - virtual void update(double); - virtual void printOn(std::ostream&) const; - - double get_hypervolume() const; - -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