MendelMin.hh

#ifndef RIVET_MendelMin_H
#define RIVET_MendelMin_H
 
#include "Rivet/Tools/Random.hh"
#include <valarray>
#include <random>
#include <functional>
#include <map>
#include <vector>
#include <cmath>
#include <iostream>
#include <iomanip>
 
namespace Rivet {
 
  using std::valarray;
 
 
  class MendelMin {
  public:
 
    using Params = std::valarray<double>;
    using FuncT = std::function<double(const Params&, const Params&)>;
    using FuncNoFixedT = std::function<double(const Params&)>;
    // /// Typedef for the [0,1] random number generator
    // using RndT = std::function<double()>;
 
 
    MendelMin(const FuncT& fin, unsigned int ndim,
              const Params& fixpar, //const RndT & rndin,
              unsigned int npop=20, //unsigned int ngen=20, //< unused
              double margin=0.1)
      : _f(fin), _q(fixpar), //_rnd(rndin),
        _NDim(ndim), _margin(margin),
        _pop(npop), _fit(npop, -1.0), showTrace(false) {}
 
 
    MendelMin(const FuncNoFixedT& fin, unsigned int ndim,
              //const RndT & rndin,
              unsigned int npop=20, unsigned int ngen=20,
              double margin=0.1)
      : MendelMin([&](const Params& ps, const Params&) -> double { return fin(ps); },
                  ndim, {}, npop, /* ngen, */ margin)
    {   }
 
 
    void guess(const Params & p) {
      _pop.push_back(p);
      limit01(_pop.back());
      _fit.push_back(f(_pop.back()));
    }
 
 
    double evolve(unsigned int nGen) {
      for ( unsigned n = 0; n < nGen; ++n ) {
        // Calculate the fitness.
        auto mm = minmax();
        // Always kill the fittest individual.
        if ( showTrace ) _debug();
        for ( unsigned int i = 1; i < _pop.size(); ++i ) {
          if ( _fit[i] > rnd()*(mm.second - mm.first) )
            // Kill all individuals that have low fitness or are just unlucky.
            _fit[i] = -1.0;
          else
            // Improve This individual to be more like the fittest.
            move(_pop[i],_pop[0]);
        }
      }
      return _fit[0];
    }
 
    Params fittest() const {
      return _pop[0];
    }
 
    double fit() const {
      return _fit[0];
    }
 
    double rnd() const {
      return rand01(); //_rnd();
    }
 
    Params rndParams() const {
      Params ret(_NDim);
      for ( unsigned int i = 0; i < _NDim; ++i ) ret[i] = rnd();
      return ret;
    }
 
    void limit01(Params & p) const {
      for ( unsigned int i = 0; i < _NDim; ++i )
        p[i] = std::max(0.0, std::min(p[i], 1.0));
    }
 
    void move(Params & bad, const Params & better) const {
      bad += (better - bad)*(rndParams()*(1.0 + 2.0*_margin) - _margin);
      limit01(bad);
    }
 
    double f(const Params & p) const {
      return _f(p, _q);
    }
 
 
    std::pair<double, double> minmax() {
      std::pair<double,double> mm(std::numeric_limits<double>::max(), 0.0);
      unsigned int iwin = 0;
      for ( unsigned int i = 0; i < _pop.size(); ++i ) {
        double & v = _fit[i];
        // negative fitness value means the individual is dead, so we
        // welocme a new immigrant.
        if ( v < 0.0 ) _pop[i] = rndParams();
 
        // The calculated fitness value cannot be negative.
        v = std::max(0.0, f(_pop[i]));
 
        // Compare to the best and worst fitness so far.
        if ( v < mm.first ) iwin = i;
        mm.first = std::min(v, mm.first);
        mm.second = std::max(mm.second, v);
      }
 
      // Move the winner to the top.
      if ( iwin != 0 ) {
        std::swap(_pop[0], _pop[iwin]);
        std::swap(_fit[0], _fit[iwin]);
      }
      return mm;
    }
 
    void _debug() {
      std::cout << "GenAlgMax population status:" << std::endl;
      for ( unsigned int i = 0; i < _pop.size(); ++i ) {
        std::cout << std::setw(10) << _fit[i] << " (" << _pop[i][0];
        for ( unsigned int ip = 1; ip < _NDim; ++ip )
          std::cout << "," << _pop[i][ip];
        std::cout << ")" << std::endl;
      }
    }
 
 
  private:
 
    const FuncT _f;
 
    Params _q;
 
    //  const double _q;
 
    //const RndT _rnd;
 
    unsigned int _NDim;
 
    double _margin;
 
    std::vector<Params> _pop;
 
 
    std::vector<double> _fit;
 
  public:
 
    bool showTrace;
 
  };
 
 
  // template <typename FuncT, typename RndT>
  // MendelMin <FuncT, RndT>
  // makeMendelMin(const FuncT & f, const RndT & rnd, unsigned int ndim,
  //               unsigned int npop = 20, double margin = 0.1) {
  //   return MendelMin<FuncT, RndT>(f, rnd, ndim, npop, margin);
  // }
 
}
 
#endif // RIVET_MendelMin_H