MendelMin.hh

#ifndef RIVET_MendelMin_H
#define RIVET_MendelMin_H
 
#include "Rivet/Tools/Random.hh"
#include <cmath>
#include <functional>
#include <iomanip>
#include <iostream>
#include <map>
#include <random>
#include <valarray>
#include <vector>
 
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