Add ex 12 solution

exercises/ex12/solutions/expression_templates/algorithmic_derivation.hpp

+/* Programming Techniques for Scientific Simulations, HS 2020
+ * Exercise 12.1
+ */
+
+#ifndef ALGO_DERIVATION
+#define ALGO_DERIVATION
+
+enum OP_enum {Add, Multiply};
+
+//********
+
+template<typename T>
+class Constant {
+public:
+    Constant(const T & v) : val_(v) {}
+    T operator()(const T &) const {
+        return val_;
+    }
+    T derivative(const T &) const {
+        return 0;
+    }
+private:
+    T val_;
+};
+
+// this function is just here so we don't need to write <int>. The function
+// can deduce it (whereas the class can't)
+template <typename T>
+Constant<T> constant(const T & x) {
+    return Constant<T>(x);
+}
+
+//********
+
+template <typename T>
+class Variable {
+public:
+    T operator()(const T & x) const {
+        return x;
+    }
+    T derivative(const T & x) const {
+        return 1;
+    }
+};
+
+//********
+
+template<typename L, typename R, OP_enum op>
+class Expression {
+public:
+    Expression(const L & l, const R & r) : l_(l), r_(r) { }
+
+    template <typename T>
+    T operator()(const T & x) const {
+        switch (op) { // we could just as well use if ... else if .... else
+            case Add:
+                return l_(x) + r_(x);
+            case Multiply:
+                return l_(x) * r_(x);
+        }
+    }
+
+    template <typename T>
+    T derivative(const T & x) const {
+        switch (op)
+        {
+            case Add:
+                return l_.derivative(x) + r_.derivative(x);
+            case Multiply:
+                return r_(x)*l_.derivative(x) + l_(x)*r_.derivative(x);
+        }
+    }
+
+private:
+    L l_;
+    R r_;
+};
+
+//********
+
+template<typename L, typename R>
+Expression<L, R, Multiply> operator*(const L & l, const R & r) {
+    return Expression<L, R, Multiply>(l, r);
+}
+
+template<typename L, typename R>
+Expression<L, R, Add> operator+(const L & l, const R & r) {
+    return Expression<L, R, Add>(l, r);
+}
+
+#endif //ALGO_DERIVATION

exercises/ex12/solutions/expression_templates/algorithmic_derivation_bonus.hpp

+/* Programming Techniques for Scientific Simulations, HS 2020
+ * Exercise 12.1
+ */
+
+#ifndef ALGO_DERIVATION
+#define ALGO_DERIVATION
+
+enum OP_enum {Add, Multiply};
+
+//********
+
+template<typename T>
+class Constant {
+public:
+    using derivative_type = Constant<T>;
+
+    Constant(const T & v) : val_(v) {}
+    T operator()(const T &) const {
+        return val_;
+    }
+    derivative_type derivative() const {
+        return derivative_type(0);
+    }
+private:
+    T val_;
+};
+
+// this function is just here so we don't need to write <int>. The function
+// can deduce it (whereas the class can't)
+template <typename T>
+Constant<T> constant(const T & x) {
+    return Constant<T>(x);
+}
+
+//********
+
+template <typename T>
+class Variable {
+public:
+    using derivative_type = Constant<T>;
+
+    T operator()(const T & x) const {
+        return x;
+    }
+    derivative_type derivative() const {
+        return derivative_type(1);
+    }
+};
+
+//********
+template<typename L, typename R, OP_enum op>
+class Expression;
+
+template<typename L, typename R, OP_enum op>
+struct expression_helper;
+
+template<typename L, typename R>
+struct expression_helper<L, R, Add> {
+    using derivative_type = Expression< typename L::derivative_type
+                                      , typename R::derivative_type
+                                      , Add
+                                      >;
+    static derivative_type derivative(const L & l, const R & r) {
+        return derivative_type(l.derivative(), r.derivative());
+    }
+};
+
+template<typename L, typename R>
+struct expression_helper<L, R, Multiply> {
+    using left_helper = Expression< typename L::derivative_type
+                                  , R
+                                  , Multiply
+                                  >;
+    using right_helper = Expression< L
+                                   , typename R::derivative_type
+                                   , Multiply
+                                   >;
+    using derivative_type = Expression< left_helper
+                                      , right_helper
+                                      , Add
+                                      >;
+    static derivative_type derivative(const L & l, const R & r) {
+        return derivative_type(left_helper(l.derivative(), r)
+                             , right_helper(l, r.derivative())
+                             );
+    }
+};
+
+
+
+template<typename L, typename R, OP_enum op>
+class Expression {
+public:
+    using derivative_type = typename expression_helper<L, R, op>::derivative_type;
+
+    Expression(const L & l, const R & r) : l_(l), r_(r) { }
+
+    template <typename T>
+    T operator()(const T & x) const {
+        switch (op) { // we could just as well use if ... else if .... else
+            case Add:
+                return l_(x) + r_(x);
+            case Multiply:
+                return l_(x) * r_(x);
+        }
+    }
+
+    derivative_type derivative() const {
+        return expression_helper<L, R, op>::derivative(l_, r_);
+    }
+
+private:
+    L l_;
+    R r_;
+};
+
+//********
+
+template<typename L, typename R>
+Expression<L, R, Multiply> operator*(const L & l, const R & r) {
+    return Expression<L, R, Multiply>(l, r);
+}
+
+template<typename L, typename R>
+Expression<L, R, Add> operator+(const L & l, const R & r) {
+    return Expression<L, R, Add>(l, r);
+}
+
+#endif //ALGO_DERIVATION

exercises/ex12/solutions/expression_templates/algorithmic_derivation_main.cpp

+/* Programming Techniques for Scientific Simulations, HS 2020
+ * Exercise 12.1
+ */
+
+#include <iostream>
+#include "algorithmic_derivation.hpp"
+
+// expression: 5*x*(x+1)+72*x+38 = 5*x^2 + 77*x + 38
+// derivative: 10*x+77
+
+int expr(const int _x)
+{
+    Variable<int> x;
+    return ( constant(5) * x * (x + constant(1))
+           + constant(72) * x + constant(38)
+           )(_x);
+}
+
+int expr_deriv(const int _x)
+{
+    Variable<int> x;
+    return ( constant(5) * x * (x + constant(1))
+           + constant(72) * x + constant(38)
+           ).derivative(_x);
+}
+
+int main()
+{
+    std::cout << expr(8) << ' ' << expr_deriv(8) << std::endl;
+    return 0;
+}

exercises/ex12/solutions/expression_templates/algorithmic_derivation_main_bonus.cpp

+/* Programming Techniques for Scientific Simulations, HS 2020
+ * Exercise 12.1
+ */
+
+#include <iostream>
+#include "algorithmic_derivation_bonus.hpp"
+
+// expression: 5*x*(x+1)+72*x+38 = 5*x^2 + 77*x + 38
+// derivative: 10*x+77
+
+int expr(const int _x)
+{
+    Variable<int> x;
+    return ( constant(5) * x * (x + constant(1))
+           + constant(72) * x + constant(38)
+           )(_x);
+}
+
+int expr_deriv(const int _x)
+{
+    Variable<int> x;
+    return ( constant(5) * x * (x + constant(1))
+           + constant(72) * x + constant(38)
+           ).derivative()(_x);
+}
+
+int main()
+{
+    std::cout << expr(8) << ' ' << expr_deriv(8) << std::endl;
+    return 0;
+}

exercises/ex12/solutions/hello_world/src/hello.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Programming Techniques for Scientific Simulations, HS 2020
+# Exercise 12.3
+
+def hello():
+    print('Hello', end=' ')
+
+
+if __name__ == '__main__':
+    print('hello.py executed as main program.')

exercises/ex12/solutions/hello_world/src/include/world.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Programming Techniques for Scientific Simulations, HS 2020
+# Exercise 12.3
+
+def world():
+    print('World!')

exercises/ex12/solutions/hello_world/src/main.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Programming Techniques for Scientific Simulations, HS 2020
+# Exercise 12.3
+
+# Import module 'hello'.
+import hello
+
+# Import object 'world' from submodule 'world' of package 'include'.
+from include.world import world
+
+if __name__ == '__main__':
+    hello.hello()
+    world()

exercises/ex12/solutions/pypenna/main.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Programming Techniques for Scientific Simulations, HS 2020
+# Exercise 12.4
+
+# For Python 2.x compatibility.
+from __future__ import division, print_function
+
+import random
+
+from pypenna import Population, Fishing
+
+
+def run(pop, time, of):
+    while pop.time < time:
+        pop.step()
+        print(pop.time)
+        of.write('{0.time}\t{0.N}\n'.format(pop))
+
+
+def main(N, M1, p1, M2, p2):
+    random.seed(0)
+
+    pop = Population(
+        mutation_rate=1,
+        bad_threshold=3,
+        repr_age=7,
+        N_max=200000,
+        N_0=20000,
+    )
+
+    with open('population.dat', 'w') as of:
+        of.write('# time\tN\n')
+        run(pop, M1, of)
+        pop.fishing = Fishing(probability=p1)
+        run(pop, M2, of)
+        pop.fishing = Fishing(probability=p2)
+        run(pop, N, of)
+
+
+if __name__ == "__main__":
+    main(N=1600, M1=1000, p1=0.17, M2=1200, p2=0.22)

exercises/ex12/solutions/pypenna/plot.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Programming Techniques for Scientific Simulations, HS 2020
+# Exercise 12.4
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+pop = np.loadtxt('population.dat')
+
+plt.figure()
+plt.plot(pop[:, 0], pop[:, 1], 'x')
+plt.yscale('log')
+plt.xlabel('year')
+plt.ylabel('population size')
+plt.savefig('population.pdf', bbox_inches='tight')
+
+plt.show()

exercises/ex12/solutions/pypenna/pypenna/__init__.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Programming Techniques for Scientific Simulations, HS 2020
+# Exercise 12.4
+
+from .population import *
+from .animal import *
+from .genome import *
+from .fishing import *

exercises/ex12/solutions/pypenna/pypenna/animal.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Programming Techniques for Scientific Simulations, HS 2020
+# Exercise 12.4
+
+from __future__ import division, print_function
+
+from .genome import Genome
+
+import random
+
+
+class Animal(object):
+    def __init__(self, genome=None, age=0):
+        self.genome = Genome(genome)
+        self.age = age
+
+    # property decorator for getter method
+    @property
+    def is_pregnant(self):
+        """
+        Returns whether a fish is pregnant or not.
+        """
+        if self.is_adult:
+            return random.random() < self.__class__.pregnancy_prob
+        return False
+
+    @property
+    def is_adult(self):
+        """
+        Returns whether the fish has reached the reproductive age.
+        """
+        return self.age >= self.__class__.repr_age
+
+    def make_child(self):
+        """
+        Returns a child fish with age 0 and a genome that is a mutated
+        version of the parent's genome.
+        """
+        child = Animal(Genome(self.genome))
+        child.genome.mutate()
+
+        return child
+
+    def grow(self):
+        """
+        Makes the fish one year older.
+        """
+        self.age += 1
+
+    @property
+    def is_sick(self):
+        """
+        Returns whether the fish is sick due to too many bad genes.
+        """
+        return (
+            self.genome.num_bad(self.age) >= self.genome.bad_threshold or
+            self.age >= self.genome.gene_size
+        )

exercises/ex12/solutions/pypenna/pypenna/fishing.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Programming Techniques for Scientific Simulations, HS 2020
+# Exercise 12.4
+
+from __future__ import division, print_function
+
+import random
+
+
+class Fishing(object):
+    """
+    Class for determining whether a given fish is being fished, given
+    a certain fishing probability and minimum fishing age.
+    """
+    def __init__(self, probability=0, min_age=0):
+        self.probability = probability
+        self.min_age = min_age
+
+    def is_fished(self, fish):
+        """
+        Determines whether the given fish is being fished. Returns True
+        if the fish gets fished, and False otherwise.
+        """
+        if fish.age >= self.min_age:
+            return random.random() < self.probability
+        return False

exercises/ex12/solutions/pypenna/pypenna/genome.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Programming Techniques for Scientific Simulations, HS 2020
+# Exercise 12.4
+
+# For Python 2.x compatibility.
+from __future__ import division, print_function
+
+import copy
+import random
+import numpy as np
+
+
+# Inherit from object to make Genome a new-style class (only required in Python 2.x).
+class Genome(object):
+    gene_size = 64 # Class variable shared by all instances.
+
+    def __init__(self, genome=None):
+        if genome is None:
+            # np.array with dtype=bool is a "bitset".
+            self.genes = np.zeros(self.__class__.gene_size, dtype=bool)
+        else:
+            # Instance variable unique to each instance.
+            self.genes = copy.deepcopy(genome.genes)
+
+    def mutate(self):
+        """
+        Mutates the Genome in-place.
+        """
+        for _ in range(self.__class__.mutation_rate):
+            self._flip(random.randint(0, self.__class__.gene_size - 1))
+
+    def _flip(self, idx):
+        """
+        Flips a single gene.
+        """
+        self.genes[idx] = not self.genes[idx]
+
+    def num_bad(self, age):
+        """
+        Returns the number of bad genes up to a certain age.
+        """
+        return sum(self.genes[:age])

exercises/ex12/solutions/pypenna/pypenna/population.py

+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+# Programming Techniques for Scientific Simulations, HS 2020
+# Exercise 12.4
+
+from __future__ import division, print_function
+
+from .animal import Animal
+from .genome import Genome
+from .fishing import Fishing
+
+import copy
+import random
+
+
+class Population(object):
+    """
+    Simulation class for the Penna model.
+    """
+    def __init__(
+            self,
+            mutation_rate,
+            bad_threshold,
+            repr_age,
+            N_max,
+            N_0,
+            fishing=None,
+            pregnancy_prob=1
+        ):
+        Genome.mutation_rate = mutation_rate
+        Genome.bad_threshold = bad_threshold
+        Animal.repr_age = repr_age
+
+        self.N_max = N_max
+        self.time = 0
+
+        if fishing is not None:
+            self.fishing = fishing
+        else:
+            self.fishing = Fishing()
+
+        Animal.pregnancy_prob = pregnancy_prob
+
+        self.population = [Animal() for _ in range(N_0)]
+
+    @property
+    def N(self):
+        """
+        Number of fish in the population
+        """
+        return len(self.population)
+
+    def step(self):
+        """
+        Performs the iteration for one year
+        """
+        for fish in self.population:
+            fish.grow()
+
+        self.population = [
+            fish for fish in self.population
+            if not fish.is_sick
+            if not random.random() < self.N / self.N_max
+            if not self.fishing.is_fished(fish)
+        ]
+
+        # Copy the population to create the "parents" list.
+        for fish in copy.copy(self.population):
+            if fish.is_pregnant:
+                self.population.append(fish.make_child())
+
+        self.time += 1