class Algorithm(object):
    
    __metaclass__ = ABCMeta
    
    def __init__(self,
                 problem,
                 evaluator=None,
                 log_frequency=None,
                 **kwargs):
        super(Algorithm, self).__init__()
        self.problem = problem
        self.evaluator = evaluator
        self.log_frequency = log_frequency
        self.nfe = 0
        
        if self.evaluator is None:
            from .config import PlatypusConfig
            self.evaluator = PlatypusConfig.default_evaluator
            
        if self.log_frequency is None:
            from .config import PlatypusConfig
            self.log_frequency = PlatypusConfig.default_log_frequency
    
    @abstractmethod
    def step(self):
        raise NotImplementedError("method not implemented")
    
    def evaluate_all(self, solutions):
        unevaluated = [s for s in solutions if not s.evaluated]
        
        jobs = [_EvaluateJob(s) for s in unevaluated]
        results = self.evaluator.evaluate_all(jobs)
            
        # if needed, update the original solution with the results
        for i, result in enumerate(results):
            if unevaluated[i] != result.solution:
                unevaluated[i].variables[:] = result.solution.variables[:]
                unevaluated[i].objectives[:] = result.solution.objectives[:]
                unevaluated[i].constraints[:] = result.solution.constraints[:]
                unevaluated[i].constraint_violation = result.solution.constraint_violation
                unevaluated[i].feasible = result.solution.feasible
                unevaluated[i].evaluated = result.solution.evaluated
        
        self.nfe += len(unevaluated)
    
    def run(self, condition, callback=None):
        if isinstance(condition, int):
            condition = MaxEvaluations(condition)
            
        if isinstance(condition, TerminationCondition):
            condition.initialize(self)
            
        last_log = self.nfe
        start_time = time.time()
        
        LOGGER.log(logging.INFO, "%s starting", type(self).__name__)

        while not condition(self):
            self.step()
            
            if self.log_frequency is not None and self.nfe >= last_log + self.log_frequency:
                LOGGER.log(logging.INFO,
                           "%s running; NFE Complete: %d, Elapsed Time: %s",
                           type(self).__name__,
                           self.nfe,
                           datetime.timedelta(seconds=time.time()-start_time))

            if callback is not None:
                callback(self)
                
        LOGGER.log(logging.INFO,
                   "%s finished; Total NFE: %d, Elapsed Time: %s",
                   type(self).__name__,
                   self.nfe,
                   datetime.timedelta(seconds=time.time()-start_time)) class ParetoDominance(Dominance):
    """Pareto dominance with constraints.
    
    If either solution violates constraints, then the solution with a smaller
    constraint violation is preferred.  If both solutions are feasible, then
    Pareto dominance is used to select the preferred solution.
    """
    
    def __init__(self):
        super(ParetoDominance, self).__init__()
    
    def compare(self, solution1, solution2):
        problem = solution1.problem
        
        # first check constraint violation
        if problem.nconstrs > 0 and solution1.constraint_violation != solution2.constraint_violation:
            if solution1.constraint_violation == 0.0:
                return -1
            elif solution2.constraint_violation == 0.0:
                return 1
            elif solution1.constraint_violation < solution2.constraint_violation:
                return -1
            elif solution2.constraint_violation < solution1.constraint_violation:
                return 1
        
        # then use Pareto dominance on the objectives
        dominate1 = False
        dominate2 = False
        
        for i in range(problem.nobjs):
            o1 = solution1.objectives[i]
            o2 = solution2.objectives[i]
            
            if problem.directions[i] == Problem.MAXIMIZE:
                o1 = -o1
                o2 = -o2

            if o1 < o2:
                dominate1 = True
                    
                if dominate2:
                    return 0
            elif o1 > o2:
                dominate2 = True
                
                if dominate1:
                    return 0
            
        if dominate1 == dominate2:
            return 0
        elif dominate1:
            return -1
        else:
            return 1 class AttributeDominance(Dominance):
    
    def __init__(self, getter, larger_preferred=True):
        super(AttributeDominance, self).__init__()
        self.larger_preferred = larger_preferred
        
        if hasattr(getter, "__call__"):
            self.getter = getter
        else:
            self.getter = operator.attrgetter(getter)
        
    def compare(self, solution1, solution2):
        a = self.getter(solution1)
        b = self.getter(solution2)
        
        if self.larger_preferred:
            a = -a
            b = -b
        
        if a < b:
            return -1
        elif a > b:
            return 1
        else:
            return 0 class AttributeDominance(Dominance):
    
    def __init__(self, getter, larger_preferred=True):
        super(AttributeDominance, self).__init__()
        self.larger_preferred = larger_preferred
        
        if hasattr(getter, "__call__"):
            self.getter = getter
        else:
            self.getter = operator.attrgetter(getter)
        
    def compare(self, solution1, solution2):
        a = self.getter(solution1)
        b = self.getter(solution2)
        
        if self.larger_preferred:
            a = -a
            b = -b
        
        if a < b:
            return -1
        elif a > b:
            return 1
        else:
            return 0 def nondominated_sort(solutions):
    """Fast non-dominated sorting.
    
    Performs fast non-dominated sorting on a collection of solutions.  The
    solutions will be assigned the following attributes:
    
    1. :code:`rank` - The index of the non-dominated front containing the
       solution.  Rank 0 stores all non-dominated solutions.
       
    2. :code:`crowding_distance` - The crowding distance of the given solution.
       Larger values indicate less crowding near the solution.
       
    Parameters
    ----------
    solutions : iterable
        The collection of solutions
    """
    rank = 0
    
    while len(solutions) > 0:
        archive = Archive()
        archive += solutions
        
        for solution in archive:
            solution.rank = rank
            
        crowding_distance(archive)
            
        solutions = [x for x in solutions if x not in archive]
        rank += 1 def nondominated_prune(solutions, size):
    """Prunes a population using non-dominated sorting.
    
    Similar to :code:`nondominated_truncate`, except the crowding distance
    is recomputed after each solution is removed.  
    
    Parameters
    ----------
    solutions : iterable
        The collection of solutions that have been non-domination sorted
    size: int
        The size of the truncated result
    """
    result, remaining = nondominated_split(solutions, size)
    
    while len(result) + len(remaining) > size:
        crowding_distance(remaining)
        remaining = sorted(remaining, key=crowding_distance_key)
        del remaining[0]
        
    return result + remaining def nondominated_truncate(solutions, size):
    """Truncates a population using non-dominated sorting.
    
    Truncates a population to the given size.  The resulting population is
    filled with the first N-1 fronts.  The Nth front is too large and must be
    split using crowding distance.
    
    Parameters
    ----------
    solutions : iterable
        The collection of solutions that have been non-domination sorted
    size: int
        The size of the truncated result
    """
    result = sorted(solutions, key=functools.cmp_to_key(nondominated_cmp)) 
    return result[:size] def nondominated_split(solutions, size):
    """Identify the front that must be truncated.
    
    When truncating a population to a fixed size using non-dominated sorting,
    identify the front N that can not completely fit within the truncated
    result.  Returns a tuple :code:`(first, last)`, where :code:`first`
    contains all solutions in the first N fronts (those solutions that will
    definitely remain after truncation), and the :code:`last` is the front
    that must be truncated.  :code:`last` can be empty.
    
    Parameters
    ----------
    solutions : iterable
        The collection of solutions that have been non-dominated sorted
    size : int
        The size of the truncated result
    """
    result = []
    rank = 0
    
    while len(result) < size:
        front = [x for x in solutions if x.rank==rank]
        
        if len(front) == 0:
            break
        
        if len(result)+len(front) <= size:
            result.extend(front)
        else:
            return (result, front)
        
        rank += 1
        
    return result, [] def crowding_distance(solutions):
    """Calculates crowding distance for a non-dominated front.
    
    Computes the crowding distance for a single non-dominated front.  It is
    assumed all solutions are non-dominated.  This method assigns the attribute
    :code:`crowding_distance` to all solutions.
    
    Parameters
    ----------
    solutions : iterable
        The collection of solutions
    """
    for solution in solutions:
        solution.crowding_distance = 0.0
        
    solutions = unique(solutions)
    
    if len(solutions) < 3:
        for solution in solutions:
            solution.crowding_distance = POSITIVE_INFINITY
    else:
        nobjs = solutions[0].problem.nobjs
            
        for i in range(nobjs):
            sorted_solutions = sorted(solutions, key=functools.partial(objective_key, index=i))
            min_value = sorted_solutions[0].objectives[i]
            max_value = sorted_solutions[-1].objectives[i]
            
            sorted_solutions[0].crowding_distance += POSITIVE_INFINITY
            sorted_solutions[-1].crowding_distance += POSITIVE_INFINITY
            
            for j in range(1, len(sorted_solutions)-1):
                if max_value - min_value < EPSILON:
                    sorted_solutions[j].crowding_distance = POSITIVE_INFINITY
                else:
                    diff = sorted_solutions[j+1].objectives[i] - sorted_solutions[j-1].objectives[i]
                    sorted_solutions[j].crowding_distance += diff / (max_value - min_value) EPSILON = sys.float_info.epsilon POSITIVE_INFINITY = float("inf") class Archive(object):
    """An archive only containing non-dominated solutions."""
    
    def __init__(self, dominance = ParetoDominance()):
        super(Archive, self).__init__()
        self._dominance = dominance
        self._contents = []
        
    def add(self, solution):
        flags = [self._dominance.compare(solution, s) for s in self._contents]
        dominates = [x > 0 for x in flags]
        nondominated = [x == 0 for x in flags]
        
        if any(dominates):
            return False
        else:
            self._contents = list(itertools.compress(self._contents, nondominated)) + [solution]
            return True
        
    def append(self, solution):
        self.add(solution)
        
    def extend(self, solutions):
        for solution in solutions:
            self.append(solution)
            
    def remove(self, solution):
        try:
            self._contents.remove(solution)
            return True
        except ValueError:
            return False
    
    def __len__(self):
        return len(self._contents)
    
    def __getitem__(self, key):
        return self._contents[key]
            
    def __iadd__(self, other):
        if hasattr(other, "__iter__"):
            for o in other:
                self.add(o)
        else:
            self.add(other)
            
        return self
    
    def __iter__(self):
        return iter(self._contents) class EpsilonDominance(Dominance):
    """Epsilon dominance.
    
    Similar to Pareto dominance except if the two solutions are contained
    within the same epsilon-box, the solution closer to the optimal corner
    or the box is preferred.
    """
    
    def __init__(self, epsilons):
        super(EpsilonDominance, self).__init__()
        
        if hasattr(epsilons, "__getitem__"):
            self.epsilons = epsilons
        else:
            self.epsilons = [epsilons]
        
    def same_box(self, solution1, solution2):
        problem = solution1.problem
        
        # first check constraint violation
        if problem.nconstrs > 0 and solution1.constraint_violation != solution2.constraint_violation:
            if solution1.constraint_violation == 0:
                return False
            elif solution2.constraint_violation == 0:
                return False
            elif solution1.constraint_violation < solution2.constraint_violation:
                return False
            elif solution2.constraint_violation < solution1.constraint_violation:
                return False
        
        # then use epsilon dominance on the objectives
        dominate1 = False
        dominate2 = False
        
        for i in range(problem.nobjs):
            o1 = solution1.objectives[i]
            o2 = solution2.objectives[i]
            
            if problem.directions[i] == Problem.MAXIMIZE:
                o1 = -o1
                o2 = -o2
                
            epsilon = float(self.epsilons[i % len(self.epsilons)])
            i1 = math.floor(o1 / epsilon)
            i2 = math.floor(o2 / epsilon)

            if i1 < i2:
                dominate1 = True
                    
                if dominate2:
                    return False
            elif i1 > i2:
                dominate2 = True
                
                if dominate1:
                    return False
        
        if not dominate1 and not dominate2:
            return True
        else:
            return False
    
    def compare(self, solution1, solution2):
        problem = solution1.problem
        
        # first check constraint violation
        if problem.nconstrs > 0 and solution1.constraint_violation != solution2.constraint_violation:
            if solution1.constraint_violation == 0:
                return -1
            elif solution2.constraint_violation == 0:
                return 1
            elif solution1.constraint_violation < solution2.constraint_violation:
                return -1
            elif solution2.constraint_violation < solution1.constraint_violation:
                return 1
        
        # then use epsilon dominance on the objectives
        dominate1 = False
        dominate2 = False
        
        for i in range(problem.nobjs):
            o1 = solution1.objectives[i]
            o2 = solution2.objectives[i]
            
            if problem.directions[i] == Problem.MAXIMIZE:
                o1 = -o1
                o2 = -o2
                
            epsilon = float(self.epsilons[i % len(self.epsilons)])
            i1 = math.floor(o1 / epsilon)
            i2 = math.floor(o2 / epsilon)

            if i1 < i2:
                dominate1 = True
                    
                if dominate2:
                    return 0
            elif i1 > i2:
                dominate2 = True
                
                if dominate1:
                    return 0
        
        if not dominate1 and not dominate2:
            dist1 = 0.0
            dist2 = 0.0
            
            for i in range(problem.nobjs):
                o1 = solution1.objectives[i]
                o2 = solution2.objectives[i]
            
                if problem.directions[i] == Problem.MAXIMIZE:
                    o1 = -o1
                    o2 = -o2
                
                epsilon = float(self.epsilons[i % len(self.epsilons)])
                i1 = math.floor(o1 / epsilon)
                i2 = math.floor(o2 / epsilon)
                
                dist1 += math.pow(o1 - i1*epsilon, 2.0)
                dist2 += math.pow(o2 - i2*epsilon, 2.0)
            
            if dist1 < dist2:
                return -1
            else:
                return 1
        elif dominate1:
            return -1
        else:
            return 1 class FitnessArchive(Archive):
    
    def __init__(self, fitness, dominance = ParetoDominance(), larger_preferred=True, getter=fitness_key):
        super(FitnessArchive, self).__init__(dominance)
        self.fitness = fitness
        self.larger_preferred = larger_preferred
        self.getter = getter
        
    def truncate(self, size):
        self.fitness(self._contents)
        self._contents = truncate_fitness(self._contents,
                                          size,
                                          larger_preferred=self.larger_preferred,
                                          getter=self.getter) class Solution(object):
    """Class representing a solution to a problem.
    
    Attributes
    ----------
    problem: Problem
        The problem.
    variables: FixedLengthArray of objects
        The values of the variables.
    objectives: FixedLengthArray of float
        The values of the objectives.  These values will only be assigned after
        the solution is evaluated.
    constraints: FixedLengthArray of float
        The values of the constraints.  These values will only be assigned
        after the solution is evaluated.
    constraint_violation: float
        The magnitude of the constraint violation.
    feasible: bool
        True if the solution does not violate any constraints; False otherwise.
    evaluated: bool
        True if the solution is evaluated; False otherwise.
    """
    
    def __init__(self, problem):
        """Creates a new solution for the given problem."""
        super(Solution, self).__init__()
        self.problem = problem
        self.variables = FixedLengthArray(problem.nvars)
        self.objectives = FixedLengthArray(problem.nobjs)
        self.constraints = FixedLengthArray(problem.nconstrs)
        self.constraint_violation = 0.0
        self.evaluated = False
        
    def evaluate(self):
        """Evaluates this solution."""
        self.problem(self)
        
    def __repr__(self):
        return self.__str__()
        
    def __str__(self):
        return "Solution[" + ",".join(list(map(str, self.variables))) + "|" + ",".join(list(map(str, self.objectives))) + "|" + str(self.constraint_violation) + "]"
    
    def __deepcopy__(self, memo):
        """Overridden to avoid cloning the problem definition."""
        result = Solution(self.problem)
        memo[id(self)] = result
        
        for k, v in self.__dict__.items():
            if k != "problem":
                setattr(result, k, copy.deepcopy(v, memo))
                
        return result class HypervolumeFitnessEvaluator(FitnessEvaluator):
    
    def __init__(self,
                 kappa = 0.05,
                 rho = 2.0,
                 dominance = ParetoDominance()):
        super(HypervolumeFitnessEvaluator, self).__init__(kappa = kappa)
        self.rho = rho
        self.dominance = dominance
    
    def calculate_indicator(self, solution1, solution2):
        problem = solution1.problem
        
        if self.dominance.compare(solution1, solution2) < 0:
            return -self.hypervolume(solution1, solution2, problem.nobjs)
        else:
            return self.hypervolume(solution2, solution1, problem.nobjs)
    
    def hypervolume(self, solution1, solution2, d):
        a = solution1.normalized_objectives[d-1]
        
        if solution2 is None:
            b = self.rho
        else:
            b = solution2.normalized_objectives[d-1]
            
        if solution1.problem.directions[d-1] == Problem.MAXIMIZE:
            a = 1.0 - a
            b = 1.0 - b

        if d == 1:
            if a < b:
                return (b - a) / self.rho
            else:
                return 0.0
        else:
            if a < b:
                return self.hypervolume(solution1, None, d-1)*(b-a)/self.rho + \
                       self.hypervolume(solution1, solution2, d-1)*(self.rho-b)/self.rho
            else:
                return self.hypervolume(solution1, solution2, d-1)*(self.rho-a)/self.rho def nondominated_cmp(x, y):
    if x.rank == y.rank:
        if -x.crowding_distance < -y.crowding_distance:
            return -1
        elif -x.crowding_distance > -y.crowding_distance:
            return 1
        else:
            return 0
    else:
        if x.rank < y.rank:
            return -1
        elif x.rank > y.rank:
            return 1
        else:
            return 0 def fitness_key(x):
    return x.fitness def crowding_distance_key(x):
    return x.crowding_distance class AdaptiveGridArchive(Archive):
    
    def __init__(self, capacity, nobjs, divisions, dominance = ParetoDominance()):
        super(AdaptiveGridArchive, self).__init__(dominance)
        self.capacity = capacity
        self.nobjs = nobjs
        self.divisions = divisions
        
        self.adapt_grid()
        
    def add(self, solution):
        # check if the candidate solution dominates or is dominated
        flags = [self._dominance.compare(solution, s) for s in self._contents]
        dominates = [x > 0 for x in flags]
        nondominated = [x == 0 for x in flags]
        
        if any(dominates):
            return False

        self._contents = list(itertools.compress(self._contents, nondominated))

        # archive is empty, add the candidate
        if len(self) == 0:
            self._contents.append(solution)
            self.adapt_grid()
            return True
        
        # temporarily add the candidate solution
        self._contents.append(solution)
        index = self.find_index(solution)
        
        if index < 0:
            self.adapt_grid()
            index = self.find_index(solution)
        else:
            self.density[index] += 1
            
        if len(self) <= self.capacity:
            # keep the candidate if size is less than capacity
            return True
        elif self.density[index] == self.density[self.find_densest()]:
            # reject candidate if in most dense cell
            self.remove(solution)
            return False
        else:
            # keep candidate and remove one from densest cell
            self.remove(self.pick_from_densest())
            return True
        
    def remove(self, solution):
        removed = super(AdaptiveGridArchive, self).remove(solution)
        
        if removed:
            index = self.find_index(solution)
            
            if self.density[index] > 1:
                self.density[index] -= 1
            else:
                self.adapt_grid()
                
        return removed
        
    def adapt_grid(self):
        self.minimum = [POSITIVE_INFINITY]*self.nobjs
        self.maximum = [-POSITIVE_INFINITY]*self.nobjs
        self.density = [0.0]*(self.divisions**self.nobjs)
        
        for solution in self:
            for i in range(self.nobjs):
                self.minimum[i] = min(self.minimum[i], solution.objectives[i])
                self.maximum[i] = max(self.maximum[i], solution.objectives[i])
                
        for solution in self:
            self.density[self.find_index(solution)] += 1
            
    def find_index(self, solution):
        index = 0
        
        for i in range(self.nobjs):
            value = solution.objectives[i]
            
            if value < self.minimum[i] or value > self.maximum[i]:
                return -1
            
            if self.maximum[i] > self.minimum[i]:
                value = (value - self.minimum[i]) / (self.maximum[i] - self.minimum[i])
            else:
                value = 0
            
            temp_index = int(self.divisions * value)
            
            if temp_index == self.divisions:
                temp_index -= 1
                
            index += temp_index * pow(self.divisions, i)
            
        return index
    
    def find_densest(self):
        index = -1
        value = -1
        
        for i in range(len(self)):
            temp_index = self.find_index(self[i])
            temp_value = self.density[temp_index]
            
            if temp_value > value:
                value = temp_value
                index = temp_index
                
        return index
    
    def pick_from_densest(self):
        solution = None
        value = -1
        
        for i in range(len(self)):
            temp_value = self.density[self.find_index(self[i])]
            
            if temp_value > value:
                solution = self[i]
                value = temp_value
                
        return solution class Selector(object):
    
    __metaclass__ = ABCMeta
    
    def __init__(self):
        super(Selector, self).__init__()
        
    def select(self, n, population):
        return list(map(self.select_one, itertools.repeat(population, n)))
        
    @abstractmethod
    def select_one(self, population):
        raise NotImplementedError("method not implemented") class EpsilonBoxArchive(Archive):
    
    def __init__(self, epsilons):
        super(EpsilonBoxArchive, self).__init__(EpsilonDominance(epsilons))
        self.improvements = 0

    def add(self, solution):
        flags = [self._dominance.compare(solution, s) for s in self._contents]
        dominates = [x > 0 for x in flags]
        nondominated = [x == 0 for x in flags]
        dominated = [x < 0 for x in flags]
        not_same_box = [not self._dominance.same_box(solution, s) for s in self._contents]
        
        if any(dominates):
            return False
        else:
            self._contents = list(itertools.compress(self._contents, nondominated)) + [solution]
            
            if dominated and not_same_box:
                self.improvements += 1 class PlatypusError(Exception):
    pass class Problem(object):
    """Class representing a problem.
    
    Attributes
    ----------
    nvars: int
        The number of decision variables
    nobjs: int
        The number of objectives.
    nconstrs: int
        The number of constraints.
    function: callable
        The function used to evaluate the problem.  If no function is given,
        it is expected that the evaluate method is overridden.
    types: FixedLengthArray of Type
        The type of each decision variable.  The type describes the bounds and
        encoding/decoding required for each decision variable.
    directions: FixedLengthArray of int
        The optimization direction of each objective, either MINIMIZE (-1) or
        MAXIMIZE (1)
    constraints: FixedLengthArray of Constraint
        Describes the types of constraints as an equality or inequality.  The
        default requires each constraint value to be 0.
    """
    
    MINIMIZE = -1
    MAXIMIZE = 1
    
    def __init__(self, nvars, nobjs, nconstrs = 0, function=None):
        """Create a new problem.
    
        Problems can be constructed by either subclassing and overriding the
        evaluate method or passing in a function of the form::
        
            def func_name(vars):
                # vars is a list of the decision variable values
                return (objs, constrs)
        
        Parameters
        ----------
        nvars: int
            The number of decision variables.
        nobjs: int
            The number of objectives.
        nconstrs: int (default 0)
            The number of constraints.
        function: callable (default None)
            The function that is used to evaluate the problem.
        """
        super(Problem, self).__init__()
        self.nvars = nvars
        self.nobjs = nobjs
        self.nconstrs = nconstrs
        self.function = function
        self.types = FixedLengthArray(nvars)
        self.directions = FixedLengthArray(nobjs, self.MINIMIZE)
        self.constraints = FixedLengthArray(nconstrs, "==0", _convert_constraint)
        
    def __call__(self, solution):
        """Evaluate the solution.
        
        This method is responsible for decoding the decision variables,
        invoking the evaluate method, and updating the solution.
        
        Parameters
        ----------
        solution: Solution
            The solution to evaluate.
        """
        problem = solution.problem
        solution.variables[:] = [problem.types[i].decode(solution.variables[i]) for i in range(problem.nvars)]
        
        self.evaluate(solution)
        
        solution.variables[:] = [problem.types[i].encode(solution.variables[i]) for i in range(problem.nvars)]
        solution.constraint_violation = sum([abs(f(x)) for (f, x) in zip(solution.problem.constraints, solution.constraints)])
        solution.feasible = solution.constraint_violation == 0.0
        solution.evaluated = True
        
    def evaluate(self, solution):
        """Evaluates the problem.
        
        By default, this method calls the function passed to the constructor.
        Alternatively, a problem can subclass and override this method.  When
        overriding, this method is responsible for updating the objectives
        and constraints stored in the solution.
        
        Parameters
        ----------
        solution: Solution
            The solution to evaluate.
        """
        if self.function is None:
            raise PlatypusError("function not defined")
        
        if self.nconstrs > 0:
            (objs, constrs) = self.function(solution.variables)
        else:
            objs = self.function(solution.variables)
            constrs = []
            
        if not hasattr(objs, "__getitem__"):
            objs = [objs]
            
        if not hasattr(constrs, "__getitem__"):
            constrs = [constrs]
            
        if len(objs) != self.nobjs:
            raise PlatypusError("incorrect number of objectives: expected %d, received %d" % (self.nobjs, len(objs)))
        
        if len(constrs) != self.nconstrs:
            raise PlatypusError("incorrect number of constraints: expected %d, received %d" % (self.nconstrs, len(constrs)))
        
        solution.objectives[:] = objs
        solution.constraints[:] = constrs class TournamentSelector(Selector):
    
    def __init__(self, tournament_size = 2, dominance = ParetoDominance()):
        super(TournamentSelector, self).__init__()
        self.tournament_size = tournament_size
        self.dominance = dominance
    
    def select_one(self, population):
        winner = random.choice(population)
        
        for _ in range(self.tournament_size-1):
            candidate = random.choice(population)
            flag = self.dominance.compare(winner, candidate)
            
            if flag > 0:
                winner = candidate
                
        return winner class RandomGenerator(Generator):
    
    def __init__(self):
        super(RandomGenerator, self).__init__()
    
    def generate(self, problem):
        solution = Solution(problem)
        solution.variables = [x.rand() for x in problem.types]
        return solution class DifferentialEvolution(Variator):
    
    def __init__(self, crossover_rate=0.1, step_size=0.5):
        super(DifferentialEvolution, self).__init__(4)
        self.crossover_rate = crossover_rate
        self.step_size = step_size
        
    def evolve(self, parents):
        result = copy.deepcopy(parents[0])
        problem = result.problem
        jrand = random.randrange(problem.nvars)
        
        for j in range(problem.nvars):
            if random.uniform(0.0, 1.0) <= self.crossover_rate or j == jrand:
                v1 = float(parents[1].variables[j])
                v2 = float(parents[2].variables[j])
                v3 = float(parents[3].variables[j])
                
                y = v3 + self.step_size*(v1 - v2)
                y = clip(y, problem.types[j].min_value, problem.types[j].max_value)
                
                result.variables[j] = y
                result.evaluated = False
                
        return [result] def clip(value, min_value, max_value):
    return max(min_value, min(value, max_value)) class UniformMutation(Mutation):
        
    def __init__(self, probability, perturbation):
        super(UniformMutation, self).__init__()
        self.probability = probability
        self.perturbation = perturbation
        
    def mutate(self, parent):
        result = copy.deepcopy(parent)
        problem = result.problem
        
        for i in range(problem.nvars):
            if random.uniform(0.0, 1.0) <= self.probability:
                type = problem.types[i]
                value = result.variables[i] + (random.uniform(0.0, 1.0) - 0.5) * self.perturbation
                result.variables[i] = clip(value, type.min_value, type.max_value)
                result.evaluated = False
                
        return result class NonUniformMutation(Mutation):
        
    def __init__(self, probability, perturbation, max_iterations, algorithm):
        super(NonUniformMutation, self).__init__()
        self.probability = probability
        self.perturbation = perturbation
        self.max_iterations = max_iterations
        self.algorithm = algorithm
        
    def _delta(self, difference):
        current_iteration = self.algorithm.nfe / self.algorithm.swarm_size
        fraction = min(1.0, current_iteration / float(self.max_iterations))
        return difference * (1.0 - math.pow(random.uniform(0.0, 1.0), math.pow(1.0 - fraction, self.perturbation)))
        
    def mutate(self, parent):
        result = copy.deepcopy(parent)
        problem = result.problem
        
        for i in range(problem.nvars):
            if random.uniform(0.0, 1.0) <= self.probability:
                type = problem.types[i]
                value = result.variables[i]
                
                if bool(random.getrandbits(1)):
                    value += self._delta(type.max_value - value)
                else:
                    value += self._delta(type.min_value - value)
                    
                result.variables[i] = clip(value, type.min_value, type.max_value)
                result.evaluated = False
                
        return result class GAOperator(Variator):
    
    def __init__(self, variation, mutation):
        super(GAOperator, self).__init__(variation.arity)
        self.variation = variation
        self.mutation = mutation
        
    def evolve(self, parents):
        return list(map(self.mutation.evolve, self.variation.evolve(parents))) class SBX(Variator):
     
    def __init__(self, probability = 1.0, distribution_index = 15.0):
        super(SBX, self).__init__(2)
        self.probability = probability
        self.distribution_index = distribution_index
         
    def evolve(self, parents):
        child1 = copy.deepcopy(parents[0])
        child2 = copy.deepcopy(parents[1])
        
        if random.uniform(0.0, 1.0) <= self.probability:
            problem = child1.problem
            nvars = problem.nvars
            
            for i in range(nvars):
                if isinstance(problem.types[i], Real):
                    if random.uniform(0.0, 1.0) <= 0.5:
                        x1 = float(child1.variables[i])
                        x2 = float(child2.variables[i])
                        lb = problem.types[i].min_value
                        ub = problem.types[i].max_value
                        
                        x1, x2 = self.sbx_crossover(x1, x2, lb, ub)
                        
                        child1.variables[i] = x1
                        child2.variables[i] = x2
                        child1.evaluated = False
                        child2.evaluated = False
                    
        return [child1, child2]
                    
    def sbx_crossover(self, x1, x2, lb, ub):
        dx = x2 - x1
        
        if dx > EPSILON:
            if x2 > x1:
                y2 = x2
                y1 = x1
            else:
                y2 = x1
                y1 = x2

            beta = 1.0 / (1.0 + (2.0 * (y1 - lb) / (y2 - y1)))
            alpha = 2.0 - pow(beta, self.distribution_index + 1.0)
            rand = random.uniform(0.0, 1.0)

            if rand <= 1.0 / alpha:
                alpha = alpha * rand
                betaq = pow(alpha, 1.0 / (self.distribution_index + 1.0))
            else:
                alpha = alpha * rand;
                alpha = 1.0 / (2.0 - alpha)
                betaq = pow(alpha, 1.0 / (self.distribution_index + 1.0))
            
            x1 = 0.5 * ((y1 + y2) - betaq * (y2 - y1))
            beta = 1.0 / (1.0 + (2.0 * (ub - y2) / (y2 - y1)));
            alpha = 2.0 - pow(beta, self.distribution_index + 1.0);
            
            if rand <= 1.0 / alpha:
                alpha = alpha * rand;
                betaq = pow(alpha, 1.0 / (self.distribution_index + 1.0));
            else:
                alpha = alpha * rand;
                alpha = 1.0 / (2.0 - alpha);
                betaq = pow(alpha, 1.0 / (self.distribution_index + 1.0));
            
            x2 = 0.5 * ((y1 + y2) + betaq * (y2 - y1));
            
            # randomly swap the values
            if bool(random.getrandbits(1)):
                x1, x2 = x2, x1
            
            x1 = clip(x1, lb, ub)
            x2 = clip(x2, lb, ub)
            
        return x1, x2 class PM(Mutation):
    
    def __init__(self, probability = 1, distribution_index = 20.0):
        super(PM, self).__init__()
        self.probability = probability
        self.distribution_index = distribution_index
        
    def mutate(self, parent):
        child = copy.deepcopy(parent)
        problem = child.problem
        probability = self.probability
        
        if isinstance(probability, int):
            probability /= float(len([t for t in problem.types if isinstance(t, Real)]))
            
        for i in range(len(child.variables)):
            if isinstance(problem.types[i], Real):
                if random.uniform(0.0, 1.0) <= probability:
                    child.variables[i] = self.pm_mutation(float(child.variables[i]),
                                                          problem.types[i].min_value,
                                                          problem.types[i].max_value)

                    child.evaluated = False
        
        return child
    
    def pm_mutation(self, x, lb, ub):
        u = random.uniform(0, 1)
        dx = ub - lb
        
        if u < 0.5:
            bl = (x - lb) / dx
            b = 2.0*u + (1.0 - 2.0*u)*pow(1.0 - bl, self.distribution_index + 1.0)
            delta = pow(b, 1.0 / (self.distribution_index + 1.0)) - 1.0
        else:
            bu = (ub - x) / dx
            b = 2.0*(1.0 - u) + 2.0*(u - 0.5)*pow(1.0 - bu, self.distribution_index + 1.0)
            delta = 1.0 - pow(b, 1.0 / (self.distribution_index + 1.0))
            
        x = x + delta*dx
        x = clip(x, lb, ub)
            
        return x class UM(Mutation):
    """Uniform mutation."""
    
    def __init__(self, probability = 1):
        super(UM, self).__init__()
        self.probability = probability
        
    def mutate(self, parent):
        child = copy.deepcopy(parent)
        problem = child.problem
        probability = self.probability
        
        if isinstance(probability, int):
            probability /= float(len([t for t in problem.types if isinstance(t, Real)]))
        
        for i in range(len(child.variables)):
            if isinstance(problem.types[i], Real):
                if random.uniform(0.0, 1.0) <= self.probability:
                    child.variables[i] = self.um_mutation(float(child.variables[i]),
                                                          problem.types[i].min_value,
                                                          problem.types[i].max_value)
                    child.evaluated = False
        
        return child
    
    def um_mutation(self, x, lb, ub):
        return random.uniform(lb, ub) class PCX(Variator):
    
    def __init__(self, nparents = 10, noffspring = 2, eta = 0.1, zeta = 0.1):
        super(PCX, self).__init__(nparents)
        self.nparents = nparents
        self.noffspring = noffspring
        self.eta = eta
        self.zeta = zeta
        
    def evolve(self, parents):
        result = []
        
        for _ in range(self.noffspring):
            index = random.randrange(len(parents))
            parents[index], parents[-1] = parents[-1], parents[index]
            
            result.append(self.pcx(parents))
            
        return result
    
    def pcx(self, parents):
        k = len(parents)
        n = parents[0].problem.nvars
        x = []
        
        for i in range(k):
            x.append(parents[i].variables[:])
            
        g = [sum([x[i][j] for i in range(k)]) / k for j in range(n)]
        D = 0.0
        
        # basis vectors defined by parents
        e_eta = []
        e_eta.append(subtract(x[k-1], g))
        
        for i in range(k-1):
            d = subtract(x[i], g)
            
            if not is_zero(d):
                e = orthogonalize(d, e_eta)
                
                if not is_zero(e):
                    D += magnitude(e)
                    e_eta.append(normalize(e))
        
        D /= k-1
        
        # construct the offspring
        variables = x[k-1]
        variables = add(variables, multiply(random.gauss(0.0, self.zeta), e_eta[0]))
        
        eta = random.gauss(0.0, self.eta)
        
        for i in range(1, len(e_eta)):
            variables = add(variables, multiply(eta*D, e_eta[i]))
            
        result = copy.deepcopy(parents[k-1])
        
        for j in range(n):
            type = result.problem.types[j]
            result.variables[j] = clip(variables[j], type.min_value, type.max_value)
            
        result.evaluated = False
        return result class UNDX(Variator):
    
    def __init__(self, nparents = 10, noffspring = 2, zeta = 0.5, eta = 0.35):
        super(UNDX, self).__init__(nparents)
        self.nparents = nparents
        self.noffspring = noffspring
        self.zeta = zeta
        self.eta = eta
        
    def evolve(self, parents):
        result = []
        
        for _ in range(self.noffspring):
            result.append(self.undx(parents))
            
        return result
    
    def undx(self, parents):
        k = len(parents)
        n = parents[0].problem.nvars
        x = []
        
        for i in range(k):
            x.append(parents[i].variables[:])
            
        g = [sum([x[i][j] for i in range(k)]) / k for j in range(n)]
        
        
        # basis vectors defined by parents
        e_zeta = []
        e_eta = []
        
        for i in range(k-1):
            d = subtract(x[i], g)
            
            if not is_zero(d):
                dbar = magnitude(d)
                e = orthogonalize(d, e_zeta)
                
                if not is_zero(e):
                    e_zeta.append(multiply(dbar, normalize(e)))
        
        D = magnitude(subtract(x[k-1], g))
        
        # create the remaining basis vectors
        for i in range(n-len(e_zeta)):
            d = random_vector(n)
            
            if not is_zero(d):
                e = orthogonalize(d, e_eta)
                
                if not is_zero(e):
                    e_eta.append(multiply(D, normalize(e)))
        
        # construct the offspring
        variables = g
        
        for i in range(len(e_zeta)):
            variables = add(variables, multiply(random.gauss(0.0, self.zeta), e_zeta[i]))
        
        for i in range(1, len(e_eta)):
            variables = add(variables, multiply(random.gauss(0.0, self.eta / math.sqrt(n)), e_eta[i]))
            
        result = copy.deepcopy(parents[k-1])
        
        for j in range(n):
            type = result.problem.types[j]
            result.variables[j] = clip(variables[j], type.min_value, type.max_value)
            
        result.evaluated = False
        return result class SPX(Variator):
    
    def __init__(self, nparents = 10, noffspring = 2, expansion = None):
        super(SPX, self).__init__(nparents)
        self.nparents = nparents
        self.noffspring = noffspring
        
        if expansion is None:
            self.expansion = math.sqrt(nparents+1)
        else:
            self.expansion = expansion
        
    def evolve(self, parents):
        n = len(parents)
        m = parents[0].problem.nvars
        x = []
        
        for i in range(n):
            x.append(parents[i].variables[:])
        
        # compute center of mass
        G = [sum([x[i][j] for i in range(n)]) / n for j in range(m)]
        
        # compute expanded simplex vertices
        for i in range(n):
            x[i] = add(G, multiply(self.expansion, subtract(x[i], G)))
                
        # generate offspring
        result = []
        
        for _ in range(self.noffspring):
            child = copy.deepcopy(parents[n-1])
            r = [math.pow(random.uniform(0.0, 1.0), 1.0 / (i + 1.0)) for i in range(n-1)]
            C = zeros(n, m)
            
            for i in range(n):
                for j in range(m):
                    if i == 0:
                        C[i][j] = 0.0
                    else:
                        C[i][j] = r[i-1] * (x[i-1][j] - x[i][j] + C[i-1][j])
            
            for j in range(m):
                type = child.problem.types[j]
                child.variables[j] = clip(x[n-1][j] + C[n-1][j], type.min_value, type.max_value)
            
            child.evaluated = False
            result.append(child)
            
        return result class Multimethod(Variator):
    
    def __init__(self, algorithm, variators, update_frequency=100):
        super(Multimethod, self).__init__(max([v.arity for v in variators]))
        self.algorithm = algorithm
        self.variators = variators
        self.update_frequency = update_frequency
        self.last_update = 0
        self.probabilities = [1.0 / len(variators) for _ in range(len(variators))]
        
        self.select()
        
    def select(self):
        self.last_update += 1
        
        if self.last_update >= self.update_frequency:
            self.last_update = 0
            counts = [1 for _ in range(len(self.variators))]
            
            if hasattr(self.algorithm, "archive"):
                for solution in self.algorithm.archive:
                    if hasattr(solution, "operator"):
                        counts[solution.operator] += 1
                    
            if hasattr(self.algorithm, "recency_list"):
                for solution in self.algorithm.recency_list:
                    if hasattr(solution, "operator"):
                        counts[solution.operator] += 1
                    
            self.probabilities = [counts[i] / float(sum(counts)) for i in range(len(self.variators))]
        
        self.next_variator = roulette(self.probabilities)
        self.arity = self.variators[self.next_variator].arity
        
    def evolve(self, parents):
        variator = self.variators[self.next_variator]
        result = variator.evolve(parents)
        
        for solution in result:
            solution.operator = self.next_variator
        
        self.select()
        return result class DistanceMatrix(object):
    """Maintains pairwise distances between solutions.
    
    The distance matrix, used by SPEA2, maintains the pairwise distances
    between solutions.  It also provides convenient routines to lookup the
    distance between any two solutions, find the most crowded solution, and
    remove a solution.
    """
    
    def __init__(self, solutions, distance_fun=euclidean_dist):
        super(DistanceMatrix, self).__init__()
        self.distances = []
        
        for i in range(len(solutions)):
            distances_i = []
            for j in range(len(solutions)):
                if i != j:
                    distances_i.append((j, distance_fun(solutions[i], solutions[j])))
                      
            self.distances.append(sorted(distances_i, key=lambda x : x[1]))                
    
    def find_most_crowded(self):
        """Finds the most crowded solution.
        
        Returns the index of the most crowded solution, which is the solution
        with the smallest distance to the nearest neighbor.  Any ties are
        broken by looking at the next distant neighbor.
        """
        minimum_distance = POSITIVE_INFINITY
        minimum_index = -1
        
        for i in range(len(self.distances)):
            distances_i = self.distances[i]
            
            if distances_i[0][1] < minimum_distance:
                minimum_distance = distances_i[0][1]
                minimum_index = i
            elif distances_i[0][1] == minimum_distance:
                for j in range(len(distances_i)):
                    dist1 = distances_i[j][1]
                    dist2 = self.distances[minimum_index][j][1]
                    
                    if dist1 < dist2:
                        minimum_index = i
                        break
                    if dist2 < dist1:
                        break
        
        return minimum_index
    
    def remove_point(self, index):
        """Removes the distance entries for the given solution.
        
        Parameters
        ----------
        index : int
            The index of the solution
        """
        del self.distances[index]
        
        for i in range(len(self.distances)):
            self.distances[i] = [(x if x < index else x-1, y) for (x, y) in self.distances[i] if x != index]
    
    def kth_distance(self, i, k):
        """Returns the distance to the k-th nearest neighbor.
        
        Parameters
        ----------
        i : int
            The index of the solution
        k : int
            Finds the k-th nearest neightbor distance
        """
        return self.distances[i][k][1]
    
    def __getitem__(self, key):
        if isinstance(key, tuple) and len(key) == 2:
            if key[0] == key[1]:
                return 0.0
            else:
                for i, d in self.distances[key[0]]:
                    if i == key[1]:
                        return d
                
                raise ValueError("key not found")
        else:
            raise ValueError("key must be a tuple") def choose(n, k):
    if 0 <= k <= n:
        ntok = 1
        ktok = 1
        for t in range(1, min(k, n - k) + 1):
            ntok *= n
            ktok *= t
            n -= 1
        return ntok / ktok
    else:
        return 0 def point_line_dist(point, line):
    return magnitude(subtract(multiply(float(dot(line, point))/float(dot(line, line)), line), point)) def lsolve(A, b):
    """Gaussian elimination with partial pivoting.
    
    This is implemented here to avoid a dependency on numpy.  This could be
    replaced by :code:`(x, _, _, _) = lstsq(A, b)`, but we prefer the pure
    Python implementation here.
    """
    N = len(b)
     
    for p in range(N):
        # find pivot row and swap
        max = p
         
        for i in range(p+1, N):
            if abs(A[i][p]) > abs(A[max][p]):
                max = i
                 
        A[p], A[max] = A[max], A[p]
        b[p], b[max] = b[max], b[p]
        
        # singular or nearly singular
        if abs(A[p][p]) <= EPSILON:
            raise SingularError("matrix is singular or nearly singular")
        
        # pivot within A and b
        for i in range(p+1, N):
            alpha = A[i][p] / A[p][p]
            b[i] -= alpha * b[p]
            
            for j in range(p, N):
                A[i][j] -= alpha * A[p][j]

    # back substitution
    x = [0.0]*N
    
    for i in range(N-1, -1, -1):
        sum = 0.0
        
        for j in range(i+1, N):
            sum += A[i][j] * x[j]
            
        x[i] = (b[i] - sum) / A[i][i]

    return x def tred2(n, V, d, e):
    """Symmetric Householder reduction to tridiagonal form.
    
    This is derived from the Algol procedures tred2 by Bowdler, Martin,
    Reinsch, and Wilkinson, Handbook for Auto. Comp., Vol.ii-Linear Algebra,
    and the corresponding Fortran subroutine in EISPACK.
    """
    for j in range(n):
        d[j] = V[n-1][j]
        
    for i in range(n-1, 0, -1):
        scale = 0.0
        h = 0.0
        
        for k in range(i):
            scale += abs(d[k])
            
        if scale == 0.0:
            e[i] = d[i-1]
            
            for j in range(i):
                d[j] = V[i-1][j]
                V[i][j] = V[j][i] = 0.0
                
        else:
            for k in range(i):
                d[k] /= scale
                h += d[k]**2
                
            f = d[i-1]
            g = math.sqrt(h)
            
            if f > 0.0:
                g = -g
                
            e[i] = scale*g
            h -= f*g
            d[i-1] = f-g
            
            for j in range(i):
                e[j] = 0.0
                
            for j in range(i):
                f = d[j]
                V[j][i] = f
                g = e[j] + V[j][j]*f
                
                for k in range(j+1, i):
                    g += V[k][j]*d[k]
                    e[k] += V[k][j]*f
                    
                e[j] = g
            
            f = 0.0
            
            for j in range(i):
                e[j] /= h
                f += e[j]*d[j]
                
            hh = f / (2*h)
            
            for j in range(i):
                e[j] -= hh*d[j]
                
            for j in range(i):
                f = d[j]
                g = e[j]
                
                for k in range(j, i):
                    V[k][j] -= f*e[k] + g*d[k]
                    
                d[j] = V[i-1][j]
                V[i][j] = 0.0
                
        d[i] = h
    
    for i in range(n-1):
        V[n-1][i] = V[i][i]
        V[i][i] = 1.0
        h = d[i+1]
        
        if h != 0.0:
            for k in range(i+1):
                d[k] = V[k][i+1] / h
                
            for j in range(i+1):
                g = 0.0
                
                for k in range(i+1):
                    g += V[k][i+1] * V[k][j]
                    
                for k in range(i+1):
                    V[k][j] -= g*d[k]
                    
        for k in range(i+1):
            V[k][i+1] = 0.0
                    
    for j in range(n):
        d[j] = V[n-1][j]
        V[n-1][j] = 0.0
            
    V[n-1][n-1] = 1.0
    e[0] = 0.0 def tql2(n, d, e, V):
    """Symmetric tridiagonal QL algorithm.
    
    This is derived from the Algol procedures tql2, by Bowdler, Martin,
    Reinsch, and Wilkinson, Handbook for Auto. Comp., Vol.ii-Linear Algebra,
    and the corresponding Fortran subroutine in EISPACK.
    """
    for i in range(1, n):
        e[i-1] = e[i]
        
    e[n-1] = 0.0
    
    f = 0.0
    tst1 = 0.0
    eps = math.pow(2.0, -52.0)
    
    for l in range(n):
        tst1 = max(tst1, abs(d[l]) + abs(e[l]))
        m = 1
        
        while m < n:
            if abs(e[m]) <= eps*tst1:
                break
            m += 1
        
        if m > l:
            iter = 0
            
            while True:
                iter += 1
                g = d[l]
                p = (d[l+1] - g) / (2.0*e[l])
                r = hypot(p, 1.0)
                
                if p < 0:
                    r = -r
                    
                d[l] = e[l] / (p + r)
                d[l+1] = e[l] * (p + r)
                dl1 = d[l+1]
                h = g - d[l]
                
                for i in range(l+2, n):
                    d[i] -= h
                    
                f += h
                p = d[m]
                c = 1.0
                c2 = c
                c3 = c
                el1 = e[l+1]
                s = 0.0
                s2 = 0.0
                
                for i in range(m-1, l-1, -1):
                    c3 = c2
                    c2 = c
                    s2 = s
                    g = c*e[i]
                    h = c*p
                    r = hypot(p, e[i])
                    e[i+1] = s*r
                    s = e[i] / r
                    c = p / r
                    p = c*d[i] - s*g
                    d[i+1] = h + s*(c*g + s*d[i])
                    
                    for k in range(n):
                        h = V[k][i+1]
                        V[k][i+1] = s*V[k][i] + c*h
                        V[k][i] = c*V[k][i] - s*h
                    
                p = -s*s2*c3*el1*e[l] / dl1
                e[l] = s*p
                d[l] = c*p
                
                if abs(e[l]) <= eps*tst1:
                    break
            
        d[l] = d[l] + f
        e[l] = 0.0
        
    for i in range(n-1):
        k = i
        p = d[i]
        
        for j in range(i+1, n):
            if d[j] < p:
                k = j
                p = d[j]
        
        if k != i:
            d[k] = d[i]
            d[i] = p
            
            for j in range(n):
                p = V[j][i]
                V[j][i] = V[j][k]
                V[j][k] = p def check_eigensystem(n, C, diag, Q):
    res = 0
    
    for i in range(n):
        for j in range(n):
            cc = 0.0
            dd = 0.0
            
            for k in range(n):
                cc += diag[k] * Q[i][k] * Q[j][k]
                dd += Q[i][k] * Q[j][k]
                
            if abs(cc - C[i if i>j else j][j if i>j else i])/math.sqrt(C[i][i]*C[j][j]) > 1e-10 and abs(cc - C[i if i>j else j][j if i>j else i]) > 1e-9:
                print >> sys.stderr, "imprecise result detected", i, j, cc, C[i if i>j else j][j if i>j else i], (cc-C[i if i>j else j][j if i>j else i])
                res += 1
            
            if abs(dd - (1 if i==j else 0)) > 1e-10:
                print >> sys.stderr, "imprecise result detected (Q not orthog.)", i, j, dd
                res += 1
                
    return res def remove_keys(d, *keys):
    """Returns a new dictionary with the given keys removed.
    
    Parameters
    ----------
    d : dict
        The original dictionary.
    keys : list of keys
        The keys to remove.  If the key is not found in the dictionary, it is
        ignored.
    """
    result = dict(d)
    for key in keys:
        result.pop(key, None)
    return result def only_keys_for(d, func):
    """Returns a new dictionary containing only keys matching function arguments.
    
    Parameters
    ----------
    d : dict
        The original dictionary.
    func: callable
        The function.
    """
    return only_keys(d, *inspect.getargspec(func)[0]) def random_weights(nobjs, population_size):
    """Returns a set of randomly-generated but uniformly distributed weights.
    
    Simply producing N randomly-generated weights does not necessarily produce
    uniformly-distributed weights.  To help produce more uniformly-distributed
    weights, this method picks weights from a large collection of randomly-
    generated weights such that the distances between weights is maximized.
    
    Parameters
    ----------
    nobjs : int
        The number of objectives.
    population_size : int
        The number of weights to generate.
    """
    
    weights = []
    
    if nobjs == 2:
        weights = [[1, 0], [0, 1]]
        weights.extend([(i/(population_size-1.0), 1.0-i/(population_size-1.0)) for i in range(1, population_size-1)])
    else:
        # generate candidate weights
        candidate_weights = []
        
        for i in range(population_size*50):
            random_values = [random.uniform(0.0, 1.0) for _ in range(nobjs)]
            candidate_weights.append([x/sum(random_values) for x in random_values])
        
        # add weights for the corners
        for i in range(nobjs):
            weights.append([0]*i + [1] + [0]*(nobjs-i-1))
            
        # iteratively fill in the remaining weights by finding the candidate
        # weight with the largest distance from the assigned weights
        while len(weights) < population_size:
            max_index = -1
            max_distance = -POSITIVE_INFINITY
            
            for i in range(len(candidate_weights)):
                distance = POSITIVE_INFINITY
                
                for j in range(len(weights)):
                    temp = math.sqrt(sum([math.pow(candidate_weights[i][k]-weights[j][k], 2.0) for k in range(nobjs)]))
                    distance = min(distance, temp)
                    
                if distance > max_distance:
                    max_index = i
                    max_distance = distance
                    
            weights.append(candidate_weights[max_index])
            del candidate_weights[max_index]
            
    return weights def chebyshev(solution, ideal_point, weights, min_weight=0.0001):
    """Chebyshev (Tchebycheff) fitness of a solution with multiple objectives.
    
    This function is designed to only work with minimized objectives.
    
    Parameters
    ----------
    solution : Solution
        The solution.
    ideal_point : list of float
        The ideal point.
    weights : list of float
        The weights.
    min_weight : float
        The minimum weight allowed.
    """
    nobjs = solution.problem.nobjs
    objs = solution.objectives
    return max([max(weights[i], min_weight) * (objs[i]-ideal_point[i]) for i in range(nobjs)]) def normal_boundary_weights(nobjs, divisions_outer, divisions_inner=0):
    """Returns weights generated by the normal boundary method.
    
    The weights produced by this method are uniformly distributed on the
    hyperplane intersecting
    
        [(1, 0, ..., 0), (0, 1, ..., 0), ..., (0, 0, ..., 1)].
        
    Parameters
    ----------
    nobjs : int
        The number of objectives.
    divisions_outer : int
        The number of divisions along the outer set of weights.
    divisions_inner : int (optional)
        The number of divisions along the inner set of weights.
    """
    
    def generate_recursive(weights, weight, left, total, index):
        if index == nobjs - 1:
            weight[index] = float(left) / float(total)
            weights.append(copy.copy(weight))
        else:
            for i in range(left+1):
                weight[index] = float(i) / float(total)
                generate_recursive(weights, weight, left-i, total, index+1)
    
    def generate_weights(divisions):
        weights = []
        generate_recursive(weights, [0.0]*nobjs, divisions, divisions, 0)
        return weights
        
    weights = generate_weights(divisions_outer)
    
    if divisions_inner > 0:
        inner_weights = generate_weights(divisions_inner)
        
        for i in range(len(inner_weights)):
            weight = inner_weights[i]
            
            for j in range(len(weight)):
                weight[j] = (1.0 / nobjs + weight[j]) / 2.0
                
            weights.append(weight)
        
    return weights def default_variator(problem):
    if len(problem.types) == 0:
        raise PlatypusError("problem has no decision variables")
    
    base_type = problem.types[0].__class__
    
    if all([isinstance(t, base_type) for t in problem.types]):
        if base_type in PlatypusConfig.default_variator:
            return PlatypusConfig.default_variator[base_type]
        else:
            for default_type in six.iterkeys(PlatypusConfig.default_variator):
                if issubclass(base_type, default_type):
                    return PlatypusConfig.default_variator[default_type]
                
            raise PlatypusError("no default variator for %s" % base_type)
    else:
        raise PlatypusError("must define variator for mixed types") def default_mutator(problem):
    if len(problem.types) == 0:
        raise PlatypusError("problem has no decision variables")
    
    base_type = problem.types[0].__class__
    
    if all([isinstance(t, base_type) for t in problem.types]):
        if base_type in PlatypusConfig.default_mutator:
            return PlatypusConfig.default_mutator[base_type]
        else:
            for default_type in six.iterkeys(PlatypusConfig.default_mutator):
                if issubclass(base_type, default_type):
                    return PlatypusConfig.default_mutator[default_type]
                
            raise PlatypusError("no default mutator for %s" % base_type)
    else:
        raise PlatypusError("must define mutator for mixed types")
import sys
import copy
import math
import random
import operator
import itertools
import functools
from abc import ABCMeta, abstractmethod
from .core import Algorithm, ParetoDominance, AttributeDominance,\
    AttributeDominance, nondominated_sort, nondominated_prune,\
    nondominated_truncate, nondominated_split, crowding_distance,\
    EPSILON, POSITIVE_INFINITY, Archive, EpsilonDominance, FitnessArchive,\
    Solution, HypervolumeFitnessEvaluator, nondominated_cmp, fitness_key,\
    crowding_distance_key, AdaptiveGridArchive, Selector, EpsilonBoxArchive,\
    PlatypusError, Problem
from .operators import TournamentSelector, RandomGenerator,\
    DifferentialEvolution, clip, UniformMutation, NonUniformMutation,\
    GAOperator, SBX, PM, UM, PCX, UNDX, SPX, Multimethod
from .tools import DistanceMatrix, choose, point_line_dist, lsolve,\
    tred2, tql2, check_eigensystem, remove_keys, only_keys_for
from .weights import random_weights, chebyshev, normal_boundary_weights
from .config import default_variator, default_mutator
    from sets import Set as set
    def eigendecomposition(self):
        self.last_eigenupdate = self.iteration
        
        if self.diagonal_iterations >= self.iteration:
            for i in range(self.problem.nvars):
                self.diag_D[i] = math.sqrt(self.C[i][i])
        else:
            for i in range(self.problem.nvars):
                for j in range(i+1):
                    self.B[i][j] = self.B[j][i] = self.C[i][j]
                    
            offdiag = [0.0]*self.problem.nvars
            tred2(self.problem.nvars, self.B, self.diag_D, offdiag)
            tql2(self.problem.nvars, self.diag_D, offdiag, self.B)
            
            if self.check_consistency:
                check_eigensystem(self.problem.nvars, self.C, self.diag_D, self.B)
                
            for i in range(self.problem.nvars):
                if self.diag_D[i] < 0.0:
                    print("an eigenvalue has become negative", file=sys.stderr)
                    self.diag_D[i] = 0.0
                    
                self.diag_D[i] = math.sqrt(self.diag_D[i])
            
    def test_and_correct_numerics(self):
        # flat fitness, test if function values are identical
        if len(self.population) > 0:
            self.population = sorted(self.population, key=lambda x : x.objectives[0])
            
            if self.population[0].objectives[0] == self.population[min(self.offspring_size-1, self.offspring_size/2 + 1) - 1].objectives[0]:
                print("flat fitness landscape, consider reformulation of fitness, step size increased", file=sys.stderr)
                self.sigma *= math.exp(0.2 + self.cs / self.damps)
                
        # align (renormalize) scale C (and consequently sigma)
        fac = 1.0
        
        if max(self.diag_D) < 1e-6:
            fac = 1.0 / max(self.diag_D)
        elif min(self.diag_D) > 1e4:
            fac = 1.0 / min(self.diag_D)
            
        if fac != 1.0:
            self.sigma /= fac
        
            for i in range(self.problem.nvars):
                self.pc[i] *= fac
                self.diag_D[i] *= fac
                
                for j in range(i+1):
                    self.C[i][j] *= fac**2
            
    def sample(self):
        if (self.iteration - self.last_eigenupdate) > 1.0 / self.ccov / self.problem.nvars / 5.0:
            self.eigendecomposition()
            
        if self.check_consistency:
            self.test_and_correct_numerics()
            
        samples = []
        
        for _ in range(self.offspring_size):
            solution = Solution(self.problem)
            
            if self.diagonal_iterations >= self.iteration:
                while True:
                    feasible = True
                    
                    for j in range(self.problem.nvars):
                        type = self.problem.types[j]
                        value = self.xmean[j] + self.sigma * self.diag_D[j] * random.gauss(0.0, 1.0)
                        if value < type.min_value or value > type.max_value:
                            feasible = False
                            break
                        
                        solution.variables[j] = value
                        
                    if feasible:
                        break
            else:
                artmp = [0.0]*self.problem.nvars
                
                while True:
                    feasible = True
                    
                    for j in range(self.problem.nvars):
                        artmp[j] = self.diag_D[j] * random.gauss(0.0, 1.0)
                        
                    for j in range(self.problem.nvars):
                        type = self.problem.types[j]
                        mutation = 0.0
                        
                        for k in range(self.problem.nvars):
                            mutation += self.B[j][k] * artmp[k]
                            
                        value = self.xmean[j] + self.sigma * mutation
                        
                        if value < type.min_value or value > type.max_value:
                            feasible = False
                            break
                        
                        solution.variables[j] = value
                        
                    if feasible:
                        break
                    
            samples.append(solution)
        
        self.iteration += 1
        return samples
    
    def update_distribution(self):
        xold = self.xmean[:]
        BDz = [0.0]*self.problem.nvars
        artmp = [0.0]*self.problem.nvars
        
        if self.problem.nobjs == 1:
            self.population = sorted(self.population, key=lambda x : x.objectives[0])
        else:
            if self.fitness_evaluator is None:

                self.population = sorted(self.population, key=functools.cmp_to_key(nondominated_cmp)) 