"""
Evolutionary optimization of something
"""

import random

import numpy as np
import numpy.random as npr

from automata import SnowDrift


class EvolutionaryOptimizer(object):
    """ Optimize!
    """
    def __init__(self):
        """ Set some parameters
        """
        self.mutation_probability = 0.02

    def init(self, size):
        """ Generate initial population
        """
        raise NotImplementedError

    def get_fitness(self, obj):
        """ Compute fitness of individual of population
        """
        raise NotImplementedError

    def mutate(self, obj):
        """ Mutate single individual
        """
        raise NotImplementedError

    def crossover(self, mom, dad):
        """ Generate offspring from parents
        """
        raise NotImplementedError

    def run(self, size, max_iter=100):
        """ Let life begin
        """
        population = self.init(size)

        for _ in range(max_iter):
            #print(population)
            pop_fitness = [self.get_fitness(o) for o in population]
            #print(pop_fitness)

            # crossover best individuals and replace worst with child
            best_indiv = np.argpartition(pop_fitness, -2)[-2:]
            mom, dad = population[best_indiv]
            child = self.crossover(mom, dad)
            #print(mom, dad, '->', child)

            worst_indiv = np.argmin(pop_fitness)
            population[worst_indiv] = child

            # apply mutations
            mut = lambda o: \
                self.mutate(o) if random.random() < self.mutation_probability \
                else o
            population = np.array([mut(o) for o in population])

            print('Mean individual:', np.mean(population, axis=0))
            #print()
            yield np.mean(population, axis=0)

class SnowdriftOptimizer(EvolutionaryOptimizer):
    """ Optimize snowdrift game by assuming each individual to be the pair of
        benefit and cost floats
    """
    def init(self, size):
        pop = []
        for _ in range(size):
            pop.append((random.uniform(0, 1), random.uniform(0, 1)))
        return np.array(pop)

    def crossover(self, mom, dad):
        return np.mean([mom, dad], axis=0)

    def mutate(self, obj):
        sigma = 0.05
        return (obj[0] * random.gauss(1, sigma), obj[1] * random.gauss(1, sigma))

    def get_fitness(self, obj):
        # setup system
        lattice = npr.random_integers(0, 1, size=(15, 15))
        model = SnowDrift(lattice)

        # generate dynamics
        iter_num = 100

        benefit, cost = obj
        res = list(model.iterate(iter_num, benefit=benefit, cost=cost))

        # cut off transient
        ss = res[-int(iter_num/10):]

        # compute fitness
        fit = -np.sum(ss)
        return fit

def main():
    """ Setup environment
    """
    opti = SnowdriftOptimizer()
    opti.run(4)

if __name__ == '__main__':
    main()

@@ -54,16 +54,39 @@
 t random
+%0Aimport multiprocessing
 %0A%0Aimport
@@ -126,16 +126,71 @@
 as npr%0A%0A
+import matplotlib.pylab as plt%0A%0Afrom tqdm import tqdm%0A%0A
 from aut
@@ -984,24 +984,41 @@
 init(size)%0A%0A
+        res = %5B%5D%0A
         for 
@@ -1018,24 +1018,29 @@
    for _ in 
+tqdm(
 range(max_it
@@ -1042,16 +1042,17 @@
 ax_iter)
+)
 :%0A      
@@ -1778,16 +1778,17 @@
         
+#
 print('M
@@ -1871,14 +1871,19 @@
     
-yield 
+res.append(
 np.m
@@ -1905,16 +1905,36 @@
  axis=0)
+)%0A        return res
 %0A%0Aclass 
@@ -2585,14 +2585,12 @@
 ze=(
-15
+2
 , 1
-5
 ))%0A 
@@ -2923,41 +2923,630 @@
 def 
-main():%0A    %22%22%22 Setup environment
+plot_runs(runs):%0A    %22%22%22 Plot population evolutions%0A    %22%22%22%0A    cmap = plt.get_cmap('viridis')%0A    for i, r in enumerate(runs):%0A        bs, cs = zip(*r)%0A        color = cmap(float(i)/len(runs))%0A%0A        plt.plot(bs, '-', c=color)%0A        plt.plot(cs, '--', c=color)%0A%0A    plt.title('population evolution overview')%0A    plt.xlabel('time')%0A    plt.ylabel('value')%0A%0A    plt.ylim((0, 1))%0A%0A    plt.plot(0, 0, '-', c='black', label='benefit value')%0A    plt.plot(0, 0, '--', c='black', label='cost value')%0A    plt.legend(loc='best')%0A%0A    plt.savefig('result.pdf')%0A    plt.show()%0A%0A%0Adef work(i):%0A    %22%22%22 Handle one optimization case
 %0A   
@@ -3582,24 +3582,31 @@
 mizer()%0A    
+return 
 opti.run(4)%0A
@@ -3606,9 +3606,280 @@
 run(
-4
+20)%0A%0Adef main():%0A    %22%22%22 Setup environment%0A    %22%22%22%0A    core_num = int(multiprocessing.cpu_count() * 4/5)%0A    print('Using %25d cores' %25 core_num)%0A%0A    with multiprocessing.Pool(core_num) as p:%0A        runs = %5Bi for i in p.imap_unordered(work, range(10))%5D%0A%0A    plot_runs(runs
 )%0A%0Ai

Add plotting and multiprocessing