[ex09] Solutions

exercises/ex09/solutions/cache_effects/.gitignore

+cache
+out.dat

exercises/ex09/solutions/cache_effects/Makefile

+#Makefile
+
+CXX = g++
+CXXFLAGS = -std=c++11 -O0 -DNDEBUG
+CXXFLAGS += -Wall -Wextra -Wpedantic -fno-tree-vectorize
+
+.PHONY: all
+all: result/plot_a.pdf
+
+result/plot_a.pdf: plot.py result/out.dat result
+	python plot.py result/out.dat
+
+result/out.dat: cache result
+	./cache | tee result/out.dat
+
+result:
+	mkdir result
+
+cache: cache.cpp
+
+.PHONY: clean
+clean:
+	rm -rf result
+	rm -f cache

exercises/ex09/solutions/cache_effects/cache.cpp

+/* Programming Techniques for Scientific Simulations, ETH Zürich
+ */
+
+#include <vector>
+#include <iostream>
+#include <iomanip>
+#include <chrono>
+
+using array_t = int;
+namespace timer = std::chrono;
+
+int main() {
+  constexpr size_t elem_size = sizeof(array_t);
+  constexpr size_t MINSIZE = 2;
+  constexpr size_t MAXSIZE = 2*2*2 * 1024 * 1024;  // 2^23
+  constexpr size_t MAXSTEP = 2*2 *1024;            // 2^12
+  constexpr size_t NUMOPS  = 20 * MAXSIZE;
+
+  array_t* arr = new array_t[MAXSIZE];
+  for(size_t N = MINSIZE; N < MAXSIZE; N *= 2) {
+    for(size_t step = 1; step <= std::min(MAXSTEP, N); step *= 2) {
+
+      const size_t num_steps = N / step;
+      const size_t num_sweeps = NUMOPS / num_steps;
+
+      for(size_t i = 0; i < N; ++i)
+        arr[i] = 0;
+
+      timer::time_point<timer::high_resolution_clock> start , stop;
+      start = timer::high_resolution_clock::now();
+
+      for(size_t sweep = 0; sweep < num_sweeps; ++sweep)
+          for(size_t i = 0; i < N; i += step)
+            arr[i]+= 1;
+
+
+      stop = timer::high_resolution_clock::now();
+      const double time = static_cast<timer::duration<double> >(stop-start).count();
+      const double numops = num_steps * num_sweeps;
+
+      std::cout << std::setprecision(12)
+                << std::setw(16) << N * elem_size      // array size in bytes
+                << std::setw(16) << step * elem_size   // touched bytes
+                << std::setw(16) << time
+                << std::setw(16) << numops
+                << std::endl;
+    }
+  }
+
+  delete[] arr;
+
+  return 0;
+}

exercises/ex09/solutions/cache_effects/plot.py

+#!/usr/bin/env python
+"""
+Usage: ./plot.py data_file
+"""
+
+import sys
+import os
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+if len(sys.argv) < 2:
+    sys.exit('Please provide a filename!')
+
+RESULTS = np.sort(
+    np.loadtxt(
+        sys.argv[1],
+        dtype={
+            'names': ['size', 'step', 'time', 'nops'],
+            'formats': ['u4', 'u4', 'f8', 'u4']
+        }
+    ),
+    order=['step', 'size']
+)
+
+STEPS = np.unique(RESULTS['step'])
+STEP_SLICES = [RESULTS['step'] == s for s in STEPS]
+
+plt.figure()
+for s, s_slice in zip(STEPS, STEP_SLICES):
+    if s > 512:
+        break
+    plt.semilogx(
+        RESULTS['size'][s_slice],
+        1e9 * RESULTS['time'][s_slice] / RESULTS['nops'][s_slice],
+        marker='+',
+        label='{} bytes'.format(s),
+        basex=2
+    )
+plt.title('Part a) Cache and Cacheline Size.')
+plt.legend(title='Step Size')
+plt.ylabel('time [ns]')
+plt.xlabel('array size [bytes]')
+plt.savefig(os.path.join('result', 'plot_a.pdf'), bbox_inches='tight')
+
+plt.figure()
+for s, s_slice in zip(STEPS, STEP_SLICES):
+    if s < 1024:
+        continue
+    plt.semilogx(
+        RESULTS['size'][s_slice] / s,
+        1e9 * RESULTS['time'][s_slice] / RESULTS['nops'][s_slice],
+        label='{} bytes'.format(s),
+        basex=2
+    )
+plt.title('Part b) Associativity.')
+plt.legend(title='Step Size')
+plt.ylabel('time [ns]')
+plt.xlabel('array size / step size')
+plt.savefig(os.path.join('result', 'plot_b.pdf'), bbox_inches='tight')

exercises/ex09/solutions/penna/.gitignore

+result
+build
+*.dat
+*.png

exercises/ex09/solutions/penna/CMakeLists.txt

+# Require 3.1 for CMAKE_CXX_STANDARD property
+cmake_minimum_required(VERSION 3.1)
+
+project(PennaFishPopulation)
+
+# C++11
+set(CMAKE_CXX_STANDARD 11)
+
+# Make it a Release build to add -O3 and -DNDEBUG flags by default
+# Comment it out to revert to standard behaviour
+set(CMAKE_BUILD_TYPE Release)
+# Also add -march-native (and -funroll-loops if you switch lines)
+add_compile_options($<$<CONFIG:RELEASE>:-march=native>)
+# add_compile_options("$<$<CONFIG:RELEASE>:-march=native;-funroll-loops>")
+
+# Setting warning compiler flags
+if(CMAKE_CXX_COMPILER_ID MATCHES "(C|c?)lang")
+  add_compile_options(-Weverything)
+else()
+  add_compile_options(-Wall -Wextra -Wpedantic)
+endif()
+
+set(DIRECTORIES src test)
+
+foreach(directory ${DIRECTORIES})
+  add_subdirectory(${directory})
+endforeach(directory)
+
+# This needs to be in both the top-level and the test CMake file.
+enable_testing()
+
+include_directories(${CMAKE_SOURCE_DIR}/src)
+
+add_executable(penna_sim main.cpp)
+target_link_libraries(penna_sim penna)
+
+add_executable(penna_sim_vector main.cpp)
+set_target_properties(penna_sim_vector PROPERTIES COMPILE_FLAGS "-DPENNA_VECTOR")
+target_link_libraries(penna_sim_vector penna_vector)

exercises/ex09/solutions/penna/main.cpp

+#include "src/fish_population.hpp"
+#include "src/animal.hpp"
+
+#include <string>
+#include <random>
+#include <vector>
+#include <cstdlib>
+#include <fstream>
+#include <iostream>
+#include <iterator>
+
+std::vector<double> gene_histogram(const Penna::FishPopulation& pop) {
+    std::vector<double> histo(Penna::Genome::number_of_genes);
+    double n = pop.size();
+    for (const auto a : pop) {
+      for (Penna::age_t i = 0; i < histo.size(); ++i)
+        histo[i] += a.genome().is_bad(i) / n;
+    }
+    return histo;
+}
+
+
+int main(int argc, const char** argv) {
+    if (!(argc==12 || argc==15)) {
+      std::cout << "Usage: " << argv[0] << " [simulational time] [mutation rate M] "
+                << "[bad threshold T] [reproduction age R] [pregnancy probability P] "
+                << "[population limit Nmax] [initial population size N0] [num measurements] "
+                << "[time at start of fishing 1] [fishing rate 1] [fishing age 1] "
+                << "optional: [time at start of fishing 2] [fishing rate 2] [fishing age 2]"
+                << std::endl;
+      return -1;
+    }
+    std::size_t sim_time        = std::stoul(argv[1]),
+                mut_rate        = std::stoul(argv[2]),
+                bad_threshold   = std::stoul(argv[3]),
+                repr_age        = std::stoul(argv[4]),
+                nmax            = std::stoul(argv[6]),
+                n0              = std::stoul(argv[7]),
+                nmeas           = std::stoul(argv[8]),
+                fishing_1_start = std::stoul(argv[9]),
+                fishing_age_1   = std::stoul(argv[11]),
+                fishing_2_start = argc > 12 ? std::stoul(argv[12]) : sim_time,
+                fishing_age_2   = argc > 12 ? std::stoul(argv[14]) : 0;
+
+    double pregnancy_prob = std::stod(argv[5]),
+           fishing_rate_1 = std::stod(argv[10]),
+           fishing_rate_2 = argc > 12 ? std::stod(argv[13]) : 0.;
+
+    srand(42);  // Setting the seed for our random number generator
+    Penna::Genome::set_mutation_rate               (mut_rate);
+    Penna::Animal::set_bad_threshold               (bad_threshold);
+    Penna::Animal::set_maturity_age                (repr_age);
+    Penna::Animal::set_probability_to_get_pregnant (pregnancy_prob);
+
+    std::cout << "## sim_time=" << sim_time << ", M=" << mut_rate << ", T=" << bad_threshold
+              << ", R=" << repr_age << ", pregnancy_prob=" << pregnancy_prob << ", Nmax=" << nmax
+              << ", N0=" << n0 << std::endl;
+    std::cout << "## Penna::Fishing period 1: starts=" << fishing_1_start
+              << ", rate=" << fishing_rate_1 << ", age>=" << fishing_age_1 << std::endl;
+    if (argc > 12)
+      std::cout << "## Penna::Fishing period 2: starts=" << fishing_2_start
+                << ", rate=" << fishing_rate_2 << ", age>=" << fishing_age_2 << std::endl;
+
+    Penna::FishPopulation pop(nmax, n0);
+    std::cout << "# time\tN" << std::endl;
+    std::size_t time = 0;
+
+    // Generating a stable population
+    for(; time < fishing_1_start; ++time) {
+      pop.simulate(1);
+      if (time % (sim_time/nmeas) == 0)
+        std::cout << time << "\t" << pop.size() << std::endl;
+    }
+
+    // Fishing period 1
+    pop.change_fishing(fishing_rate_1,fishing_age_1);
+    for(; time < fishing_2_start; ++time) {
+      pop.simulate(1);
+      if (time % (sim_time/nmeas) == 0)
+        std::cout << time << "\t" << pop.size() << std::endl;
+    }
+
+    // Fishing period 2
+    pop.change_fishing(fishing_rate_2,fishing_age_2);
+    for(; time < sim_time; ++time) {
+      pop.simulate(1);
+      if (time % (sim_time/nmeas) == 0)
+        std::cout << time << "\t" << pop.size() << std::endl;
+    }
+
+    std::vector<double> genes = gene_histogram(pop);
+    std::ofstream genefile("result/gene_histogram.dat");
+    std::copy(genes.begin(),genes.end(),std::ostream_iterator<double>(genefile,"\n"));
+
+  return 0;
+}

exercises/ex09/solutions/penna/plot.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+from numpy import *
+import matplotlib
+matplotlib.use("agg")
+from matplotlib.pyplot import *
+import sys
+import os
+
+for pop_name in ['population_list', 'population_vec']:
+    pop = loadtxt(os.path.join('result', '{}.dat'.format(pop_name)))
+    figure()
+    plot(pop[:,0],pop[:,1],'x')
+    xscale('linear')
+    yscale('log')
+    xlabel('year')
+    ylabel('population size')
+    savefig(os.path.join('result', '{}.png'.format(pop_name)))
+
+figure()
+gen = loadtxt(os.path.join('result', 'gene_histogram.dat'))
+plot(gen)
+ylim(ymin=0)
+xlabel('genome position')
+ylabel('bad fraction')
+savefig(os.path.join('result', 'gene_histogram.png'))

exercises/ex09/solutions/penna/run.sh

+#!/bin/bash
+
+TIME=1600
+M=1
+T=3
+R=7
+P=1.0
+NMAX=200000
+N0=$((NMAX/10))
+NMEAS=$TIME
+M1=1000
+M2=1200
+r1=0.17
+a1=0
+r2=0.22
+a2=0
+
+DATFILE_LIST=population_list.dat
+DATFILE_VEC=population_vec.dat
+
+cmake -Bbuild -H./
+make -C build
+mkdir -p result
+echo "Timing with TIME=$TIME M=$M T=$T R=$R P=$P NMAX=$NMAX, N0=$N0, M1=$M1, r1=$r1, a1=$a1, M2=$M2, r2=$r2, a2=$a2 ..."
+echo "List version:"
+time ./build/penna_sim $TIME $M $T $R $P $NMAX $N0 $NMEAS $M1 $r1 $a1 $M2 $r2 $a2 > result/$DATFILE_LIST
+echo ""
+echo "Penna vector version:"
+time ./build/penna_sim_vector $TIME $M $T $R $P $NMAX $N0 $NMEAS $M1 $r1 $a1 $M2 $r2 $a2 > result/$DATFILE_VEC
+./plot.py

exercises/ex09/solutions/penna/src/CMakeLists.txt

+include_directories(${PROJECT_SOURCE_DIR}/src)
+
+add_library(penna STATIC genome.cpp animal.cpp population.cpp fish_population.cpp)
+
+add_library(penna_vector STATIC genome.cpp animal.cpp population.cpp fish_population.cpp)
+set_target_properties(penna_vector PROPERTIES COMPILE_FLAGS "-DPENNA_VECTOR")

exercises/ex09/solutions/penna/src/animal.cpp

+/**
+ * Implementation of the Penna animal class.
+ * Programming Techniques for Scientific Simulations, ETH Zürich
+ */
+
+#include "animal.hpp"
+#include <cassert>
+
+namespace Penna {
+
+// Definition of static data members
+age_t Animal::bad_threshold_;
+age_t Animal::maturity_age_;
+double Animal::probability_to_get_pregnant_;
+
+void Animal::set_maturity_age( age_t r ) {
+    maturity_age_ = r;
+}
+
+void Animal::set_probability_to_get_pregnant( double p ) {
+    probability_to_get_pregnant_ = p;
+}
+
+void Animal::set_bad_threshold( age_t t ) {
+    bad_threshold_ = t;
+}
+
+age_t Animal::get_bad_threshold() {
+    return bad_threshold_;
+}
+
+Animal::Animal() : gen_(), age_(0), pregnant_(false) {
+}
+
+Animal::Animal( const Genome& gen ) : gen_(gen), age_(0), pregnant_(false) {
+}
+
+bool Animal::is_dead() const {
+    return age_ > maximum_age || gen_.count_bad(age_) > bad_threshold_;
+}
+
+bool Animal::is_pregnant() const {
+    return pregnant_;
+}
+
+age_t Animal::age() const {
+    return age_;
+}
+
+void Animal::grow() {
+    assert( !is_dead() );
+    ++age_;
+    if (age_ > maturity_age_ && !pregnant_)
+      if (probability_to_get_pregnant_ > drand48())
+        pregnant_ = true;
+}
+
+Animal Animal::give_birth() {
+    assert( pregnant_ == true );
+    pregnant_ = false;
+    Genome child_genome = gen_;
+    // Note: this "should" return a copy such that we don't expose a freely callable mutate
+    //       method but not doing it in-place is slow -> ideal use for C++11 move semantics
+    child_genome.mutate();
+    return Animal( child_genome );
+}
+
+} // end namespace Penna

exercises/ex09/solutions/penna/src/animal.hpp

+/**
+ * Header for the Penna animal class.
+ * Programming Techniques for Scientific Simulations, ETH Zürich
+ */
+
+#ifndef ANIMAL_HPP
+#define ANIMAL_HPP
+
+#include "genome.hpp"
+
+/// Penna namespace.
+namespace Penna {
+
+/// Animal with a Genome and age.
+class Animal {
+public:
+    static const age_t maximum_age = Genome::number_of_genes; ///< upper age limit
+
+    static void set_maturity_age( age_t );   ///< set maturity age
+    static void set_probability_to_get_pregnant ( double ); ///< set birth rate modifier
+
+    static void  set_bad_threshold( age_t ); ///< set threshold for deletrious mutation
+    static age_t get_bad_threshold();        ///< get deletrious mutation threshold
+
+    Animal();                    ///< default constructor: Uses all good genome
+    Animal( const Genome& );     ///< constructor using a given genome
+
+    bool is_dead() const;        ///< get death status
+    bool is_pregnant() const;    ///< get pregnancy status
+    age_t age() const;           ///< get age
+
+    void grow();                 ///< grow the animal older by one year
+    Animal give_birth();         ///< create baby with inherited, randomly mutated genes
+
+    Genome const & genome() const { ///< const getter to gather statistics about the genes
+        return gen_;
+    }
+private:
+    static age_t bad_threshold_; ///< number T of bad genes an animal can tolerate
+    static age_t maturity_age_;  ///< maturity age parameter R
+    static double probability_to_get_pregnant_; ///< birth rate modifier b
+
+    // We remove the const-qualifier here because we need the Copy-Assignable concept
+    #ifdef PENNA_VECTOR
+    Genome gen_;
+    #else
+    const Genome gen_;
+    #endif
+    age_t age_;
+    bool pregnant_;
+};
+
+} // namespace Penna
+
+#endif // ANIMAL_HPP

exercises/ex09/solutions/penna/src/fish_population.cpp

+#include "fish_population.hpp"
+#include <vector>
+#include <algorithm>
+#include <functional>
+#include <cstdlib>
+
+namespace Penna  {
+
+FishPopulation::FishPopulation(const std::size_t nmax,
+                               const std::size_t n0,
+                               const double f_rate,
+                               const std::size_t f_age)
+    : Population(nmax, n0), f_rate_(f_rate), f_age_(f_age) {}
+
+void FishPopulation::change_fishing(const double f_rate, const std::size_t f_age) {
+    f_rate_ = f_rate;
+    f_age_ = f_age;
+}
+
+class FishingPredicate {
+public:
+    FishingPredicate(const double probability, const std::size_t minage)
+        : probability_(probability), minage_(minage) {};
+
+    bool operator()(const Animal& a) const {
+        return a.age() > minage_ && (probability_ >= 1. || (rand() / double(RAND_MAX)) < probability_);
+    };
+
+private:
+    const double probability_;
+    const std::size_t minage_;
+};
+
+void FishPopulation::step() {
+    // Do normal aging.
+    Population::step();
+
+    // Fishing
+    if(f_rate_ > 0) {
+      population_.remove_if( FishingPredicate( f_rate_, f_age_ ) );
+    }
+}
+
+} // end namespace Penna

exercises/ex09/solutions/penna/src/fish_population.hpp

+#ifndef FISH_POPULATION_HPP
+#define FISH_POPULATION_HPP
+
+#include "population.hpp"
+
+namespace Penna {
+/**
+ * Population of animals with fishing.
+ */
+class FishPopulation : public Population {
+public:
+    /**
+     * Constructor.
+     * @param nmax Maximum population size. Parameter N_{max} in Penna's paper.
+     * @param n0 Initial population size.
+     * @param frate: Percentage of fish caught.
+     * @param fage: Minimum age for removal by fishing.
+     */
+    FishPopulation(const std::size_t nmax, const std::size_t n0, const double f_rate=0, const std::size_t f_age=0);
+
+    /// Change fishing rate and minium age of fishable fish.
+    void change_fishing(const double f_rate, const std::size_t f_age=0);
+    
+    /// Simulate one time step (year).
+    void step();
+    
+private:
+    double f_rate_;
+    std::size_t f_age_;
+};
+
+} // end namespace Penna
+
+#endif // !defined FISH_POPULATION_HPP

exercises/ex09/solutions/penna/src/genome.cpp

+/**
+ * Implementation of Penna genome class.
+ * Programming Techniques for Scientific Simulations, ETH Zürich
+ */
+
+#include "genome.hpp"
+#include <cstdlib> // drand48()
+#include <array>
+#include <numeric>
+
+namespace Penna {
+
+// Declaration of the static variables.
+age_t Genome::mutation_rate_;
+
+void Genome::set_mutation_rate( age_t m ) {
+    mutation_rate_ = m;
+}
+
+age_t Genome::count_bad( age_t n ) const {
+    // Intricate but fast implementation for builtin type sized bitsets:
+    // ignore the _first_ (max_age - age) genes and count defective in the rest.
+    return (genes_ << (number_of_genes-n-1)).count();
+}
+
+/// Pre: mutation_rate <= number_of_genes
+void Genome::mutate() {
+    // Fisher-Yates shuffle on random access data structure
+    // Mutate a random selection of M genes (without replacement)
+    std::array<unsigned char,number_of_genes> ages;
+    std::iota(std::begin(ages), std::end(ages), 0);
+    for(size_t i = 0; i < mutation_rate_; ++i) {
+        const size_t pivot = i + drand48() * (number_of_genes - i);
+        genes_.flip(ages[pivot]);
+        std::swap(ages[pivot], ages[i]);
+    }
+
+    // Branchless retry-mask with constant access and good low-density behaviour
+    //~ std::bitset<number_of_genes> used(0);
+    //~ size_t used_cnt = 0;
+    //~ while(used_cnt < mutation_rate_) {
+        //~ const size_t pivot = drand48() * (number_of_genes);
+        //~ genes_[pivot] = !(used[pivot] ^ genes_[pivot]);
+        //~ used_cnt += !used[pivot] * 1;
+        //~ used[pivot] = 1;
+    //~ }
+
+    // "Incorrect" implementation for future reference
+    //~ for(size_t i = 0; i < mutation_rate_; ++i)
+        //~ genes_.flip(int(drand48() * number_of_genes));
+}
+
+} // namespace Penna

exercises/ex09/solutions/penna/src/genome.hpp

+/**