cmaes_mixed_integers.py

import itertools

import numpy

from deap import cma

class CovarianceConditionError(Exception):
    pass

class StrategyMixedInteger(cma.Strategy):
    def __init__(self, centroid, sigma, S, **params):
        super(StrategyMixedInteger, self).__init__(centroid, sigma, **params)
        self.S_int = numpy.array(S)
        self.i_I_R = numpy.flatnonzero(2 * self.sigma * numpy.diag(self.C)**0.5 < self.S_int)
        print(self.S_int)
    
    def generate(self, ind_init):
        # print(self.i_I_R)
        n_I_R = self.i_I_R.shape[0]
        lambda_int = None
        
        if n_I_R == 0:
            lambda_int = 0
        elif n_I_R >= self.dim:
            lambda_int = int(numpy.floor(self.lambda_ / 2))
        else:
            lambda_int = int(min(numpy.floor(self.lambda_ / 10) + n_I_R + 1,
                                 numpy.floor(self.lambda_ / 2) - 1))
        
        indices_int = numpy.arange(lambda_int, dtype=int)
        numpy.random.shuffle(indices_int)

        Rp = numpy.zeros((self.lambda_, self.dim))
        Rpp = numpy.zeros((self.lambda_, self.dim))
        
        # Ri' has exactly one of its components set to one.
        # The Ri' are dependent in that the number of mutations for each coordinate
        # differs at most by one
        for i, j in zip(indices_int, itertools.cycle(self.i_I_R)):
            Rp[i, j] = 1
            Rpp[i, j] = numpy.random.geometric(p=0.7**(1.0/n_I_R)) - 1

        I_pm1 = (-1)**numpy.random.randint(0, 2, (self.lambda_, self.dim))
        R_int = I_pm1 * (Rp + Rpp)

        if self.update_count > 0:
            R_int[-1, :] = numpy.floor(-self.S_int - self.last_best) - numpy.floor(-self.S_int - self.centroid)

        y = numpy.random.standard_normal((self.lambda_, self.dim))
        arz = self.centroid + self.sigma * numpy.dot(y, self.BD.T) + self.S_int * R_int

        # The update method requires to remember the y of each individual
        population = list(map(ind_init, arz))
        for ind, yi in zip(population, y):
            ind.__y = yi

        return population
    
    def update(self, population):
        population.sort(key=lambda ind: ind.fitness, reverse=True)

        self.last_best = numpy.array(population[0])

        old_centroid = self.centroid
        self.centroid = numpy.dot(self.weights, population[0:self.mu])

        z = numpy.array([ind.__y for ind in population[:self.mu]])
        zmean = numpy.dot(self.weights, z)

        # Cumulation: update evolution path
        self.ps = (1 - self.cs) * self.ps \
            + numpy.sqrt(self.cs * (2 - self.cs) * self.mueff) \
            * numpy.dot(self.B, zmean)

        self.update_count += 1
        hsig = float(numpy.sum(self.ps**2) / \
                     (1 - (1 - self.cs)**(2 * self.update_count)) / self.dim \
                     < 2 + 4 / (self.dim + 1))

        self.pc = (1 - self.cc) * self.pc \
            + hsig * numpy.sqrt(self.cc * (2 - self.cc) * self.mueff) \
            * numpy.dot(self.BD, zmean)

        if self.i_I_R.shape[0] == 0:
            artmp = population[0:self.mu] - old_centroid
        else:
            artmp = z

        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

        if self.i_I_R.shape[0] == 0:
            self.sigma *= numpy.exp(min(1, (numpy.linalg.norm(self.ps) / self.chiN - 1.))
                                    * self.cs / self.damps)
        elif self.i_I_R.shape[0] < self.dim:
            sig_ix = numpy.flatnonzero(self.sigma * numpy.diag(self.C)**0.5 * numpy.sqrt(1. / self.cs) > 0.2 * self.S_int)
            M = sig_ix.shape[0]
            self.sigma *= numpy.exp(min(1, numpy.linalg.norm(self.ps[sig_ix]) / \
                                           M**0.5 * (1. - 1. / (4. * M) + 1. / (21.*M**2))) \
                                    * 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

        self.i_I_R = numpy.flatnonzero(2 * self.sigma * numpy.diag(self.C)**0.5 < self.S_int)

        if self.cond < 0:
            raise CovarianceConditionError