Tuner: Moved tuning logic into the python wrapper - draft of Android tuning app using kivy

cmake/python/setup.py

@@ -143,7 +143,7 @@ def main():
                     libraries=libraries)]
     
     #External
-    extensions += [Extension('external._tree',
+    extensions += [Extension('autotuning.external.sklearn._tree',
                              ['external/sklearn/_tree.c'],
                              include_dirs = [numpy_include])]
 
@@ -155,7 +155,7 @@ def main():
                 author='Philippe Tillet',
                 author_email='ptillet@g.harvard.edu',
                 license='MPL 2.0',
-                packages=['isaac','isaac.external'],
+                packages=['isaac','isaac.autotuning', 'isaac.autotuning.external', 'isaac.autotuning.external.deap', 'isaac.autotuning.external.deap.tools', 'isaac.autotuning.external.sklearn'],
                 ext_package="isaac",
                 ext_modules=extensions,
                 cmdclass={'build_py': build_py, 'build_ext': build_ext_subclass},

python/external/sklearn/_tree.pyx

@@ -23,11 +23,8 @@ from cpython cimport Py_INCREF, PyObject
 import numpy as np
 cimport numpy as np
 np.import_array()
-
 from _tree cimport StackRecord
 
-from scipy.sparse import issparse, csc_matrix, csr_matrix
-
 cdef extern from "numpy/arrayobject.h":
     object PyArray_NewFromDescr(object subtype, np.dtype descr,
                                 int nd, np.npy_intp* dims,
@@ -2063,9 +2060,6 @@ cdef class BaseSparseSplitter(Splitter):
         # Call parent init
         Splitter.init(self, X, y, sample_weight)
 
-        if not isinstance(X, csc_matrix):
-            raise ValueError("X should be in csc format")
-
         cdef SIZE_t* samples = self.samples
         cdef SIZE_t n_samples = self.n_samples
 
@@ -2790,18 +2784,7 @@ cdef class TreeBuilder:
     cdef inline _check_input(self, object X, np.ndarray y,
                              np.ndarray sample_weight):
         """Check input dtype, layout and format"""
-        if issparse(X):
+        if X.dtype != DTYPE:
-            X = X.tocsc()
-            X.sort_indices()
-
-            if X.data.dtype != DTYPE:
-                X.data = np.ascontiguousarray(X.data, dtype=DTYPE)
-
-            if X.indices.dtype != np.int32 or X.indptr.dtype != np.int32:
-                raise ValueError("No support for np.int64 index based "
-                                 "sparse matrices")
-
-        elif X.dtype != DTYPE:
             # since we have to copy we will make it fortran for efficiency
             X = np.asfortranarray(X, dtype=DTYPE)
 
@@ -3430,10 +3413,7 @@ cdef class Tree:
 
     cpdef np.ndarray apply(self, object X):
         """Finds the terminal region (=leaf node) for each sample in X."""
-        if issparse(X):
+        return self._apply_dense(X)
-            return self._apply_sparse_csr(X)
-        else:
-            return self._apply_dense(X)
 
 
     cdef inline np.ndarray _apply_dense(self, object X):
@@ -3482,10 +3462,6 @@ cdef class Tree:
         """Finds the terminal region (=leaf node) for each sample in sparse X.
 
         """
-        # Check input
-        if not isinstance(X, csr_matrix):
-            raise ValueError("X should be in csr_matrix format, got %s"
-                             % type(X))
 
         if X.dtype != DTYPE:
             raise ValueError("X.dtype should be np.float32, got %s" % X.dtype)

python/isaac/autotuning/external/deap/__init__.py

@@ -0,0 +1,17 @@
+#    This file is part of DEAP.
+#
+#    DEAP is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU Lesser General Public License as
+#    published by the Free Software Foundation, either version 3 of
+#    the License, or (at your option) any later version.
+#
+#    DEAP is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#    GNU Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.
+__author__ = "DEAP Team"
+__version__ = "1.0"
+__revision__ = "1.0.2"

python/isaac/autotuning/external/deap/algorithms.py

@@ -0,0 +1,472 @@
+#    This file is part of DEAP.
+#
+#    DEAP is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU Lesser General Public License as
+#    published by the Free Software Foundation, either version 3 of
+#    the License, or (at your option) any later version.
+#
+#    DEAP is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#    GNU Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.
+
+"""The :mod:`algorithms` module is intended to contain some specific algorithms
+in order to execute very common evolutionary algorithms. The method used here
+are more for convenience than reference as the implementation of every 
+evolutionary algorithm may vary infinitely. Most of the algorithms in this 
+module use operators registered in the toolbox. Generaly, the keyword used are
+:meth:`mate` for crossover, :meth:`mutate` for mutation, :meth:`~deap.select`
+for selection and :meth:`evaluate` for evaluation.
+
+You are encouraged to write your own algorithms in order to make them do what
+you really want them to do.
+"""
+
+import random
+
+import tools
+
+def varAnd(population, toolbox, cxpb, mutpb):
+    """Part of an evolutionary algorithm applying only the variation part
+    (crossover **and** mutation). The modified individuals have their
+    fitness invalidated. The individuals are cloned so returned population is
+    independent of the input population.
+    
+    :param population: A list of individuals to vary.
+    :param toolbox: A :class:`~deap.base.Toolbox` that contains the evolution
+                    operators.
+    :param cxpb: The probability of mating two individuals.
+    :param mutpb: The probability of mutating an individual.
+    :returns: A list of varied individuals that are independent of their
+              parents.
+    
+    The variation goes as follow. First, the parental population
+    :math:`P_\mathrm{p}` is duplicated using the :meth:`toolbox.clone` method
+    and the result is put into the offspring population :math:`P_\mathrm{o}`.
+    A first loop over :math:`P_\mathrm{o}` is executed to mate pairs of consecutive
+    individuals. According to the crossover probability *cxpb*, the
+    individuals :math:`\mathbf{x}_i` and :math:`\mathbf{x}_{i+1}` are mated
+    using the :meth:`toolbox.mate` method. The resulting children
+    :math:`\mathbf{y}_i` and :math:`\mathbf{y}_{i+1}` replace their respective
+    parents in :math:`P_\mathrm{o}`. A second loop over the resulting
+    :math:`P_\mathrm{o}` is executed to mutate every individual with a
+    probability *mutpb*. When an individual is mutated it replaces its not
+    mutated version in :math:`P_\mathrm{o}`. The resulting
+    :math:`P_\mathrm{o}` is returned.
+    
+    This variation is named *And* beceause of its propention to apply both
+    crossover and mutation on the individuals. Note that both operators are
+    not applied systematicaly, the resulting individuals can be generated from
+    crossover only, mutation only, crossover and mutation, and reproduction
+    according to the given probabilities. Both probabilities should be in
+    :math:`[0, 1]`.
+    """
+    offspring = [toolbox.clone(ind) for ind in population]
+    
+    # Apply crossover and mutation on the offspring
+    for i in range(1, len(offspring), 2):
+        if random.random() < cxpb:
+            offspring[i-1], offspring[i] = toolbox.mate(offspring[i-1], offspring[i])
+            del offspring[i-1].fitness.values, offspring[i].fitness.values
+    
+    for i in range(len(offspring)):
+        if random.random() < mutpb:
+            offspring[i], = toolbox.mutate(offspring[i])
+            del offspring[i].fitness.values
+    
+    return offspring
+
+def eaSimple(population, toolbox, cxpb, mutpb, ngen, stats=None,
+             halloffame=None, verbose=__debug__):
+    """This algorithm reproduce the simplest evolutionary algorithm as
+    presented in chapter 7 of [Back2000]_.
+    
+    :param population: A list of individuals.
+    :param toolbox: A :class:`~deap.base.Toolbox` that contains the evolution
+                    operators.
+    :param cxpb: The probability of mating two individuals.
+    :param mutpb: The probability of mutating an individual.
+    :param ngen: The number of generation.
+    :param stats: A :class:`~deap.tools.Statistics` object that is updated
+                  inplace, optional.
+    :param halloffame: A :class:`~deap.tools.HallOfFame` object that will
+                       contain the best individuals, optional.
+    :param verbose: Whether or not to log the statistics.
+    :returns: The final population and a :class:`~deap.tools.Logbook`
+              with the statistics of the evolution.
+    
+    The algorithm takes in a population and evolves it in place using the
+    :meth:`varAnd` method. It returns the optimized population and a
+    :class:`~deap.tools.Logbook` with the statistics of the evolution (if
+    any). The logbook will contain the generation number, the number of
+    evalutions for each generation and the statistics if a
+    :class:`~deap.tools.Statistics` if any. The *cxpb* and *mutpb* arguments
+    are passed to the :func:`varAnd` function. The pseudocode goes as follow
+    ::
+
+        evaluate(population)
+        for g in range(ngen):
+            population = select(population, len(population))
+            offspring = varAnd(population, toolbox, cxpb, mutpb)
+            evaluate(offspring)
+            population = offspring
+
+    As stated in the pseudocode above, the algorithm goes as follow. First, it
+    evaluates the individuals with an invalid fitness. Second, it enters the
+    generational loop where the selection procedure is applied to entirely
+    replace the parental population. The 1:1 replacement ratio of this
+    algorithm **requires** the selection procedure to be stochastic and to
+    select multiple times the same individual, for example,
+    :func:`~deap.tools.selTournament` and :func:`~deap.tools.selRoulette`.
+    Third, it applies the :func:`varAnd` function to produce the next
+    generation population. Fourth, it evaluates the new individuals and
+    compute the statistics on this population. Finally, when *ngen*
+    generations are done, the algorithm returns a tuple with the final
+    population and a :class:`~deap.tools.Logbook` of the evolution.
+
+    .. note::
+
+        Using a non-stochastic selection method will result in no selection as
+        the operator selects *n* individuals from a pool of *n*.
+    
+    This function expects the :meth:`toolbox.mate`, :meth:`toolbox.mutate`,
+    :meth:`toolbox.select` and :meth:`toolbox.evaluate` aliases to be
+    registered in the toolbox.
+    
+    .. [Back2000] Back, Fogel and Michalewicz, "Evolutionary Computation 1 :
+       Basic Algorithms and Operators", 2000.
+    """
+    logbook = tools.Logbook()
+    logbook.header = ['gen', 'nevals'] + (stats.fields if stats else [])
+
+    # Evaluate the individuals with an invalid fitness
+    invalid_ind = [ind for ind in population if not ind.fitness.valid]
+    fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
+    for ind, fit in zip(invalid_ind, fitnesses):
+        ind.fitness.values = fit
+
+    if halloffame is not None:
+        halloffame.update(population)
+
+    record = stats.compile(population) if stats else {}
+    logbook.record(gen=0, nevals=len(invalid_ind), **record)
+    if verbose:
+        print logbook.stream
+
+    # Begin the generational process
+    for gen in range(1, ngen+1):
+        # Select the next generation individuals
+        offspring = toolbox.select(population, len(population))
+        
+        # Vary the pool of individuals
+        offspring = varAnd(offspring, toolbox, cxpb, mutpb)
+        
+        # Evaluate the individuals with an invalid fitness
+        invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
+        fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
+        for ind, fit in zip(invalid_ind, fitnesses):
+            ind.fitness.values = fit
+        
+        # Update the hall of fame with the generated individuals
+        if halloffame is not None:
+            halloffame.update(offspring)
+            
+        # Replace the current population by the offspring
+        population[:] = offspring
+        
+        # Append the current generation statistics to the logbook
+        record = stats.compile(population) if stats else {}
+        logbook.record(gen=gen, nevals=len(invalid_ind), **record)
+        if verbose:
+            print logbook.stream        
+
+    return population, logbook
+
+def varOr(population, toolbox, lambda_, cxpb, mutpb):
+    """Part of an evolutionary algorithm applying only the variation part
+    (crossover, mutation **or** reproduction). The modified individuals have
+    their fitness invalidated. The individuals are cloned so returned
+    population is independent of the input population.
+    
+    :param population: A list of individuals to vary.
+    :param toolbox: A :class:`~deap.base.Toolbox` that contains the evolution
+                    operators.
+    :param lambda\_: The number of children to produce
+    :param cxpb: The probability of mating two individuals.
+    :param mutpb: The probability of mutating an individual.
+    :returns: A list of varied individuals that are independent of their
+              parents.
+    
+    The variation goes as follow. On each of the *lambda_* iteration, it
+    selects one of the three operations; crossover, mutation or reproduction.
+    In the case of a crossover, two individuals are selected at random from
+    the parental population :math:`P_\mathrm{p}`, those individuals are cloned
+    using the :meth:`toolbox.clone` method and then mated using the
+    :meth:`toolbox.mate` method. Only the first child is appended to the
+    offspring population :math:`P_\mathrm{o}`, the second child is discarded.
+    In the case of a mutation, one individual is selected at random from
+    :math:`P_\mathrm{p}`, it is cloned and then mutated using using the
+    :meth:`toolbox.mutate` method. The resulting mutant is appended to
+    :math:`P_\mathrm{o}`. In the case of a reproduction, one individual is
+    selected at random from :math:`P_\mathrm{p}`, cloned and appended to
+    :math:`P_\mathrm{o}`.
+    
+    This variation is named *Or* beceause an offspring will never result from
+    both operations crossover and mutation. The sum of both probabilities
+    shall be in :math:`[0, 1]`, the reproduction probability is
+    1 - *cxpb* - *mutpb*.
+    """
+    assert (cxpb + mutpb) <= 1.0, ("The sum of the crossover and mutation "
+        "probabilities must be smaller or equal to 1.0.")
+    
+    offspring = []
+    for _ in xrange(lambda_):
+        op_choice = random.random()
+        if op_choice < cxpb:            # Apply crossover
+            ind1, ind2 = map(toolbox.clone, random.sample(population, 2))
+            ind1, ind2 = toolbox.mate(ind1, ind2)
+            del ind1.fitness.values
+            offspring.append(ind1)
+        elif op_choice < cxpb + mutpb:  # Apply mutation
+            ind = toolbox.clone(random.choice(population))
+            ind, = toolbox.mutate(ind)
+            del ind.fitness.values
+            offspring.append(ind)
+        else:                           # Apply reproduction
+            offspring.append(random.choice(population))
+    
+    return offspring
+
+def eaMuPlusLambda(population, toolbox, mu, lambda_, cxpb, mutpb, ngen,
+                   stats=None, halloffame=None, verbose=__debug__):
+    """This is the :math:`(\mu + \lambda)` evolutionary algorithm.
+    
+    :param population: A list of individuals.
+    :param toolbox: A :class:`~deap.base.Toolbox` that contains the evolution
+                    operators.
+    :param mu: The number of individuals to select for the next generation.
+    :param lambda\_: The number of children to produce at each generation.
+    :param cxpb: The probability that an offspring is produced by crossover.
+    :param mutpb: The probability that an offspring is produced by mutation.
+    :param ngen: The number of generation.
+    :param stats: A :class:`~deap.tools.Statistics` object that is updated
+                  inplace, optional.
+    :param halloffame: A :class:`~deap.tools.HallOfFame` object that will
+                       contain the best individuals, optional.
+    :param verbose: Whether or not to log the statistics.
+    :returns: The final population and a :class:`~deap.tools.Logbook`
+              with the statistics of the evolution.
+    
+    The algorithm takes in a population and evolves it in place using the
+    :meth:`varOr` method. It returns the optimized population and a
+    :class:`~deap.tools.Logbook` with the statistics of the evolution (if
+    any). The logbook will contain the generation number, the number of
+    evalutions for each generation and the statistics if a
+    :class:`~deap.tools.Statistics` if any. The *cxpb* and *mutpb* arguments
+    are passed to the :func:`varAnd` function. The pseudocode goes as follow
+    ::
+
+        evaluate(population)
+        for g in range(ngen):
+            offspring = varOr(population, toolbox, lambda_, cxpb, mutpb)
+            evaluate(offspring)
+            population = select(population + offspring, mu)
+
+    First, the individuals having an invalid fitness are evaluated. Second,
+    the evolutionary loop begins by producing *lambda_* offspring from the
+    population, the offspring are generated by the :func:`varOr` function. The
+    offspring are then evaluated and the next generation population is
+    selected from both the offspring **and** the population. Finally, when
+    *ngen* generations are done, the algorithm returns a tuple with the final
+    population and a :class:`~deap.tools.Logbook` of the evolution.
+
+    .. note::
+
+        Care must be taken when the lambda:mu ratio is 1 to 1 as a non-stochastic
+        selection will result in no selection at all as
+        the operator selects *lambda* individuals from a pool of *mu*.
+
+    This function expects :meth:`toolbox.mate`, :meth:`toolbox.mutate`,
+    :meth:`toolbox.select` and :meth:`toolbox.evaluate` aliases to be
+    registered in the toolbox. This algorithm uses the :func:`varOr`
+    variation.
+    """
+    logbook = tools.Logbook()
+    logbook.header = ['gen', 'nevals'] + (stats.fields if stats else [])
+
+    # Evaluate the individuals with an invalid fitness
+    invalid_ind = [ind for ind in population if not ind.fitness.valid]
+    fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
+    for ind, fit in zip(invalid_ind, fitnesses):
+        ind.fitness.values = fit
+
+    if halloffame is not None:
+        halloffame.update(population)
+
+    record = stats.compile(population) if stats is not None else {}
+    logbook.record(gen=0, nevals=len(invalid_ind), **record)
+    if verbose:
+        print logbook.stream
+
+    # Begin the generational process
+    for gen in range(1, ngen+1):
+        # Vary the population
+        offspring = varOr(population, toolbox, lambda_, cxpb, mutpb)
+        
+        # Evaluate the individuals with an invalid fitness
+        invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
+        fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
+        for ind, fit in zip(invalid_ind, fitnesses):
+            ind.fitness.values = fit
+        
+        # Update the hall of fame with the generated individuals
+        if halloffame is not None:
+            halloffame.update(offspring)
+
+        # Select the next generation population
+        population[:] = toolbox.select(population + offspring, mu)
+
+        # Update the statistics with the new population
+        record = stats.compile(population) if stats is not None else {}
+        logbook.record(gen=gen, nevals=len(invalid_ind), **record)
+        if verbose:
+            print logbook.stream
+
+    return population, logbook
+    
+def eaMuCommaLambda(population, toolbox, mu, lambda_, cxpb, mutpb, ngen,
+                    stats=None, halloffame=None, verbose=__debug__):
+    """This is the :math:`(\mu~,~\lambda)` evolutionary algorithm.
+    
+    :param population: A list of individuals.    
+    :param toolbox: A :class:`~deap.base.Toolbox` that contains the evolution
+                    operators.
+    :param mu: The number of individuals to select for the next generation.
+    :param lambda\_: The number of children to produce at each generation.
+    :param cxpb: The probability that an offspring is produced by crossover.
+    :param mutpb: The probability that an offspring is produced by mutation.
+    :param ngen: The number of generation.
+    :param stats: A :class:`~deap.tools.Statistics` object that is updated
+                  inplace, optional.
+    :param halloffame: A :class:`~deap.tools.HallOfFame` object that will
+                       contain the best individuals, optional.
+    :param verbose: Whether or not to log the statistics.
+    :returns: The final population.
+    
+    First, the individuals having an invalid fitness are evaluated. Then, the
+    evolutionary loop begins by producing *lambda_* offspring from the
+    population, the offspring are generated by a crossover, a mutation or a
+    reproduction proportionally to the probabilities *cxpb*, *mutpb* and 1 -
+    (cxpb + mutpb). The offspring are then evaluated and the next generation
+    population is selected **only** from the offspring. Briefly, the operators
+    are applied as following ::
+    
+        evaluate(population)
+        for i in range(ngen):
+            offspring = varOr(population, toolbox, lambda_, cxpb, mutpb)
+            evaluate(offspring)
+            population = select(offspring, mu)
+    
+    This function expects :meth:`toolbox.mate`, :meth:`toolbox.mutate`,
+    :meth:`toolbox.select` and :meth:`toolbox.evaluate` aliases to be
+    registered in the toolbox. This algorithm uses the :func:`varOr`
+    variation.
+    """
+    assert lambda_ >= mu, "lambda must be greater or equal to mu."
+
+    # Evaluate the individuals with an invalid fitness
+    invalid_ind = [ind for ind in population if not ind.fitness.valid]
+    fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
+    for ind, fit in zip(invalid_ind, fitnesses):
+        ind.fitness.values = fit
+
+    if halloffame is not None:
+        halloffame.update(population)
+
+    logbook = tools.Logbook()
+    logbook.header = ['gen', 'nevals'] + (stats.fields if stats else [])
+
+    record = stats.compile(population) if stats is not None else {}
+    logbook.record(gen=0, nevals=len(invalid_ind), **record)
+    if verbose:
+        print logbook.stream
+
+    # Begin the generational process
+    for gen in range(1, ngen+1):
+        # Vary the population
+        offspring = varOr(population, toolbox, lambda_, cxpb, mutpb)
+
+        # Evaluate the individuals with an invalid fitness
+        invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
+        fitnesses = toolbox.map(toolbox.evaluate, invalid_ind)
+        for ind, fit in zip(invalid_ind, fitnesses):
+            ind.fitness.values = fit
+
+        # Update the hall of fame with the generated individuals
+        if halloffame is not None:
+            halloffame.update(offspring)
+
+        # Select the next generation population
+        population[:] = toolbox.select(offspring, mu)
+
+        # Update the statistics with the new population
+        record = stats.compile(population) if stats is not None else {}
+        logbook.record(gen=gen, nevals=len(invalid_ind), **record)
+        if verbose:
+            print logbook.stream
+
+
+    return population, logbook
+
+def eaGenerateUpdate(toolbox, ngen, halloffame=None, stats=None, 
+                     verbose=__debug__):
+    """This is algorithm implements the ask-tell model proposed in 
+    [Colette2010]_, where ask is called `generate` and tell is called `update`.
+    
+    :param toolbox: A :class:`~deap.base.Toolbox` that contains the evolution
+                    operators.
+    :param ngen: The number of generation.
+    :param stats: A :class:`~deap.tools.Statistics` object that is updated
+                  inplace, optional.
+    :param halloffame: A :class:`~deap.tools.HallOfFame` object that will
+                       contain the best individuals, optional.
+    :param verbose: Whether or not to log the statistics.
+
+    :returns: The final population.
+    
+    The toolbox should contain a reference to the generate and the update method 
+    of the chosen strategy.
+
+    .. [Colette2010] Collette, Y., N. Hansen, G. Pujol, D. Salazar Aponte and
+       R. Le Riche (2010). On Object-Oriented Programming of Optimizers -
+       Examples in Scilab. In P. Breitkopf and R. F. Coelho, eds.:
+       Multidisciplinary Design Optimization in Computational Mechanics,
+       Wiley, pp. 527-565;
+
+    """
+    logbook = tools.Logbook()
+    logbook.header = ['gen', 'nevals'] + (stats.fields if stats else [])
+
+    for gen in xrange(ngen):
+        # Generate a new population
+        population = toolbox.generate()
+        # Evaluate the individuals
+        fitnesses = toolbox.map(toolbox.evaluate, population)
+        for ind, fit in zip(population, fitnesses):
+            ind.fitness.values = fit
+        
+        if halloffame is not None:
+            halloffame.update(population)
+        
+        # Update the strategy with the evaluated individuals
+        toolbox.update(population)
+        
+        record = stats.compile(population) if stats is not None else {}
+        logbook.record(gen=gen, nevals=len(population), **record)
+        if verbose:
+            print logbook.stream
+
+    return population, logbook

python/isaac/autotuning/external/deap/base.py

@@ -0,0 +1,264 @@
+#    This file is part of DEAP.
+#
+#    DEAP is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU Lesser General Public License as
+#    published by the Free Software Foundation, either version 3 of
+#    the License, or (at your option) any later version.
+#
+#    DEAP is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#    GNU Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.
+
+"""The :mod:`~deap.base` module provides basic structures to build
+evolutionary algorithms. It contains the :class:`~deap.base.Toolbox`, useful
+to store evolutionary operators, and a virtual :class:`~deap.base.Fitness`
+class used as base class, for the fitness member of any individual. """
+
+import sys
+
+from collections import Sequence
+from copy import deepcopy
+from functools import partial 
+from operator import mul, truediv
+
+class Toolbox(object):
+    """A toolbox for evolution that contains the evolutionary operators. At
+    first the toolbox contains a :meth:`~deap.toolbox.clone` method that
+    duplicates any element it is passed as argument, this method defaults to
+    the :func:`copy.deepcopy` function. and a :meth:`~deap.toolbox.map`
+    method that applies the function given as first argument to every items
+    of the iterables given as next arguments, this method defaults to the
+    :func:`map` function. You may populate the toolbox with any other
+    function by using the :meth:`~deap.base.Toolbox.register` method.
+
+    Concrete usages of the toolbox are shown for initialization in the
+    :ref:`creating-types` tutorial and for tools container in the
+    :ref:`next-step` tutorial.
+    """
+
+    def __init__(self):
+        self.register("clone", deepcopy)
+        self.register("map", map)
+
+    def register(self, alias, function, *args, **kargs):
+        """Register a *function* in the toolbox under the name *alias*. You 
+        may provide default arguments that will be passed automatically when 
+        calling the registered function. Fixed arguments can then be overriden 
+        at function call time.
+        
+        :param alias: The name the operator will take in the toolbox. If the
+                      alias already exist it will overwrite the the operator
+                      already present.
+        :param function: The function to which refer the alias.
+        :param argument: One or more argument (and keyword argument) to pass
+                         automatically to the registered function when called,
+                         optional.
+        
+        The following code block is an example of how the toolbox is used. ::
+
+            >>> def func(a, b, c=3):
+            ...     print a, b, c
+            ... 
+            >>> tools = Toolbox()
+            >>> tools.register("myFunc", func, 2, c=4)
+            >>> tools.myFunc(3)
+            2 3 4
+        
+        The registered function will be given the attributes :attr:`__name__`
+        set to the alias and :attr:`__doc__` set to the original function's
+        documentation. The :attr:`__dict__` attribute will also be updated
+        with the original function's instance dictionnary, if any.
+        """
+        pfunc = partial(function, *args, **kargs)
+        pfunc.__name__ = alias
+        pfunc.__doc__ = function.__doc__
+        
+        if hasattr(function, "__dict__") and not isinstance(function, type):
+            # Some functions don't have a dictionary, in these cases 
+            # simply don't copy it. Moreover, if the function is actually
+            # a class, we do not want to copy the dictionary.
+            pfunc.__dict__.update(function.__dict__.copy())
+        
+        setattr(self, alias, pfunc)
+
+    def unregister(self, alias):
+        """Unregister *alias* from the toolbox.
+        
+        :param alias: The name of the operator to remove from the toolbox.
+        """
+        delattr(self, alias)
+
+    def decorate(self, alias, *decorators):
+        """Decorate *alias* with the specified *decorators*, *alias*
+        has to be a registered function in the current toolbox.
+        
+        :param alias: The name of the operator to decorate.
+        :param decorator: One or more function decorator. If multiple
+                          decorators are provided they will be applied in
+                          order, with the last decorator decorating all the
+                          others.
+        
+        .. note::
+            Decorate a function using the toolbox makes it unpicklable, and
+            will produce an error on pickling. Although this limitation is not
+            relevant in most cases, it may have an impact on distributed
+            environments like multiprocessing.
+            A function can still be decorated manually before it is added to
+            the toolbox (using the @ notation) in order to be picklable.
+        """
+        pfunc = getattr(self, alias)
+        function, args, kargs = pfunc.func, pfunc.args, pfunc.keywords
+        for decorator in decorators:
+            function = decorator(function)
+        self.register(alias, function, *args, **kargs)
+
+
+class Fitness(object):
+    """The fitness is a measure of quality of a solution. If *values* are
+    provided as a tuple, the fitness is initalized using those values,
+    otherwise it is empty (or invalid).
+    
+    :param values: The initial values of the fitness as a tuple, optional.
+
+    Fitnesses may be compared using the ``>``, ``<``, ``>=``, ``<=``, ``==``,
+    ``!=``. The comparison of those operators is made lexicographically.
+    Maximization and minimization are taken care off by a multiplication
+    between the :attr:`weights` and the fitness :attr:`values`. The comparison
+    can be made between fitnesses of different size, if the fitnesses are
+    equal until the extra elements, the longer fitness will be superior to the
+    shorter.
+
+    Different types of fitnesses are created in the :ref:`creating-types`
+    tutorial.
+
+    .. note::
+       When comparing fitness values that are **minimized**, ``a > b`` will
+       return :data:`True` if *a* is **smaller** than *b*.
+    """
+    
+    weights = None
+    """The weights are used in the fitness comparison. They are shared among
+    all fitnesses of the same type. When subclassing :class:`Fitness`, the
+    weights must be defined as a tuple where each element is associated to an
+    objective. A negative weight element corresponds to the minimization of
+    the associated objective and positive weight to the maximization.
+
+    .. note::
+        If weights is not defined during subclassing, the following error will 
+        occur at instantiation of a subclass fitness object: 
+        
+        ``TypeError: Can't instantiate abstract <class Fitness[...]> with
+        abstract attribute weights.``
+    """
+    
+    wvalues = ()
+    """Contains the weighted values of the fitness, the multiplication with the
+    weights is made when the values are set via the property :attr:`values`.
+    Multiplication is made on setting of the values for efficiency.
+    
+    Generally it is unnecessary to manipulate wvalues as it is an internal
+    attribute of the fitness used in the comparison operators.
+    """
+    
+    def __init__(self, values=()):
+        if self.weights is None:
+            raise TypeError("Can't instantiate abstract %r with abstract "
+                "attribute weights." % (self.__class__))
+        
+        if not isinstance(self.weights, Sequence):
+            raise TypeError("Attribute weights of %r must be a sequence." 
+                % self.__class__)
+        
+        if len(values) > 0:
+            self.values = values
+        
+    def getValues(self):
+        return tuple(map(truediv, self.wvalues, self.weights))
+            
+    def setValues(self, values):
+        try:
+            self.wvalues = tuple(map(mul, values, self.weights))
+        except TypeError:
+            _, _, traceback = sys.exc_info()
+            raise TypeError, ("Both weights and assigned values must be a "
+            "sequence of numbers when assigning to values of "
+            "%r. Currently assigning value(s) %r of %r to a fitness with "
+            "weights %s."
+            % (self.__class__, values, type(values), self.weights)), traceback
+            
+    def delValues(self):
+        self.wvalues = ()
+
+    values = property(getValues, setValues, delValues,
+        ("Fitness values. Use directly ``individual.fitness.values = values`` "
+         "in order to set the fitness and ``del individual.fitness.values`` "
+         "in order to clear (invalidate) the fitness. The (unweighted) fitness "
+         "can be directly accessed via ``individual.fitness.values``."))
+    
+    def dominates(self, other, obj=slice(None)):
+        """Return true if each objective of *self* is not strictly worse than 
+        the corresponding objective of *other* and at least one objective is 
+        strictly better.
+
+        :param obj: Slice indicating on which objectives the domination is 
+                    tested. The default value is `slice(None)`, representing
+                    every objectives.
+        """
+        not_equal = False
+        for self_wvalue, other_wvalue in zip(self.wvalues[obj], other.wvalues[obj]):
+            if self_wvalue > other_wvalue:
+                not_equal = True
+            elif self_wvalue < other_wvalue:
+                return False                
+        return not_equal
+
+    @property
+    def valid(self):
+        """Assess if a fitness is valid or not."""
+        return len(self.wvalues) != 0
+        
+    def __hash__(self):
+        return hash(self.wvalues)
+
+    def __gt__(self, other):
+        return not self.__le__(other)
+        
+    def __ge__(self, other):
+        return not self.__lt__(other)
+
+    def __le__(self, other):
+        return self.wvalues <= other.wvalues
+
+    def __lt__(self, other):
+        return self.wvalues < other.wvalues
+
+    def __eq__(self, other):
+        return self.wvalues == other.wvalues
+    
+    def __ne__(self, other):
+        return not self.__eq__(other)
+
+    def __deepcopy__(self, memo):
+        """Replace the basic deepcopy function with a faster one.
+        
+        It assumes that the elements in the :attr:`values` tuple are 
+        immutable and the fitness does not contain any other object 
+        than :attr:`values` and :attr:`weights`.
+        """
+        copy_ = self.__class__()
+        copy_.wvalues = self.wvalues
+        return copy_
+
+    def __str__(self):
+        """Return the values of the Fitness object."""
+        return str(self.values if self.valid else tuple())
+
+    def __repr__(self):
+        """Return the Python code to build a copy of the object."""
+        return "%s.%s(%r)" % (self.__module__, self.__class__.__name__,
+            self.values if self.valid else tuple())
+

python/isaac/autotuning/external/deap/benchmarks/__init__.py

@@ -0,0 +1,620 @@
+#    This file is part of DEAP.
+#
+#    DEAP is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU Lesser General Public License as
+#    published by the Free Software Foundation, either version 3 of
+#    the License, or (at your option) any later version.
+#
+#    DEAP is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#    GNU Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.
+"""
+Regroup typical EC benchmarks functions to import easily and benchmark
+examples.
+"""
+
+import random
+from math import sin, cos, pi, exp, e, sqrt
+from operator import mul
+from functools import reduce
+
+# Unimodal
+def rand(individual):
+    """Random test objective function.
+
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Type
+         - minimization or maximization
+       * - Range
+         - none
+       * - Global optima
+         - none
+       * - Function
+         - :math:`f(\mathbf{x}) = \\text{\\texttt{random}}(0,1)`
+    """
+    return random.random(),
+    
+def plane(individual):
+    """Plane test objective function.
+
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Type
+         - minimization
+       * - Range
+         - none
+       * - Global optima
+         - :math:`x_i = 0, \\forall i \in \\lbrace 1 \\ldots N\\rbrace`, :math:`f(\mathbf{x}) = 0`
+       * - Function
+         - :math:`f(\mathbf{x}) = x_0`
+    """
+    return individual[0],
+
+def sphere(individual):
+    """Sphere test objective function.
+
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Type
+         - minimization
+       * - Range
+         - none
+       * - Global optima
+         - :math:`x_i = 0, \\forall i \in \\lbrace 1 \\ldots N\\rbrace`, :math:`f(\mathbf{x}) = 0`
+       * - Function
+         - :math:`f(\mathbf{x}) = \sum_{i=1}^Nx_i^2`
+    """
+    return sum(gene * gene for gene in individual),
+
+def cigar(individual):
+    """Cigar test objective function.
+    
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Type
+         - minimization
+       * - Range
+         - none
+       * - Global optima
+         - :math:`x_i = 0, \\forall i \in \\lbrace 1 \\ldots N\\rbrace`, :math:`f(\mathbf{x}) = 0`
+       * - Function
+         - :math:`f(\mathbf{x}) = x_0^2 + 10^6\\sum_{i=1}^N\,x_i^2`
+    """
+    return individual[0]**2 + 1e6 * sum(gene * gene for gene in individual),
+
+def rosenbrock(individual):  
+    """Rosenbrock test objective function.
+
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Type
+         - minimization
+       * - Range
+         - none
+       * - Global optima
+         - :math:`x_i = 1, \\forall i \in \\lbrace 1 \\ldots N\\rbrace`, :math:`f(\mathbf{x}) = 0`
+       * - Function
+         - :math:`f(\\mathbf{x}) = \\sum_{i=1}^{N-1} (1-x_i)^2 + 100 (x_{i+1} - x_i^2 )^2`
+    
+    .. plot:: code/benchmarks/rosenbrock.py
+       :width: 67 %
+    """
+    return sum(100 * (x * x - y)**2 + (1. - x)**2 \
+                   for x, y in zip(individual[:-1], individual[1:])),
+
+def h1(individual):
+    """ Simple two-dimensional function containing several local maxima.
+    From: The Merits of a Parallel Genetic Algorithm in Solving Hard 
+    Optimization Problems, A. J. Knoek van Soest and L. J. R. Richard 
+    Casius, J. Biomech. Eng. 125, 141 (2003)
+
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Type
+         - maximization
+       * - Range
+         - :math:`x_i \in [-100, 100]`
+       * - Global optima
+         - :math:`\mathbf{x} = (8.6998, 6.7665)`, :math:`f(\mathbf{x}) = 2`\n
+       * - Function
+         - :math:`f(\mathbf{x}) = \\frac{\sin(x_1 - \\frac{x_2}{8})^2 + \
+            \\sin(x_2 + \\frac{x_1}{8})^2}{\\sqrt{(x_1 - 8.6998)^2 + \
+            (x_2 - 6.7665)^2} + 1}`
+
+    .. plot:: code/benchmarks/h1.py
+       :width: 67 %
+    """
+    num = (sin(individual[0] - individual[1] / 8))**2 + (sin(individual[1] + individual[0] / 8))**2
+    denum = ((individual[0] - 8.6998)**2 + (individual[1] - 6.7665)**2)**0.5 + 1
+    return num / denum,
+
+ # 
+# Multimodal
+def ackley(individual):
+    """Ackley test objective function.
+
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Type
+         - minimization
+       * - Range
+         - :math:`x_i \in [-15, 30]`
+       * - Global optima
+         - :math:`x_i = 0, \\forall i \in \\lbrace 1 \\ldots N\\rbrace`, :math:`f(\mathbf{x}) = 0`
+       * - Function
+         - :math:`f(\\mathbf{x}) = 20 - 20\exp\left(-0.2\sqrt{\\frac{1}{N} \
+            \\sum_{i=1}^N x_i^2} \\right) + e - \\exp\\left(\\frac{1}{N}\sum_{i=1}^N \\cos(2\pi x_i) \\right)`
+    
+    .. plot:: code/benchmarks/ackley.py
+       :width: 67 %
+    """
+    N = len(individual)
+    return 20 - 20 * exp(-0.2*sqrt(1.0/N * sum(x**2 for x in individual))) \
+            + e - exp(1.0/N * sum(cos(2*pi*x) for x in individual)),
+            
+def bohachevsky(individual):
+    """Bohachevsky test objective function.
+
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Type
+         - minimization
+       * - Range
+         - :math:`x_i \in [-100, 100]`
+       * - Global optima
+         - :math:`x_i = 0, \\forall i \in \\lbrace 1 \\ldots N\\rbrace`, :math:`f(\mathbf{x}) = 0`
+       * - Function
+         -  :math:`f(\mathbf{x}) = \sum_{i=1}^{N-1}(x_i^2 + 2x_{i+1}^2 - \
+                   0.3\cos(3\pi x_i) - 0.4\cos(4\pi x_{i+1}) + 0.7)`
+    
+    .. plot:: code/benchmarks/bohachevsky.py
+       :width: 67 %
+    """
+    return sum(x**2 + 2*x1**2 - 0.3*cos(3*pi*x) - 0.4*cos(4*pi*x1) + 0.7 
+                for x, x1 in zip(individual[:-1], individual[1:])),
+
+def griewank(individual):
+    """Griewank test objective function.
+    
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Type
+         - minimization
+       * - Range
+         - :math:`x_i \in [-600, 600]`
+       * - Global optima
+         - :math:`x_i = 0, \\forall i \in \\lbrace 1 \\ldots N\\rbrace`, :math:`f(\mathbf{x}) = 0`
+       * - Function
+         - :math:`f(\\mathbf{x}) = \\frac{1}{4000}\\sum_{i=1}^N\,x_i^2 - \
+                  \prod_{i=1}^N\\cos\\left(\\frac{x_i}{\sqrt{i}}\\right) + 1`
+
+    .. plot:: code/benchmarks/griewank.py
+       :width: 67 %
+    """
+    return 1.0/4000.0 * sum(x**2 for x in individual) - \
+        reduce(mul, (cos(x/sqrt(i+1.0)) for i, x in enumerate(individual)), 1) + 1,
+            
+def rastrigin(individual):
+    """Rastrigin test objective function.
+
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Type
+         - minimization
+       * - Range
+         - :math:`x_i \in [-5.12, 5.12]`
+       * - Global optima
+         - :math:`x_i = 0, \\forall i \in \\lbrace 1 \\ldots N\\rbrace`, :math:`f(\mathbf{x}) = 0`
+       * - Function
+         - :math:`f(\\mathbf{x}) = 10N \sum_{i=1}^N x_i^2 - 10 \\cos(2\\pi x_i)`
+
+    .. plot:: code/benchmarks/rastrigin.py
+       :width: 67 %
+    """     
+    return 10 * len(individual) + sum(gene * gene - 10 * \
+                        cos(2 * pi * gene) for gene in individual),
+
+def rastrigin_scaled(individual):
+    """Scaled Rastrigin test objective function.
+    
+    :math:`f_{\\text{RastScaled}}(\mathbf{x}) = 10N + \sum_{i=1}^N \
+        \left(10^{\left(\\frac{i-1}{N-1}\\right)} x_i \\right)^2 x_i)^2 - \
+        10\cos\\left(2\\pi 10^{\left(\\frac{i-1}{N-1}\\right)} x_i \\right)`
+    """
+    N = len(individual)
+    return 10*N + sum((10**(i/(N-1))*x)**2 - 
+                      10*cos(2*pi*10**(i/(N-1))*x) for i, x in enumerate(individual)),
+
+def rastrigin_skew(individual):
+    """Skewed Rastrigin test objective function.
+    
+     :math:`f_{\\text{RastSkew}}(\mathbf{x}) = 10N \sum_{i=1}^N \left(y_i^2 - 10 \\cos(2\\pi x_i)\\right)`
+        
+     :math:`\\text{with } y_i = \
+                            \\begin{cases} \
+                                10\\cdot x_i & \\text{ if } x_i > 0,\\\ \
+                                x_i & \\text{ otherwise } \
+                            \\end{cases}`
+    """
+    N = len(individual)
+    return 10*N + sum((10*x if x > 0 else x)**2 
+                    - 10*cos(2*pi*(10*x if x > 0 else x)) for x in individual),
+def schaffer(individual):
+    """Schaffer test objective function.
+    
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Type
+         - minimization
+       * - Range
+         - :math:`x_i \in [-100, 100]`
+       * - Global optima
+         - :math:`x_i = 0, \\forall i \in \\lbrace 1 \\ldots N\\rbrace`, :math:`f(\mathbf{x}) = 0`
+       * - Function
+         -  :math:`f(\mathbf{x}) = \sum_{i=1}^{N-1} (x_i^2+x_{i+1}^2)^{0.25} \cdot \
+                  \\left[ \sin^2(50\cdot(x_i^2+x_{i+1}^2)^{0.10}) + 1.0 \
+                  \\right]`
+
+    .. plot:: code/benchmarks/schaffer.py
+        :width: 67 %
+    """
+    return sum((x**2+x1**2)**0.25 * ((sin(50*(x**2+x1**2)**0.1))**2+1.0) 
+                for x, x1 in zip(individual[:-1], individual[1:])),
+
+def schwefel(individual):
+    """Schwefel test objective function.
+
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Type
+         - minimization
+       * - Range
+         - :math:`x_i \in [-500, 500]`
+       * - Global optima
+         - :math:`x_i = 420.96874636, \\forall i \in \\lbrace 1 \\ldots N\\rbrace`, :math:`f(\mathbf{x}) = 0`
+       * - Function
+         - :math:`f(\mathbf{x}) = 418.9828872724339\cdot N - \
+            \sum_{i=1}^N\,x_i\sin\\left(\sqrt{|x_i|}\\right)`
+
+
+    .. plot:: code/benchmarks/schwefel.py
+        :width: 67 %
+    """    
+    N = len(individual)
+    return 418.9828872724339*N-sum(x*sin(sqrt(abs(x))) for x in individual),
+
+def himmelblau(individual):
+    """The Himmelblau's function is multimodal with 4 defined minimums in 
+    :math:`[-6, 6]^2`.
+
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Type
+         - minimization
+       * - Range
+         - :math:`x_i \in [-6, 6]`
+       * - Global optima
+         - :math:`\mathbf{x}_1 = (3.0, 2.0)`, :math:`f(\mathbf{x}_1) = 0`\n
+           :math:`\mathbf{x}_2 = (-2.805118, 3.131312)`, :math:`f(\mathbf{x}_2) = 0`\n
+           :math:`\mathbf{x}_3 = (-3.779310, -3.283186)`, :math:`f(\mathbf{x}_3) = 0`\n
+           :math:`\mathbf{x}_4 = (3.584428, -1.848126)`, :math:`f(\mathbf{x}_4) = 0`\n
+       * - Function
+         - :math:`f(x_1, x_2) = (x_1^2 + x_2 - 11)^2 + (x_1 + x_2^2 -7)^2`
+
+    .. plot:: code/benchmarks/himmelblau.py
+        :width: 67 %
+    """
+    return (individual[0] * individual[0] + individual[1] - 11)**2 + \
+        (individual[0] + individual[1] * individual[1] - 7)**2,
+
+def shekel(individual, a, c):
+    """The Shekel multimodal function can have any number of maxima. The number
+    of maxima is given by the length of any of the arguments *a* or *c*, *a*
+    is a matrix of size :math:`M\\times N`, where *M* is the number of maxima
+    and *N* the number of dimensions and *c* is a :math:`M\\times 1` vector.
+    The matrix :math:`\\mathcal{A}` can be seen as the position of the maxima
+    and the vector :math:`\\mathbf{c}`, the width of the maxima.
+    
+    :math:`f_\\text{Shekel}(\mathbf{x}) = \\sum_{i = 1}^{M} \\frac{1}{c_{i} + 
+    \\sum_{j = 1}^{N} (x_{j} - a_{ij})^2 }`
+    
+    The following figure uses
+    
+    :math:`\\mathcal{A} = \\begin{bmatrix} 0.5 & 0.5 \\\\ 0.25 & 0.25 \\\\ 
+    0.25 & 0.75 \\\\ 0.75 & 0.25 \\\\ 0.75 & 0.75 \\end{bmatrix}` and
+    :math:`\\mathbf{c} = \\begin{bmatrix} 0.002 \\\\ 0.005 \\\\ 0.005
+    \\\\ 0.005 \\\\ 0.005 \\end{bmatrix}`, thus defining 5 maximums in
+    :math:`\\mathbb{R}^2`.
+    
+    .. plot:: code/benchmarks/shekel.py
+        :width: 67 %
+    """
+    return sum((1. / (c[i] + sum((x - a[i][j])**2 for j, x in enumerate(individual)))) for i in range(len(c))),
+
+# Multiobjectives
+def kursawe(individual):
+    """Kursawe multiobjective function.
+    
+    :math:`f_{\\text{Kursawe}1}(\\mathbf{x}) = \\sum_{i=1}^{N-1} -10 e^{-0.2 \\sqrt{x_i^2 + x_{i+1}^2} }`
+    
+    :math:`f_{\\text{Kursawe}2}(\\mathbf{x}) = \\sum_{i=1}^{N} |x_i|^{0.8} + 5 \\sin(x_i^3)`
+
+    .. plot:: code/benchmarks/kursawe.py
+       :width: 100 %
+    """
+    f1 = sum(-10 * exp(-0.2 * sqrt(x * x + y * y)) for x, y in zip(individual[:-1], individual[1:]))
+    f2 = sum(abs(x)**0.8 + 5 * sin(x * x * x) for x in individual)
+    return f1, f2
+
+
+def schaffer_mo(individual):
+    """Schaffer's multiobjective function on a one attribute *individual*.
+    From: J. D. Schaffer, "Multiple objective optimization with vector
+    evaluated genetic algorithms", in Proceedings of the First International
+    Conference on Genetic Algorithms, 1987. 
+    
+    :math:`f_{\\text{Schaffer}1}(\\mathbf{x}) = x_1^2`
+    
+    :math:`f_{\\text{Schaffer}2}(\\mathbf{x}) = (x_1-2)^2`
+    """
+    return individual[0] ** 2, (individual[0] - 2) ** 2
+    
+def zdt1(individual):
+    """ZDT1 multiobjective function.
+    
+    :math:`g(\\mathbf{x}) = 1 + \\frac{9}{n-1}\\sum_{i=2}^n x_i`
+    
+    :math:`f_{\\text{ZDT1}1}(\\mathbf{x}) = x_1`
+    
+    :math:`f_{\\text{ZDT1}2}(\\mathbf{x}) = g(\\mathbf{x})\\left[1 - \\sqrt{\\frac{x_1}{g(\\mathbf{x})}}\\right]`
+    """
+    g  = 1.0 + 9.0*sum(individual[1:])/(len(individual)-1)
+    f1 = individual[0]
+    f2 = g * (1 - sqrt(f1/g))
+    return f1, f2
+
+def zdt2(individual):
+    """ZDT2 multiobjective function.
+    
+    :math:`g(\\mathbf{x}) = 1 + \\frac{9}{n-1}\\sum_{i=2}^n x_i`
+    
+    :math:`f_{\\text{ZDT2}1}(\\mathbf{x}) = x_1`
+    
+    :math:`f_{\\text{ZDT2}2}(\\mathbf{x}) = g(\\mathbf{x})\\left[1 - \\left(\\frac{x_1}{g(\\mathbf{x})}\\right)^2\\right]`
+    
+    """
+
+    g  = 1.0 + 9.0*sum(individual[1:])/(len(individual)-1)
+    f1 = individual[0]
+    f2 = g * (1 - (f1/g)**2)
+    return f1, f2
+    
+def zdt3(individual):
+    """ZDT3 multiobjective function.
+
+    :math:`g(\\mathbf{x}) = 1 + \\frac{9}{n-1}\\sum_{i=2}^n x_i`
+
+    :math:`f_{\\text{ZDT3}1}(\\mathbf{x}) = x_1`
+
+    :math:`f_{\\text{ZDT3}2}(\\mathbf{x}) = g(\\mathbf{x})\\left[1 - \\sqrt{\\frac{x_1}{g(\\mathbf{x})}} - \\frac{x_1}{g(\\mathbf{x})}\\sin(10\\pi x_1)\\right]`
+
+    """
+
+    g  = 1.0 + 9.0*sum(individual[1:])/(len(individual)-1)
+    f1 = individual[0]
+    f2 = g * (1 - sqrt(f1/g) - f1/g * sin(10*pi*f1))
+    return f1, f2
+
+def zdt4(individual):
+    """ZDT4 multiobjective function.
+    
+    :math:`g(\\mathbf{x}) = 1 + 10(n-1) + \\sum_{i=2}^n \\left[ x_i^2 - 10\\cos(4\\pi x_i) \\right]`
+
+    :math:`f_{\\text{ZDT4}1}(\\mathbf{x}) = x_1`
+    
+    :math:`f_{\\text{ZDT4}2}(\\mathbf{x}) = g(\\mathbf{x})\\left[ 1 - \\sqrt{x_1/g(\\mathbf{x})} \\right]`
+    
+    """
+    g  = 1 + 10*(len(individual)-1) + sum(xi**2 - 10*cos(4*pi*xi) for xi in individual[1:])
+    f1 = individual[0]
+    f2 = g * (1 - sqrt(f1/g))
+    return f1, f2
+    
+def zdt6(individual):
+    """ZDT6 multiobjective function.
+    
+    :math:`g(\\mathbf{x}) = 1 + 9 \\left[ \\left(\\sum_{i=2}^n x_i\\right)/(n-1) \\right]^{0.25}`
+    
+    :math:`f_{\\text{ZDT6}1}(\\mathbf{x}) = 1 - \\exp(-4x_1)\\sin^6(6\\pi x_1)`
+    
+    :math:`f_{\\text{ZDT6}2}(\\mathbf{x}) = g(\\mathbf{x}) \left[ 1 - (f_{\\text{ZDT6}1}(\\mathbf{x})/g(\\mathbf{x}))^2 \\right]`
+    
+    """
+    g  = 1 + 9 * (sum(individual[1:]) / (len(individual)-1))**0.25
+    f1 = 1 - exp(-4*individual[0]) * sin(6*pi*individual[0])**6
+    f2 = g * (1 - (f1/g)**2)
+    return f1, f2
+
+def dtlz1(individual, obj):
+    """DTLZ1 mutliobjective function. It returns a tuple of *obj* values. 
+    The individual must have at least *obj* elements.
+    From: K. Deb, L. Thiele, M. Laumanns and E. Zitzler. Scalable Multi-Objective 
+    Optimization Test Problems. CEC 2002, p. 825 - 830, IEEE Press, 2002.
+
+    :math:`g(\\mathbf{x}_m) = 100\\left(|\\mathbf{x}_m| + \sum_{x_i \in \\mathbf{x}_m}\\left((x_i - 0.5)^2 - \\cos(20\pi(x_i - 0.5))\\right)\\right)`
+
+    :math:`f_{\\text{DTLZ1}1}(\\mathbf{x}) = \\frac{1}{2} (1 + g(\\mathbf{x}_m)) \\prod_{i=1}^{m-1}x_i`
+    
+    :math:`f_{\\text{DTLZ1}2}(\\mathbf{x}) = \\frac{1}{2} (1 + g(\\mathbf{x}_m)) (1-x_{m-1}) \\prod_{i=1}^{m-2}x_i`
+    
+    :math:`\\ldots`
+    
+    :math:`f_{\\text{DTLZ1}m-1}(\\mathbf{x}) = \\frac{1}{2} (1 + g(\\mathbf{x}_m)) (1 - x_2) x_1`
+    
+    :math:`f_{\\text{DTLZ1}m}(\\mathbf{x}) = \\frac{1}{2} (1 - x_1)(1 + g(\\mathbf{x}_m))`
+    
+    Where :math:`m` is the number of objectives and :math:`\\mathbf{x}_m` is a
+    vector of the remaining attributes :math:`[x_m~\\ldots~x_n]` of the
+    individual in :math:`n > m` dimensions.
+    
+    """
+    g = 100 * (len(individual[obj-1:]) + sum((xi-0.5)**2 - cos(20*pi*(xi-0.5)) for xi in individual[obj-1:]))
+    f = [0.5 * reduce(mul, individual[:obj-1], 1) * (1 + g)]
+    f.extend(0.5 * reduce(mul, individual[:m], 1) * (1 - individual[m]) * (1 + g) for m in reversed(xrange(obj-1)))
+    return f
+
+def dtlz2(individual, obj):
+    """DTLZ2 mutliobjective function. It returns a tuple of *obj* values. 
+    The individual must have at least *obj* elements.
+    From: K. Deb, L. Thiele, M. Laumanns and E. Zitzler. Scalable Multi-Objective 
+    Optimization Test Problems. CEC 2002, p. 825 - 830, IEEE Press, 2002.
+    
+    :math:`g(\\mathbf{x}_m) = \\sum_{x_i \in \\mathbf{x}_m} (x_i - 0.5)^2`
+    
+    :math:`f_{\\text{DTLZ2}1}(\\mathbf{x}) = (1 + g(\\mathbf{x}_m)) \\prod_{i=1}^{m-1} \\cos(0.5x_i\pi)`
+    
+    :math:`f_{\\text{DTLZ2}2}(\\mathbf{x}) = (1 + g(\\mathbf{x}_m)) \\sin(0.5x_{m-1}\pi ) \\prod_{i=1}^{m-2} \\cos(0.5x_i\pi)`
+    
+    :math:`\\ldots`
+    
+    :math:`f_{\\text{DTLZ2}m}(\\mathbf{x}) = (1 + g(\\mathbf{x}_m)) \\sin(0.5x_{1}\pi )`
+    
+    Where :math:`m` is the number of objectives and :math:`\\mathbf{x}_m` is a
+    vector of the remaining attributes :math:`[x_m~\\ldots~x_n]` of the
+    individual in :math:`n > m` dimensions.
+    """
+    xc = individual[:obj-1]
+    xm = individual[obj-1:]
+    g = sum((xi-0.5)**2 for xi in xm)
+    f = [(1.0+g) *  reduce(mul, (cos(0.5*xi*pi) for xi in xc), 1.0)]
+    f.extend((1.0+g) * reduce(mul, (cos(0.5*xi*pi) for xi in xc[:m]), 1) * sin(0.5*xc[m]*pi) for m in range(obj-2, -1, -1))
+
+    return f
+
+def dtlz3(individual, obj):
+    """DTLZ3 mutliobjective function. It returns a tuple of *obj* values. 
+    The individual must have at least *obj* elements.
+    From: K. Deb, L. Thiele, M. Laumanns and E. Zitzler. Scalable Multi-Objective 
+    Optimization Test Problems. CEC 2002, p. 825 - 830, IEEE Press, 2002.
+    
+    :math:`g(\\mathbf{x}_m) = 100\\left(|\\mathbf{x}_m| + \sum_{x_i \in \\mathbf{x}_m}\\left((x_i - 0.5)^2 - \\cos(20\pi(x_i - 0.5))\\right)\\right)`
+    
+    :math:`f_{\\text{DTLZ3}1}(\\mathbf{x}) = (1 + g(\\mathbf{x}_m)) \\prod_{i=1}^{m-1} \\cos(0.5x_i\pi)`
+    
+    :math:`f_{\\text{DTLZ3}2}(\\mathbf{x}) = (1 + g(\\mathbf{x}_m)) \\sin(0.5x_{m-1}\pi ) \\prod_{i=1}^{m-2} \\cos(0.5x_i\pi)`
+    
+    :math:`\\ldots`
+    
+    :math:`f_{\\text{DTLZ3}m}(\\mathbf{x}) = (1 + g(\\mathbf{x}_m)) \\sin(0.5x_{1}\pi )`
+    
+    Where :math:`m` is the number of objectives and :math:`\\mathbf{x}_m` is a
+    vector of the remaining attributes :math:`[x_m~\\ldots~x_n]` of the
+    individual in :math:`n > m` dimensions.
+    """
+    xc = individual[:obj-1]
+    xm = individual[obj-1:]
+    g = 100 * (len(xm) + sum((xi-0.5)**2 - cos(20*pi*(xi-0.5)) for xi in xm))
+    f = [(1.0+g) *  reduce(mul, (cos(0.5*xi*pi) for xi in xc), 1.0)]
+    f.extend((1.0+g) * reduce(mul, (cos(0.5*xi*pi) for xi in xc[:m]), 1) * sin(0.5*xc[m]*pi) for m in range(obj-2, -1, -1))
+    return f
+
+def dtlz4(individual, obj, alpha):
+    """DTLZ4 mutliobjective function. It returns a tuple of *obj* values. The
+    individual must have at least *obj* elements. The *alpha* parameter allows
+    for a meta-variable mapping in :func:`dtlz2` :math:`x_i \\rightarrow
+    x_i^\\alpha`, the authors suggest :math:`\\alpha = 100`.
+    From: K. Deb, L. Thiele, M. Laumanns and E. Zitzler. Scalable Multi-Objective 
+    Optimization Test Problems. CEC 2002, p. 825 - 830, IEEE Press, 2002.
+    
+    :math:`g(\\mathbf{x}_m) = \\sum_{x_i \in \\mathbf{x}_m} (x_i - 0.5)^2`
+    
+    :math:`f_{\\text{DTLZ4}1}(\\mathbf{x}) = (1 + g(\\mathbf{x}_m)) \\prod_{i=1}^{m-1} \\cos(0.5x_i^\\alpha\pi)`
+    
+    :math:`f_{\\text{DTLZ4}2}(\\mathbf{x}) = (1 + g(\\mathbf{x}_m)) \\sin(0.5x_{m-1}^\\alpha\pi ) \\prod_{i=1}^{m-2} \\cos(0.5x_i^\\alpha\pi)`
+    
+    :math:`\\ldots`
+    
+    :math:`f_{\\text{DTLZ4}m}(\\mathbf{x}) = (1 + g(\\mathbf{x}_m)) \\sin(0.5x_{1}^\\alpha\pi )`
+    
+    Where :math:`m` is the number of objectives and :math:`\\mathbf{x}_m` is a
+    vector of the remaining attributes :math:`[x_m~\\ldots~x_n]` of the
+    individual in :math:`n > m` dimensions.
+    """
+    xc = individual[:obj-1]
+    xm = individual[obj-1:]
+    g = sum((xi-0.5)**2 for xi in xm)
+    f = [(1.0+g) *  reduce(mul, (cos(0.5*xi**alpha*pi) for xi in xc), 1.0)]
+    f.extend((1.0+g) * reduce(mul, (cos(0.5*xi**alpha*pi) for xi in xc[:m]), 1) * sin(0.5*xc[m]**alpha*pi) for m in range(obj-2, -1, -1))
+    return f
+
+def fonseca(individual):
+    """Fonseca and Fleming's multiobjective function.
+    From: C. M. Fonseca and P. J. Fleming, "Multiobjective optimization and
+    multiple constraint handling with evolutionary algorithms -- Part II:
+    Application example", IEEE Transactions on Systems, Man and Cybernetics,
+    1998.
+    
+    :math:`f_{\\text{Fonseca}1}(\\mathbf{x}) = 1 - e^{-\\sum_{i=1}^{3}(x_i - \\frac{1}{\\sqrt{3}})^2}`
+    
+    :math:`f_{\\text{Fonseca}2}(\\mathbf{x}) = 1 - e^{-\\sum_{i=1}^{3}(x_i + \\frac{1}{\\sqrt{3}})^2}`
+    """
+    f_1 = 1 - exp(-sum((xi - 1/sqrt(3))**2 for xi in individual[:3]))
+    f_2 = 1 - exp(-sum((xi + 1/sqrt(3))**2 for xi in individual[:3]))
+    return f_1, f_2
+
+def poloni(individual):
+    """Poloni's multiobjective function on a two attribute *individual*. From:
+    C. Poloni, "Hybrid GA for multi objective aerodynamic shape optimization",
+    in Genetic Algorithms in Engineering and Computer Science, 1997.
+    
+    :math:`A_1 = 0.5 \\sin (1) - 2 \\cos (1) + \\sin (2) - 1.5 \\cos (2)`
+
+    :math:`A_2 = 1.5 \\sin (1) - \\cos (1) + 2 \\sin (2) - 0.5 \\cos (2)`
+
+    :math:`B_1 = 0.5 \\sin (x_1) - 2 \\cos (x_1) + \\sin (x_2) - 1.5 \\cos (x_2)`
+
+    :math:`B_2 = 1.5 \\sin (x_1) - cos(x_1) + 2 \\sin (x_2) - 0.5 \\cos (x_2)`
+    
+    :math:`f_{\\text{Poloni}1}(\\mathbf{x}) = 1 + (A_1 - B_1)^2 + (A_2 - B_2)^2`
+    
+    :math:`f_{\\text{Poloni}2}(\\mathbf{x}) = (x_1 + 3)^2 + (x_2 + 1)^2`
+    """
+    x_1 = individual[0]
+    x_2 = individual[1]
+    A_1 = 0.5 * sin(1) - 2 * cos(1) + sin(2) - 1.5 * cos(2)
+    A_2 = 1.5 * sin(1) - cos(1) + 2 * sin(2) - 0.5 * cos(2)
+    B_1 = 0.5 * sin(x_1) - 2 * cos(x_1) + sin(x_2) - 1.5 * cos(x_2)
+    B_2 = 1.5 * sin(x_1) - cos(x_1) + 2 * sin(x_2) - 0.5 * cos(x_2)
+    return 1 + (A_1 - B_1)**2 + (A_2 - B_2)**2, (x_1 + 3)**2 + (x_2 + 1)**2
+

python/isaac/autotuning/external/deap/benchmarks/binary.py

@@ -0,0 +1,134 @@
+#    This file is part of DEAP.
+#
+#    DEAP is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU Lesser General Public License as
+#    published by the Free Software Foundation, either version 3 of
+#    the License, or (at your option) any later version.
+#
+#    DEAP is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#    GNU Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.
+
+import math
+
+def bin2float(min_, max_, nbits):
+    """Convert a binary array into an array of float where each
+    float is composed of *nbits* and is between *min_* and *max_*
+    and return the result of the decorated function.
+
+    .. note::
+        This decorator requires the first argument of
+        the evaluation function to be named *individual*.
+    """
+    def wrap(function):
+        def wrapped_function(individual, *args, **kargs):
+            nelem = len(individual)/nbits
+            decoded = [0] * nelem
+            for i in xrange(nelem):
+                gene = int("".join(map(str, individual[i*nbits:i*nbits+nbits])), 2)
+                div = 2**nbits - 1
+                temp = float(gene)/float(div)
+                decoded[i] = min_ + (temp * (max_ - min_))
+            return function(decoded, *args, **kargs)
+        return wrapped_function
+    return wrap
+
+def trap(individual):
+    u = sum(individual)
+    k = len(individual)
+    if u == k:
+        return k
+    else:
+        return k - 1 - u
+
+def inv_trap(individual):
+    u = sum(individual)
+    k = len(individual)
+    if u == 0:
+        return k
+    else:
+        return u - 1
+
+def chuang_f1(individual):
+    """Binary deceptive function from : Multivariate Multi-Model Approach for
+    Globally Multimodal Problems by Chung-Yao Chuang and Wen-Lian Hsu.
+    
+    The function takes individual of 40+1 dimensions and has two global optima
+    in [1,1,...,1] and [0,0,...,0].
+    """    
+    total = 0
+    if individual[-1] == 0:
+        for i in xrange(0,len(individual)-1,4):
+            total += inv_trap(individual[i:i+4])
+    else:
+        for i in xrange(0,len(individual)-1,4):
+            total += trap(individual[i:i+4])
+    return total,
+
+def chuang_f2(individual):
+    """Binary deceptive function from : Multivariate Multi-Model Approach for
+    Globally Multimodal Problems by Chung-Yao Chuang and Wen-Lian Hsu.
+    
+    The function takes individual of 40+1 dimensions and has four global optima
+    in [1,1,...,0,0], [0,0,...,1,1], [1,1,...,1] and [0,0,...,0].    
+    """    
+    total = 0
+    if individual[-2] == 0 and individual[-1] == 0:
+        for i in xrange(0,len(individual)-2,8):
+            total += inv_trap(individual[i:i+4]) + inv_trap(individual[i+4:i+8])
+    elif individual[-2] == 0 and individual[-1] == 1:
+        for i in xrange(0,len(individual)-2,8):
+            total += inv_trap(individual[i:i+4]) + trap(individual[i+4:i+8])
+    elif individual[-2] == 1 and individual[-1] == 0:
+        for i in xrange(0,len(individual)-2,8):
+            total += trap(individual[i:i+4]) + inv_trap(individual[i+4:i+8])
+    else:
+        for i in xrange(0,len(individual)-2,8):
+            total += trap(individual[i:i+4]) + trap(individual[i+4:i+8])
+    return total,
+
+def chuang_f3(individual):
+    """Binary deceptive function from : Multivariate Multi-Model Approach for
+    Globally Multimodal Problems by Chung-Yao Chuang and Wen-Lian Hsu.
+
+    The function takes individual of 40+1 dimensions and has two global optima
+    in [1,1,...,1] and [0,0,...,0].
+    """
+    total = 0
+    if individual[-1] == 0:
+        for i in xrange(0,len(individual)-1,4):
+            total += inv_trap(individual[i:i+4])
+    else:
+        for i in xrange(2,len(individual)-3,4):
+            total += inv_trap(individual[i:i+4])
+        total += trap(individual[-2:]+individual[:2])
+    return total,
+
+# Royal Road Functions
+def royal_road1(individual, order):
+    """Royal Road Function R1 as presented by Melanie Mitchell in : 
+    "An introduction to Genetic Algorithms".
+    """
+    nelem = len(individual) / order
+    max_value = int(2**order - 1)
+    total = 0
+    for i in xrange(nelem):
+        value = int("".join(map(str, individual[i*order:i*order+order])), 2)
+        total += int(order) * int(value/max_value)
+    return total,
+
+def royal_road2(individual, order):
+    """Royal Road Function R2 as presented by Melanie Mitchell in : 
+    "An introduction to Genetic Algorithms".
+    """
+    total = 0
+    norder = order
+    while norder < order**2:
+        total += royal_road1(norder, individual)[0]
+        norder *= 2
+    return total,
+

python/isaac/autotuning/external/deap/benchmarks/gp.py

@@ -0,0 +1,128 @@
+#    This file is part of DEAP.
+#
+#    DEAP is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU Lesser General Public License as
+#    published by the Free Software Foundation, either version 3 of
+#    the License, or (at your option) any later version.
+#
+#    DEAP is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#    GNU Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.
+
+from math import exp, sin, cos
+
+def kotanchek(data):
+    """Kotanchek benchmark function.
+    
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+       
+       * - Range
+         - :math:`\mathbf{x} \in [-1, 7]^2`
+       * - Function
+         - :math:`f(\mathbf{x}) = \\frac{e^{-(x_1 - 1)^2}}{3.2 + (x_2 - 2.5)^2}`
+    """
+    return exp(-(data[0] - 1)**2) / (3.2 + (data[1] - 2.5)**2)
+
+def salustowicz_1d(data):
+    """Salustowicz benchmark function.
+    
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Range
+         - :math:`x \in [0, 10]`
+       * - Function
+         - :math:`f(x) = e^x x^3 \cos(x) \sin(x) (\cos(x) \sin^2(x) - 1)`
+    """
+    return exp(-data[0]) * data[0]**3 * cos(data[0]) * sin(data[0]) * (cos(data[0]) * sin(data[0])**2 - 1)
+    
+def salustowicz_2d(data):
+    """Salustowicz benchmark function.
+    
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Range
+         - :math:`\mathbf{x} \in [0, 7]^2`
+       * - Function
+         - :math:`f(\mathbf{x}) = e^{x_1} x_1^3 \cos(x_1) \sin(x_1) (\cos(x_1) \sin^2(x_1) - 1) (x_2 -5)`
+    """
+    return exp(-data[0]) * data[0]**3 * cos(data[0]) * sin(data[0]) * (cos(data[0]) * sin(data[0])**2 - 1) * (data[1] - 5)
+
+def unwrapped_ball(data):
+    """Unwrapped ball benchmark function.
+    
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Range
+         - :math:`\mathbf{x} \in [-2, 8]^n`
+       * - Function
+         - :math:`f(\mathbf{x}) = \\frac{10}{5 + \sum_{i=1}^n (x_i - 3)^2}`
+    """
+    return 10. / (5. + sum((d - 3)**2 for d in data))
+
+def rational_polynomial(data):
+    """Rational polynomial ball benchmark function.
+    
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Range
+         - :math:`\mathbf{x} \in [0, 2]^3`
+       * - Function
+         - :math:`f(\mathbf{x}) = \\frac{30 * (x_1 - 1) (x_3 - 1)}{x_2^2 (x_1 - 10)}`
+    """
+    return 30. * (data[0] - 1) * (data[2] - 1) / (data[1]**2 * (data[0] - 10))
+
+def sin_cos(data):
+    """Sine cosine benchmark function.
+    
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Range
+         - :math:`\mathbf{x} \in [0, 6]^2`
+       * - Function
+         - :math:`f(\mathbf{x}) = 6\sin(x_1)\cos(x_2)`
+    """
+    6 * sin(data[0]) * cos(data[1])
+
+def ripple(data):
+    """Ripple benchmark function.
+    
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Range
+         - :math:`\mathbf{x} \in [-5, 5]^2`
+       * - Function
+         - :math:`f(\mathbf{x}) = (x_1 - 3) (x_2 - 3) + 2 \sin((x_1 - 4) (x_2 -4))`
+    """
+    return (data[0] - 3) * (data[1] - 3) + 2 * sin((data[0] - 4) * (data[1] - 4))
+
+def rational_polynomial2(data):
+    """Rational polynomial benchmark function.
+    
+    .. list-table:: 
+       :widths: 10 50
+       :stub-columns: 1
+
+       * - Range
+         - :math:`\mathbf{x} \in [0, 6]^2`
+       * - Function
+         - :math:`f(\mathbf{x}) = \\frac{(x_1 - 3)^4 + (x_2 - 3)^3 - (x_2 - 3)}{(x_2 - 2)^4 + 10}`
+    """
+    return ((data[0] - 3)**4 + (data[1] - 3)**3 - (data[1] - 3)) / ((data[1] - 2)**4 + 10)

python/isaac/autotuning/external/deap/benchmarks/movingpeaks.py

@@ -0,0 +1,394 @@
+#    This file is part of DEAP.
+#
+#    DEAP is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU Lesser General Public License as
+#    published by the Free Software Foundation, either version 3 of
+#    the License, or (at your option) any later version.
+#
+#    DEAP is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#    GNU Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.
+
+"""
+Re-implementation of the `Moving Peaks Benchmark 
+<http://people.aifb.kit.edu/jbr/MovPeaks/>`_ by Jurgen Branke. With the
+addition of the fluctuating number of peaks presented in *du Plessis and
+Engelbrecht, 2013, Self-Adaptive Environment with Fluctuating Number of
+Optima.*
+"""
+
+import math
+import itertools
+import random
+from collections import Sequence
+
+def cone(individual, position, height, width):
+    """The cone peak function to be used with scenario 2 and 3.
+    
+    :math:`f(\mathbf{x}) = h - w \sqrt{\sum_{i=1}^N (x_i - p_i)^2}`
+    
+    """
+    value = 0.0
+    for x, p in zip(individual, position):
+        value += (x - p)**2
+    return height - width * math.sqrt(value)
+
+def sphere(individual, position, height, width):
+    value = 0.0
+    for x, p in zip(individual, position):
+        value += (x - p)**2
+    return height * value
+
+def function1(individual, position, height, width):
+    """The function1 peak function to be used with scenario 1.
+    
+    :math:`f(\mathbf{x}) = \\frac{h}{1 + w \sqrt{\sum_{i=1}^N (x_i - p_i)^2}}`
+    
+    """
+    value = 0.0
+    for x, p in zip(individual, position):
+        value += (x - p)**2
+    return height / (1 + width * value)
+
+class MovingPeaks:
+    """The Moving Peaks Benchmark is a fitness function changing over time. It
+    consists of a number of peaks, changing in height, width and location. The
+    peaks function is given by *pfunc*, wich is either a function object or a
+    list of function objects (the default is :func:`function1`). The number of
+    peaks is determined by *npeaks* (which defaults to 5). This parameter can
+    be either a integer or a sequence. If it is set to an integer the number
+    of peaks won't change, while if set to a sequence of 3 elements, the
+    number of peaks will fluctuate between the first and third element of that
+    sequence, the second element is the inital number of peaks. When
+    fluctuating the number of peaks, the parameter *number_severity* must be
+    included, it represents the number of peak fraction that is allowed to
+    change. The dimensionality of the search domain is *dim*. A basis function
+    *bfunc* can also be given to act as static landscape (the default is no
+    basis function). The argument *random* serves to grant an independent
+    random number generator to the moving peaks so that the evolution is not
+    influenced by number drawn by this object (the default uses random
+    functions from the Python module :mod:`random`). Various other keyword
+    parameters listed in the table below are required to setup the benchmark,
+    default parameters are based on scenario 1 of this benchmark.
+    
+    =================== ============================= =================== =================== ======================================================================================================================
+    Parameter           :data:`SCENARIO_1` (Default)  :data:`SCENARIO_2`  :data:`SCENARIO_3`    Details
+    =================== ============================= =================== =================== ======================================================================================================================
+    ``pfunc``           :func:`function1`             :func:`cone`        :func:`cone`        The peak function or a list of peak function.
+    ``npeaks``          5                             10                  50                  Number of peaks. If an integer, the number of peaks won't change, if a sequence it will fluctuate [min, current, max].
+    ``bfunc``           :obj:`None`                   :obj:`None`         ``lambda x: 10``    Basis static function.
+    ``min_coord``       0.0                           0.0                 0.0                 Minimum coordinate for the centre of the peaks.
+    ``max_coord``       100.0                         100.0               100.0               Maximum coordinate for the centre of the peaks.
+    ``min_height``      30.0                          30.0                30.0                Minimum height of the peaks.
+    ``max_height``      70.0                          70.0                70.0                Maximum height of the peaks.
+    ``uniform_height``  50.0                          50.0                0                   Starting height for all peaks, if ``uniform_height <= 0`` the initial height is set randomly for each peak.
+    ``min_width``       0.0001                        1.0                 1.0                 Minimum width of the peaks.
+    ``max_width``       0.2                           12.0                12.0                Maximum width of the peaks
+    ``uniform_width``   0.1                           0                   0                   Starting width for all peaks, if ``uniform_width <= 0`` the initial width is set randomly for each peak.
+    ``lambda_``         0.0                           0.5                 0.5                 Correlation between changes.
+    ``move_severity``   1.0                           1.5                 1.0                 The distance a single peak moves when peaks change.
+    ``height_severity`` 7.0                           7.0                 1.0                 The standard deviation of the change made to the height of a peak when peaks change.
+    ``width_severity``  0.01                          1.0                 0.5                 The standard deviation of the change made to the width of a peak when peaks change.
+    ``period``          5000                          5000                1000                Period between two changes.
+    =================== ============================= =================== =================== ======================================================================================================================
+    
+    Dictionnaries :data:`SCENARIO_1`, :data:`SCENARIO_2` and
+    :data:`SCENARIO_3` of this module define the defaults for these
+    parameters. The scenario 3 requires a constant basis function
+    which can be given as a lambda function ``lambda x: constant``.
+    
+    The following shows an example of scenario 1 with non uniform heights and
+    widths.
+    
+    .. plot:: code/benchmarks/movingsc1.py
+       :width: 67 %
+    """
+    def __init__(self, dim, random=random, **kargs):
+        # Scenario 1 is the default
+        sc = SCENARIO_1.copy()
+        sc.update(kargs)
+        
+        pfunc = sc.get("pfunc")
+        npeaks = sc.get("npeaks")
+        self.dim = dim
+
+        self.minpeaks, self.maxpeaks = None, None
+        if hasattr(npeaks, "__getitem__"):
+            self.minpeaks, npeaks, self.maxpeaks = npeaks
+            self.number_severity = sc.get("number_severity")
+            
+        try:
+            if len(pfunc) == npeaks:
+                self.peaks_function = pfunc
+            else:
+                self.peaks_function = self.random.sample(pfunc, npeaks)
+            self.pfunc_pool = tuple(pfunc)
+        except TypeError:
+            self.peaks_function = list(itertools.repeat(pfunc, npeaks))
+            self.pfunc_pool = (pfunc,)
+        
+        self.random = random
+        self.basis_function = sc.get("bfunc")
+        
+        self.min_coord = sc.get("min_coord")
+        self.max_coord = sc.get("max_coord")
+        
+        self.min_height = sc.get("min_height")
+        self.max_height = sc.get("max_height")
+        uniform_height = sc.get("uniform_height")
+        
+        self.min_width = sc.get("min_width")
+        self.max_width = sc.get("max_width")
+        uniform_width = sc.get("uniform_width")
+        
+        self.lambda_ = sc.get("lambda_")
+        self.move_severity = sc.get("move_severity")
+        self.height_severity = sc.get("height_severity")
+        self.width_severity = sc.get("width_severity")
+        
+        self.peaks_position = [[self.random.uniform(self.min_coord, self.max_coord) for _ in range(dim)] for _ in range(npeaks)]
+        
+        if uniform_height != 0:
+            self.peaks_height = [uniform_height for _ in range(npeaks)]
+        else:
+            self.peaks_height = [self.random.uniform(self.min_height, self.max_height) for _ in range(npeaks)]
+
+
+        if uniform_width != 0:
+            self.peaks_width = [uniform_width for _ in range(npeaks)]
+        else:
+            self.peaks_width = [self.random.uniform(self.min_width, self.max_width) for _ in range(npeaks)]
+        
+        self.last_change_vector = [[self.random.random() - 0.5 for _ in range(dim)] for _ in range(npeaks)]
+        
+        self.period = sc.get("period")
+
+        # Used by the Offline Error calculation
+        self._optimum = None
+        self._error = None
+        self._offline_error = 0
+
+        # Also used for auto change
+        self.nevals = 0
+
+    def globalMaximum(self):
+        """Returns the global maximum value and position."""
+        # The global maximum is at one peak's position
+        potential_max = list()
+        for func, pos, height, width in zip(self.peaks_function,
+                                            self.peaks_position,
+                                            self.peaks_height,
+                                            self.peaks_width):
+            potential_max.append((func(pos, pos, height, width), pos))
+        return max(potential_max)
+    
+    def maximums(self):
+        """Returns all visible maximums value and position sorted with the
+        global maximum first.
+        """
+        # The maximums are at the peaks position but might be swallowed by 
+        # other peaks
+        maximums = list()
+        for func, pos, height, width in zip(self.peaks_function,
+                                            self.peaks_position,
+                                            self.peaks_height,
+                                            self.peaks_width):
+            val = func(pos, pos, height, width)
+            if val >= self.__call__(pos, count=False)[0]:
+                maximums.append((val, pos))
+        return sorted(maximums, reverse=True)
+    
+    def __call__(self, individual, count=True):
+        """Evaluate a given *individual* with the current benchmark
+        configuration.
+
+        :param indidivudal: The individual to evaluate.
+        :param count: Wether or not to count this evaluation in
+                      the total evaluation count. (Defaults to
+                      :data:`True`)
+        """
+        possible_values = []
+        
+        for func, pos, height, width in zip(self.peaks_function,
+                                            self.peaks_position,
+                                            self.peaks_height,
+                                            self.peaks_width):
+            possible_values.append(func(individual, pos, height, width))
+        
+        if self.basis_function:
+            possible_values.append(self.basis_function(individual))
+
+        fitness = max(possible_values)
+
+        if count:
+            # Compute the offline error
+            self.nevals += 1
+            if self._optimum is None:
+                self._optimum = self.globalMaximum()[0]
+                self._error = abs(fitness - self._optimum)
+            self._error = min(self._error, abs(fitness - self._optimum))
+            self._offline_error += self._error
+
+            # We exausted the number of evaluation, change peaks for the next one.
+            if self.period > 0 and self.nevals % self.period == 0:
+                self.changePeaks()
+        
+        return fitness,
+    
+    def offlineError(self):
+        return self._offline_error / self.nevals
+
+    def currentError(self):
+        return self._error
+
+    def changePeaks(self):
+        """Order the peaks to change position, height, width and number."""
+        # Change the number of peaks
+        if self.minpeaks is not None and self.maxpeaks is not None:
+            npeaks = len(self.peaks_function)
+            u = self.random.random()
+            r = self.maxpeaks - self.minpeaks
+            if u < 0.5:
+                # Remove n peaks or less depending on the minimum number of peaks
+                u = self.random.random()
+                n = min(npeaks - self.minpeaks, int(round(r * u * self.number_severity)))
+                for i in range(n):
+                    idx = self.random.randrange(len(self.peaks_function))
+                    self.peaks_function.pop(idx)
+                    self.peaks_position.pop(idx)
+                    self.peaks_height.pop(idx)
+                    self.peaks_width.pop(idx)
+                    self.last_change_vector.pop(idx)
+            else:
+                # Add n peaks or less depending on the maximum number of peaks
+                u = self.random.random()
+                n = min(self.maxpeaks - npeaks, int(round(r * u * self.number_severity)))
+                for i in range(n):
+                    self.peaks_function.append(self.random.choice(self.pfunc_pool))
+                    self.peaks_position.append([self.random.uniform(self.min_coord, self.max_coord) for _ in range(self.dim)])
+                    self.peaks_height.append(self.random.uniform(self.min_height, self.max_height))
+                    self.peaks_width.append(self.random.uniform(self.min_width, self.max_width))
+                    self.last_change_vector.append([self.random.random() - 0.5 for _ in range(self.dim)])
+
+        for i in range(len(self.peaks_function)):
+            # Change peak position
+            shift = [self.random.random() - 0.5 for _ in range(len(self.peaks_position[i]))]
+            shift_length = sum(s**2 for s in shift)
+            shift_length = self.move_severity / math.sqrt(shift_length) if shift_length > 0 else 0
+            
+            shift = [shift_length * (1.0 - self.lambda_) * s \
+                     + self.lambda_ * c for s, c in zip(shift, self.last_change_vector[i])]
+            
+            shift_length = sum(s**2 for s in shift)
+            shift_length = self.move_severity / math.sqrt(shift_length) if shift_length > 0 else 0
+
+            shift = [s*shift_length for s in shift]
+            
+            new_position = []
+            final_shift = []
+            for pp, s in zip(self.peaks_position[i], shift):
+                new_coord = pp + s
+                if new_coord < self.min_coord:
+                    new_position.append(2.0 * self.min_coord - pp - s)
+                    final_shift.append(-1.0 * s)
+                elif new_coord > self.max_coord:
+                    new_position.append(2.0 * self.max_coord - pp - s)
+                    final_shift.append(-1.0 * s)
+                else:
+                    new_position.append(new_coord)
+                    final_shift.append(s)
+
+            self.peaks_position[i] = new_position
+            self.last_change_vector[i] = final_shift
+
+            # Change peak height
+            change = self.random.gauss(0, 1) * self.height_severity
+            new_value = change + self.peaks_height[i]
+            if new_value < self.min_height:
+                self.peaks_height[i] = 2.0 * self.min_height - self.peaks_height[i] - change
+            elif new_value > self.max_height:
+                self.peaks_height[i] = 2.0 * self.max_height - self.peaks_height[i] - change
+            else:
+                self.peaks_height[i] = new_value
+
+            # Change peak width
+            change = self.random.gauss(0, 1) * self.width_severity
+            new_value = change + self.peaks_width[i]
+            if new_value < self.min_width:
+                self.peaks_width[i] = 2.0 * self.min_width - self.peaks_width[i] - change
+            elif new_value > self.max_width:
+                self.peaks_width[i] = 2.0 * self.max_width - self.peaks_width[i] - change
+            else:
+                self.peaks_width[i] = new_value
+
+        self._optimum = None
+
+SCENARIO_1 = {"pfunc" : function1,
+              "npeaks" : 5,
+              "bfunc": None,
+              "min_coord": 0.0,
+              "max_coord": 100.0,
+              "min_height": 30.0,
+              "max_height": 70.0,
+              "uniform_height": 50.0,
+              "min_width": 0.0001,
+              "max_width": 0.2,
+              "uniform_width": 0.1,
+              "lambda_": 0.0,
+              "move_severity": 1.0,
+              "height_severity": 7.0,
+              "width_severity": 0.01,
+              "period": 5000}
+
+SCENARIO_2 = {"pfunc" : cone,
+              "npeaks" : 10,
+              "bfunc" : None,
+              "min_coord": 0.0,
+              "max_coord": 100.0,
+              "min_height": 30.0,
+              "max_height": 70.0,
+              "uniform_height": 50.0,
+              "min_width": 1.0,
+              "max_width": 12.0,
+              "uniform_width": 0,
+              "lambda_": 0.5,
+              "move_severity": 1.0,
+              "height_severity": 7.0,
+              "width_severity": 1.0,
+              "period": 5000}
+
+SCENARIO_3 = {"pfunc" : cone,
+              "npeaks" : 50,
+              "bfunc" : lambda x: 10,
+              "min_coord": 0.0,
+              "max_coord": 100.0,
+              "min_height": 30.0,
+              "max_height": 70.0,
+              "uniform_height": 0,
+              "min_width": 1.0,
+              "max_width": 12.0,
+              "uniform_width": 0,
+              "lambda_": 0.5,
+              "move_severity": 1.0,
+              "height_severity": 1.0,
+              "width_severity": 0.5,
+              "period": 1000}
+
+def diversity(population):
+    nind = len(population)
+    ndim = len(population[0])
+    d = [0.0] * ndim
+    for x in population:
+        d = [di + xi for di, xi in zip(d, x)]
+    d = [di / nind for di in d]
+    return math.sqrt(sum((di - xi)**2 for x in population for di, xi in zip(d, x)))
+
+if __name__ == "__main__":
+    mpb = MovingPeaks(dim=2, npeaks=[1,1,10], number_severity=0.1)
+    print mpb.maximums()
+    mpb.changePeaks()
+    print mpb.maximums()

python/isaac/autotuning/external/deap/benchmarks/tools.py

@@ -0,0 +1,282 @@
+"""Module containing tools that are useful when benchmarking algorithms
+"""
+from math import hypot, sqrt
+from functools import wraps
+from itertools import repeat
+try:
+    import numpy
+except ImportError:
+    numpy = False
+
+class translate(object):
+    """Decorator for evaluation functions, it translates the objective
+    function by *vector* which should be the same length as the individual
+    size. When called the decorated function should take as first argument the
+    individual to be evaluated. The inverse translation vector is actually
+    applied to the individual and the resulting list is given to the
+    evaluation function. Thus, the evaluation function shall not be expecting
+    an individual as it will receive a plain list.
+    
+    This decorator adds a :func:`translate` method to the decorated function.
+    """
+    def __init__(self, vector):
+        self.vector = vector
+    
+    def __call__(self, func):
+        # wraps is used to combine stacked decorators that would add functions
+        @wraps(func)
+        def wrapper(individual, *args, **kargs):
+            # A subtraction is applied since the translation is applied to the
+            # individual and not the function
+            return func([v - t for v, t in zip(individual, self.vector)], 
+                *args, **kargs)
+        wrapper.translate = self.translate
+        return wrapper
+    
+    def translate(self, vector):
+        """Set the current translation to *vector*. After decorating the
+        evaluation function, this function will be available directly from
+        the function object. ::
+            
+            @translate([0.25, 0.5, ..., 0.1])
+            def evaluate(individual):
+                return sum(individual),
+            
+            # This will cancel the translation
+            evaluate.translate([0.0, 0.0, ..., 0.0])
+        """
+        self.vector = vector
+
+class rotate(object):
+    """Decorator for evaluation functions, it rotates the objective function
+    by *matrix* which should be a valid orthogonal NxN rotation matrix, with N
+    the length of an individual. When called the decorated function should
+    take as first argument the individual to be evaluated. The inverse
+    rotation matrix is actually applied to the individual and the resulting
+    list is given to the evaluation function. Thus, the evaluation function
+    shall not be expecting an individual as it will receive a plain list
+    (numpy.array). The multiplication is done using numpy.
+    
+    This decorator adds a :func:`rotate` method to the decorated function.
+    
+    .. note::
+    
+       A random orthogonal matrix Q can be created via QR decomposition. ::
+           
+           A = numpy.random.random((n,n))
+           Q, _ = numpy.linalg.qr(A)
+    """
+    def __init__(self, matrix):
+        if not numpy:
+            raise RuntimeError("Numpy is required for using the rotation "
+                "decorator")
+        # The inverse is taken since the rotation is applied to the individual
+        # and not the function which is the inverse
+        self.matrix = numpy.linalg.inv(matrix)
+
+    def __call__(self, func):
+        # wraps is used to combine stacked decorators that would add functions
+        @wraps(func)
+        def wrapper(individual, *args, **kargs):
+            return func(numpy.dot(self.matrix, individual), *args, **kargs)
+        wrapper.rotate = self.rotate
+        return wrapper
+
+    def rotate(self, matrix):
+        """Set the current rotation to *matrix*. After decorating the
+        evaluation function, this function will be available directly from
+        the function object. ::
+            
+            # Create a random orthogonal matrix
+            A = numpy.random.random((n,n))
+            Q, _ = numpy.linalg.qr(A)
+            
+            @rotate(Q)
+            def evaluate(individual):
+                return sum(individual),
+            
+            # This will reset rotation to identity
+            evaluate.rotate(numpy.identity(n))
+        """
+        self.matrix = numpy.linalg.inv(matrix)
+
+class noise(object):
+    """Decorator for evaluation functions, it evaluates the objective function
+    and adds noise by calling the function(s) provided in the *noise*
+    argument. The noise functions are called without any argument, consider
+    using the :class:`~deap.base.Toolbox` or Python's
+    :func:`functools.partial` to provide any required argument. If a single
+    function is provided it is applied to all objectives of the evaluation
+    function. If a list of noise functions is provided, it must be of length
+    equal to the number of objectives. The noise argument also accept
+    :obj:`None`, which will leave the objective without noise.
+
+    This decorator adds a :func:`noise` method to the decorated
+    function.
+    """
+    def __init__(self, noise):
+        try:
+            self.rand_funcs = tuple(noise)
+        except TypeError:
+            self.rand_funcs = repeat(noise)
+
+    def __call__(self, func):
+        # wraps is used to combine stacked decorators that would add functions
+        @wraps(func)
+        def wrapper(individual, *args, **kargs):
+            result = func(individual, *args, **kargs)
+            noisy = list()
+            for r, f in zip(result, self.rand_funcs):
+                if f is None:
+                    noisy.append(r)
+                else:
+                    noisy.append(r + f())
+            return tuple(noisy)
+        wrapper.noise = self.noise
+        return wrapper
+    
+    def noise(self, noise):
+        """Set the current noise to *noise*. After decorating the
+        evaluation function, this function will be available directly from
+        the function object. ::
+        
+            prand = functools.partial(random.gauss, mu=0.0, sigma=1.0)
+        
+            @noise(prand)
+            def evaluate(individual):
+                return sum(individual),
+        
+            # This will remove noise from the evaluation function
+            evaluate.noise(None)
+        """
+        try:
+            self.rand_funcs = tuple(noise)
+        except TypeError:
+            self.rand_funcs = repeat(noise)
+            
+class scale(object):
+    """Decorator for evaluation functions, it scales the objective function by
+    *factor* which should be the same length as the individual size. When
+    called the decorated function should take as first argument the individual
+    to be evaluated. The inverse factor vector is actually applied to the
+    individual and the resulting list is given to the evaluation function.
+    Thus, the evaluation function shall not be expecting an individual as it
+    will receive a plain list.
+    
+    This decorator adds a :func:`scale` method to the decorated function.
+    """
+    def __init__(self, factor):
+        # Factor is inverted since it is aplied to the individual and not the
+        # objective function
+        self.factor = tuple(1.0/f for f in factor)
+
+    def __call__(self, func):
+        # wraps is used to combine stacked decorators that would add functions
+        @wraps(func)
+        def wrapper(individual, *args, **kargs):
+            return func([v * f for v, f in zip(individual, self.factor)], 
+                *args, **kargs)
+        wrapper.scale = self.scale
+        return wrapper
+
+    def scale(self, factor):
+        """Set the current scale to *factor*. After decorating the
+        evaluation function, this function will be available directly from
+        the function object. ::
+            
+            @scale([0.25, 2.0, ..., 0.1])
+            def evaluate(individual):
+                return sum(individual),
+            
+            # This will cancel the scaling
+            evaluate.scale([1.0, 1.0, ..., 1.0])
+        """
+        # Factor is inverted since it is aplied to the individual and not the
+        # objective function
+        self.factor = tuple(1.0/f for f in factor)
+
+class bound(object):
+    """Decorator for crossover and mutation functions, it changes the
+    individuals after the modification is done to bring it back in the allowed
+    *bounds*. The *bounds* are functions taking individual and returning
+    wheter of not the variable is allowed. You can provide one or multiple such
+    functions. In the former case, the function is used on all dimensions and
+    in the latter case, the number of functions must be greater or equal to
+    the number of dimension of the individuals.
+
+    The *type* determines how the attributes are brought back into the valid
+    range
+    
+    This decorator adds a :func:`bound` method to the decorated function.
+    """
+    def _clip(self, individual):
+        return individual
+
+    def _wrap(self, individual):
+        return individual
+
+    def _mirror(self, individual):
+        return individual
+
+    def __call__(self, func):
+        @wraps(func)
+        def wrapper(*args, **kargs):
+            individuals = func(*args, **kargs)
+            return self.bound(individuals)
+        wrapper.bound = self.bound
+        return wrapper
+
+    def __init__(self, bounds, type):
+        try:
+            self.bounds = tuple(bounds)
+        except TypeError:
+            self.bounds = itertools.repeat(bounds)
+
+        if type == "mirror":
+            self.bound = self._mirror
+        elif type == "wrap":
+            self.bound = self._wrap
+        elif type == "clip":
+            self.bound = self._clip
+        
+def diversity(first_front, first, last):
+    """Given a Pareto front `first_front` and the two extreme points of the 
+    optimal Pareto front, this function returns a metric of the diversity 
+    of the front as explained in the original NSGA-II article by K. Deb.
+    The smaller the value is, the better the front is.
+    """
+    df = hypot(first_front[0].fitness.values[0] - first[0],
+               first_front[0].fitness.values[1] - first[1])
+    dl = hypot(first_front[-1].fitness.values[0] - last[0],
+               first_front[-1].fitness.values[1] - last[1])
+    dt = [hypot(first.fitness.values[0] - second.fitness.values[0],
+                first.fitness.values[1] - second.fitness.values[1])
+          for first, second in zip(first_front[:-1], first_front[1:])]
+
+    if len(first_front) == 1:
+        return df + dl
+
+    dm = sum(dt)/len(dt)
+    di = sum(abs(d_i - dm) for d_i in dt)
+    delta = (df + dl + di)/(df + dl + len(dt) * dm )
+    return delta
+
+def convergence(first_front, optimal_front):
+    """Given a Pareto front `first_front` and the optimal Pareto front, 
+    this function returns a metric of convergence
+    of the front as explained in the original NSGA-II article by K. Deb.
+    The smaller the value is, the closer the front is to the optimal one.
+    """
+    distances = []
+    
+    for ind in first_front:
+        distances.append(float("inf"))
+        for opt_ind in optimal_front:
+            dist = 0.
+            for i in xrange(len(opt_ind)):
+                dist += (ind.fitness.values[i] - opt_ind[i])**2
+            if dist < distances[-1]:
+                distances[-1] = dist
+        distances[-1] = sqrt(distances[-1])
+        
+    return sum(distances) / len(distances)

python/isaac/autotuning/external/deap/cma.py

@@ -0,0 +1,292 @@
+#    This file is part of DEAP.
+#
+#    DEAP is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU Lesser General Public License as
+#    published by the Free Software Foundation, either version 3 of
+#    the License, or (at your option) any later version.
+#
+#    DEAP is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#    GNU Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.
+
+#    Special thanks to Nikolaus Hansen for providing major part of 
+#    this code. The CMA-ES algorithm is provided in many other languages
+#    and advanced versions at http://www.lri.fr/~hansen/cmaesintro.html.
+
+"""A module that provides support for the Covariance Matrix Adaptation 
+Evolution Strategy.
+"""
+import numpy
+import copy
+from math import sqrt, log, exp
+
+class Strategy(object):
+    """
+    A strategy that will keep track of the basic parameters of the CMA-ES
+    algorithm.
+    
+    :param centroid: An iterable object that indicates where to start the
+                     evolution.
+    :param sigma: The initial standard deviation of the distribution.
+    :param parameter: One or more parameter to pass to the strategy as
+                      described in the following table, optional.
+    
+    +----------------+---------------------------+----------------------------+
+    | Parameter      | Default                   | Details                    |
+    +================+===========================+============================+
+    | ``lambda_``    | ``int(4 + 3 * log(N))``   | Number of children to      |
+    |                |                           | produce at each generation,|
+    |                |                           | ``N`` is the individual's  |
+    |                |                           | size (integer).            |
+    +----------------+---------------------------+----------------------------+
+    | ``mu``         | ``int(lambda_ / 2)``      | The number of parents to   | 
+    |                |                           | keep from the              |
+    |                |                           | lambda children (integer). |
+    +----------------+---------------------------+----------------------------+
+    | ``cmatrix``    | ``identity(N)``           | The initial covariance     |
+    |                |                           | matrix of the distribution |
+    |                |                           | that will be sampled.      |
+    +----------------+---------------------------+----------------------------+
+    | ``weights``    | ``"superlinear"``         | Decrease speed, can be     |
+    |                |                           | ``"superlinear"``,         |
+    |                |                           | ``"linear"`` or            |
+    |                |                           | ``"equal"``.               |
+    +----------------+---------------------------+----------------------------+
+    | ``cs``         | ``(mueff + 2) /           | Cumulation constant for    |
+    |                | (N + mueff + 3)``         | step-size.                 |
+    +----------------+---------------------------+----------------------------+
+    | ``damps``      | ``1 + 2 * max(0, sqrt((   | Damping for step-size.     |
+    |                | mueff - 1) / (N + 1)) - 1)|                            |
+    |                | + cs``                    |                            |
+    +----------------+---------------------------+----------------------------+
+    | ``ccum``       | ``4 / (N + 4)``           | Cumulation constant for    |
+    |                |                           | covariance matrix.         |
+    +----------------+---------------------------+----------------------------+
+    | ``ccov1``      | ``2 / ((N + 1.3)^2 +      | Learning rate for rank-one |
+    |                | mueff)``                  | update.                    |
+    +----------------+---------------------------+----------------------------+
+    | ``ccovmu``     | ``2 * (mueff - 2 + 1 /    | Learning rate for rank-mu  |
+    |                | mueff) / ((N + 2)^2 +     | update.                    |
+    |                | mueff)``                  |                            |
+    +----------------+---------------------------+----------------------------+
+
+    """
+    def __init__(self, centroid, sigma, **kargs):
+        self.params = kargs
+        
+        # Create a centroid as a numpy array
+        self.centroid = numpy.array(centroid)
+        
+        self.dim = len(self.centroid)
+        self.sigma = sigma
+        self.pc = numpy.zeros(self.dim)
+        self.ps = numpy.zeros(self.dim)
+        self.chiN = sqrt(self.dim) * (1 - 1. / (4. * self.dim) + \
+                                      1. / (21. * self.dim**2))
+        
+        self.C = self.params.get("cmatrix", numpy.identity(self.dim))
+        self.diagD, self.B = numpy.linalg.eigh(self.C)
+
+        indx = numpy.argsort(self.diagD)
+        self.diagD = self.diagD[indx]**0.5
+        self.B = self.B[:, indx]
+        self.BD = self.B * self.diagD
+        
+        self.cond = self.diagD[indx[-1]]/self.diagD[indx[0]]
+        
+        self.lambda_ = self.params.get("lambda_", int(4 + 3 * log(self.dim)))
+        self.update_count = 0
+        self.computeParams(self.params)
+        
+    def generate(self, ind_init):
+        """Generate a population of :math:`\lambda` individuals of type
+        *ind_init* from the current strategy.
+        
+        :param ind_init: A function object that is able to initialize an
+                         individual from a list.
+        :returns: A list of individuals.
+        """
+        arz = numpy.random.standard_normal((self.lambda_, self.dim))
+        arz = self.centroid + self.sigma * numpy.dot(arz, self.BD.T)
+        return map(ind_init, arz)
+        
+    def update(self, population):
+        """Update the current covariance matrix strategy from the
+        *population*.
+        
+        :param population: A list of individuals from which to update the
+                           parameters.
+        """
+        population.sort(key=lambda ind: ind.fitness, reverse=True)
+        
+        old_centroid = self.centroid
+        self.centroid = numpy.dot(self.weights, population[0:self.mu])
+        
+        c_diff = self.centroid - old_centroid
+        
+        # Cumulation : update evolution path
+        self.ps = (1 - self.cs) * self.ps \
+             + sqrt(self.cs * (2 - self.cs) * self.mueff) / self.sigma \
+             * numpy.dot(self.B, (1. / self.diagD) \
+                          * numpy.dot(self.B.T, c_diff))
+        
+        hsig = float((numpy.linalg.norm(self.ps) / 
+                sqrt(1. - (1. - self.cs)**(2. * (self.update_count + 1.))) / self.chiN
+                < (1.4 + 2. / (self.dim + 1.))))
+        
+        self.update_count += 1
+        
+        self.pc = (1 - self.cc) * self.pc + hsig \
+                  * sqrt(self.cc * (2 - self.cc) * self.mueff) / self.sigma \
+                  * c_diff
+        
+        # Update covariance matrix
+        artmp = population[0:self.mu] - old_centroid
+        self.C = (1 - self.ccov1 - self.ccovmu + (1 - hsig) \
+                   * self.ccov1 * self.cc * (2 - self.cc)) * self.C \
+                + self.ccov1 * numpy.outer(self.pc, self.pc) \
+                + self.ccovmu * numpy.dot((self.weights * artmp.T), artmp) \
+                / self.sigma**2
+        
+        
+        self.sigma *= numpy.exp((numpy.linalg.norm(self.ps) / self.chiN - 1.) \
+                                * self.cs / self.damps)
+        
+        self.diagD, self.B = numpy.linalg.eigh(self.C)
+        indx = numpy.argsort(self.diagD)
+        
+        self.cond = self.diagD[indx[-1]]/self.diagD[indx[0]]
+        
+        self.diagD = self.diagD[indx]**0.5
+        self.B = self.B[:, indx]
+        self.BD = self.B * self.diagD
+
+    def computeParams(self, params):
+        """Computes the parameters depending on :math:`\lambda`. It needs to
+        be called again if :math:`\lambda` changes during evolution.
+        
+        :param params: A dictionary of the manually set parameters.
+        """
+        self.mu = params.get("mu", int(self.lambda_ / 2))
+        rweights = params.get("weights", "superlinear")
+        if rweights == "superlinear":
+            self.weights = log(self.mu + 0.5) - \
+                        numpy.log(numpy.arange(1, self.mu + 1))
+        elif rweights == "linear":
+            self.weights = self.mu + 0.5 - numpy.arange(1, self.mu + 1)
+        elif rweights == "equal":
+            self.weights = numpy.ones(self.mu)
+        else:
+            raise RuntimeError("Unknown weights : %s" % rweights)
+        
+        self.weights /= sum(self.weights)
+        self.mueff = 1. / sum(self.weights**2)
+        
+        self.cc = params.get("ccum", 4. / (self.dim + 4.))
+        self.cs = params.get("cs", (self.mueff + 2.) / 
+                                   (self.dim + self.mueff + 3.))
+        self.ccov1 = params.get("ccov1", 2. / ((self.dim + 1.3)**2 + \
+                                         self.mueff))
+        self.ccovmu = params.get("ccovmu", 2. * (self.mueff - 2. +  \
+                                                 1. / self.mueff) / \
+                                           ((self.dim + 2.)**2 + self.mueff))
+        self.ccovmu = min(1 - self.ccov1, self.ccovmu)
+        self.damps = 1. + 2. * max(0, sqrt((self.mueff - 1.) / \
+                                            (self.dim + 1.)) - 1.) + self.cs
+        self.damps = params.get("damps", self.damps)
+        
+class StrategyOnePlusLambda(object):
+    """
+    A CMA-ES strategy that uses the :math:`1 + \lambda` paradigme.
+    
+    :param parent: An iterable object that indicates where to start the
+                   evolution. The parent requires a fitness attribute.
+    :param sigma: The initial standard deviation of the distribution.
+    :param parameter: One or more parameter to pass to the strategy as
+                      described in the following table, optional.
+    """
+    def __init__(self, parent, sigma, **kargs):
+        self.parent = parent
+        self.sigma = sigma
+        self.dim = len(self.parent)
+
+        self.C = numpy.identity(self.dim)
+        self.A = numpy.identity(self.dim)
+        
+        self.pc = numpy.zeros(self.dim)
+        
+        self.computeParams(kargs)
+        self.psucc = self.ptarg
+        
+    def computeParams(self, params):
+        """Computes the parameters depending on :math:`\lambda`. It needs to
+        be called again if :math:`\lambda` changes during evolution.
+        
+        :param params: A dictionary of the manually set parameters.
+        """
+        # Selection :
+        self.lambda_ = params.get("lambda_", 1)
+        
+        # Step size control :
+        self.d = params.get("d", 1.0 + self.dim/(2.0*self.lambda_))
+        self.ptarg = params.get("ptarg", 1.0/(5+sqrt(self.lambda_)/2.0))
+        self.cp = params.get("cp", self.ptarg*self.lambda_/(2+self.ptarg*self.lambda_))
+        
+        # Covariance matrix adaptation
+        self.cc = params.get("cc", 2.0/(self.dim+2.0))
+        self.ccov = params.get("ccov", 2.0/(self.dim**2 + 6.0))
+        self.pthresh = params.get("pthresh", 0.44)
+    
+    def generate(self, ind_init):
+        """Generate a population of :math:`\lambda` individuals of type
+        *ind_init* from the current strategy.
+        
+        :param ind_init: A function object that is able to initialize an
+                         individual from a list.
+        :returns: A list of individuals.
+        """
+        # self.y = numpy.dot(self.A, numpy.random.standard_normal(self.dim))
+        arz = numpy.random.standard_normal((self.lambda_, self.dim))
+        arz = self.parent + self.sigma * numpy.dot(arz, self.A.T)        
+        return map(ind_init, arz)
+    
+    def update(self, population):
+        """Update the current covariance matrix strategy from the
+        *population*.
+        
+        :param population: A list of individuals from which to update the
+                           parameters.
+        """
+        population.sort(key=lambda ind: ind.fitness, reverse=True)
+        lambda_succ = sum(self.parent.fitness <= ind.fitness for ind in population)
+        p_succ = float(lambda_succ) / self.lambda_
+        self.psucc = (1-self.cp)*self.psucc + self.cp*p_succ
+        
+        if self.parent.fitness <= population[0].fitness:
+            x_step = (population[0] - numpy.array(self.parent)) / self.sigma
+            self.parent = copy.deepcopy(population[0])
+            if self.psucc < self.pthresh:
+                self.pc = (1 - self.cc)*self.pc + sqrt(self.cc * (2 - self.cc)) * x_step
+                self.C = (1-self.ccov)*self.C + self.ccov * numpy.outer(self.pc, self.pc)
+            else:
+                self.pc = (1 - self.cc)*self.pc
+                self.C = (1-self.ccov)*self.C + self.ccov * (numpy.outer(self.pc, self.pc) + self.cc*(2-self.cc)*self.C)
+
+        self.sigma = self.sigma * exp(1.0/self.d * (self.psucc - self.ptarg)/(1.0-self.ptarg))
+        
+        # We use Cholesky since for now we have no use of eigen decomposition
+        # Basically, Cholesky returns a matrix A as C = A*A.T
+        # Eigen decomposition returns two matrix B and D^2 as C = B*D^2*B.T = B*D*D*B.T
+        # So A == B*D
+        # To compute the new individual we need to multiply each vector z by A
+        # as y = centroid + sigma * A*z
+        # So the Cholesky is more straightforward as we don't need to compute 
+        # the squareroot of D^2, and multiply B and D in order to get A, we directly get A.
+        # This can't be done (without cost) with the standard CMA-ES as the eigen decomposition is used
+        # to compute covariance matrix inverse in the step-size evolutionary path computation.
+        self.A = numpy.linalg.cholesky(self.C)
+        

python/isaac/autotuning/external/deap/creator.py

@@ -0,0 +1,154 @@
+#    This file is part of DEAP.
+#
+#    DEAP is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU Lesser General Public License as
+#    published by the Free Software Foundation, either version 3 of
+#    the License, or (at your option) any later version.
+#
+#    DEAP is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#    GNU Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.
+
+"""The :mod:`~deap.creator` is a meta-factory allowing to create classes that
+will fulfill the needs of your evolutionary algorithms. In effect, new
+classes can be built from any imaginable type, from :class:`list` to
+:class:`set`, :class:`dict`, :class:`~deap.gp.PrimitiveTree` and more,
+providing the possibility to implement genetic algorithms, genetic
+programming, evolution strategies, particle swarm optimizers, and many more.
+"""
+
+import array
+import copy
+
+class_replacers = {}
+"""Some classes in Python's standard library as well as third party library
+may be in part incompatible with the logic used in DEAP. To palliate
+this problem, the method :func:`create` uses the dictionary
+`class_replacers` to identify if the base type provided is problematic, and if
+so  the new class inherits from the replacement class instead of the
+original base class.
+
+`class_replacers` keys are classes to be replaced and the values are the
+replacing classes.
+"""
+
+try:
+    import numpy
+    (numpy.ndarray, numpy.array)
+except ImportError:
+    # Numpy is not present, skip the definition of the replacement class.
+    pass
+except AttributeError:
+    # Numpy is present, but there is either no ndarray or array in numpy,
+    # also skip the definition of the replacement class.
+    pass
+else:
+    class _numpy_array(numpy.ndarray):
+        def __deepcopy__(self, memo):
+            """Overrides the deepcopy from numpy.ndarray that does not copy
+            the object's attributes. This one will deepcopy the array and its
+            :attr:`__dict__` attribute.
+            """
+            copy_ = numpy.ndarray.copy(self)
+            copy_.__dict__.update(copy.deepcopy(self.__dict__, memo))
+            return copy_
+
+        @staticmethod
+        def __new__(cls, iterable):
+            """Creates a new instance of a numpy.ndarray from a function call.
+            Adds the possibility to instanciate from an iterable."""
+            return numpy.array(list(iterable)).view(cls)
+            
+        def __setstate__(self, state):
+            self.__dict__.update(state)
+
+        def __reduce__(self):
+            return (self.__class__, (list(self),), self.__dict__)
+
+    class_replacers[numpy.ndarray] = _numpy_array
+
+class _array(array.array):
+    @staticmethod
+    def __new__(cls, seq=()):
+        return super(_array, cls).__new__(cls, cls.typecode, seq)
+    
+    def __deepcopy__(self, memo):
+        """Overrides the deepcopy from array.array that does not copy
+        the object's attributes and class type.
+        """
+        cls = self.__class__
+        copy_ = cls.__new__(cls, self)
+        memo[id(self)] = copy_
+        copy_.__dict__.update(copy.deepcopy(self.__dict__, memo))
+        return copy_
+
+    def __reduce__(self):
+        return (self.__class__, (list(self),), self.__dict__)
+class_replacers[array.array] = _array
+
+def create(name, base, **kargs):
+    """Creates a new class named *name* inheriting from *base* in the
+    :mod:`~deap.creator` module. The new class can have attributes defined by
+    the subsequent keyword arguments passed to the function create. If the
+    argument is a class (without the parenthesis), the __init__ function is
+    called in the initialization of an instance of the new object and the
+    returned instance is added as an attribute of the class' instance.
+    Otherwise, if the argument is not a class, (for example an :class:`int`),
+    it is added as a "static" attribute of the class.
+    
+    :param name: The name of the class to create.
+    :param base: A base class from which to inherit.
+    :param attribute: One or more attributes to add on instanciation of this
+                      class, optional.
+    
+    The following is used to create a class :class:`Foo` inheriting from the
+    standard :class:`list` and having an attribute :attr:`bar` being an empty
+    dictionary and a static attribute :attr:`spam` initialized to 1. ::
+    
+        create("Foo", list, bar=dict, spam=1)
+        
+    This above line is exactly the same as defining in the :mod:`creator`
+    module something like the following. ::
+    
+        class Foo(list):
+            spam = 1
+            
+            def __init__(self):
+                self.bar = dict()
+
+    The :ref:`creating-types` tutorial gives more examples of the creator
+    usage.
+    """
+    dict_inst = {}
+    dict_cls = {}
+    for obj_name, obj in kargs.iteritems():
+        if isinstance(obj, type):
+            dict_inst[obj_name] = obj
+        else:
+            dict_cls[obj_name] = obj
+
+    # Check if the base class has to be replaced
+    if base in class_replacers:
+        base = class_replacers[base]
+
+    # A DeprecationWarning is raised when the object inherits from the 
+    # class "object" which leave the option of passing arguments, but
+    # raise a warning stating that it will eventually stop permitting
+    # this option. Usually this happens when the base class does not
+    # override the __init__ method from object.
+    def initType(self, *args, **kargs):
+        """Replace the __init__ function of the new type, in order to
+        add attributes that were defined with **kargs to the instance.
+        """
+        for obj_name, obj in dict_inst.iteritems():
+            setattr(self, obj_name, obj())
+        if base.__init__ is not object.__init__:
+            base.__init__(self, *args, **kargs)
+
+    objtype = type(str(name), (base,), dict_cls)
+    objtype.__init__ = initType
+    globals()[name] = objtype

python/isaac/autotuning/external/deap/gp.py

@@ -0,0 +1,952 @@
+#    This file is part of DEAP.
+#
+#    DEAP is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU Lesser General Public License as
+#    published by the Free Software Foundation, either version 3 of
+#    the License, or (at your option) any later version.
+#
+#    DEAP is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#    GNU Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.
+
+"""The :mod:`gp` module provides the methods and classes to perform
+Genetic Programming with DEAP. It essentially contains the classes to
+build a Genetic Program Tree, and the functions to evaluate it.
+
+This module support both strongly and loosely typed GP.
+"""
+import copy
+import random
+import re
+import sys
+import warnings
+
+from collections import defaultdict, deque
+from functools import partial, wraps
+from inspect import isclass
+from operator import eq, lt
+
+######################################
+# GP Data structure                  #
+######################################
+
+# Define the name of type for any types.
+__type__ = object
+
+class PrimitiveTree(list):
+    """Tree spefically formated for optimization of genetic
+    programming operations. The tree is represented with a
+    list where the nodes are appended in a depth-first order.
+    The nodes appended to the tree are required to
+    have an attribute *arity* which defines the arity of the
+    primitive. An arity of 0 is expected from terminals nodes.
+    """
+    def __init__(self, content):
+        list.__init__(self, content)
+
+    def __deepcopy__(self, memo):
+        new = self.__class__(self)
+        new.__dict__.update(copy.deepcopy(self.__dict__, memo))
+        return new
+
+    def __setitem__(self, key, val):
+        # Check for most common errors
+        # Does NOT check for STGP constraints
+        if isinstance(key, slice):
+            if key.start >= len(self):
+                raise IndexError("Invalid slice object (try to assign a %s"
+                    " in a tree of size %d). Even if this is allowed by the"
+                    " list object slice setter, this should not be done in"
+                    " the PrimitiveTree context, as this may lead to an"
+                    " unpredictable behavior for searchSubtree or evaluate."
+                     % (key, len(self)))
+            total = val[0].arity
+            for node in val[1:]:
+                total += node.arity - 1
+            if total != 0:
+                raise ValueError("Invalid slice assignation : insertion of"
+                    " an incomplete subtree is not allowed in PrimitiveTree."
+                    " A tree is defined as incomplete when some nodes cannot"
+                    " be mapped to any position in the tree, considering the"
+                    " primitives' arity. For instance, the tree [sub, 4, 5,"
+                    " 6] is incomplete if the arity of sub is 2, because it"
+                    " would produce an orphan node (the 6).")
+        elif val.arity != self[key].arity:
+            raise ValueError("Invalid node replacement with a node of a"
+                             " different arity.")
+        list.__setitem__(self, key, val)
+
+    def __str__(self):
+        """Return the expression in a human readable string.
+        """
+        string = ""
+        stack = []
+        for node in self:
+            stack.append((node, []))
+            while len(stack[-1][1]) == stack[-1][0].arity:
+                prim, args = stack.pop()
+                string = prim.format(*args)
+                if len(stack) == 0:
+                    break   # If stack is empty, all nodes should have been seen
+                stack[-1][1].append(string)
+
+        return string
+
+    @classmethod
+    def from_string(cls, string, pset):
+        """Try to convert a string expression into a PrimitiveTree given a
+        PrimitiveSet *pset*. The primitive set needs to contain every primitive
+        present in the expression.
+
+        :param string: String representation of a Python expression.
+        :param pset: Primitive set from which primitives are selected.
+        :returns: PrimitiveTree populated with the deserialized primitives.
+        """
+        tokens = re.split("[ \t\n\r\f\v(),]", string)
+        expr = []
+        ret_types = deque()
+        for token in tokens:
+            if token == '':
+                continue
+            if len(ret_types) != 0:
+                type_ = ret_types.popleft()
+            else:
+                type_ = None
+
+            if token in pset.mapping:
+                primitive = pset.mapping[token]
+
+                if type_ is not None and not issubclass(primitive.ret, type_):
+                    raise TypeError("Primitive {} return type {} does not "
+                                    "match the expected one: {}."
+                                    .format(primitive, primitive.ret, type_))
+
+                expr.append(primitive)
+                if isinstance(primitive, Primitive):
+                    ret_types.extendleft(reversed(primitive.args))
+            else:
+                try:
+                    token = eval(token)
+                except NameError:
+                    raise TypeError("Unable to evaluate terminal: {}.".format(token))
+
+                if type_ is None:
+                    type_ = type(token)
+
+                if not issubclass(type(token), type_):
+                    raise TypeError("Terminal {} type {} does not "
+                                    "match the expected one: {}."
+                                    .format(token, type(token), type_))
+
+                expr.append(Terminal(token, False, type_))
+        return cls(expr)
+
+    @property
+    def height(self):
+        """Return the height of the tree, or the depth of the
+        deepest node.
+        """
+        stack = [0]
+        max_depth = 0
+        for elem in self:
+            depth = stack.pop()
+            max_depth = max(max_depth, depth)
+            stack.extend([depth+1] * elem.arity)
+        return max_depth
+
+    @property
+    def root(self):
+        """Root of the tree, the element 0 of the list.
+        """
+        return self[0]
+
+    def searchSubtree(self, begin):
+        """Return a slice object that corresponds to the
+        range of values that defines the subtree which has the
+        element with index *begin* as its root.
+        """
+        end = begin + 1
+        total = self[begin].arity
+        while total > 0:
+            total += self[end].arity - 1
+            end += 1
+        return slice(begin, end)
+
+
+class Primitive(object):
+    """Class that encapsulates a primitive and when called with arguments it
+    returns the Python code to call the primitive with the arguments.
+
+        >>> pr = Primitive("mul", (int, int), int)
+        >>> pr.format(1, 2)
+        'mul(1, 2)'
+    """
+    __slots__ = ('name', 'arity', 'args', 'ret', 'seq')
+    def __init__(self, name, args, ret):
+        self.name = name
+        self.arity = len(args)
+        self.args = args
+        self.ret = ret
+        args = ", ".join(map("{{{0}}}".format, range(self.arity)))
+        self.seq = "{name}({args})".format(name=self.name, args=args)
+
+    def format(self, *args):
+        return self.seq.format(*args)
+
+class Terminal(object):
+    """Class that encapsulates terminal primitive in expression. Terminals can
+    be values or 0-arity functions.
+    """
+    __slots__ = ('name', 'value', 'ret', 'conv_fct')
+    def __init__(self, terminal, symbolic, ret):
+        self.ret = ret
+        self.value = terminal
+        self.name = str(terminal)
+        self.conv_fct = str if symbolic else repr
+
+    @property
+    def arity(self):
+        return 0
+
+    def format(self):
+        return self.conv_fct(self.value)
+
+class Ephemeral(Terminal):
+    """Class that encapsulates a terminal which value is set when the
+    object is created. To mutate the value, a new object has to be
+    generated. This is an abstract base class. When subclassing, a
+    staticmethod 'func' must be defined.
+    """
+    def __init__(self):
+        Terminal.__init__(self, self.func(), symbolic=False, ret=self.ret)
+
+    @staticmethod
+    def func():
+        """Return a random value used to define the ephemeral state.
+        """
+        raise NotImplementedError
+
+class PrimitiveSetTyped(object):
+    """Class that contains the primitives that can be used to solve a
+    Strongly Typed GP problem. The set also defined the researched
+    function return type, and input arguments type and number.
+    """
+    def __init__(self, name, in_types, ret_type, prefix="ARG"):
+        self.terminals = defaultdict(list)
+        self.primitives = defaultdict(list)
+        self.arguments = []
+        # setting "__builtins__" to None avoid the context
+        # being polluted by builtins function when evaluating
+        # GP expression.
+        self.context = {"__builtins__" : None}
+        self.mapping = dict()
+        self.terms_count = 0
+        self.prims_count = 0
+
+        self.name = name
+        self.ret = ret_type
+        self.ins = in_types
+        for i, type_ in enumerate(in_types):
+            arg_str = "{prefix}{index}".format(prefix=prefix, index=i)
+            self.arguments.append(arg_str)
+            term = Terminal(arg_str, True, type_)
+            self._add(term)
+            self.terms_count += 1
+
+    def renameArguments(self, **kargs):
+        """Rename function arguments with new names from *kargs*.
+        """
+        for i, old_name in enumerate(self.arguments):
+            if old_name in kargs:
+                new_name = kargs[old_name]
+                self.arguments[i] = new_name
+                self.mapping[new_name] = self.mapping[old_name]
+                self.mapping[new_name].value = new_name
+                del self.mapping[old_name]
+
+    def _add(self, prim):
+        def addType(dict_, ret_type):
+            if not ret_type in dict_:
+                new_list = []
+                for type_, list_ in dict_.items():
+                    if issubclass(type_, ret_type):
+                        for item in list_:
+                            if not item in new_list:
+                                new_list.append(item)
+                dict_[ret_type] = new_list
+
+        addType(self.primitives, prim.ret)
+        addType(self.terminals, prim.ret)
+
+        self.mapping[prim.name] = prim
+        if isinstance(prim, Primitive):
+            for type_ in prim.args:
+                addType(self.primitives, type_)
+                addType(self.terminals, type_)
+            dict_ = self.primitives
+        else:
+            dict_ = self.terminals
+
+        for type_ in dict_:
+            if issubclass(prim.ret, type_):
+                dict_[type_].append(prim)
+
+    def addPrimitive(self, primitive, in_types, ret_type, name=None):
+        """Add a primitive to the set.
+
+        :param primitive: callable object or a function.
+        :parma in_types: list of primitives arguments' type
+        :param ret_type: type returned by the primitive.
+        :param name: alternative name for the primitive instead
+                     of its __name__ attribute.
+        """
+        if name is None:
+            name = primitive.__name__
+        prim = Primitive(name, in_types, ret_type)
+
+        assert name not in self.context or \
+               self.context[name] is primitive, \
+               "Primitives are required to have a unique name. " \
+               "Consider using the argument 'name' to rename your "\
+               "second '%s' primitive." % (name,)
+
+        self._add(prim)
+        self.context[prim.name] = primitive
+        self.prims_count += 1
+
+    def addTerminal(self, terminal, ret_type, name=None):
+        """Add a terminal to the set. Terminals can be named
+        using the optional *name* argument. This should be
+        used : to define named constant (i.e.: pi); to speed the
+        evaluation time when the object is long to build; when
+        the object does not have a __repr__ functions that returns
+        the code to build the object; when the object class is
+        not a Python built-in.
+
+        :param terminal: Object, or a function with no arguments.
+        :param ret_type: Type of the terminal.
+        :param name: defines the name of the terminal in the expression.
+        """
+        symbolic = False
+        if name is None and callable(terminal):
+            name = terminal.__name__
+
+        assert name not in self.context, \
+               "Terminals are required to have a unique name. " \
+               "Consider using the argument 'name' to rename your "\
+               "second %s terminal." % (name,)
+
+        if name is not None:
+            self.context[name] = terminal
+            terminal = name
+            symbolic = True
+        elif terminal in (True, False):
+            # To support True and False terminals with Python 2.
+            self.context[str(terminal)] = terminal
+
+        prim = Terminal(terminal, symbolic, ret_type)
+        self._add(prim)
+        self.terms_count += 1
+
+    def addEphemeralConstant(self, name, ephemeral, ret_type):
+        """Add an ephemeral constant to the set. An ephemeral constant
+        is a no argument function that returns a random value. The value
+        of the constant is constant for a Tree, but may differ from one
+        Tree to another.
+
+        :param name: name used to refers to this ephemeral type.
+        :param ephemeral: function with no arguments returning a random value.
+        :param ret_type: type of the object returned by *ephemeral*.
+        """
+        module_gp = globals()
+        if not name in module_gp:
+            class_ = type(name, (Ephemeral,), {'func' : staticmethod(ephemeral),
+                                               'ret' : ret_type})
+            module_gp[name] = class_
+        else:
+            class_ = module_gp[name]
+            if issubclass(class_, Ephemeral):
+                if class_.func is not ephemeral:
+                    raise Exception("Ephemerals with different functions should "
+                                    "be named differently, even between psets.")
+                elif class_.ret is not ret_type:
+                    raise Exception("Ephemerals with the same name and function "
+                                    "should have the same type, even between psets.")
+            else:
+                raise Exception("Ephemerals should be named differently "
+                                "than classes defined in the gp module.")
+        
+        self._add(class_)
+        self.terms_count += 1
+
+    def addADF(self, adfset):
+        """Add an Automatically Defined Function (ADF) to the set.
+
+        :param adfset: PrimitiveSetTyped containing the primitives with which
+                       the ADF can be built.
+        """
+        prim = Primitive(adfset.name, adfset.ins, adfset.ret)
+        self._add(prim)
+        self.prims_count += 1
+
+    @property
+    def terminalRatio(self):
+        """Return the ratio of the number of terminals on the number of all
+        kind of primitives.
+        """
+        return self.terms_count / float(self.terms_count + self.prims_count)
+
+class PrimitiveSet(PrimitiveSetTyped):
+    """Class same as :class:`~deap.gp.PrimitiveSetTyped`, except there is no
+    definition of type.
+    """
+    def __init__(self, name, arity, prefix="ARG"):
+        args = [__type__]*arity
+        PrimitiveSetTyped.__init__(self, name, args, __type__, prefix)
+
+    def addPrimitive(self, primitive, arity, name=None):
+        """Add primitive *primitive* with arity *arity* to the set.
+        If a name *name* is provided, it will replace the attribute __name__
+        attribute to represent/identify the primitive.
+        """
+        assert arity > 0, "arity should be >= 1"
+        args = [__type__] * arity
+        PrimitiveSetTyped.addPrimitive(self, primitive, args, __type__, name)
+
+    def addTerminal(self, terminal, name=None):
+        """Add a terminal to the set."""
+        PrimitiveSetTyped.addTerminal(self, terminal, __type__, name)
+
+    def addEphemeralConstant(self, name, ephemeral):
+        """Add an ephemeral constant to the set."""
+        PrimitiveSetTyped.addEphemeralConstant(self, name, ephemeral, __type__)
+
+
+######################################
+# GP Tree compilation functions      #
+######################################
+def compile(expr, pset):
+    """Compile the expression *expr*.
+
+    :param expr: Expression to compile. It can either be a PrimitiveTree,
+                 a string of Python code or any object that when
+                 converted into string produced a valid Python code
+                 expression.
+    :param pset: Primitive set against which the expression is compile.
+    :returns: a function if the primitive set has 1 or more arguments,
+              or return the results produced by evaluating the tree.
+    """
+    code = str(expr)
+    if len(pset.arguments) > 0:
+        # This section is a stripped version of the lambdify
+        # function of SymPy 0.6.6.
+        args = ",".join(arg for arg in pset.arguments)
+        code = "lambda {args}: {code}".format(args=args, code=code)
+    try:
+        return eval(code, pset.context, {})
+    except MemoryError:
+        _, _, traceback = sys.exc_info()
+        raise MemoryError, ("DEAP : Error in tree evaluation :"
+        " Python cannot evaluate a tree higher than 90. "
+        "To avoid this problem, you should use bloat control on your "
+        "operators. See the DEAP documentation for more information. "
+        "DEAP will now abort."), traceback
+
+def compileADF(expr, psets):
+    """Compile the expression represented by a list of trees. The first
+    element of the list is the main tree, and the following elements are
+    automatically defined functions (ADF) that can be called by the first
+    tree.
+
+
+    :param expr: Expression to compile. It can either be a PrimitiveTree,
+                 a string of Python code or any object that when
+                 converted into string produced a valid Python code
+                 expression.
+    :param psets: List of primitive sets. Each set corresponds to an ADF
+                  while the last set is associated with the expression
+                  and should contain reference to the preceding ADFs.
+    :returns: a function if the main primitive set has 1 or more arguments,
+              or return the results produced by evaluating the tree.
+    """
+    adfdict = {}
+    func = None
+    for pset, subexpr in reversed(zip(psets, expr)):
+        pset.context.update(adfdict)
+        func = compile(subexpr, pset)
+        adfdict.update({pset.name : func})
+    return func
+
+######################################
+# GP Program generation functions    #
+######################################
+def genFull(pset, min_, max_, type_=None):
+    """Generate an expression where each leaf has a the same depth
+    between *min* and *max*.
+
+    :param pset: Primitive set from which primitives are selected.
+    :param min_: Minimum height of the produced trees.
+    :param max_: Maximum Height of the produced trees.
+    :param type_: The type that should return the tree when called, when
+                  :obj:`None` (default) no return type is enforced.
+    :returns: A full tree with all leaves at the same depth.
+    """
+    def condition(height, depth):
+        """Expression generation stops when the depth is equal to height."""
+        return depth == height
+    return generate(pset, min_, max_, condition, type_)
+
+def genGrow(pset, min_, max_, type_=None):
+    """Generate an expression where each leaf might have a different depth
+    between *min* and *max*.
+
+    :param pset: Primitive set from which primitives are selected.
+    :param min_: Minimum height of the produced trees.
+    :param max_: Maximum Height of the produced trees.
+    :param type_: The type that should return the tree when called, when
+                  :obj:`None` (default) no return type is enforced.
+    :returns: A grown tree with leaves at possibly different depths.
+    """
+    def condition(height, depth):
+        """Expression generation stops when the depth is equal to height
+        or when it is randomly determined that a a node should be a terminal.
+        """
+        return depth == height or \
+               (depth >= min_ and random.random() < pset.terminalRatio)
+    return generate(pset, min_, max_, condition, type_)
+
+def genHalfAndHalf(pset, min_, max_, type_=None):
+    """Generate an expression with a PrimitiveSet *pset*.
+    Half the time, the expression is generated with :func:`~deap.gp.genGrow`,
+    the other half, the expression is generated with :func:`~deap.gp.genFull`.
+
+    :param pset: Primitive set from which primitives are selected.
+    :param min_: Minimum height of the produced trees.
+    :param max_: Maximum Height of the produced trees.
+    :param type_: The type that should return the tree when called, when
+                  :obj:`None` (default) no return type is enforced.
+    :returns: Either, a full or a grown tree.
+    """
+    method = random.choice((genGrow, genFull))
+    return method(pset, min_, max_, type_)
+
+def genRamped(pset, min_, max_, type_=None):
+    """
+    .. deprecated:: 1.0
+        The function has been renamed. Use :func:`~deap.gp.genHalfAndHalf` instead.
+    """
+    warnings.warn("gp.genRamped has been renamed. Use genHalfAndHalf instead.",
+                  FutureWarning)
+    return genHalfAndHalf(pset, min_, max_, type_)
+
+def generate(pset, min_, max_, condition, type_=None):
+    """Generate a Tree as a list of list. The tree is build
+    from the root to the leaves, and it stop growing when the
+    condition is fulfilled.
+
+    :param pset: Primitive set from which primitives are selected.
+    :param min_: Minimum height of the produced trees.
+    :param max_: Maximum Height of the produced trees.
+    :param condition: The condition is a function that takes two arguments,
+                      the height of the tree to build and the current
+                      depth in the tree.
+    :param type_: The type that should return the tree when called, when
+                  :obj:`None` (default) no return type is enforced.
+    :returns: A grown tree with leaves at possibly different depths
+              dependending on the condition function.
+    """
+    if type_ is None:
+        type_ = pset.ret
+    expr = []
+    height = random.randint(min_, max_)
+    stack = [(0, type_)]
+    while len(stack) != 0:
+        depth, type_ = stack.pop()
+        if condition(height, depth):
+            try:
+                term = random.choice(pset.terminals[type_])
+            except IndexError:
+                _, _, traceback = sys.exc_info()
+                raise IndexError, "The gp.generate function tried to add "\
+                                  "a terminal of type '%s', but there is "\
+                                  "none available." % (type_,), traceback
+            if isclass(term):
+                term = term()
+            expr.append(term)
+        else:
+            try:
+                prim = random.choice(pset.primitives[type_])
+            except IndexError:
+                _, _, traceback = sys.exc_info()
+                raise IndexError, "The gp.generate function tried to add "\
+                                  "a primitive of type '%s', but there is "\
+                                  "none available." % (type_,), traceback
+            expr.append(prim)
+            for arg in reversed(prim.args):
+                stack.append((depth+1, arg))
+    return expr
+
+
+######################################
+# GP Crossovers                      #
+######################################
+
+def cxOnePoint(ind1, ind2):
+    """Randomly select in each individual and exchange each subtree with the
+    point as root between each individual.
+
+    :param ind1: First tree participating in the crossover.
+    :param ind2: Second tree participating in the crossover.
+    :returns: A tuple of two trees.
+    """
+    if len(ind1) < 2 or len(ind2) < 2:
+        # No crossover on single node tree
+        return ind1, ind2
+
+    # List all available primitive types in each individual
+    types1 = defaultdict(list)
+    types2 = defaultdict(list)
+    if ind1.root.ret == __type__:
+        # Not STGP optimization
+        types1[__type__] = xrange(1, len(ind1))
+        types2[__type__] = xrange(1, len(ind2))
+        common_types = [__type__]
+    else:
+        for idx, node in enumerate(ind1[1:], 1):
+            types1[node.ret].append(idx)
+        for idx, node in enumerate(ind2[1:], 1):
+            types2[node.ret].append(idx)
+        common_types = set(types1.keys()).intersection(set(types2.keys()))
+
+    if len(common_types) > 0:
+        type_ = random.choice(list(common_types))
+
+        index1 = random.choice(types1[type_])
+        index2 = random.choice(types2[type_])
+
+        slice1 = ind1.searchSubtree(index1)
+        slice2 = ind2.searchSubtree(index2)
+        ind1[slice1], ind2[slice2] = ind2[slice2], ind1[slice1]
+
+    return ind1, ind2
+
+
+def cxOnePointLeafBiased(ind1, ind2, termpb):
+    """Randomly select crossover point in each individual and exchange each
+    subtree with the point as root between each individual.
+
+    :param ind1: First typed tree participating in the crossover.
+    :param ind2: Second typed tree participating in the crossover.
+    :param termpb: The probability of chosing a terminal node (leaf).
+    :returns: A tuple of two typed trees.
+
+    When the nodes are strongly typed, the operator makes sure the
+    second node type corresponds to the first node type.
+
+    The parameter *termpb* sets the probability to choose between a terminal
+    or non-terminal crossover point. For instance, as defined by Koza, non-
+    terminal primitives are selected for 90% of the crossover points, and
+    terminals for 10%, so *termpb* should be set to 0.1.
+    """
+
+    if len(ind1) < 2 or len(ind2) < 2:
+        # No crossover on single node tree
+        return ind1, ind2
+
+    # Determine wether we keep terminals or primitives for each individual
+    terminal_op = partial(eq, 0)
+    primitive_op = partial(lt, 0)
+    arity_op1 = terminal_op if random.random() < termpb else primitive_op
+    arity_op2 = terminal_op if random.random() < termpb else primitive_op
+
+    # List all available primitive or terminal types in each individual
+    types1 = defaultdict(list)
+    types2 = defaultdict(list)
+
+    for idx, node in enumerate(ind1[1:], 1):
+        if arity_op1(node.arity):
+            types1[node.ret].append(idx)
+
+    for idx, node in enumerate(ind2[1:], 1):
+        if arity_op2(node.arity):
+            types2[node.ret].append(idx)
+
+    common_types = set(types1.keys()).intersection(set(types2.keys()))
+
+    if len(common_types) > 0:
+        # Set does not support indexing
+        type_ = random.sample(common_types, 1)[0]
+        index1 = random.choice(types1[type_])
+        index2 = random.choice(types2[type_])
+
+        slice1 = ind1.searchSubtree(index1)
+        slice2 = ind2.searchSubtree(index2)
+        ind1[slice1], ind2[slice2] = ind2[slice2], ind1[slice1]
+
+    return ind1, ind2
+
+
+######################################
+# GP Mutations                       #
+######################################
+def mutUniform(individual, expr, pset):
+    """Randomly select a point in the tree *individual*, then replace the
+    subtree at that point as a root by the expression generated using method
+    :func:`expr`.
+
+    :param individual: The tree to be mutated.
+    :param expr: A function object that can generate an expression when
+                 called.
+    :returns: A tuple of one tree.
+    """
+    index = random.randrange(len(individual))
+    slice_ = individual.searchSubtree(index)
+    type_ = individual[index].ret
+    individual[slice_] = expr(pset=pset, type_=type_)
+    return individual,
+
+
+def mutNodeReplacement(individual, pset):
+    """Replaces a randomly chosen primitive from *individual* by a randomly
+    chosen primitive with the same number of arguments from the :attr:`pset`
+    attribute of the individual.
+
+    :param individual: The normal or typed tree to be mutated.
+    :returns: A tuple of one tree.
+    """
+    if len(individual) < 2:
+        return individual,
+
+    index = random.randrange(1, len(individual))
+    node = individual[index]
+
+    if node.arity == 0: # Terminal
+        term = random.choice(pset.terminals[node.ret])
+        if isclass(term):
+            term = term()
+        individual[index] = term
+    else:   # Primitive
+        prims = [p for p in pset.primitives[node.ret] if p.args == node.args]
+        individual[index] = random.choice(prims)
+
+    return individual,
+
+def mutEphemeral(individual, mode):
+    """This operator works on the constants of the tree *individual*. In
+    *mode* ``"one"``, it will change the value of one of the individual
+    ephemeral constants by calling its generator function. In *mode*
+    ``"all"``, it will change the value of **all** the ephemeral constants.
+
+    :param individual: The normal or typed tree to be mutated.
+    :param mode: A string to indicate to change ``"one"`` or ``"all"``
+                 ephemeral constants.
+    :returns: A tuple of one tree.
+    """
+    if mode not in ["one", "all"]:
+        raise ValueError("Mode must be one of \"one\" or \"all\"")
+
+    ephemerals_idx = [index
+                      for index, node in enumerate(individual)
+                      if isinstance(node, Ephemeral)]
+
+    if len(ephemerals_idx) > 0:
+        if mode == "one":
+            ephemerals_idx = (random.choice(ephemerals_idx),)
+
+        for i in ephemerals_idx:
+            individual[i] = type(individual[i])()
+
+    return individual,
+
+def mutInsert(individual, pset):
+    """Inserts a new branch at a random position in *individual*. The subtree
+    at the chosen position is used as child node of the created subtree, in
+    that way, it is really an insertion rather than a replacement. Note that
+    the original subtree will become one of the children of the new primitive
+    inserted, but not perforce the first (its position is randomly selected if
+    the new primitive has more than one child).
+
+    :param individual: The normal or typed tree to be mutated.
+    :returns: A tuple of one tree.
+    """
+    index = random.randrange(len(individual))
+    node = individual[index]
+    slice_ = individual.searchSubtree(index)
+    choice = random.choice
+
+    # As we want to keep the current node as children of the new one,
+    # it must accept the return value of the current node
+    primitives = [p for p in pset.primitives[node.ret] if node.ret in p.args]
+
+    if len(primitives) == 0:
+        return individual,
+
+    new_node = choice(primitives)
+    new_subtree = [None] * len(new_node.args)
+    position = choice([i for i, a in enumerate(new_node.args) if a == node.ret])
+
+    for i, arg_type in enumerate(new_node.args):
+        if i != position:
+            term = choice(pset.terminals[arg_type])
+            if isclass(term):
+                term = term()
+            new_subtree[i] = term
+
+    new_subtree[position:position+1] = individual[slice_]
+    new_subtree.insert(0, new_node)
+    individual[slice_] = new_subtree
+    return individual,
+
+def mutShrink(individual):
+    """This operator shrinks the *individual* by chosing randomly a branch and
+    replacing it with one of the branch's arguments (also randomly chosen).
+
+    :param individual: The tree to be shrinked.
+    :returns: A tuple of one tree.
+    """
+    # We don't want to "shrink" the root
+    if len(individual) < 3 or individual.height <= 1:
+        return individual,
+
+    iprims = []
+    for i, node in enumerate(individual[1:], 1):
+        if isinstance(node, Primitive) and node.ret in node.args:
+            iprims.append((i, node))
+
+    if len(iprims) != 0:
+        index, prim = random.choice(iprims)
+        arg_idx = random.choice([i for i, type_ in enumerate(prim.args) if type_ == prim.ret])
+        rindex = index+1
+        for _ in range(arg_idx+1):
+            rslice = individual.searchSubtree(rindex)
+            subtree = individual[rslice]
+            rindex += len(subtree)
+
+        slice_ = individual.searchSubtree(index)
+        individual[slice_] = subtree
+
+    return individual,
+
+######################################
+# GP bloat control decorators        #
+######################################
+
+def staticLimit(key, max_value):
+    """Implement a static limit on some measurement on a GP tree, as defined
+    by Koza in [Koza1989]. It may be used to decorate both crossover and
+    mutation operators. When an invalid (over the limit) child is generated,
+    it is simply replaced by one of its parents, randomly selected.
+
+    This operator can be used to avoid memory errors occuring when the tree
+    gets higher than 90 levels (as Python puts a limit on the call stack
+    depth), because it can ensure that no tree higher than this limit will ever
+    be accepted in the population, except if it was generated at initialization
+    time.
+    
+    :param key: The function to use in order the get the wanted value. For 
+                instance, on a GP tree, ``operator.attrgetter('height')`` may
+                be used to set a depth limit, and ``len`` to set a size limit.
+    :param max_value: The maximum value allowed for the given measurement.
+    :returns: A decorator that can be applied to a GP operator using \
+    :func:`~deap.base.Toolbox.decorate`
+
+    .. note::
+       If you want to reproduce the exact behavior intended by Koza, set
+       *key* to ``operator.attrgetter('height')`` and *max_value* to 17.
+
+    .. [Koza1989] J.R. Koza, Genetic Programming - On the Programming of
+        Computers by Means of Natural Selection (MIT Press,
+        Cambridge, MA, 1992)
+
+    """
+    def decorator(func):
+        @wraps(func)
+        def wrapper(*args, **kwargs):
+            keep_inds = [copy.deepcopy(ind) for ind in args]
+            new_inds = list(func(*args, **kwargs))
+            for i, ind in enumerate(new_inds):
+                if key(ind) > max_value:
+                    new_inds[i] = random.choice(keep_inds)
+            return new_inds
+        return wrapper
+    return decorator
+
+def graph(expr):
+    """Construct the graph of a tree expression. The tree expression must be
+    valid. It returns in order a node list, an edge list, and a dictionary of
+    the per node labels. The node are represented by numbers, the edges are
+    tuples connecting two nodes (number), and the labels are values of a
+    dictionary for which keys are the node numbers.
+
+    :param expr: A tree expression to convert into a graph.
+    :returns: A node list, an edge list, and a dictionary of labels.
+
+    The returned objects can be used directly to populate a
+    `pygraphviz <http://networkx.lanl.gov/pygraphviz/>`_ graph::
+
+        import pygraphviz as pgv
+
+        # [...] Execution of code that produce a tree expression
+
+        nodes, edges, labels = graph(expr)
+
+        g = pgv.AGraph()
+        g.add_nodes_from(nodes)
+        g.add_edges_from(edges)
+        g.layout(prog="dot")
+
+        for i in nodes:
+            n = g.get_node(i)
+            n.attr["label"] = labels[i]
+
+        g.draw("tree.pdf")
+
+    or a `NetworX <http://networkx.github.com/>`_ graph::
+
+        import matplotlib.pyplot as plt
+        import networkx as nx
+
+        # [...] Execution of code that produce a tree expression
+
+        nodes, edges, labels = graph(expr)
+
+        g = nx.Graph()
+        g.add_nodes_from(nodes)
+        g.add_edges_from(edges)
+        pos = nx.graphviz_layout(g, prog="dot")
+
+        nx.draw_networkx_nodes(g, pos)
+        nx.draw_networkx_edges(g, pos)
+        nx.draw_networkx_labels(g, pos, labels)
+        plt.show()
+
+
+    .. note::
+
+       We encourage you to use `pygraphviz
+       <http://networkx.lanl.gov/pygraphviz/>`_ as the nodes might be plotted
+       out of order when using `NetworX <http://networkx.github.com/>`_.
+    """
+    nodes = range(len(expr))
+    edges = list()
+    labels = dict()
+
+    stack = []
+    for i, node in enumerate(expr):
+        if stack:
+            edges.append((stack[-1][0], i))
+            stack[-1][1] -= 1
+        labels[i] = node.name if isinstance(node, Primitive) else node.value
+        stack.append([i, node.arity])
+        while stack and stack[-1][1] == 0:
+            stack.pop()
+
+    return nodes, edges, labels
+
+if __name__ == "__main__":
+    import doctest
+    doctest.testmod()

python/isaac/autotuning/external/deap/tests/__init__.py

@@ -0,0 +1,14 @@
+#    This file is part of DEAP.
+#
+#    DEAP is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU Lesser General Public License as
+#    published by the Free Software Foundation, either version 3 of
+#    the License, or (at your option) any later version.
+#
+#    DEAP is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#    GNU Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.

python/isaac/autotuning/external/deap/tests/test_logbook.py

@@ -0,0 +1,46 @@
+import sys
+import unittest
+
+sys.path.append("..")
+
+import tools
+
+class LogbookTest(unittest.TestCase):
+
+    def setUp(self):
+        self.logbook = tools.Logbook()
+        print
+
+    def test_multi_chapters(self):
+        self.logbook.record(gen=0, evals=100, fitness={'obj 1' : {'avg' : 1.0, 'max' : 10},
+                                          'avg' : 1.0, 'max' : 10},
+                            length={'avg' : 1.0, 'max' : 30},
+                            test={'avg' : 1.0, 'max' : 20})
+        self.logbook.record(gen=0, evals=100, fitness={'obj 1' : {'avg' : 1.0, 'max' : 10},
+                                          'avg' : 1.0, 'max' : 10},
+                            length={'avg' : 1.0, 'max' : 30},
+                            test={'avg' : 1.0, 'max' : 20})
+        print(self.logbook.stream)
+
+
+    def test_one_chapter(self):
+        self.logbook.record(gen=0, evals=100, fitness={'avg' : 1.0, 'max' : 10})
+        self.logbook.record(gen=0, evals=100, fitness={'avg' : 1.0, 'max' : 10})
+        print(self.logbook.stream)
+
+    def test_one_big_chapter(self):
+        self.logbook.record(gen=0, evals=100, fitness={'obj 1' : {'avg' : 1.0, 'max' : 10}, 'obj 2' : {'avg' : 1.0, 'max' : 10}})
+        self.logbook.record(gen=0, evals=100, fitness={'obj 1' : {'avg' : 1.0, 'max' : 10}, 'obj 2' : {'avg' : 1.0, 'max' : 10}})
+        print(self.logbook.stream)
+
+    def test_no_chapters(self):
+        self.logbook.record(gen=0, evals=100, **{'avg' : 1.0, 'max' : 10})
+        self.logbook.record(gen=0, evals=100, **{'avg' : 1.0, 'max' : 10})        
+        print(self.logbook.stream)
+
+
+
+if __name__ == "__main__":
+    suite = unittest.TestLoader().loadTestsFromTestCase(LogbookTest)
+    unittest.TextTestRunner(verbosity=2).run(suite)
+    

python/isaac/autotuning/external/deap/tests/test_pickle.py

@@ -0,0 +1,115 @@
+
+import sys
+import unittest
+import array
+import pickle
+import operator
+from test import test_support
+
+sys.path.append("..")
+
+import numpy
+
+import creator
+import base
+import gp
+import tools
+
+def func():
+    return "True"
+
+class Pickling(unittest.TestCase):
+
+    def setUp(self):
+        creator.create("FitnessMax", base.Fitness, weights=(1.0,))
+        creator.create("IndList", list, fitness=creator.FitnessMax)
+        creator.create("IndArray", array.array,  typecode='f', fitness=creator.FitnessMax)
+        creator.create("IndTree", gp.PrimitiveTree, fitness=creator.FitnessMax)
+        self.toolbox = base.Toolbox()
+        self.toolbox.register("func", func)
+        self.toolbox.register("lambda_func", lambda: "True")
+    
+    def test_pickle_fitness(self):
+        fitness = creator.FitnessMax()
+        fitness.values = (1.0,)
+        fitness_s = pickle.dumps(fitness)
+        fitness_l = pickle.loads(fitness_s)
+        self.failUnlessEqual(fitness, fitness_l, "Unpickled fitness != pickled fitness")
+
+    def test_pickle_ind_list(self):
+        ind = creator.IndList([1.0, 2.0, 3.0])
+        ind.fitness.values = (4.0,)
+        ind_s = pickle.dumps(ind)
+        ind_l = pickle.loads(ind_s)
+        self.failUnlessEqual(ind, ind_l, "Unpickled individual list != pickled individual list")
+        self.failUnlessEqual(ind.fitness, ind_l.fitness, "Unpickled individual fitness != pickled individual fitness")
+         
+    def test_pickle_ind_array(self):
+        ind = creator.IndArray([1.0, 2.0, 3.0])
+        ind.fitness.values = (4.0,)
+        ind_s = pickle.dumps(ind)
+        ind_l = pickle.loads(ind_s)
+        self.failUnlessEqual(ind, ind_l, "Unpickled individual array != pickled individual array")
+        self.failUnlessEqual(ind.fitness, ind_l.fitness, "Unpickled individual fitness != pickled individual fitness")
+    
+    def test_pickle_tree(self):
+        ind = creator.IndTree([operator.add, 1, 2])
+        ind.fitness.values = (1.0,)
+        ind_s = pickle.dumps(ind)
+        ind_l = pickle.loads(ind_s)
+        msg =  "Unpickled individual %s != pickled individual %s" % (str(ind), str(ind_l))
+        self.failUnlessEqual(ind, ind_l, msg)
+        msg =  "Unpickled fitness %s != pickled fitness %s" % (str(ind.fitness), str(ind_l.fitness))
+        self.failUnlessEqual(ind.fitness, ind_l.fitness, msg)
+    
+    def test_pickle_population(self):
+        ind1 = creator.IndList([1.0,2.0,3.0])
+        ind1.fitness.values = (1.0,)
+        ind2 = creator.IndList([4.0,5.0,6.0])
+        ind2.fitness.values = (2.0,)
+        ind3 = creator.IndList([7.0,8.0,9.0])
+        ind3.fitness.values = (3.0,)
+        
+        pop = [ind1, ind2, ind3]
+        
+        pop_s = pickle.dumps(pop)
+        pop_l = pickle.loads(pop_s)
+        
+        self.failUnlessEqual(pop[0], pop_l[0], "Unpickled individual list != pickled individual list")
+        self.failUnlessEqual(pop[0].fitness, pop_l[0].fitness, "Unpickled individual fitness != pickled individual fitness")
+        self.failUnlessEqual(pop[1], pop_l[1], "Unpickled individual list != pickled individual list")
+        self.failUnlessEqual(pop[1].fitness, pop_l[1].fitness, "Unpickled individual fitness != pickled individual fitness")
+        self.failUnlessEqual(pop[2], pop_l[2], "Unpickled individual list != pickled individual list")
+        self.failUnlessEqual(pop[2].fitness, pop_l[2].fitness, "Unpickled individual fitness != pickled individual fitness")
+    
+    def test_pickle_logbook(self):
+        stats = tools.Statistics()
+        logbook = tools.Logbook()
+
+        stats.register("mean", numpy.mean)
+        record = stats.compile([1,2,3,4,5,6,8,9,10])
+        logbook.record(**record)
+
+        stats_s = pickle.dumps(logbook)
+        logbook_r = pickle.loads(stats_s)
+
+        self.failUnlessEqual(logbook, logbook_r, "Unpickled logbook != pickled logbook")
+
+
+    if not sys.version_info < (2, 7):
+        def test_pickle_partial(self):
+            func_s = pickle.dumps(self.toolbox.func)
+            func_l = pickle.loads(func_s)
+
+            self.failUnlessEqual(self.toolbox.func(), func_l())
+        
+        @unittest.expectedFailure
+        def test_pickle_lambda(self):
+            func_s = pickle.dumps(self.toolbox.lambda_func)
+            func_l = pickle.loads(func_s)
+        
+            self.failUnlessEqual(self.toolbox.lambda_func(), func_l())
+
+if __name__ == "__main__":
+    suite = unittest.TestLoader().loadTestsFromTestCase(Pickling)
+    unittest.TextTestRunner(verbosity=2).run(suite)

python/isaac/autotuning/external/deap/tools/__init__.py

@@ -0,0 +1,30 @@
+#    This file is part of DEAP.
+#
+#    DEAP is free software: you can redistribute it and/or modify
+#    it under the terms of the GNU Lesser General Public License as
+#    published by the Free Software Foundation, either version 3 of
+#    the License, or (at your option) any later version.
+#
+#    DEAP is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+#    GNU Lesser General Public License for more details.
+#
+#    You should have received a copy of the GNU Lesser General Public
+#    License along with DEAP. If not, see <http://www.gnu.org/licenses/>.
+"""The :mod:`~deap.tools` module contains the operators for evolutionary
+algorithms. They are used to modify, select and move the individuals in their
+environment. The set of operators it contains are readily usable in the
+:class:`~deap.base.Toolbox`. In addition to the basic operators this module
+also contains utility tools to enhance the basic algorithms with
+:class:`Statistics`, :class:`HallOfFame`, :class:`Checkpoint`, and
+:class:`History`.
+"""
+
+from .crossover import *
+from .emo import *
+from .init import *
+from .migration import *
+from .mutation import *
+from .selection import *
+from .support import *

python/isaac/autotuning/external/deap/tools/crossover.py

@@ -0,0 +1,445 @@
+from __future__ import division
+import random
+import warnings
+
+from collections import Sequence
+from itertools import repeat
+
+######################################
+# GA Crossovers                      #
+######################################
+
+def cxOnePoint(ind1, ind2):
+    """Executes a one point crossover on the input :term:`sequence` individuals.
+    The two individuals are modified in place. The resulting individuals will
+    respectively have the length of the other.
+    
+    :param ind1: The first individual participating in the crossover.
+    :param ind2: The second individual participating in the crossover.
+    :returns: A tuple of two individuals.
+
+    This function uses the :func:`~random.randint` function from the
+    python base :mod:`random` module.
+    """
+    size = min(len(ind1), len(ind2))
+    cxpoint = random.randint(1, size - 1)
+    ind1[cxpoint:], ind2[cxpoint:] = ind2[cxpoint:], ind1[cxpoint:]
+    
+    return ind1, ind2
+
+def cxTwoPoint(ind1, ind2):
+    """Executes a two-point crossover on the input :term:`sequence`
+    individuals. The two individuals are modified in place and both keep
+    their original length. 
+    
+    :param ind1: The first individual participating in the crossover.
+    :param ind2: The second individual participating in the crossover.
+    :returns: A tuple of two individuals.
+
+    This function uses the :func:`~random.randint` function from the Python 
+    base :mod:`random` module.
+    """
+    size = min(len(ind1), len(ind2))
+    cxpoint1 = random.randint(1, size)
+    cxpoint2 = random.randint(1, size - 1)
+    if cxpoint2 >= cxpoint1:
+        cxpoint2 += 1
+    else: # Swap the two cx points
+        cxpoint1, cxpoint2 = cxpoint2, cxpoint1
+   
+    ind1[cxpoint1:cxpoint2], ind2[cxpoint1:cxpoint2] \
+        = ind2[cxpoint1:cxpoint2], ind1[cxpoint1:cxpoint2]
+        
+    return ind1, ind2
+
+def cxTwoPoints(ind1, ind2):
+    """
+    .. deprecated:: 1.0
+       The function has been renamed.  Use :func:`~deap.tools.cxTwoPoint` instead.
+    """
+    warnings.warn("tools.cxTwoPoints has been renamed. Use cxTwoPoint instead.",
+                  FutureWarning)
+    return cxTwoPoint(ind1, ind2)
+
+def cxUniform(ind1, ind2, indpb):
+    """Executes a uniform crossover that modify in place the two
+    :term:`sequence` individuals. The attributes are swapped accordingto the
+    *indpb* probability.
+    
+    :param ind1: The first individual participating in the crossover.
+    :param ind2: The second individual participating in the crossover.
+    :param indpb: Independent probabily for each attribute to be exchanged.
+    :returns: A tuple of two individuals.
+    
+    This function uses the :func:`~random.random` function from the python base
+    :mod:`random` module.
+    """
+    size = min(len(ind1), len(ind2))    
+    for i in xrange(size):
+        if random.random() < indpb:
+            ind1[i], ind2[i] = ind2[i], ind1[i]
+    
+    return ind1, ind2
+    
+def cxPartialyMatched(ind1, ind2):
+    """Executes a partially matched crossover (PMX) on the input individuals.
+    The two individuals are modified in place. This crossover expects
+    :term:`sequence` individuals of indices, the result for any other type of
+    individuals is unpredictable.
+    
+    :param ind1: The first individual participating in the crossover.
+    :param ind2: The second individual participating in the crossover.
+    :returns: A tuple of two individuals.
+
+    Moreover, this crossover generates two children by matching
+    pairs of values in a certain range of the two parents and swapping the values
+    of those indexes. For more details see [Goldberg1985]_.
+
+    This function uses the :func:`~random.randint` function from the python base
+    :mod:`random` module.
+    
+    .. [Goldberg1985] Goldberg and Lingel, "Alleles, loci, and the traveling
+       salesman problem", 1985.
+    """
+    size = min(len(ind1), len(ind2))
+    p1, p2 = [0]*size, [0]*size
+
+    # Initialize the position of each indices in the individuals
+    for i in xrange(size):
+        p1[ind1[i]] = i
+        p2[ind2[i]] = i
+    # Choose crossover points
+    cxpoint1 = random.randint(0, size)
+    cxpoint2 = random.randint(0, size - 1)
+    if cxpoint2 >= cxpoint1:
+        cxpoint2 += 1
+    else: # Swap the two cx points
+        cxpoint1, cxpoint2 = cxpoint2, cxpoint1
+    
+    # Apply crossover between cx points
+    for i in xrange(cxpoint1, cxpoint2):
+        # Keep track of the selected values
+        temp1 = ind1[i]
+        temp2 = ind2[i]
+        # Swap the matched value
+        ind1[i], ind1[p1[temp2]] = temp2, temp1
+        ind2[i], ind2[p2[temp1]] = temp1, temp2
+        # Position bookkeeping
+        p1[temp1], p1[temp2] = p1[temp2], p1[temp1]
+        p2[temp1], p2[temp2] = p2[temp2], p2[temp1]
+    
+    return ind1, ind2
+
+def cxUniformPartialyMatched(ind1, ind2, indpb):
+    """Executes a uniform partially matched crossover (UPMX) on the input
+    individuals. The two individuals are modified in place. This crossover
+    expects :term:`sequence` individuals of indices, the result for any other
+    type of individuals is unpredictable.
+    
+    :param ind1: The first individual participating in the crossover.
+    :param ind2: The second individual participating in the crossover.
+    :returns: A tuple of two individuals.
+
+    Moreover, this crossover generates two children by matching
+    pairs of values chosen at random with a probability of *indpb* in the two
+    parents and swapping the values of those indexes. For more details see
+    [Cicirello2000]_.
+
+    This function uses the :func:`~random.random` and :func:`~random.randint`
+    functions from the python base :mod:`random` module.
+    
+    .. [Cicirello2000] Cicirello and Smith, "Modeling GA performance for
+       control parameter optimization", 2000.
+    """
+    size = min(len(ind1), len(ind2))
+    p1, p2 = [0]*size, [0]*size
+
+    # Initialize the position of each indices in the individuals
+    for i in xrange(size):
+        p1[ind1[i]] = i
+        p2[ind2[i]] = i
+    
+    for i in xrange(size):
+        if random.random() < indpb:
+            # Keep track of the selected values
+            temp1 = ind1[i]
+            temp2 = ind2[i]
+            # Swap the matched value
+            ind1[i], ind1[p1[temp2]] = temp2, temp1
+            ind2[i], ind2[p2[temp1]] = temp1, temp2
+            # Position bookkeeping
+            p1[temp1], p1[temp2] = p1[temp2], p1[temp1]
+            p2[temp1], p2[temp2] = p2[temp2], p2[temp1]
+    
+    return ind1, ind2
+
+def cxOrdered(ind1, ind2):
+    """Executes an ordered crossover (OX) on the input
+    individuals. The two individuals are modified in place. This crossover
+    expects :term:`sequence` individuals of indices, the result for any other
+    type of individuals is unpredictable.
+    
+    :param ind1: The first individual participating in the crossover.
+    :param ind2: The second individual participating in the crossover.
+    :returns: A tuple of two individuals.
+
+    Moreover, this crossover generates holes in the input
+    individuals. A hole is created when an attribute of an individual is
+    between the two crossover points of the other individual. Then it rotates
+    the element so that all holes are between the crossover points and fills
+    them with the removed elements in order. For more details see
+    [Goldberg1989]_.
+    
+    This function uses the :func:`~random.sample` function from the python base
+    :mod:`random` module.
+    
+    .. [Goldberg1989] Goldberg. Genetic algorithms in search, 
+       optimization and machine learning. Addison Wesley, 1989
+    """
+    size = min(len(ind1), len(ind2))
+    a, b = random.sample(xrange(size), 2)
+    if a > b:
+        a, b = b, a
+
+    holes1, holes2 = [True]*size, [True]*size
+    for i in range(size):
+        if i < a or i > b:
+            holes1[ind2[i]] = False
+            holes2[ind1[i]] = False
+    
+    # We must keep the original values somewhere before scrambling everything
+    temp1, temp2 = ind1, ind2
+    k1 , k2 = b + 1, b + 1
+    for i in range(size):
+        if not holes1[temp1[(i + b + 1) % size]]:
+            ind1[k1 % size] = temp1[(i + b + 1) % size]
+            k1 += 1
+        
+        if not holes2[temp2[(i + b + 1) % size]]:
+            ind2[k2 % size] = temp2[(i + b + 1) % size]
+            k2 += 1
+    
+    # Swap the content between a and b (included)
+    for i in range(a, b + 1):
+        ind1[i], ind2[i] = ind2[i], ind1[i]
+    
+    return ind1, ind2
+
+def cxBlend(ind1, ind2, alpha):
+    """Executes a blend crossover that modify in-place the input individuals.
+    The blend crossover expects :term:`sequence` individuals of floating point
+    numbers.
+    
+    :param ind1: The first individual participating in the crossover.
+    :param ind2: The second individual participating in the crossover.
+    :param alpha: Extent of the interval in which the new values can be drawn
+                  for each attribute on both side of the parents' attributes.
+    :returns: A tuple of two individuals.
+    
+    This function uses the :func:`~random.random` function from the python base
+    :mod:`random` module.
+    """
+    for i, (x1, x2) in enumerate(zip(ind1, ind2)):
+        gamma = (1. + 2. * alpha) * random.random() - alpha
+        ind1[i] = (1. - gamma) * x1 + gamma * x2
+        ind2[i] = gamma * x1 + (1. - gamma) * x2
+
+    return ind1, ind2
+
+def cxSimulatedBinary(ind1, ind2, eta):
+    """Executes a simulated binary crossover that modify in-place the input
+    individuals. The simulated binary crossover expects :term:`sequence`
+    individuals of floating point numbers.
+    
+    :param ind1: The first individual participating in the crossover.
+    :param ind2: The second individual participating in the crossover.
+    :param eta: Crowding degree of the crossover. A high eta will produce
+                children resembling to their parents, while a small eta will
+                produce solutions much more different.
+    :returns: A tuple of two individuals.
+    
+    This function uses the :func:`~random.random` function from the python base
+    :mod:`random` module.
+    """
+    for i, (x1, x2) in enumerate(zip(ind1, ind2)):
+        rand = random.random()
+        if rand <= 0.5:
+            beta = 2. * rand
+        else:
+            beta = 1. / (2. * (1. - rand))
+        beta **= 1. / (eta + 1.)
+        ind1[i] = 0.5 * (((1 + beta) * x1) + ((1 - beta) * x2))
+        ind2[i] = 0.5 * (((1 - beta) * x1) + ((1 + beta) * x2))
+    
+    return ind1, ind2
+
+
+def cxSimulatedBinaryBounded(ind1, ind2, eta, low, up):
+    """Executes a simulated binary crossover that modify in-place the input
+    individuals. The simulated binary crossover expects :term:`sequence`
+    individuals of floating point numbers.
+    
+    :param ind1: The first individual participating in the crossover.
+    :param ind2: The second individual participating in the crossover.
+    :param eta: Crowding degree of the crossover. A high eta will produce
+                children resembling to their parents, while a small eta will
+                produce solutions much more different.
+    :param low: A value or an :term:`python:sequence` of values that is the lower
+                bound of the search space.
+    :param up: A value or an :term:`python:sequence` of values that is the upper
+               bound of the search space.
+    :returns: A tuple of two individuals.
+
+    This function uses the :func:`~random.random` function from the python base
+    :mod:`random` module.
+
+    .. note::
+       This implementation is similar to the one implemented in the 
+       original NSGA-II C code presented by Deb.
+    """
+    size = min(len(ind1), len(ind2))
+    if not isinstance(low, Sequence):
+        low = repeat(low, size)
+    elif len(low) < size:
+        raise IndexError("low must be at least the size of the shorter individual: %d < %d" % (len(low), size))
+    if not isinstance(up, Sequence):
+        up = repeat(up, size)
+    elif len(up) < size:
+        raise IndexError("up must be at least the size of the shorter individual: %d < %d" % (len(up), size))
+    
+    for i, xl, xu in zip(xrange(size), low, up):
+        if random.random() <= 0.5:
+            # This epsilon should probably be changed for 0 since 
+            # floating point arithmetic in Python is safer
+            if abs(ind1[i] - ind2[i]) > 1e-14:
+                x1 = min(ind1[i], ind2[i])
+                x2 = max(ind1[i], ind2[i])
+                rand = random.random()
+                
+                beta = 1.0 + (2.0 * (x1 - xl) / (x2 - x1))
+                alpha = 2.0 - beta**-(eta + 1)
+                if rand <= 1.0 / alpha:
+                    beta_q = (rand * alpha)**(1.0 / (eta + 1))
+                else:
+                    beta_q = (1.0 / (2.0 - rand * alpha))**(1.0 / (eta + 1))
+                
+                c1 = 0.5 * (x1 + x2 - beta_q * (x2 - x1))
+                
+                beta = 1.0 + (2.0 * (xu - x2) / (x2 - x1))
+                alpha = 2.0 - beta**-(eta + 1)
+                if rand <= 1.0 / alpha:
+                    beta_q = (rand * alpha)**(1.0 / (eta + 1))
+                else:
+                    beta_q = (1.0 / (2.0 - rand * alpha))**(1.0 / (eta + 1))
+                c2 = 0.5 * (x1 + x2 + beta_q * (x2 - x1))
+                
+                c1 = min(max(c1, xl), xu)
+                c2 = min(max(c2, xl), xu)
+                
+                if random.random() <= 0.5:
+                    ind1[i] = c2
+                    ind2[i] = c1
+                else:
+                    ind1[i] = c1
+                    ind2[i] = c2
+    
+    return ind1, ind2   
+
+
+######################################
+# Messy Crossovers                   #
+######################################
+
+def cxMessyOnePoint(ind1, ind2):
+    """Executes a one point crossover on :term:`sequence` individual.
+    The crossover will in most cases change the individuals size. The two
+    individuals are modified in place.
+    
+    :param ind1: The first individual participating in the crossover.
+    :param ind2: The second individual participating in the crossover.
+    :returns: A tuple of two individuals.
+    
+    This function uses the :func:`~random.randint` function from the python base
+    :mod:`random` module.        
+    """
+    cxpoint1 = random.randint(0, len(ind1))
+    cxpoint2 = random.randint(0, len(ind2))
+    ind1[cxpoint1:], ind2[cxpoint2:] = ind2[cxpoint2:], ind1[cxpoint1:]
+    
+    return ind1, ind2
+    
+######################################
+# ES Crossovers                      #
+######################################
+
+def cxESBlend(ind1, ind2, alpha):
+    """Executes a blend crossover on both, the individual and the strategy. The
+    individuals shall be a :term:`sequence` and must have a :term:`sequence`
+    :attr:`strategy` attribute. Adjustement of the minimal strategy shall be done
+    after the call to this function, consider using a decorator.
+    
+    :param ind1: The first evolution strategy participating in the crossover.
+    :param ind2: The second evolution strategy participating in the crossover.
+    :param alpha: Extent of the interval in which the new values can be drawn
+                  for each attribute on both side of the parents' attributes.
+    :returns: A tuple of two evolution strategies.
+
+    This function uses the :func:`~random.random` function from the python base
+    :mod:`random` module.  
+    """
+    for i, (x1, s1, x2, s2) in enumerate(zip(ind1, ind1.strategy, 
+                                             ind2, ind2.strategy)):
+        # Blend the values
+        gamma = (1. + 2. * alpha) * random.random() - alpha
+        ind1[i] = (1. - gamma) * x1 + gamma * x2
+        ind2[i] = gamma * x1 + (1. - gamma) * x2
+        # Blend the strategies
+        gamma = (1. + 2. * alpha) * random.random() - alpha
+        ind1.strategy[i] = (1. - gamma) * s1 + gamma * s2
+        ind2.strategy[i] = gamma * s1 + (1. - gamma) * s2
+    
+    return ind1, ind2
+
+def cxESTwoPoint(ind1, ind2):
+    """Executes a classical two points crossover on both the individuals and their
+    strategy. The individuals shall be a :term:`sequence` and must have a
+    :term:`sequence` :attr:`strategy` attribute. The crossover points for the
+    individual and the strategy are the same.
+    
+    :param ind1: The first evolution strategy participating in the crossover.
+    :param ind2: The second evolution strategy participating in the crossover.
+    :returns: A tuple of two evolution strategies.
+    
+    This function uses the :func:`~random.randint` function from the python base
+    :mod:`random` module.  
+    """
+    size = min(len(ind1), len(ind2))
+    
+    pt1 = random.randint(1, size)
+    pt2 = random.randint(1, size - 1)
+    if pt2 >= pt1:
+        pt2 += 1
+    else: # Swap the two cx points
+        pt1, pt2 = pt2, pt1
+   
+    ind1[pt1:pt2], ind2[pt1:pt2] = ind2[pt1:pt2], ind1[pt1:pt2]     
+    ind1.strategy[pt1:pt2], ind2.strategy[pt1:pt2] = \
+        ind2.strategy[pt1:pt2], ind1.strategy[pt1:pt2]
+    
+    return ind1, ind2
+
+def cxESTwoPoints(ind1, ind2):
+    """
+    .. deprecated:: 1.0
+        The function has been renamed. Use :func:`cxESTwoPoint` instead.
+    """
+    return cxESTwoPoints(ind1, ind2)
+
+# List of exported function names.
+__all__ = ['cxOnePoint', 'cxTwoPoint', 'cxUniform', 'cxPartialyMatched',
+           'cxUniformPartialyMatched', 'cxOrdered', 'cxBlend',
+           'cxSimulatedBinary','cxSimulatedBinaryBounded', 'cxMessyOnePoint', 
+           'cxESBlend', 'cxESTwoPoint']
+
+# Deprecated functions
+__all__.extend(['cxTwoPoints', 'cxESTwoPoints'])

python/isaac/autotuning/external/deap/tools/emo.py

@@ -0,0 +1,582 @@
+from __future__ import division
+import bisect
+import math
+import random
+
+from itertools import chain
+from operator import attrgetter, itemgetter
+from collections import defaultdict
+
+######################################
+# Non-Dominated Sorting   (NSGA-II)  #
+######################################
+
+def selNSGA2(individuals, k):
+    """Apply NSGA-II selection operator on the *individuals*. Usually, the
+    size of *individuals* will be larger than *k* because any individual
+    present in *individuals* will appear in the returned list at most once.
+    Having the size of *individuals* equals to *k* will have no effect other
+    than sorting the population according to their front rank. The
+    list returned contains references to the input *individuals*. For more
+    details on the NSGA-II operator see [Deb2002]_.
+    
+    :param individuals: A list of individuals to select from.
+    :param k: The number of individuals to select.
+    :returns: A list of selected individuals.
+    
+    .. [Deb2002] Deb, Pratab, Agarwal, and Meyarivan, "A fast elitist
+       non-dominated sorting genetic algorithm for multi-objective
+       optimization: NSGA-II", 2002.
+    """
+    pareto_fronts = sortNondominated(individuals, k)
+    for front in pareto_fronts:
+        assignCrowdingDist(front)
+    
+    chosen = list(chain(*pareto_fronts[:-1]))
+    k = k - len(chosen)
+    if k > 0:
+        sorted_front = sorted(pareto_fronts[-1], key=attrgetter("fitness.crowding_dist"), reverse=True)
+        chosen.extend(sorted_front[:k])
+        
+    return chosen
+
+def sortNondominated(individuals, k, first_front_only=False):
+    """Sort the first *k* *individuals* into different nondomination levels 
+    using the "Fast Nondominated Sorting Approach" proposed by Deb et al.,
+    see [Deb2002]_. This algorithm has a time complexity of :math:`O(MN^2)`, 
+    where :math:`M` is the number of objectives and :math:`N` the number of 
+    individuals.
+    
+    :param individuals: A list of individuals to select from.
+    :param k: The number of individuals to select.
+    :param first_front_only: If :obj:`True` sort only the first front and
+                             exit.
+    :returns: A list of Pareto fronts (lists), the first list includes 
+              nondominated individuals.
+
+    .. [Deb2002] Deb, Pratab, Agarwal, and Meyarivan, "A fast elitist
+       non-dominated sorting genetic algorithm for multi-objective
+       optimization: NSGA-II", 2002.
+    """
+    if k == 0:
+        return []
+
+    map_fit_ind = defaultdict(list)
+    for ind in individuals:
+        map_fit_ind[ind.fitness].append(ind)
+    fits = map_fit_ind.keys()
+    
+    current_front = []
+    next_front = []
+    dominating_fits = defaultdict(int)
+    dominated_fits = defaultdict(list)
+    
+    # Rank first Pareto front
+    for i, fit_i in enumerate(fits):
+        for fit_j in fits[i+1:]:
+            if fit_i.dominates(fit_j):
+                dominating_fits[fit_j] += 1
+                dominated_fits[fit_i].append(fit_j)
+            elif fit_j.dominates(fit_i):
+                dominating_fits[fit_i] += 1
+                dominated_fits[fit_j].append(fit_i)
+        if dominating_fits[fit_i] == 0:
+            current_front.append(fit_i)
+    
+    fronts = [[]]
+    for fit in current_front:
+        fronts[-1].extend(map_fit_ind[fit])
+    pareto_sorted = len(fronts[-1])
+
+    # Rank the next front until all individuals are sorted or 
+    # the given number of individual are sorted.
+    if not first_front_only:
+        N = min(len(individuals), k)
+        while pareto_sorted < N:
+            fronts.append([])
+            for fit_p in current_front:
+                for fit_d in dominated_fits[fit_p]:
+                    dominating_fits[fit_d] -= 1
+                    if dominating_fits[fit_d] == 0:
+                        next_front.append(fit_d)
+                        pareto_sorted += len(map_fit_ind[fit_d])
+                        fronts[-1].extend(map_fit_ind[fit_d])
+            current_front = next_front
+            next_front = []
+    
+    return fronts
+
+def assignCrowdingDist(individuals):
+    """Assign a crowding distance to each individual's fitness. The 
+    crowding distance can be retrieve via the :attr:`crowding_dist` 
+    attribute of each individual's fitness.
+    """
+    if len(individuals) == 0:
+        return
+    
+    distances = [0.0] * len(individuals)
+    crowd = [(ind.fitness.values, i) for i, ind in enumerate(individuals)]
+    
+    nobj = len(individuals[0].fitness.values)
+    
+    for i in xrange(nobj):
+        crowd.sort(key=lambda element: element[0][i])
+        distances[crowd[0][1]] = float("inf")
+        distances[crowd[-1][1]] = float("inf")
+        if crowd[-1][0][i] == crowd[0][0][i]:
+            continue
+        norm = nobj * float(crowd[-1][0][i] - crowd[0][0][i])
+        for prev, cur, next in zip(crowd[:-2], crowd[1:-1], crowd[2:]):
+            distances[cur[1]] += (next[0][i] - prev[0][i]) / norm
+
+    for i, dist in enumerate(distances):
+        individuals[i].fitness.crowding_dist = dist
+
+def selTournamentDCD(individuals, k):
+    """Tournament selection based on dominance (D) between two individuals, if
+    the two individuals do not interdominate the selection is made
+    based on crowding distance (CD). The *individuals* sequence length has to
+    be a multiple of 4. Starting from the beginning of the selected
+    individuals, two consecutive individuals will be different (assuming all
+    individuals in the input list are unique). Each individual from the input
+    list won't be selected more than twice.
+    
+    This selection requires the individuals to have a :attr:`crowding_dist`
+    attribute, which can be set by the :func:`assignCrowdingDist` function.
+    
+    :param individuals: A list of individuals to select from.
+    :param k: The number of individuals to select.
+    :returns: A list of selected individuals.
+    """
+    def tourn(ind1, ind2):
+        if ind1.fitness.dominates(ind2.fitness):
+            return ind1
+        elif ind2.fitness.dominates(ind1.fitness):
+            return ind2
+
+        if ind1.fitness.crowding_dist < ind2.fitness.crowding_dist:
+            return ind2
+        elif ind1.fitness.crowding_dist > ind2.fitness.crowding_dist:
+            return ind1
+
+        if random.random() <= 0.5:
+            return ind1
+        return ind2
+
+    individuals_1 = random.sample(individuals, len(individuals))
+    individuals_2 = random.sample(individuals, len(individuals))
+
+    chosen = []
+    for i in xrange(0, k, 4):
+        chosen.append(tourn(individuals_1[i],   individuals_1[i+1]))
+        chosen.append(tourn(individuals_1[i+2], individuals_1[i+3]))
+        chosen.append(tourn(individuals_2[i],   individuals_2[i+1]))
+        chosen.append(tourn(individuals_2[i+2], individuals_2[i+3]))
+
+    return chosen
+
+#######################################
+# Generalized Reduced runtime ND sort #
+#######################################
+
+def identity(obj):
+    """Returns directly the argument *obj*.
+    """
+    return obj
+
+def isDominated(wvalues1, wvalues2):
+    """Returns whether or not *wvalues1* dominates *wvalues2*.
+    
+    :param wvalues1: The weighted fitness values that would be dominated.
+    :param wvalues2: The weighted fitness values of the dominant.
+    :returns: :obj:`True` if wvalues2 dominates wvalues1, :obj:`False`
+              otherwise.
+    """
+    not_equal = False
+    for self_wvalue, other_wvalue in zip(wvalues1, wvalues2):
+        if self_wvalue > other_wvalue:
+            return False
+        elif self_wvalue < other_wvalue:
+            not_equal = True
+    return not_equal
+
+def median(seq, key=identity):
+    """Returns the median of *seq* - the numeric value separating the higher 
+    half of a sample from the lower half. If there is an even number of 
+    elements in *seq*, it returns the mean of the two middle values.
+    """
+    sseq = sorted(seq, key=key)
+    length = len(seq)
+    if length % 2 == 1:
+        return key(sseq[(length - 1) // 2])
+    else:
+        return (key(sseq[(length - 1) // 2]) + key(sseq[length // 2])) / 2.0
+
+def sortLogNondominated(individuals, k, first_front_only=False):
+    """Sort *individuals* in pareto non-dominated fronts using the Generalized
+    Reduced Run-Time Complexity Non-Dominated Sorting Algorithm presented by
+    Fortin et al. (2013).
+    
+    :param individuals: A list of individuals to select from.
+    :returns: A list of Pareto fronts (lists), with the first list being the
+              true Pareto front.
+    """
+    if k == 0:
+        return []
+    
+    #Separate individuals according to unique fitnesses
+    unique_fits = defaultdict(list)
+    for i, ind in enumerate(individuals):
+        unique_fits[ind.fitness.wvalues].append(ind)
+    
+    #Launch the sorting algorithm
+    obj = len(individuals[0].fitness.wvalues)-1
+    fitnesses = unique_fits.keys()
+    front = dict.fromkeys(fitnesses, 0)
+
+    # Sort the fitnesses lexicographically.
+    fitnesses.sort(reverse=True)
+    sortNDHelperA(fitnesses, obj, front)    
+    
+    #Extract individuals from front list here
+    nbfronts = max(front.values())+1
+    pareto_fronts = [[] for i in range(nbfronts)]
+    for fit in fitnesses:
+        index = front[fit]
+        pareto_fronts[index].extend(unique_fits[fit])
+
+    # Keep only the fronts required to have k individuals.
+    if not first_front_only:
+        count = 0
+        for i, front in enumerate(pareto_fronts):
+            count += len(front)
+            if count >= k:
+                return pareto_fronts[:i+1]
+        return pareto_fronts
+    else:
+        return pareto_fronts[0]
+
+def sortNDHelperA(fitnesses, obj, front):
+    """Create a non-dominated sorting of S on the first M objectives"""
+    if len(fitnesses) < 2:
+        return
+    elif len(fitnesses) == 2:
+        # Only two individuals, compare them and adjust front number
+        s1, s2 = fitnesses[0], fitnesses[1]
+        if isDominated(s2[:obj+1], s1[:obj+1]):
+            front[s2] = max(front[s2], front[s1] + 1)
+    elif obj == 1:
+        sweepA(fitnesses, front)
+    elif len(frozenset(map(itemgetter(obj), fitnesses))) == 1:
+        #All individuals for objective M are equal: go to objective M-1
+        sortNDHelperA(fitnesses, obj-1, front)
+    else:
+        # More than two individuals, split list and then apply recursion
+        best, worst = splitA(fitnesses, obj)
+        sortNDHelperA(best, obj, front)
+        sortNDHelperB(best, worst, obj-1, front)
+        sortNDHelperA(worst, obj, front)
+
+def splitA(fitnesses, obj):
+    """Partition the set of fitnesses in two according to the median of
+    the objective index *obj*. The values equal to the median are put in
+    the set containing the least elements.
+    """
+    median_ = median(fitnesses, itemgetter(obj))
+    best_a, worst_a = [], []
+    best_b, worst_b = [], []
+
+    for fit in fitnesses:
+        if fit[obj] > median_:
+            best_a.append(fit)
+            best_b.append(fit)
+        elif fit[obj] < median_:
+            worst_a.append(fit)
+            worst_b.append(fit)
+        else:
+            best_a.append(fit)
+            worst_b.append(fit)
+
+    balance_a = abs(len(best_a) - len(worst_a))
+    balance_b = abs(len(best_b) - len(worst_b))
+
+    if balance_a <= balance_b:
+        return best_a, worst_a
+    else:
+        return best_b, worst_b
+
+def sweepA(fitnesses, front):
+    """Update rank number associated to the fitnesses according
+    to the first two objectives using a geometric sweep procedure.
+    """
+    stairs = [-fitnesses[0][1]]
+    fstairs = [fitnesses[0]]  
+    for fit in fitnesses[1:]:
+        idx = bisect.bisect_right(stairs, -fit[1])
+        if 0 < idx <= len(stairs):
+            fstair = max(fstairs[:idx], key=front.__getitem__)
+            front[fit] = max(front[fit], front[fstair]+1)
+        for i, fstair in enumerate(fstairs[idx:], idx):
+            if front[fstair] == front[fit]:
+                del stairs[i]
+                del fstairs[i]
+                break
+        stairs.insert(idx, -fit[1])
+        fstairs.insert(idx, fit)
+
+def sortNDHelperB(best, worst, obj, front):
+    """Assign front numbers to the solutions in H according to the solutions 
+    in L. The solutions in L are assumed to have correct front numbers and the 
+    solutions in H are not compared with each other, as this is supposed to 
+    happen after sortNDHelperB is called."""
+    key = itemgetter(obj)    
+    if len(worst) == 0 or len(best) == 0:
+        #One of the lists is empty: nothing to do
+        return
+    elif len(best) == 1 or len(worst) == 1:
+        #One of the lists has one individual: compare directly
+        for hi in worst:
+            for li in best:
+                if isDominated(hi[:obj+1], li[:obj+1]) or hi[:obj+1] == li[:obj+1]:
+                    front[hi] = max(front[hi], front[li] + 1)
+    elif obj == 1:
+        sweepB(best, worst, front)
+    elif key(min(best, key=key)) >= key(max(worst, key=key)):
+        #All individuals from L dominate H for objective M:
+        #Also supports the case where every individuals in L and H 
+        #has the same value for the current objective
+        #Skip to objective M-1
+        sortNDHelperB(best, worst, obj-1, front)
+    elif key(max(best, key=key)) >= key(min(worst, key=key)):
+        best1, best2, worst1, worst2 = splitB(best, worst, obj)
+        sortNDHelperB(best1, worst1, obj, front)
+        sortNDHelperB(best1, worst2, obj-1, front)
+        sortNDHelperB(best2, worst2, obj, front)
+
+def splitB(best, worst, obj):
+    """Split both best individual and worst sets of fitnesses according
+    to the median of objective *obj* computed on the set containing the
+    most elements. The values equal to the median are attributed so as 
+    to balance the four resulting sets as much as possible.
+    """
+    median_ = median(best if len(best) > len(worst) else worst, itemgetter(obj))
+    best1_a, best2_a, best1_b, best2_b = [], [], [], []
+    for fit in best:
+        if fit[obj] > median_:
+            best1_a.append(fit)
+            best1_b.append(fit)
+        elif fit[obj] < median_:
+            best2_a.append(fit)
+            best2_b.append(fit)
+        else:
+            best1_a.append(fit)
+            best2_b.append(fit)
+    
+    worst1_a, worst2_a, worst1_b, worst2_b = [], [], [], []        
+    for fit in worst:
+        if fit[obj] > median_:
+            worst1_a.append(fit)
+            worst1_b.append(fit)
+        elif fit[obj] < median_:
+            worst2_a.append(fit)
+            worst2_b.append(fit)
+        else:
+            worst1_a.append(fit)
+            worst2_b.append(fit)
+    
+    balance_a = abs(len(best1_a) - len(best2_a) + len(worst1_a) - len(worst2_a))
+    balance_b = abs(len(best1_b) - len(best2_b) + len(worst1_b) - len(worst2_b))
+
+    if balance_a <= balance_b:
+        return best1_a, best2_a, worst1_a, worst2_a
+    else:
+        return best1_b, best2_b, worst1_b, worst2_b
+
+def sweepB(best, worst, front):
+    """Adjust the rank number of the worst fitnesses according to
+    the best fitnesses on the first two objectives using a sweep
+    procedure.
+    """
+    stairs, fstairs = [], []
+    iter_best = iter(best)
+    next_best = next(iter_best, False)
+    for h in worst:
+        while next_best and h[:2] <= next_best[:2]:
+            insert = True
+            for i, fstair in enumerate(fstairs):
+                if front[fstair] == front[next_best]:
+                    if fstair[1] > next_best[1]:
+                        insert = False
+                    else:
+                        del stairs[i], fstairs[i]
+                    break
+            if insert:
+                idx = bisect.bisect_right(stairs, -next_best[1])
+                stairs.insert(idx, -next_best[1])
+                fstairs.insert(idx, next_best)
+            next_best = next(iter_best, False)
+
+        idx = bisect.bisect_right(stairs, -h[1])
+        if 0 < idx <= len(stairs):
+            fstair = max(fstairs[:idx], key=front.__getitem__)
+            front[h] = max(front[h], front[fstair]+1)
+
+######################################
+# Strength Pareto         (SPEA-II)  #
+######################################
+
+def selSPEA2(individuals, k):
+    """Apply SPEA-II selection operator on the *individuals*. Usually, the
+    size of *individuals* will be larger than *n* because any individual
+    present in *individuals* will appear in the returned list at most once.
+    Having the size of *individuals* equals to *n* will have no effect other
+    than sorting the population according to a strength Pareto scheme. The
+    list returned contains references to the input *individuals*. For more
+    details on the SPEA-II operator see [Zitzler2001]_.
+    
+    :param individuals: A list of individuals to select from.
+    :param k: The number of individuals to select.
+    :returns: A list of selected individuals.
+    
+    .. [Zitzler2001] Zitzler, Laumanns and Thiele, "SPEA 2: Improving the
+       strength Pareto evolutionary algorithm", 2001.
+    """
+    N = len(individuals)
+    L = len(individuals[0].fitness.values)
+    K = math.sqrt(N)
+    strength_fits = [0] * N
+    fits = [0] * N
+    dominating_inds = [list() for i in xrange(N)]
+    
+    for i, ind_i in enumerate(individuals):
+        for j, ind_j in enumerate(individuals[i+1:], i+1):
+            if ind_i.fitness.dominates(ind_j.fitness):
+                strength_fits[i] += 1
+                dominating_inds[j].append(i)
+            elif ind_j.fitness.dominates(ind_i.fitness):
+                strength_fits[j] += 1
+                dominating_inds[i].append(j)
+    
+    for i in xrange(N):
+        for j in dominating_inds[i]:
+            fits[i] += strength_fits[j]
+    
+    # Choose all non-dominated individuals
+    chosen_indices = [i for i in xrange(N) if fits[i] < 1]
+    
+    if len(chosen_indices) < k:     # The archive is too small
+        for i in xrange(N):
+            distances = [0.0] * N
+            for j in xrange(i + 1, N):
+                dist = 0.0
+                for l in xrange(L):
+                    val = individuals[i].fitness.values[l] - \
+                          individuals[j].fitness.values[l]
+                    dist += val * val
+                distances[j] = dist
+            kth_dist = _randomizedSelect(distances, 0, N - 1, K)
+            density = 1.0 / (kth_dist + 2.0)
+            fits[i] += density
+            
+        next_indices = [(fits[i], i) for i in xrange(N)
+                        if not i in chosen_indices]
+        next_indices.sort()
+        #print next_indices
+        chosen_indices += [i for _, i in next_indices[:k - len(chosen_indices)]]
+                
+    elif len(chosen_indices) > k:   # The archive is too large
+        N = len(chosen_indices)
+        distances = [[0.0] * N for i in xrange(N)]
+        sorted_indices = [[0] * N for i in xrange(N)]
+        for i in xrange(N):
+            for j in xrange(i + 1, N):
+                dist = 0.0
+                for l in xrange(L):
+                    val = individuals[chosen_indices[i]].fitness.values[l] - \
+                          individuals[chosen_indices[j]].fitness.values[l]
+                    dist += val * val
+                distances[i][j] = dist
+                distances[j][i] = dist
+            distances[i][i] = -1
+        
+        # Insert sort is faster than quick sort for short arrays
+        for i in xrange(N):
+            for j in xrange(1, N):
+                l = j
+                while l > 0 and distances[i][j] < distances[i][sorted_indices[i][l - 1]]:
+                    sorted_indices[i][l] = sorted_indices[i][l - 1]
+                    l -= 1
+                sorted_indices[i][l] = j
+        
+        size = N
+        to_remove = []
+        while size > k:
+            # Search for minimal distance
+            min_pos = 0
+            for i in xrange(1, N):
+                for j in xrange(1, size):
+                    dist_i_sorted_j = distances[i][sorted_indices[i][j]]
+                    dist_min_sorted_j = distances[min_pos][sorted_indices[min_pos][j]]
+                    
+                    if dist_i_sorted_j < dist_min_sorted_j:
+                        min_pos = i
+                        break
+                    elif dist_i_sorted_j > dist_min_sorted_j:
+                        break
+            
+            # Remove minimal distance from sorted_indices
+            for i in xrange(N):
+                distances[i][min_pos] = float("inf")
+                distances[min_pos][i] = float("inf")
+                
+                for j in xrange(1, size - 1):
+                    if sorted_indices[i][j] == min_pos:
+                        sorted_indices[i][j] = sorted_indices[i][j + 1]
+                        sorted_indices[i][j + 1] = min_pos
+            
+            # Remove corresponding individual from chosen_indices
+            to_remove.append(min_pos)
+            size -= 1
+        
+        for index in reversed(sorted(to_remove)):
+            del chosen_indices[index]
+    
+    return [individuals[i] for i in chosen_indices]
+    
+def _randomizedSelect(array, begin, end, i):
+    """Allows to select the ith smallest element from array without sorting it.
+    Runtime is expected to be O(n).
+    """
+    if begin == end:
+        return array[begin]
+    q = _randomizedPartition(array, begin, end)
+    k = q - begin + 1
+    if i < k:
+        return _randomizedSelect(array, begin, q, i)
+    else:
+        return _randomizedSelect(array, q + 1, end, i - k)
+
+def _randomizedPartition(array, begin, end):
+    i = random.randint(begin, end)
+    array[begin], array[i] = array[i], array[begin]
+    return _partition(array, begin, end)
+    
+def _partition(array, begin, end):
+    x = array[begin]
+    i = begin - 1
+    j = end + 1
+    while True:
+        j -= 1
+        while array[j] > x:
+            j -= 1
+        i += 1
+        while array[i] < x:
+            i += 1
+        if i < j:
+            array[i], array[j] = array[j], array[i]
+        else:
+            return j
+
+
+__all__ = ['selNSGA2', 'selSPEA2', 'sortNondominated', 'sortLogNondominated',
+           'selTournamentDCD']

python/isaac/autotuning/external/deap/tools/init.py

@@ -0,0 +1,86 @@
+from __future__ import division
+
+def initRepeat(container, func, n):
+    """Call the function *container* with a generator function corresponding
+    to the calling *n* times the function *func*.
+    
+    :param container: The type to put in the data from func.
+    :param func: The function that will be called n times to fill the
+                 container.
+    :param n: The number of times to repeat func.
+    :returns: An instance of the container filled with data from func.
+    
+    This helper function can can be used in conjunction with a Toolbox 
+    to register a generator of filled containers, as individuals or 
+    population.
+    
+        >>> initRepeat(list, random.random, 2) # doctest: +ELLIPSIS, 
+        ...                                    # doctest: +NORMALIZE_WHITESPACE
+        [0.4761..., 0.6302...]
+
+    See the :ref:`list-of-floats` and :ref:`population` tutorials for more examples.
+    """
+    return container(func() for _ in xrange(n))
+
+def initIterate(container, generator):
+    """Call the function *container* with an iterable as 
+    its only argument. The iterable must be returned by 
+    the method or the object *generator*.
+    
+    :param container: The type to put in the data from func.
+    :param generator: A function returning an iterable (list, tuple, ...),
+                      the content of this iterable will fill the container.
+    :returns: An instance of the container filled with data from the
+              generator.
+    
+    This helper function can can be used in conjunction with a Toolbox 
+    to register a generator of filled containers, as individuals or 
+    population.
+
+        >>> from random import sample
+        >>> from functools import partial
+        >>> gen_idx = partial(sample, range(10), 10)
+        >>> initIterate(list, gen_idx)
+        [4, 5, 3, 6, 0, 9, 2, 7, 1, 8]
+
+    See the :ref:`permutation` and :ref:`arithmetic-expr` tutorials for
+    more examples.
+    """
+    return container(generator())
+
+def initCycle(container, seq_func, n=1):
+    """Call the function *container* with a generator function corresponding
+    to the calling *n* times the functions present in *seq_func*.
+    
+    :param container: The type to put in the data from func.
+    :param seq_func: A list of function objects to be called in order to
+                     fill the container.
+    :param n: Number of times to iterate through the list of functions.
+    :returns: An instance of the container filled with data from the
+              returned by the functions.
+    
+    This helper function can can be used in conjunction with a Toolbox 
+    to register a generator of filled containers, as individuals or 
+    population.
+    
+        >>> func_seq = [lambda:1 , lambda:'a', lambda:3]
+        >>> initCycle(list, func_seq, n=2)
+        [1, 'a', 3, 1, 'a', 3]
+
+    See the :ref:`funky` tutorial for an example.
+    """
+    return container(func() for _ in xrange(n) for func in seq_func)
+
+__all__ = ['initRepeat', 'initIterate', 'initCycle']
+
+
+if __name__ == "__main__":
+    import doctest
+    import random
+    random.seed(64)
+    doctest.run_docstring_examples(initRepeat, globals())
+    
+    random.seed(64)
+    doctest.run_docstring_examples(initIterate, globals())
+    doctest.run_docstring_examples(initCycle, globals())
+    

python/isaac/autotuning/external/deap/tools/migration.py

@@ -0,0 +1,53 @@
+from __future__ import division
+
+
+def migRing(populations, k, selection, replacement=None, migarray=None):
+    """Perform a ring migration between the *populations*. The migration first
+    select *k* emigrants from each population using the specified *selection*
+    operator and then replace *k* individuals from the associated population
+    in the *migarray* by the emigrants. If no *replacement* operator is
+    specified, the immigrants will replace the emigrants of the population,
+    otherwise, the immigrants will replace the individuals selected by the
+    *replacement* operator. The migration array, if provided, shall contain
+    each population's index once and only once. If no migration array is
+    provided, it defaults to a serial ring migration (1 -- 2 -- ... -- n --
+    1). Selection and replacement function are called using the signature
+    ``selection(populations[i], k)`` and ``replacement(populations[i], k)``.
+    It is important to note that the replacement strategy must select *k*
+    **different** individuals. For example, using a traditional tournament for
+    replacement strategy will thus give undesirable effects, two individuals
+    will most likely try to enter the same slot.
+    
+    :param populations: A list of (sub-)populations on which to operate
+                        migration.
+    :param k: The number of individuals to migrate.
+    :param selection: The function to use for selection.
+    :param replacement: The function to use to select which individuals will
+                        be replaced. If :obj:`None` (default) the individuals
+                        that leave the population are directly replaced.
+    :param migarray: A list of indices indicating where the individuals from 
+                     a particular position in the list goes. This defaults
+                     to a ring migration.
+    """
+    nbr_demes = len(populations)
+    if migarray is None:
+        migarray = range(1, nbr_demes) + [0]
+    
+    immigrants = [[] for i in xrange(nbr_demes)]
+    emigrants = [[] for i in xrange(nbr_demes)]
+
+    for from_deme in xrange(nbr_demes):
+        emigrants[from_deme].extend(selection(populations[from_deme], k))
+        if replacement is None:
+            # If no replacement strategy is selected, replace those who migrate
+            immigrants[from_deme] = emigrants[from_deme]
+        else:
+            # Else select those who will be replaced
+            immigrants[from_deme].extend(replacement(populations[from_deme], k))
+
+    for from_deme, to_deme in enumerate(migarray):
+        for i, immigrant in enumerate(immigrants[to_deme]):
+            indx = populations[to_deme].index(immigrant)
+            populations[to_deme][indx] = emigrants[from_deme][i]
+
+__all__ = ['migRing']

python/isaac/autotuning/external/deap/tools/mutation.py

@@ -0,0 +1,209 @@
+from __future__ import division
+import math
+import random
+
+from itertools import repeat
+from collections import Sequence
+
+######################################
+# GA Mutations                       #
+######################################
+
+def mutGaussian(individual, mu, sigma, indpb):
+    """This function applies a gaussian mutation of mean *mu* and standard
+    deviation *sigma* on the input individual. This mutation expects a
+    :term:`sequence` individual composed of real valued attributes.
+    The *indpb* argument is the probability of each attribute to be mutated.
+    
+    :param individual: Individual to be mutated.
+    :param mu: Mean or :term:`python:sequence` of means for the
+               gaussian addition mutation.
+    :param sigma: Standard deviation or :term:`python:sequence` of 
+                  standard deviations for the gaussian addition mutation.
+    :param indpb: Independent probability for each attribute to be mutated.
+    :returns: A tuple of one individual.
+    
+    This function uses the :func:`~random.random` and :func:`~random.gauss`
+    functions from the python base :mod:`random` module.
+    """
+    size = len(individual)
+    if not isinstance(mu, Sequence):
+        mu = repeat(mu, size)
+    elif len(mu) < size:
+        raise IndexError("mu must be at least the size of individual: %d < %d" % (len(mu), size))
+    if not isinstance(sigma, Sequence):
+        sigma = repeat(sigma, size)
+    elif len(sigma) < size:
+        raise IndexError("sigma must be at least the size of individual: %d < %d" % (len(sigma), size))
+    
+    for i, m, s in zip(xrange(size), mu, sigma):
+        if random.random() < indpb:
+            individual[i] += random.gauss(m, s)
+    
+    return individual,
+
+def mutPolynomialBounded(individual, eta, low, up, indpb):
+    """Polynomial mutation as implemented in original NSGA-II algorithm in
+    C by Deb.
+    
+    :param individual: :term:`Sequence <sequence>` individual to be mutated.
+    :param eta: Crowding degree of the mutation. A high eta will produce
+                a mutant resembling its parent, while a small eta will
+                produce a solution much more different.
+    :param low: A value or a :term:`python:sequence` of values that
+                is the lower bound of the search space.
+    :param up: A value or a :term:`python:sequence` of values that
+               is the upper bound of the search space.
+    :returns: A tuple of one individual.
+    """
+    size = len(individual)
+    if not isinstance(low, Sequence):
+        low = repeat(low, size)
+    elif len(low) < size:
+        raise IndexError("low must be at least the size of individual: %d < %d" % (len(low), size))
+    if not isinstance(up, Sequence):
+        up = repeat(up, size)
+    elif len(up) < size:
+        raise IndexError("up must be at least the size of individual: %d < %d" % (len(up), size))
+    
+    for i, xl, xu in zip(xrange(size), low, up):
+        if random.random() <= indpb:
+            x = individual[i]
+            delta_1 = (x - xl) / (xu - xl)
+            delta_2 = (xu - x) / (xu - xl)
+            rand = random.random()
+            mut_pow = 1.0 / (eta + 1.)
+
+            if rand < 0.5:
+                xy = 1.0 - delta_1
+                val = 2.0 * rand + (1.0 - 2.0 * rand) * xy**(eta + 1)
+                delta_q = val**mut_pow - 1.0
+            else:
+                xy = 1.0 - delta_2
+                val = 2.0 * (1.0 - rand) + 2.0 * (rand - 0.5) * xy**(eta + 1)
+                delta_q = 1.0 - val**mut_pow
+
+            x = x + delta_q * (xu - xl)
+            x = min(max(x, xl), xu)
+            individual[i] = x
+    return individual,
+
+def mutShuffleIndexes(individual, indpb):
+    """Shuffle the attributes of the input individual and return the mutant.
+    The *individual* is expected to be a :term:`sequence`. The *indpb* argument is the
+    probability of each attribute to be moved. Usually this mutation is applied on 
+    vector of indices.
+    
+    :param individual: Individual to be mutated.
+    :param indpb: Independent probability for each attribute to be exchanged to
+                  another position.
+    :returns: A tuple of one individual.
+    
+    This function uses the :func:`~random.random` and :func:`~random.randint`
+    functions from the python base :mod:`random` module.
+    """
+    size = len(individual)
+    for i in xrange(size):
+        if random.random() < indpb:
+            swap_indx = random.randint(0, size - 2)
+            if swap_indx >= i:
+                swap_indx += 1
+            individual[i], individual[swap_indx] = \
+                individual[swap_indx], individual[i]
+    
+    return individual,
+
+def mutFlipBit(individual, indpb):
+    """Flip the value of the attributes of the input individual and return the
+    mutant. The *individual* is expected to be a :term:`sequence` and the values of the
+    attributes shall stay valid after the ``not`` operator is called on them.
+    The *indpb* argument is the probability of each attribute to be
+    flipped. This mutation is usually applied on boolean individuals.
+    
+    :param individual: Individual to be mutated.
+    :param indpb: Independent probability for each attribute to be flipped.
+    :returns: A tuple of one individual.
+    
+    This function uses the :func:`~random.random` function from the python base
+    :mod:`random` module.
+    """
+    for i in xrange(len(individual)):
+        if random.random() < indpb:
+            individual[i] = type(individual[i])(not individual[i])
+    
+    return individual,
+
+def mutUniformInt(individual, low, up, indpb):
+    """Mutate an individual by replacing attributes, with probability *indpb*,
+    by a integer uniformly drawn between *low* and *up* inclusively.
+    
+    :param individual: :term:`Sequence <sequence>` individual to be mutated.
+    :param low: The lower bound or a :term:`python:sequence` of
+                of lower bounds of the range from wich to draw the new
+                integer.
+    :param up: The upper bound or a :term:`python:sequence` of
+               of upper bounds of the range from wich to draw the new
+               integer.
+    :param indpb: Independent probability for each attribute to be mutated.
+    :returns: A tuple of one individual.
+    """
+    size = len(individual)
+    if not isinstance(low, Sequence):
+        low = repeat(low, size)
+    elif len(low) < size:
+        raise IndexError("low must be at least the size of individual: %d < %d" % (len(low), size))
+    if not isinstance(up, Sequence):
+        up = repeat(up, size)
+    elif len(up) < size:
+        raise IndexError("up must be at least the size of individual: %d < %d" % (len(up), size))
+    
+    for i, xl, xu in zip(xrange(size), low, up):
+        if random.random() < indpb:
+            individual[i] = random.randint(xl, xu)
+    
+    return individual,
+
+
+######################################
+# ES Mutations                       #
+######################################
+
+def mutESLogNormal(individual, c, indpb):
+    """Mutate an evolution strategy according to its :attr:`strategy`
+    attribute as described in [Beyer2002]_. First the strategy is mutated
+    according to an extended log normal rule, :math:`\\boldsymbol{\sigma}_t =
+    \\exp(\\tau_0 \mathcal{N}_0(0, 1)) \\left[ \\sigma_{t-1, 1}\\exp(\\tau
+    \mathcal{N}_1(0, 1)), \ldots, \\sigma_{t-1, n} \\exp(\\tau
+    \mathcal{N}_n(0, 1))\\right]`, with :math:`\\tau_0 =
+    \\frac{c}{\\sqrt{2n}}` and :math:`\\tau = \\frac{c}{\\sqrt{2\\sqrt{n}}}`,
+    the the individual is mutated by a normal distribution of mean 0 and
+    standard deviation of :math:`\\boldsymbol{\sigma}_{t}` (its current
+    strategy) then . A recommended choice is ``c=1`` when using a :math:`(10,
+    100)` evolution strategy [Beyer2002]_ [Schwefel1995]_.
+    
+    :param individual: :term:`Sequence <sequence>` individual to be mutated.
+    :param c: The learning parameter.
+    :param indpb: Independent probability for each attribute to be mutated.
+    :returns: A tuple of one individual.
+    
+    .. [Beyer2002] Beyer and Schwefel, 2002, Evolution strategies - A
+       Comprehensive Introduction
+       
+    .. [Schwefel1995] Schwefel, 1995, Evolution and Optimum Seeking.
+       Wiley, New York, NY
+    """
+    size = len(individual)
+    t = c / math.sqrt(2. * math.sqrt(size))
+    t0 = c / math.sqrt(2. * size)
+    n = random.gauss(0, 1)
+    t0_n = t0 * n
+    
+    for indx in xrange(size):
+        if random.random() < indpb:
+            individual.strategy[indx] *= math.exp(t0_n + t * random.gauss(0, 1))
+            individual[indx] += individual.strategy[indx] * random.gauss(0, 1)
+    
+    return individual,
+
+__all__ = ['mutGaussian', 'mutPolynomialBounded', 'mutShuffleIndexes', 
+           'mutFlipBit', 'mutUniformInt', 'mutESLogNormal']

python/isaac/autotuning/external/deap/tools/selection.py

@@ -0,0 +1,177 @@
+from __future__ import division
+import random
+
+from functools import partial
+from operator import attrgetter
+
+######################################
+# Selections                         #
+######################################
+
+def selRandom(individuals, k):
+    """Select *k* individuals at random from the input *individuals* with
+    replacement. The list returned contains references to the input
+    *individuals*.
+    
+    :param individuals: A list of individuals to select from.
+    :param k: The number of individuals to select.
+    :returns: A list of selected individuals.
+    
+    This function uses the :func:`~random.choice` function from the
+    python base :mod:`random` module.
+    """
+    return [random.choice(individuals) for i in xrange(k)]
+
+
+def selBest(individuals, k):
+    """Select the *k* best individuals among the input *individuals*. The
+    list returned contains references to the input *individuals*.
+    
+    :param individuals: A list of individuals to select from.
+    :param k: The number of individuals to select.
+    :returns: A list containing the k best individuals.
+    """
+    return sorted(individuals, key=attrgetter("fitness"), reverse=True)[:k]
+
+
+def selWorst(individuals, k):
+    """Select the *k* worst individuals among the input *individuals*. The
+    list returned contains references to the input *individuals*.
+    
+    :param individuals: A list of individuals to select from.
+    :param k: The number of individuals to select.
+    :returns: A list containing the k worst individuals.
+    """
+    return sorted(individuals, key=attrgetter("fitness"))[:k]
+
+
+def selTournament(individuals, k, tournsize):
+    """Select *k* individuals from the input *individuals* using *k*
+    tournaments of *tournsize* individuals. The list returned contains
+    references to the input *individuals*.
+    
+    :param individuals: A list of individuals to select from.
+    :param k: The number of individuals to select.
+    :param tournsize: The number of individuals participating in each tournament.
+    :returns: A list of selected individuals.
+    
+    This function uses the :func:`~random.choice` function from the python base
+    :mod:`random` module.
+    """
+    chosen = []
+    for i in xrange(k):
+        aspirants = selRandom(individuals, tournsize)
+        chosen.append(max(aspirants, key=attrgetter("fitness")))
+    return chosen
+
+def selRoulette(individuals, k):
+    """Select *k* individuals from the input *individuals* using *k*
+    spins of a roulette. The selection is made by looking only at the first
+    objective of each individual. The list returned contains references to
+    the input *individuals*.
+    
+    :param individuals: A list of individuals to select from.
+    :param k: The number of individuals to select.
+    :returns: A list of selected individuals.
+    
+    This function uses the :func:`~random.random` function from the python base
+    :mod:`random` module.
+    
+    .. warning::
+       The roulette selection by definition cannot be used for minimization 
+       or when the fitness can be smaller or equal to 0.
+    """
+    s_inds = sorted(individuals, key=attrgetter("fitness"), reverse=True)
+    sum_fits = sum(ind.fitness.values[0] for ind in individuals)
+    
+    chosen = []
+    for i in xrange(k):
+        u = random.random() * sum_fits
+        sum_ = 0
+        for ind in s_inds:
+            sum_ += ind.fitness.values[0]
+            if sum_ > u:
+                chosen.append(ind)
+                break
+    
+    return chosen
+
+
+def selDoubleTournament(individuals, k, fitness_size, parsimony_size, fitness_first):
+    """Tournament selection which use the size of the individuals in order
+    to discriminate good solutions. This kind of tournament is obviously
+    useless with fixed-length representation, but has been shown to
+    significantly reduce excessive growth of individuals, especially in GP,
+    where it can be used as a bloat control technique (see 
+    [Luke2002fighting]_). This selection operator implements the double 
+    tournament technique presented in this paper.
+    
+    The core principle is to use a normal tournament selection, but using a
+    special sample function to select aspirants, which is another tournament
+    based on the size of the individuals. To ensure that the selection
+    pressure is not too high, the size of the size tournament (the number
+    of candidates evaluated) can be a real number between 1 and 2. In this
+    case, the smaller individual among two will be selected with a probability
+    *size_tourn_size*/2. For instance, if *size_tourn_size* is set to 1.4,
+    then the smaller individual will have a 0.7 probability to be selected.
+    
+    .. note::
+        In GP, it has been shown that this operator produces better results
+        when it is combined with some kind of a depth limit.
+    
+    :param individuals: A list of individuals to select from.
+    :param k: The number of individuals to select.
+    :param fitness_size: The number of individuals participating in each \
+    fitness tournament
+    :param parsimony_size: The number of individuals participating in each \
+    size tournament. This value has to be a real number\
+    in the range [1,2], see above for details.
+    :param fitness_first: Set this to True if the first tournament done should \
+    be the fitness one (i.e. the fitness tournament producing aspirants for \
+    the size tournament). Setting it to False will behaves as the opposite \
+    (size tournament feeding fitness tournaments with candidates). It has been \
+    shown that this parameter does not have a significant effect in most cases\
+    (see [Luke2002fighting]_).
+    :returns: A list of selected individuals.
+    
+    .. [Luke2002fighting] Luke and Panait, 2002, Fighting bloat with 
+        nonparametric parsimony pressure
+    """
+    assert (1 <= parsimony_size <= 2), "Parsimony tournament size has to be in the range [1, 2]."
+
+    def _sizeTournament(individuals, k, select):
+        chosen = []
+        for i in xrange(k):
+            # Select two individuals from the population
+            # The first individual has to be the shortest
+            prob = parsimony_size / 2.
+            ind1, ind2 = select(individuals, k=2)
+
+            if len(ind1) > len(ind2):
+                ind1, ind2 = ind2, ind1
+            elif len(ind1) == len(ind2):
+                # random selection in case of a tie
+                prob = 0.5
+
+            # Since size1 <= size2 then ind1 is selected
+            # with a probability prob
+            chosen.append(ind1 if random.random() < prob else ind2)
+
+        return chosen
+    
+    def _fitTournament(individuals, k, select):
+        chosen = []
+        for i in xrange(k):
+            aspirants = select(individuals, k=fitness_size)
+            chosen.append(max(aspirants, key=attrgetter("fitness")))
+        return chosen
+    
+    if fitness_first:
+        tfit = partial(_fitTournament, select=selRandom)
+        return _sizeTournament(individuals, k, tfit)
+    else:
+        tsize = partial(_sizeTournament, select=selRandom)
+        return _fitTournament(individuals, k, tsize)
+
+__all__ = ['selRandom', 'selBest', 'selWorst', 'selRoulette',
+           'selTournament', 'selDoubleTournament']

python/isaac/autotuning/external/deap/tools/support.py

@@ -0,0 +1,647 @@
+from __future__ import division
+
+try:
+    import cPickle as pickle
+except ImportError:
+    import pickle
+
+from bisect import bisect_right
+from collections import defaultdict
+from copy import deepcopy
+from functools import partial
+from itertools import chain
+from operator import eq
+
+
+def identity(obj):
+    """Returns directly the argument *obj*.
+    """
+    return obj
+
+class History(object):
+    """The :class:`History` class helps to build a genealogy of all the
+    individuals produced in the evolution. It contains two attributes,
+    the :attr:`genealogy_tree` that is a dictionary of lists indexed by
+    individual, the list contain the indices of the parents. The second
+    attribute :attr:`genealogy_history` contains every individual indexed
+    by their individual number as in the genealogy tree.
+    
+    The produced genealogy tree is compatible with `NetworkX
+    <http://networkx.lanl.gov/index.html>`_, here is how to plot the genealogy
+    tree ::
+    
+        history = History()
+        
+        # Decorate the variation operators
+        toolbox.decorate("mate", history.decorator)
+        toolbox.decorate("mutate", history.decorator)
+        
+        # Create the population and populate the history
+        population = toolbox.population(n=POPSIZE)
+        history.update(population)
+        
+        # Do the evolution, the decorators will take care of updating the
+        # history
+        # [...]
+        
+        import matplotlib.pyplot as plt
+        import networkx
+        
+        graph = networkx.DiGraph(history.genealogy_tree)
+        graph = graph.reverse()     # Make the grah top-down
+        colors = [toolbox.evaluate(history.genealogy_history[i])[0] for i in graph]
+        networkx.draw(graph, node_color=colors)
+        plt.show()
+    
+    Using NetworkX in combination with `pygraphviz
+    <http://networkx.lanl.gov/pygraphviz/>`_ (dot layout) this amazing
+    genealogy tree can be obtained from the OneMax example with a population
+    size of 20 and 5 generations, where the color of the nodes indicate there
+    fitness, blue is low and red is high.
+    
+    .. image:: /_images/genealogy.png
+       :width: 67%
+     
+    .. note::
+       The genealogy tree might get very big if your population and/or the 
+       number of generation is large.
+        
+    """
+    def __init__(self):
+        self.genealogy_index = 0
+        self.genealogy_history = dict()
+        self.genealogy_tree = dict()
+        
+    def update(self, individuals):
+        """Update the history with the new *individuals*. The index present in
+        their :attr:`history_index` attribute will be used to locate their
+        parents, it is then modified to a unique one to keep track of those
+        new individuals. This method should be called on the individuals after
+        each variation.
+        
+        :param individuals: The list of modified individuals that shall be
+                            inserted in the history.
+        
+        If the *individuals* do not have a :attr:`history_index` attribute, 
+        the attribute is added and this individual is considered as having no
+        parent. This method should be called with the initial population to
+        initialize the history.
+        
+        Modifying the internal :attr:`genealogy_index` of the history or the
+        :attr:`history_index` of an individual may lead to unpredictable
+        results and corruption of the history.
+        """
+        try:
+            parent_indices = tuple(ind.history_index for ind in individuals)
+        except AttributeError:
+            parent_indices = tuple()
+        
+        for ind in individuals:
+            self.genealogy_index += 1
+            ind.history_index = self.genealogy_index
+            self.genealogy_history[self.genealogy_index] = deepcopy(ind)
+            self.genealogy_tree[self.genealogy_index] = parent_indices
+    
+    @property
+    def decorator(self):
+        """Property that returns an appropriate decorator to enhance the
+        operators of the toolbox. The returned decorator assumes that the
+        individuals are returned by the operator. First the decorator calls
+        the underlying operation and then calls the :func:`update` function
+        with what has been returned by the operator. Finally, it returns the
+        individuals with their history parameters modified according to the
+        update function.
+        """
+        def decFunc(func):
+            def wrapFunc(*args, **kargs):
+                individuals = func(*args, **kargs)
+                self.update(individuals)
+                return individuals
+            return wrapFunc
+        return decFunc
+
+    def getGenealogy(self, individual, max_depth=float("inf")):
+        """Provide the genealogy tree of an *individual*. The individual must
+        have an attribute :attr:`history_index` as defined by
+        :func:`~deap.tools.History.update` in order to retrieve its associated
+        genealogy tree. The returned graph contains the parents up to
+        *max_depth* variations before this individual. If not provided
+        the maximum depth is up to the begining of the evolution.
+
+        :param individual: The individual at the root of the genealogy tree.
+        :param max_depth: The approximate maximum distance between the root
+                          (individual) and the leaves (parents), optional.
+        :returns: A dictionary where each key is an individual index and the
+                  values are a tuple corresponding to the index of the parents.
+        """
+        gtree = {}
+        visited = set()     # Adds memory to the breadth first search
+        def genealogy(index, depth):
+            if index not in self.genealogy_tree:
+                return             
+            depth += 1
+            if depth > max_depth:
+                return
+            parent_indices = self.genealogy_tree[index]
+            gtree[index] = parent_indices
+            for ind in parent_indices:
+                if ind not in visited:
+                    genealogy(ind, depth)
+                visited.add(ind)
+        genealogy(individual.history_index, 0)
+        return gtree
+
+class Statistics(object):
+    """Object that compiles statistics on a list of arbitrary objects. 
+    When created the statistics object receives a *key* argument that 
+    is used to get the values on which the function will be computed. 
+    If not provided the *key* argument defaults to the identity function.
+
+    The value returned by the key may be a multi-dimensional object, i.e.:
+    a tuple or a list, as long as the statistical function registered
+    support it. So for example, statistics can be computed directly on
+    multi-objective fitnesses when using numpy statistical function.
+
+    :param key: A function to access the values on which to compute the
+                statistics, optional.
+        
+    ::
+    
+        >>> s = Statistics()
+        >>> s.register("mean", numpy.mean)
+        >>> s.register("max", max)
+        >>> s.compile([1, 2, 3, 4])
+        {'max': 4, 'mean': 2.5}
+        >>> s.compile([5, 6, 7, 8])
+        {'max': 8, 'mean': 6.5}
+    """
+    def __init__(self, key=identity):
+        self.key = key
+        self.functions = dict()
+        self.fields = []
+
+    def register(self, name, function, *args, **kargs):
+        """Register a *function* that will be applied on the sequence each
+        time :meth:`record` is called.
+
+        :param name: The name of the statistics function as it would appear
+                     in the dictionnary of the statistics object.
+        :param function: A function that will compute the desired statistics
+                         on the data as preprocessed by the key.
+        :param argument: One or more argument (and keyword argument) to pass
+                         automatically to the registered function when called,
+                         optional.
+        """        
+        self.functions[name] = partial(function, *args, **kargs)
+        self.fields.append(name)
+
+    def compile(self, data):
+        """Apply to the input sequence *data* each registered function 
+        and return the results as a dictionnary.
+        
+        :param data: Sequence of objects on which the statistics are computed.
+        """
+        values = tuple(self.key(elem) for elem in data)
+        
+        entry = dict()
+        for key, func in self.functions.iteritems():
+            entry[key] = func(values)
+        return entry
+
+class MultiStatistics(dict):
+    """Dictionary of :class:`Statistics` object allowing to compute
+    statistics on multiple keys using a single call to :meth:`compile`. It
+    takes a set of key-value pairs associating a statistics object to a
+    unique name. This name can then be used to retrieve the statistics object.
+
+    The following code computes statistics simultaneously on the length and
+    the first value of the provided objects.
+    ::
+
+        >>> len_stats = Statistics(key=len)
+        >>> itm0_stats = Statistics(key=itemgetter(0))
+        >>> mstats = MultiStatistics(length=len_stats, item=itm0_stats)
+        >>> mstats.register("mean", numpy.mean, axis=0)
+        >>> mstats.register("max", numpy.max, axis=0)
+        >>> mstats.compile([[0.0, 1.0, 1.0, 5.0], [2.0, 5.0]])
+        {'length': {'max': 4, 'mean': 3.0}, 'item': {'max': 2.0, 'mean': 1.0}}
+    """ 
+    def compile(self, data):
+        """Calls :meth:`Statistics.compile` with *data* of each
+        :class:`Statistics` object.
+        
+        :param data: Sequence of objects on which the statistics are computed.
+        """
+        record = {}
+        for name, stats in self.items():
+            record[name] = stats.compile(data)
+        return record
+
+    @property
+    def fields(self):
+        return sorted(self.keys())
+
+    def register(self, name, function, *args, **kargs):
+        """Register a *function* in each :class:`Statistics` object.
+        
+        :param name: The name of the statistics function as it would appear
+                     in the dictionnary of the statistics object.
+        :param function: A function that will compute the desired statistics
+                         on the data as preprocessed by the key.
+        :param argument: One or more argument (and keyword argument) to pass
+                         automatically to the registered function when called,
+                         optional.
+        """
+        for stats in self.values():
+            stats.register(name, function, *args, **kargs)
+
+class Logbook(list):
+    """Evolution records as a chronological list of dictionaries.
+
+    Data can be retrieved via the :meth:`select` method given the appropriate
+    names.
+
+    The :class:`Logbook` class may also contain other logbooks refered to 
+    as chapters. Chapters are used to store information associated to a
+    specific part of the evolution. For example when computing statistics
+    on different components of individuals (namely :class:`MultiStatistics`),
+    chapters can be used to distinguish the average fitness and the average
+    size.
+    """
+    
+    def __init__(self):
+        self.buffindex = 0
+        self.chapters = defaultdict(Logbook)
+        """Dictionary containing the sub-sections of the logbook which are also
+        :class:`Logbook`. Chapters are automatically created when the right hand
+        side of a keyworded argument, provided to the *record* function, is a
+        dictionnary. The keyword determines the chapter's name. For example, the
+        following line adds a new chapter "size" that will contain the fields
+        "max" and "mean". ::
+
+            logbook.record(gen=0, size={'max' : 10.0, 'mean' : 7.5})
+
+        To access a specific chapter, use the name of the chapter as a
+        dictionnary key. For example, to access the size chapter and select
+        the mean use ::
+
+            logbook.chapters["size"].select("mean")
+
+        Compiling a :class:`MultiStatistics` object returns a dictionary
+        containing dictionnaries, therefore when recording such an object in a
+        logbook using the keyword argument unpacking operator (**), chapters
+        will be automatically added to the logbook.
+        ::
+            
+            >>> fit_stats = Statistics(key=attrgetter("fitness.values"))
+            >>> size_stats = Statistics(key=len)
+            >>> mstats = MultiStatistics(fitness=fit_stats, size=size_stats)
+            >>> # [...]
+            >>> record = mstats.compile(population)
+            >>> logbook.record(**record)
+            >>> print logbook
+              fitness          length
+            ------------    ------------
+            max     mean    max     mean
+            2       1       4       3
+
+        """
+
+        self.columns_len = None
+        self.header = None
+        """Order of the columns to print when using the :data:`stream` and
+        :meth:`__str__` methods. The syntax is a single iterable containing
+        string elements. For example, with the previously
+        defined statistics class, one can print the generation and the
+        fitness average, and maximum with
+        ::
+
+            logbook.header = ("gen", "mean", "max")
+        
+        If not set the header is built with all fields, in arbritrary order
+        on insertion of the first data. The header can be removed by setting
+        it to :data:`None`.
+        """
+        
+        self.log_header = True
+        """Tells the log book to output or not the header when streaming the
+        first line or getting its entire string representation. This defaults
+        :data:`True`.
+        """
+
+    def record(self, **infos):
+        """Enter a record of event in the logbook as a list of key-value pairs.
+        The informations are appended chronogically to a list as a dictionnary.
+        When the value part of a pair is a dictionnary, the informations contained
+        in the dictionnary are recorded in a chapter entitled as the name of the
+        key part of the pair. Chapters are also Logbook.
+        """
+        for key, value in infos.items():
+            if isinstance(value, dict):
+                self.chapters[key].record(**value)
+                del infos[key]
+        self.append(infos)
+
+    def select(self, *names):
+        """Return a list of values associated to the *names* provided
+        in argument in each dictionary of the Statistics object list.
+        One list per name is returned in order.
+        ::
+
+            >>> log = Logbook()
+            >>> log.record(gen = 0, mean = 5.4, max = 10.0)
+            >>> log.record(gen = 1, mean = 9.4, max = 15.0)
+            >>> log.select("mean")
+            [5.4, 9.4]
+            >>> log.select("gen", "max")
+            ([0, 1], [10.0, 15.0])
+
+        With a :class:`MultiStatistics` object, the statistics for each
+        measurement can be retrieved using the :data:`chapters` member :
+        ::
+
+            >>> log = Logbook()
+            >>> log.record(**{'gen' : 0, 'fit' : {'mean' : 0.8, 'max' : 1.5}, 
+            ... 'size' : {'mean' : 25.4, 'max' : 67}})
+            >>> log.record(**{'gen' : 1, 'fit' : {'mean' : 0.95, 'max' : 1.7}, 
+            ... 'size' : {'mean' : 28.1, 'max' : 71}})
+            >>> log.chapters['size'].select("mean")
+            [25.4, 28.1]
+            >>> log.chapters['fit'].select("gen", "max")
+            ([0, 1], [1.5, 1.7])
+        """
+        if len(names) == 1:
+            return [entry.get(names[0], None) for entry in self]
+        return tuple([entry.get(name, None) for entry in self] for name in names)
+
+    @property
+    def stream(self):
+        """Retrieve the formatted not streamed yet entries of the database 
+        including the headers.
+        ::
+
+            >>> log = Logbook()
+            >>> log.append({'gen' : 0})
+            >>> print log.stream
+            gen
+              0
+            >>> log.append({'gen' : 1})
+            >>> print log.stream
+              1
+        """
+        startindex, self.buffindex = self.buffindex, len(self)
+        return self.__str__(startindex)
+
+    def __delitem__(self, key):
+        if isinstance(key, slice):
+            for i, in range(*key.indices(len(self))):
+                self.pop(i)
+                for chapter in self.chapters.values():
+                    chapter.pop(i)
+        else:
+            self.pop(key)
+            for chapter in self.chapters.values():
+                chapter.pop(key)
+        
+    def pop(self, index=0):
+        """Retrieve and delete element *index*. The header and stream will be
+        adjusted to follow the modification.
+
+        :param item: The index of the element to remove, optional. It defaults
+                     to the first element.
+        
+        You can also use the following syntax to delete elements.
+        ::
+        
+            del log[0]
+            del log[1::5]
+        """
+        if index < self.buffindex:
+            self.buffindex -= 1
+        return super(self.__class__, self).pop(index)
+
+    def __txt__(self, startindex):
+        columns = self.header
+        if not columns:
+            columns = sorted(self[0].keys()) + sorted(self.chapters.keys())
+        if not self.columns_len or len(self.columns_len) != len(columns):
+            self.columns_len = map(len, columns)
+
+        chapters_txt = {}
+        offsets = defaultdict(int)
+        for name, chapter in self.chapters.items():
+            chapters_txt[name] = chapter.__txt__(startindex)
+            if startindex == 0:
+                offsets[name] = len(chapters_txt[name]) - len(self)
+
+        str_matrix = []
+        for i, line in enumerate(self[startindex:]):
+            str_line = []
+            for j, name in enumerate(columns):
+                if name in chapters_txt:
+                    column = chapters_txt[name][i+offsets[name]]
+                else:
+                    value = line.get(name, "")
+                    string = "{0:n}" if isinstance(value, float) else "{0}"
+                    column = string.format(value)
+                self.columns_len[j] = max(self.columns_len[j], len(column))
+                str_line.append(column)
+            str_matrix.append(str_line)
+
+        if startindex == 0 and self.log_header:
+            header = []
+            nlines = 1
+            if len(self.chapters) > 0:
+                nlines += max(map(len, chapters_txt.values())) - len(self) + 1
+            header = [[] for i in xrange(nlines)]
+            for j, name in enumerate(columns):
+                if name in chapters_txt:
+                    length = max(len(line.expandtabs()) for line in chapters_txt[name])
+                    blanks = nlines - 2 - offsets[name]
+                    for i in xrange(blanks):
+                        header[i].append(" " * length)
+                    header[blanks].append(name.center(length))
+                    header[blanks+1].append("-" * length)
+                    for i in xrange(offsets[name]):
+                        header[blanks+2+i].append(chapters_txt[name][i])
+                else:
+                    length = max(len(line[j].expandtabs()) for line in str_matrix)
+                    for line in header[:-1]:
+                        line.append(" " * length)
+                    header[-1].append(name)
+            str_matrix = chain(header, str_matrix)
+
+        template = "\t".join("{%i:<%i}" % (i, l) for i, l in enumerate(self.columns_len))
+        text = [template.format(*line) for line in str_matrix]
+
+        return text
+
+    def __str__(self, startindex=0):
+        text = self.__txt__(startindex)
+        return "\n".join(text)
+
+
+class HallOfFame(object):
+    """The hall of fame contains the best individual that ever lived in the
+    population during the evolution. It is lexicographically sorted at all
+    time so that the first element of the hall of fame is the individual that
+    has the best first fitness value ever seen, according to the weights
+    provided to the fitness at creation time.
+    
+    The insertion is made so that old individuals have priority on new
+    individuals. A single copy of each individual is kept at all time, the
+    equivalence between two individuals is made by the operator passed to the
+    *similar* argument.
+
+    :param maxsize: The maximum number of individual to keep in the hall of
+                    fame.
+    :param similar: An equivalence operator between two individuals, optional.
+                    It defaults to operator :func:`operator.eq`.
+    
+    The class :class:`HallOfFame` provides an interface similar to a list
+    (without being one completely). It is possible to retrieve its length, to
+    iterate on it forward and backward and to get an item or a slice from it.
+    """
+    def __init__(self, maxsize, similar=eq):
+        self.maxsize = maxsize
+        self.keys = list()
+        self.items = list()
+        self.similar = similar
+    
+    def update(self, population):
+        """Update the hall of fame with the *population* by replacing the
+        worst individuals in it by the best individuals present in
+        *population* (if they are better). The size of the hall of fame is
+        kept constant.
+        
+        :param population: A list of individual with a fitness attribute to
+                           update the hall of fame with.
+        """
+        if len(self) == 0 and self.maxsize !=0:
+            # Working on an empty hall of fame is problematic for the
+            # "for else"
+            self.insert(population[0])
+        
+        for ind in population:
+            if ind.fitness > self[-1].fitness or len(self) < self.maxsize:
+                for hofer in self:
+                    # Loop through the hall of fame to check for any
+                    # similar individual
+                    if self.similar(ind, hofer):
+                        break
+                else:
+                    # The individual is unique and strictly better than
+                    # the worst
+                    if len(self) >= self.maxsize:
+                        self.remove(-1)
+                    self.insert(ind)
+    
+    def insert(self, item):
+        """Insert a new individual in the hall of fame using the
+        :func:`~bisect.bisect_right` function. The inserted individual is
+        inserted on the right side of an equal individual. Inserting a new 
+        individual in the hall of fame also preserve the hall of fame's order.
+        This method **does not** check for the size of the hall of fame, in a
+        way that inserting a new individual in a full hall of fame will not
+        remove the worst individual to maintain a constant size.
+        
+        :param item: The individual with a fitness attribute to insert in the
+                     hall of fame.
+        """
+        item = deepcopy(item)
+        i = bisect_right(self.keys, item.fitness)
+        self.items.insert(len(self) - i, item)
+        self.keys.insert(i, item.fitness)
+    
+    def remove(self, index):
+        """Remove the specified *index* from the hall of fame.
+        
+        :param index: An integer giving which item to remove.
+        """
+        del self.keys[len(self) - (index % len(self) + 1)]
+        del self.items[index]
+    
+    def clear(self):
+        """Clear the hall of fame."""
+        del self.items[:]
+        del self.keys[:]
+
+    def __len__(self):
+        return len(self.items)
+
+    def __getitem__(self, i):
+        return self.items[i]
+
+    def __iter__(self):
+        return iter(self.items)
+
+    def __reversed__(self):
+        return reversed(self.items)
+    
+    def __str__(self):
+        return str(self.items)
+
+
+class ParetoFront(HallOfFame):
+    """The Pareto front hall of fame contains all the non-dominated individuals
+    that ever lived in the population. That means that the Pareto front hall of
+    fame can contain an infinity of different individuals.
+    
+    :param similar: A function that tels the Pareto front whether or not two
+                    individuals are similar, optional.
+    
+    The size of the front may become very large if it is used for example on
+    a continuous function with a continuous domain. In order to limit the number
+    of individuals, it is possible to specify a similarity function that will
+    return :data:`True` if the genotype of two individuals are similar. In that
+    case only one of the two individuals will be added to the hall of fame. By
+    default the similarity function is :func:`operator.eq`.
+    
+    Since, the Pareto front hall of fame inherits from the :class:`HallOfFame`, 
+    it is sorted lexicographically at every moment.
+    """
+    def __init__(self, similar=eq):
+        HallOfFame.__init__(self, None, similar)
+    
+    def update(self, population):
+        """Update the Pareto front hall of fame with the *population* by adding 
+        the individuals from the population that are not dominated by the hall
+        of fame. If any individual in the hall of fame is dominated it is
+        removed.
+        
+        :param population: A list of individual with a fitness attribute to
+                           update the hall of fame with.
+        """
+        for ind in population:
+            is_dominated = False
+            has_twin = False
+            to_remove = []
+            for i, hofer in enumerate(self):    # hofer = hall of famer
+                if hofer.fitness.dominates(ind.fitness):
+                    is_dominated = True
+                    break
+                elif ind.fitness.dominates(hofer.fitness):
+                    to_remove.append(i)
+                elif ind.fitness == hofer.fitness and self.similar(ind, hofer):
+                    has_twin = True
+                    break
+            
+            for i in reversed(to_remove):       # Remove the dominated hofer
+                self.remove(i)
+            if not is_dominated and not has_twin:
+                self.insert(ind)
+
+__all__ = ['HallOfFame', 'ParetoFront', 'History', 'Statistics', 'MultiStatistics', 'Logbook']
+
+if __name__ == "__main__":
+    import doctest
+    from operator import itemgetter
+
+    import numpy
+    doctest.run_docstring_examples(Statistics, globals())
+    doctest.run_docstring_examples(Statistics.register, globals())
+    doctest.run_docstring_examples(Statistics.compile, globals())
+
+    doctest.run_docstring_examples(MultiStatistics, globals())
+    doctest.run_docstring_examples(MultiStatistics.register, globals())
+    doctest.run_docstring_examples(MultiStatistics.compile, globals())

python/isaac/autotuning/external/sklearn/setup.py

@@ -0,0 +1,22 @@
+import os
+
+import numpy
+from numpy.distutils.misc_util import Configuration
+
+
+def configuration(parent_package="", top_path=None):
+    config = Configuration("tree", parent_package, top_path)
+    libraries = []
+    if os.name == 'posix':
+        libraries.append('m')
+    config.add_extension("_tree",
+                         sources=["_tree.c"],
+                         include_dirs=[numpy.get_include()],
+                         libraries=libraries,
+                         extra_compile_args=["-O3"])
+    return config
+
+if __name__ == "__main__":
+    from numpy.distutils.core import setup
+    setup(**configuration().todict())
+    print 'aaa'

python/isaac/autotuning/model.py

@@ -1,4 +1,4 @@
-from isaac.external.forest import RandomForestRegressor
+from external.sklearn.forest import RandomForestRegressor
 import numpy as np
 
 def gmean(a, axis=0, dtype=None):

python/isaac/autotuning/optimize.py

@@ -5,10 +5,9 @@ from copy import deepcopy
 from sys import stdout
 from itertools import product
 
-from deap import algorithms
+from external.deap import base
-from deap import base
+from external.deap import creator
-from deap import creator
+from external.deap import tools as deap_tools
-from deap import tools as deap_tools
 
 from numpy import cumsum
 
@@ -74,6 +73,7 @@ def genetic(template, sizes, context, naccept=200, niter = 1000, cxpb=0.4, mutpb
     def evaluate(genome):
         idx = tuple(genome)
         if idx not in cache:
+            print decode(genome)
             cache[idx] = tools.benchmark(template, decode(genome), tree)
         return cache[idx],
         

python/isaac/autotuning/tune.py

@@ -4,10 +4,11 @@ from itertools import chain, product
 from numpy import argsort, argmax
 from operator import mul
 import isaac as sc
-from isaac.external.forest import RandomForestRegressor
+from external.sklearn.forest import RandomForestRegressor
 import optimize, tools, model
-
 from json import encoder
+import json
+
 encoder.FLOAT_REPR = lambda o: format(o, '.2f')
 encoder.separators = (',',':')
 
@@ -21,8 +22,7 @@ def pow2range(a, b):
 
 
 def tune(device, operation, json_path): 
-    #List devices
+    #Context
-    platforms = sc.driver.get_platforms()
     context = sc.driver.context(device)
     
     #List of size tuples to use
@@ -119,6 +119,7 @@ def tune(device, operation, json_path):
 
     
     #Export to JSON
+    json_path = tools.sanitize(device.name) + '.json' if not json_path else json_path
     if os.path.isfile(json_path):
         json_data = json.load(open(json_path, 'r'))
     else:
@@ -138,38 +139,4 @@ def tune(device, operation, json_path):
                             'value': e.tree_.value[:,:,0].astype('float64').tolist()} for e in clf.estimators_]
     D['profiles'] = [map(int, x) for x in profiles]
     json.dump(json_data, open(json_path,'w'))
-    
-
-def parse_arguments():
-    platforms = sc.driver.get_platforms()
-    devices = [d for platform in platforms for d in platform.get_devices()]
-    #Command line arguments
-    parser = argparse.ArgumentParser()
-    parser.add_argument("-d", "--device", default=0, type=int, help='Device to tune for')
-    parser.add_argument("-o", "--operation", type=str, required=True, help='Operation to tune for')
-    parser.add_argument("-j", "--json", default='', type=str)
-    args = parser.parse_args()
-    
-    device = devices[int(args.device)]
-    print("----------------")
-    print("Devices available:")
-    print("----------------")
-    for (i, d) in enumerate(devices):
-        selected = '[' + ('x' if device==d else ' ') + ']'
-        print selected , '-',  sc.driver.device_type_to_string(d.type), '-', d.name, 'on', d.platform.name
-    print("----------------")
-    
-    
-    operation = {'axpy': sc.templates.axpy, 'dot': sc.templates.dot,
-                 'ger': sc.templates.ger, 'gemv_n': sc.templates.gemv_n, 'gemv_t': sc.templates.gemv_t,
-                 'gemm_nn': sc.templates.gemm_nn, 'gemm_tn': sc.templates.gemm_tn, 'gemm_nt': sc.templates.gemm_nt, 'gemm_tt':sc.templates.gemm_tt}[args.operation]
-    json = tools.sanitize(device.name) + '.json' if not args.json else args.json
-        
-    return (device, operation, json)
-        
-            
-if __name__ == "__main__":
-    sc.driver.default.queue_properties = sc.driver.PROFILING_ENABLE
-    args = parse_arguments()
-    tune(*args)
 

python/setup.py

@@ -143,7 +143,7 @@ def main():
                     libraries=libraries)]
     
     #External
-    extensions += [Extension('external._tree',
+    extensions += [Extension('autotuning.external.sklearn._tree',
                              ['external/sklearn/_tree.c'],
                              include_dirs = [numpy_include])]
 
@@ -155,7 +155,7 @@ def main():
                 author='Philippe Tillet',
                 author_email='ptillet@g.harvard.edu',
                 license='MPL 2.0',
-                packages=['isaac','isaac.external'],
+                packages=['isaac','isaac.autotuning', 'isaac.autotuning.external', 'isaac.autotuning.external.deap', 'isaac.autotuning.external.deap.tools', 'isaac.autotuning.external.sklearn'],
                 ext_package="isaac",
                 ext_modules=extensions,
                 cmdclass={'build_py': build_py, 'build_ext': build_ext_subclass},

tune/android/android.txt

@@ -0,0 +1,3 @@
+title=Isaac
+author=Philippe Tillet
+orientation=portrait

tune/android/isaac.ini

@@ -0,0 +1,8 @@
+[hardware]
+device = Intel(R) HD Graphics
+
+[autotuning]
+blas3 = intermediate
+blas2 = intermediate
+blas1 = intermediate
+

tune/android/isaac.kv

@@ -0,0 +1,42 @@
+#:kivy 1.8.0
+
+
+<TunerScreen>:
+    ScrollView:
+        do_scroll_x: False
+        do_scroll_y: False if root.fullscreen else (content.height > root.height - dp(16))
+        AnchorLayout:
+            size_hint_y: None
+            height: root.height if root.fullscreen else max(root.height, content.height)
+            GridLayout:
+                id: content
+                cols: 1
+                spacing: '8dp'
+                padding: '8dp'
+                size_hint: (1, 1) if root.fullscreen else (.8, None)
+                height: self.height if root.fullscreen else self.minimum_height
+
+BoxLayout:
+    orientation: 'vertical'
+
+
+    ActionBar:
+        
+        ActionView:
+            id: av
+            ActionPrevious:
+                with_previous: False
+
+            ActionToggleButton:
+                id: Tune
+                text: 'Tune'
+                on_release: app.show_tune()
+                group: 'menu'
+    
+            ActionToggleButton:
+                text: 'Settings'
+                on_release: app.open_settings()
+                group: 'menu'
+                
+    ScreenManager:
+        id: sm

tune/android/main.py

@@ -0,0 +1,116 @@
+from os.path import dirname, realpath, join
+
+from kivy.uix.boxlayout import BoxLayout
+from kivy.uix.label import Label
+from kivy.uix.checkbox import CheckBox
+from kivy.properties import BooleanProperty
+from kivy.app import App
+from kivy.lang import Builder
+from kivy.properties import NumericProperty, StringProperty, BooleanProperty,ListProperty
+from kivy.uix.screenmanager import Screen
+from kivy.uix.settings import SettingsWithNoMenu
+
+import isaac as sc
+import json
+
+from isaac.autotuning.tune import tune
+
+__version__ = '1.0'
+
+class IsaacScreen(Screen):
+    fullscreen = BooleanProperty(False)
+
+    def add_widget(self, *args):
+        if 'content' in self.ids:
+            return self.ids.content.add_widget(*args)
+        return super(IsaacScreen, self).add_widget(*args)
+
+class IsaacHandler:
+    
+    def __init__(self):
+        platforms = sc.driver.get_platforms()
+        self.devices = [d for platform in platforms for d in platform.get_devices()]
+
+     
+class IsaacApp(App):
+    
+    screen_names = ListProperty([])
+
+
+    def build(self):        
+        self.tuner = 'ISAAC Kernels Tuner'
+
+        #Settings
+        self.settings_cls = SettingsWithNoMenu
+        self.use_kivy_settings = False
+
+        #Screen Manager
+        self.screen_names = ['Tune']
+        self.screens = {}
+        current_directory = dirname(realpath(__file__))
+        for name in self.screen_names:
+            path = join(current_directory, 'screens',  '{}.kv'.format(name.lower()))
+            self.screens[name] = Builder.load_file(path)
+            
+        #Default view
+        self.show_tune()
+    
+    def start_tuning(self):
+        #FIXME: will be buggy if two devices from two different platforms have the same name
+        device = next(x for x in self.isaac_handler.devices if x.name==self.config.get('hardware', 'device'))
+        operation = sc.templates.axpy
+        json_path = ''
+        #FIXME: Move profiling logics into tuning
+        sc.driver.default.queue_properties = sc.driver.PROFILING_ENABLE
+        tune(device, operation, json_path)
+        
+    def show_benchmark(self):
+        pass
+        
+    def show_tune(self):
+        if self.root.ids.sm.current != 'Tune':
+            self.root.ids.sm.switch_to(self.screens['Tune'], direction='left')
+        
+    
+    def display_settings(self, settings):
+        if 'Settings' not in self.screens:
+            self.screens['Settings'] = Screen(name='Settings')
+            self.screens['Settings'].add_widget(settings)
+        self.root.ids.sm.switch_to(self.screens['Settings'], direction='left')
+
+    def build_config(self, config):
+        self.isaac_handler = IsaacHandler()
+        config.setdefaults('hardware', {'device':  self.isaac_handler.devices[0].name})
+        config.setdefaults('autotuning', {'blas1': 'Intermediate', 'blas2': 'Intermediate', 'blas3': 'Intermediate'})
+        
+    def build_settings(self, settings):
+        
+        layout = [{'type': 'title',
+                  'title': 'Hardware'},
+                  {'type': 'options',
+                  'title': 'Device',
+                  'section': 'hardware',
+                  'key': 'device',
+                  'options': [device.name for device in self.isaac_handler.devices]},
+                  {'type': 'title',
+                  'title': 'Auto-tuning'}]
+        
+        for operation in ['BLAS1', 'BLAS2', 'BLAS3']:
+            layout +=   [{'type': 'options',
+                        'desc': 'Desired level of auto-tuning for ' + operation,
+                        'title': operation,
+                        'section': 'autotuning',
+                        'key': operation.lower(),
+                        'options': ['Simple', 'Intermediate', 'Full']}]
+        
+        settings.add_json_panel('Settings',
+                                self.config,
+                                data=json.dumps(layout))
+                                
+    def close_settings(self, *args):
+        pass
+
+        
+            
+if __name__ == '__main__':
+    IsaacApp().run()

tune/android/screens/tune.kv

@@ -0,0 +1,16 @@
+IsaacScreen:
+	name: 'Tune'
+	fullscreen: True
+
+    BoxLayout:
+        orientation: 'vertical'
+        
+        Button:
+            text: 'Start tuning'
+            pos: 0, 0
+            size_hint: 1, .1
+            on_release: app.start_tuning()
+        
+        Label:
+            text: ''
+            

tune/android/settings.py

@@ -0,0 +1 @@
+

tune/main.py

@@ -0,0 +1,35 @@
+import argparse
+import isaac as sc
+from isaac.autotuning.tune import tune
+
+def parse_arguments():
+    platforms = sc.driver.get_platforms()
+    devices = [d for platform in platforms for d in platform.get_devices()]
+    #Command line arguments
+    parser = argparse.ArgumentParser()
+    parser.add_argument("-d", "--device", default=0, type=int, help='Device to tune for')
+    parser.add_argument("-o", "--operation", type=str, required=True, help='Operation to tune for')
+    parser.add_argument("-j", "--json", default='', type=str)
+    args = parser.parse_args()
+    
+    device = devices[int(args.device)]
+    print("----------------")
+    print("Devices available:")
+    print("----------------")
+    for (i, d) in enumerate(devices):
+        selected = '[' + ('x' if device==d else ' ') + ']'
+        print selected , '-',  sc.driver.device_type_to_string(d.type), '-', d.name, 'on', d.platform.name
+    print("----------------")
+    
+    
+    operation = {'axpy': sc.templates.axpy, 'dot': sc.templates.dot,
+                 'ger': sc.templates.ger, 'gemv_n': sc.templates.gemv_n, 'gemv_t': sc.templates.gemv_t,
+                 'gemm_nn': sc.templates.gemm_nn, 'gemm_tn': sc.templates.gemm_tn, 'gemm_nt': sc.templates.gemm_nt, 'gemm_tt':sc.templates.gemm_tt}[args.operation]
+    
+    return (device, operation, args.json)
+        
+            
+if __name__ == "__main__":
+    sc.driver.default.queue_properties = sc.driver.PROFILING_ENABLE
+    args = parse_arguments()
+    tune(*args)