cxSimulatedBinaryBounded generating complex numbers

[markstrefford]

I am getting the following error in cxSimulatedBinaryBounded in crossover.py (note I've added debug print statements in so the line numbers are different from a vanilla crossover.py):

~\anaconda3\envs\tf_gpu\lib\site-packages\deap\tools\crossover.py in cxSimulatedBinaryBounded(ind1, ind2, eta, low, up)
    363                 print(f'c2={c2}, xl={xl}, xu={xu}')
    364                 c1 = min(max(c1, xl), xu)
--> 365                 c2 = min(max(c2, xl), xu)
    366 
    367                 if random.random() <= 0.5:

TypeError: '>' not supported between instances of 'float' and 'complex'

I have added some logging to this function in order to determine why I am getting this error:

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 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 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.
    """
    print(f'cxSimulatedBinaryBounded(ind1={ind1}, ind2={ind2}, eta={eta}, low={low}, up={up})')

    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(range(size), low, up):
        print(f'**** cxSimulatedBinaryBounded(): i={i}, xl={xl}, xu={xu}')
        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))
                print('c1 = 0.5 * (x1 + x2 + beta_q * (x2 - x1))')
                print(f'c1=0.5*({x1}+{x2}-{beta_q} * ({x2}-{x1}))')
                c1 = 0.5 * (x1 + x2 - beta_q * (x2 - x1))

                beta = 1.0 + (2.0 * (xu - x2) / (x2 - x1))
                print(f'beta: {beta} = 1.0 + (2.0 * ({xu} - {x2}) / ({x2} - {x1}))')
                alpha = 2.0 - beta ** -(eta + 1)
                print(f'alpha: {alpha} = 2.0 - {beta} ** -({eta} + 1)')
                if rand <= 1.0 / alpha:
                    print(f'rand: {rand} <= 1.0 / {alpha}')
                    beta_q = (rand * alpha) ** (1.0 / (eta + 1))
                    print(f'beta_q: {beta_q} = ({rand} * {alpha}) ** (1.0 / ({eta} + 1))')
                else:
                    print(f'rand: {rand} > 1.0 / {alpha}')
                    beta_q = (1.0 / (2.0 - rand * alpha)) ** (1.0 / (eta + 1))
                    print(f'beta_q: {beta_q} = (1.0 / (2.0 - {rand} * {alpha})) ** (1.0 / ({eta} + 1))')
                print('c2 = 0.5 * (x1 + x2 + beta_q * (x2 - x1))')
                print(f'c2=0.5*({x1}+{x2}+{beta_q} * ({x2}-{x1}))')
                c2 = 0.5 * (x1 + x2 + beta_q * (x2 - x1))

                print(f'c1={c1}, xl={xl}, xu={xu}')
                print(f'c2={c2}, xl={xl}, xu={xu}')
                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

This is the output I get leading up to this failure:

cxSimulatedBinaryBounded(ind1=[1.0, 5.05858450116003, 1.56831225826889, 1.7195955330240773, 5.919825525499579, -2.978016152529392, 4.074854693281623, 11.340567036771533, 1.4097374370230575, 4.026871378702991, 3.0167342271558253, 1.8481551193368633, 12.841151455161754, -4.898455396546717, -1.2479643613994247, 5.703764759224595, 2.3737432351348846, 6.359158356821031, 6.6273397619700365, 1.5934576412693926, -3.9175331946264706, -1.441787369014464, -3.1477384758050926, 6.453222131328323, 0.42670885633103184, 4.238552158055388, 6.548399963626881, 2.785132111958214, 12.980953737409225, 14.22683190177186, 3.3809525951590906, 9.624375643313943, 0.9122735716343002, 5.610034076651223, 2.9057233520098418, 2.619944687794042, 3.3495744303609545, 6.955824895924831, 13.512015009080205, 13.079873156141769, 1.811140922837732, 4.103459569553342, 2.602470910052789, 3.0228951015165135, 12.861629638388603, 6.710265963470125, -6.033227235150492, 13.647603800317762, 2.464263735775151, 6.304987938410987, 6.341763029201938, 1.5255721838139178, 5.176686422221088, 1.1014298223435848, -0.24147712004553057, 14.030323527283485, 1.3076549612310093, 7.256979778162414, 7.783386652459925, 2.410031184881055, 2.687832102175344, 0.72381075902286, -2.298153928009963, 9.286598821270683, 1.234935269363381, 5.981876083857351, 5.638084758940002, 2.2013567644409693, 10.61124487129728, -5.674314567941554, -4.5186546625447574, 10.766791956595455, 0.3477648042928736, 7.090191704159231, 7.153486044048032, 3.4747816098061026, -6.474100982417475, 13.060768397214677, 3.894467762098241, 10.325224512260371, 1.5553181773821865, 6.7793700456156625, 4.48945171667532, 2.3261228671667844, -5.338735058307199, -0.4912346388590194, -3.5522023486149528, 12.440145724941644, 2.6425398520734165, 5.938847187799992, 1.2020860970916427, 2.8900668310121818, 11.484558998675588, 9.470113564755831, 0.49187811385003943, 12.054697857864243, 1.2239720618697074, 4.898802637823138, 6.889839771319249, 1.3403833537127348, 14.808866703085513, 8.07346799763457, -3.8143452910951554, 15.797959179802312], 
ind2=[1.0, 4.087081831572911, 2.781200533473287, 1.009046414773266, -1.5972205548584206, -4.85683566690221, 3.0331911575190666, 11.505773171457639, 1.6234905122397167, 4.24102252526828, 7.25952911736916, 3.9584117442468307, 4.121602019027133, -0.04045909494529809, -2.450420146421558, 3.8960041308521287, 1.924575000517069, 6.917157997068191, 7.746477968562644, 2.1233827619148027, 13.555006117416852, 3.592085484910431, 14.13508922533888, 7.907452545191288, 0.24074997788810015, 7.856577505522743, 7.121876051836475, 1.39896941900702, 12.856448075360227, -3.8307537360812987, 15.839052461745705, 8.716703581408488, 0.38204126196782573, 7.914897508105344, 6.843751281681008, 3.2685382445767965, 15.516040360280627, 3.9175877768825487, 9.67961443774417, 9.665408327067418, 1.1378294041013812, 4.918656267151595, 4.224376949840174, 1.4653641196113019, -2.259189847173296, 13.192183929033867, 10.72352340846512, 6.32217314909981, 1.980593057406856, 6.459692828835649, 7.553665512698116, 3.667124025873373, -1.0850083817777314, -4.985954684986354, 6.86374331828554, 11.017969596237421, 2.0160709106166568, 7.488849671145568, 2.713139058528747, 1.1327101901679935, 3.3766632584678646, 5.644595700417138, 10.4736784468139, 14.208508971717599, 0.4204067253447138, 7.973933764972193, 2.6391306199915983, 3.749191099364701, 6.462482550556459, 6.493183329352345, 10.882109780847784, 14.72867742363107, 0.7107121043301619, 7.505033288695229, 7.8998689250758805, 2.5637531627947703, 2.805104645053083, -0.3090367885563774, 6.0261576573677775, 13.550002733456218, 2.205636327972509, 4.700532753050123, 5.751250863003761, 2.6675519979511377, -5.445107818569763, 1.2493336698837139, 7.974611469396999, 12.905773162312975, 1.5733891796782862, 6.785233629236359, 3.857814237807119, 1.8996654833062032, -2.3627229781868895, 5.603385948793392, -1.078215747250833, 12.901378193519605, 2.572623715719062, 4.577406996660714, 4.241176116490276, 1.1728718033966261, 2.7852435020743975, 2.8714269998591533, 5.424815815069428, 14.165257146422382], 
eta=10.0, 
low=[0.0, 4.0, 1.0, 1.0, -6.5, -6.5, -6.5, -6.5, 0.0, 4.0, 1.0, 1.0, -6.5, -6.5, -6.5, -6.5, 0.0, 4.0, 1.0, 1.0, -6.5, -6.5, -6.5, -6.5, 0.0, 4.0, 1.0, 1.0, -6.5, -6.5, -6.5, -6.5, 0.0, 4.0, 1.0, 1.0, -6.5, -6.5, -6.5, -6.5, 0.0, 4.0, 1.0, 1.0, -6.5, -6.5, -6.5, -6.5, 0.0, 4.0, 1.0, 1.0, -6.5, -6.5, -6.5, -6.5, 0.0, 4.0, 1.0, 1.0, -6.5, -6.5, -6.5, -6.5, 0.0, 4.0, 1.0, 1.0, -6.5, -6.5, -6.5, -6.5, 0.0, 4.0, 1.0, 1.0, -6.5, -6.5, -6.5, -6.5, 0.0, 4.0, 1.0, 1.0, -6.5, -6.5, -6.5, -6.5, 0.0, 4.0, 1.0, 1.0, -6.5, -6.5, -6.5, -6.5, 0.0, 4.0, 1.0, 1.0, -6.5, -6.5, -6.5, -6.5], up=[2, 8.0, 8.0, 4.0, 16.0, 16.0, 16.0, 16.0, 2, 8.0, 8.0, 4.0, 16.0, 16.0, 16.0, 16.0, 2, 8.0, 8.0, 4.0, 16.0, 16.0, 16.0, 16.0, 2, 8.0, 8.0, 4.0, 16.0, 16.0, 16.0, 16.0, 2, 8.0, 8.0, 4.0, 16.0, 16.0, 16.0, 16.0, 2, 8.0, 8.0, 4.0, 16.0, 16.0, 16.0, 16.0, 2, 8.0, 8.0, 4.0, 16.0, 16.0, 16.0, 16.0, 2, 8.0, 8.0, 4.0, 16.0, 16.0, 16.0, 16.0, 2, 8.0, 8.0, 4.0, 16.0, 16.0, 16.0, 16.0, 2, 8.0, 8.0, 4.0, 16.0, 16.0, 16.0, 16.0, 2, 8.0, 8.0, 4.0, 16.0, 16.0, 16.0, 16.0, 2, 8.0, 8.0, 4.0, 16.0, 16.0, 16.0, 16.0, 2, 8.0, 8.0, 4.0, 16.0, 16.0, 16.0, 16.0])

**** cxSimulatedBinaryBounded(): i=16, xl=0.0, xu=2
c1 = 0.5 * (x1 + x2 + beta_q * (x2 - x1))
c1=0.5*(1.924575000517069+2.3737432351348846-1.0041416936268404 * (2.3737432351348846-1.924575000517069))
beta: -0.6641570188189823 = 1.0 + (2.0 * (2 - 2.3737432351348846) / (2.3737432351348846 - 1.924575000517069))
alpha: 92.16159106468214 = 2.0 - -0.6641570188189823 ** -(10.0 + 1)
rand: 0.5222232592740141 > 1.0 / 92.16159106468214
beta_q: (0.6772850578932906+0.1988688362687656j) = (1.0 / (2.0 - 0.5222232592740141 * 92.16159106468214)) ** (1.0 / (10.0 + 1))
c2 = 0.5 * (x1 + x2 + beta_q * (x2 - x1))
c2=0.5*(1.924575000517069+2.3737432351348846+(0.6772850578932906+0.1988688362687656j) * (2.3737432351348846-1.924575000517069))
c1=1.923644841909721, xl=0.0, xu=2
c2=(2.301266584719454+0.044662782053670434j), xl=0.0, xu=2

What I can see is this calculation returns a complex number:

beta_q:
>>> (1.0 / (2.0 - 0.5222232592740141 * 92.16159106468214)) ** (1.0 / (10.0 + 1))  
(0.6772850578932906+0.1988688362687656j)

However, if I calculate this manually I don't get the complex number:

>>> (1.0 / (2.0 - 0.5222232592740141 * 92.16159106468214))
-0.021678371395529073
>>> (1.0 / (10.0 + 1))
0.09090909090909091

>>> -0.021678371395529073 ** 0.09090909090909091
-0.7058780799007794

Why is this?

[markstrefford]

It seems the problem is operator precedence:

https://stackoverflow.com/a/66903341/1378071

I've changed the code and testing it, I'll submit a patch if it works.

[markstrefford]

Pull request here #567