Decimate vector

/*
 * Copyright 2013 Florian Philipp
 *
 * Licensed under the Apache License, Version 2.0 (the "License");
 * you may not use this file except in compliance with the License.
 * You may obtain a copy of the License at
 *
 * http://www.apache.org/licenses/LICENSE-2.0

 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
#include <cstdlib>
// using std::rand, RAND_MAX, EXIT_SUCCESS
#include <vector>
// using std::vector
#include <queue>
// using std::priority_queue
#include <array>
// using std::array
#include <algorithm>
// using std::sort, std::for_each, std::make_heap, std::pop_heap, std::remove_if
#include <functional>
// using std::greater
#include <cstdio>
// using stdout, stderr, std::FILE, std::fprintf
#include <chrono>
// using std::chrono::high_resolution_clock, std::chrono::duration,
// using std::chrono::duration_cast, std::chrono::milliseconds
#include <stdexcept>
// using std::logic_error
#include <utility>
// using std::swap

namespace {
  typedef std::chrono::high_resolution_clock clock_type;

  /**
   * Mockup of a candidate solution in a genetic algorithm
   */
  class Individuum
  {
    double fitness;
    static double generate_fitness()
    {
      return static_cast<double>(std::rand())/RAND_MAX;
    }
  public:
    Individuum():
      fitness(0)
    {}
    void evolve()
    {
      fitness = generate_fitness();
    }
    bool operator<(const Individuum& o) const
    {
      return fitness < o.fitness;
    }
    bool operator<=(const Individuum& o)  const
    {
      return fitness <= o.fitness;
    }
    bool operator>(const Individuum& o) const
    {
      return o < *this;
    }
    bool operator>=(const Individuum& o) const
    {
      return o <= *this;
    }
    bool operator==(const Individuum& o) const
    {
      return fitness == o.fitness;
    }
    bool operator!=(const Individuum& o) const
    {
      return !(*this == o);
    }
    double get_fitness() const
    {
      return fitness;
    }
    void swap(Individuum& o)
    {
      using std::swap;
      swap(fitness, o.fitness);
    }
  };
} // namespace
namespace std {
  template<>
  void swap<Individuum>(Individuum& a, Individuum& b)
  {
    a.swap(b);
  }
}
namespace {
  /**
   * Function for use with std::for_each
   */
  void evolve(Individuum& i) { i.evolve(); }

  /**
   * A heap of fitness values
   *
   * The top() element has the highest fitness
   */
  typedef std::priority_queue<double> heap_t;
  /**
   * Fills the heap given in the constructor with the n least fit Individuums
   */
  struct FillHeap
  {
    typedef heap_t::size_type size_type;
    heap_t& least_fit_n;
    const size_type capacity;
    FillHeap(heap_t& least_fit_n, size_type n):
      least_fit_n(least_fit_n), capacity(n)
    {}
    void operator()(Individuum& i)
    {
      double fitness = i.get_fitness();
      if(least_fit_n.size() == capacity) {
    if(fitness < least_fit_n.top()) {
      least_fit_n.pop();
      least_fit_n.push(fitness);
    }
      }
      else
    least_fit_n.push(fitness);
    }
  };
  /**
   * Unary functor that returns true for Individuums that are less or equal fit
   * compared to the value given in the constructor
   */
  struct LessEqual
  {
    double comparison;
    LessEqual(double comparison): comparison(comparison) {}
    bool operator()(const Individuum& o) const
    {
      return o.get_fitness() <= comparison;
    }
  };
  /**
   * Similar to std::remove_if but does not preserve relative order of the
   * remaining elements
   *
   * \tparam BidirectionalIterator an iterator supporting incrementing and
   * decrementing. The type pointed to has to be move assignable
   * \tparam Predicate a unary predicate which returns true for all elements
   * that have to be removed
   * \param first beginning of the range that is to be evaluated
   * \param last iterator behind the end of the evaluated range
   * \param predicate the evaluation functor
   * \return an iterator pointing to the end of the range of remaining elements.
   * Elements behind this point and last are in valid but unordered and
   * undefined state. Elements between this point and first are in valid state
   * but unsorted
   */
  template<class BidirectionalIterator, class Predicate>
  BidirectionalIterator remove_unstable_if(BidirectionalIterator first,
                       BidirectionalIterator last,
                       Predicate predicate)
  {
    while(first != last) {
      if(predicate(*first)) {
    --last;
    *first = std::move(*last);
      }
      else
    ++first;
    }
    return last;
  }
  /**
   * Storage of all candidate solutions of a genetic algorithm
   *
   * Implements various decimation strategies for benchmarking purposes
   */
  class Population
  {
    typedef std::vector<Individuum> population_t;
  public:
    typedef population_t::size_type size_type;
  private:
    population_t population;
    size_type n_to_decimate(double percentage) const
    {
      return static_cast<size_type>(percentage * population.size());
    }
    size_type n_to_survive(double percentage_killoff) const
    {
      return population.size() - n_to_decimate(percentage_killoff);
    }
  public:
    void evolve_all()
    {
      std::for_each(population.begin(), population.end(), evolve);
    }
    void grow_to(size_type n)
    {
      population.resize(n);
    }
    /**
     * Removes the given percentage of Individuums by sorting the candidate
     * solutions by descending fitness and then shrinking the container
     *
     * Usually the slowest approach unless the population is already sorted
     */
    void decimate_with_sort(double percentage)
    {
      std::sort(population.begin(), population.end(),
        std::greater<Individuum>());
      population.resize(n_to_survive(percentage));
    }
    /**
     * Removes the given percentage of Individuums by reordering the candidate
     * solutions as a heap, then popping elements from the heap and shrinking
     * the container
     *
     * Usually the fastest approach
     */
    void decimate_with_heap(double percentage)
    {
      std::make_heap(population.begin(), population.end(),
             std::greater<Individuum>());
      const size_type n = n_to_survive(percentage);
      const population_t::iterator new_end = population.begin() + n;
      for(population_t::iterator cur_end = population.end(); cur_end != new_end;
      --cur_end)
    std::pop_heap(population.begin(), cur_end, std::greater<Individuum>());
      population.erase(new_end, population.end());
    }
    /**
     * Removes the given percentage of Individuums by finding the least fit
     * percentage and then removing them in a second pass
     *
     * Slightly slower than decimate_with_heap. Can be faster with very low
     * percentages
     */
    void decimate_with_twopass_heap(double percentage)
    {
      const size_type n = n_to_decimate(percentage);
      if(n == 0)
    return;
      heap_t heap;
      std::for_each(population.begin(), population.end(), FillHeap(heap, n));
      double cutoff = heap.top();
      population.erase(std::remove_if(population.begin(), population.end(),
                      LessEqual(cutoff)),
               population.end());
    }
    /**
     * Equivalent to decimate_with_twopass_heap but uses remove_unstable_if
     *
     * Slightly faster than decimate_with_twopass_heap but probably
     * insignificant
     */
    void decimate_with_twopass_unstable(double percentage)
    {
      const size_type n = n_to_decimate(percentage);
      if(n == 0)
    return;
      heap_t heap;
      std::for_each(population.begin(), population.end(), FillHeap(heap, n));
      double cutoff = heap.top();
      population.erase(remove_unstable_if(population.begin(), population.end(),
                      LessEqual(cutoff)),
               population.end());
    }
    void print(std::FILE* f) const
    {
      population_t::const_iterator i = population.begin();
      const population_t::const_iterator end = population.end();
      if(i != end) {
    std::fprintf(f, "%4f", i->get_fitness());
    for(++i; i != end; ++i)
      std::fprintf(f, " %4f", i->get_fitness());
      }
      std::fprintf(f, "\n");
    }
    typedef void (Population::* decimation_fun_t)(double);
    clock_type::duration benchmark(std::size_t iterations, double percentile,
                   decimation_fun_t decimation_fun)
    {
      clock_type::duration total_time;
      const size_type pop_count = population.size();
      for(std::size_t i = 0; i < iterations; ++i) {
    evolve_all();
    clock_type::time_point start_time = clock_type::now();
    (this->*decimation_fun)(percentile);
    total_time += clock_type::now() - start_time;
    grow_to(pop_count);
      }
      return total_time;
    }
    bool operator==(Population& o)
    {
      if(population.size() != o.population.size())
    return false;
      std::sort(population.begin(), population.end());
      std::sort(o.population.begin(), o.population.end());
      return std::equal(population.begin(), population.end(),
            o.population.begin());
    }
    bool operator!=(Population& o)
    {
      return ! (*this == o);
    }
    void test()
    {
      grow_to(1000);
      evolve_all();
      std::array<Population, 3> copies;
      copies.fill(*this);
      const double percentage = 0.2;
      this->decimate_with_sort(percentage);
      copies[0].decimate_with_heap(percentage);
      if(copies[0] != *this)
    throw std::logic_error("Decimate with heap not equivalent");
      copies[1].decimate_with_twopass_heap(percentage);
      if(copies[1] != *this)
    throw std::logic_error("Decimate with twopass heap not equivalent");
      copies[2].decimate_with_twopass_unstable(percentage);
      if(copies[2] != *this)
    throw std::logic_error("Decimate with twopass unstable not equivalent");
}
  };
  std::size_t to_ms(const clock_type::duration time)
  {
    return std::chrono::duration_cast<std::chrono::milliseconds>(time).count();
  }
  void run_benchmark(Population& population,
             Population::decimation_fun_t decimation_fun,
             const char* fun_name)
  {
    const std::size_t iterations = 1000;
    const Population::size_type pop_count = 100000;
    const double percentage = 0.1;
    population.grow_to(pop_count);
    clock_type::duration time = population.benchmark(iterations, percentage,
                             decimation_fun);
    std::fprintf(stdout, "%s: %zd ms\n", fun_name, to_ms(time));
  }
} // namespace

int main()
{
  Population population;
  population.test();
  run_benchmark(population, &Population::decimate_with_sort, "Sort");
  run_benchmark(population, &Population::decimate_with_heap, "Heap");
  run_benchmark(population, &Population::decimate_with_twopass_heap, "Twopass");
  run_benchmark(population, &Population::decimate_with_twopass_unstable,
        "Twopass unstable");
  return EXIT_SUCCESS;
}