GA.py

# import python packages
import pickle
import copy as copy
from deap import base, creator, gp, tools
from deap import algorithms as algo
from random import random, seed, shuffle, randint, sample, choice
import numpy as numpy
import operator
import math as math
import gc as gc
from collections import defaultdict, deque
from itertools import chain
from operator import attrgetter, itemgetter

# import other pieces of our software
from simulation import paramClass, NP
from utils import findEnd, genRandBits, bitList, genInitValueList

# generates random bitstring with at least one value for each node
def genBits(model):
    # generate random bitlist
    startInd = list(genRandBits(model.size))
    counter = 0
    # make sure bitlist isn't zero
    while numpy.sum(startInd) == 0 and counter < 10000:
        startInd = list(genRandBits(model.size))
        counter += 1
    # go through nodes and make sure that there are 1-5 ones in the random list
    for node in range(0, len(model.nodeList)):
        end = findEnd(node, model)
        start = model.individualParse[node]
        if (end - start) > 1:
            counter = 0
            while numpy.sum(startInd[start:end]) > 5 and counter < 10000:
                chosen = math.floor(random() * (end - start))
                startInd[start + int(chosen)] = 0
                counter += 1
            if numpy.sum(startInd[start:end]) == 0:
                chosen = math.floor(random() * (end - start))
                startInd[start + int(chosen)] = 1
        elif (end - start) == 1:
            startInd[start] = 1
    return [copy.deepcopy(model), startInd]


# finds the lowest error individual in a population
def findPopBest(population):
    saveVal = -1
    minny = 100000
    for i in range(len(population)):
        if numpy.sum(population[i].fitness.values) < minny:
            minny = numpy.sum(population[i].fitness.values)
            saveVal = i
    ultimate = population[saveVal]
    minvals = population[saveVal].fitness.values
    return minvals, ultimate[1], ultimate[0]


# executes two point crossover at node junctions
def cxTwoPointNode(ind1, ind2):
    # copied and modified from deap
    # needed to account for bistring only being one of two components of individual
    """Executes a two-point crossover on the input :term:`sequence`
    individuals. The two individuals are modified in place and both keep
    their original length. 
    :returns: A tuple of two individuals.
    This function uses the :func:`~random.randint` function from the Python 
    base :mod:`random` module.

    Modified to cross over between rules
    """
    size = len(ind1[0].nodeList)
    cxpointer1 = randint(1, size)
    cxpointer2 = randint(1, size - 1)
    # make sure pointers are in right order
    if cxpointer2 >= cxpointer1:
        cxpointer2 += 1
    else:  # Swap the two cx points
        cxpointer1, cxpointer2 = cxpointer2, cxpointer1
    cxpoint1 = ind1[0].individualParse[cxpointer1]
    cxpoint2 = ind1[0].individualParse[cxpointer2]
    # cross over both bitlists and the andNodeLists (as well as andNodeInvertLists)
    ind1[1][cxpoint1:cxpoint2], ind2[1][cxpoint1:cxpoint2] = (
        ind2[1][cxpoint1:cxpoint2],
        ind1[1][cxpoint1:cxpoint2],
    )
    (
        ind1[0].andNodeList[cxpointer1:cxpointer2],
        ind2[0].andNodeList[cxpointer1:cxpointer2],
    ) = (
        ind2[0].andNodeList[cxpointer1:cxpointer2],
        ind1[0].andNodeList[cxpointer1:cxpointer2],
    )
    (
        ind1[0].andNodeInvertList[cxpointer1:cxpointer2],
        ind2[0].andNodeInvertList[cxpointer1:cxpointer2],
    ) = (
        ind2[0].andNodeInvertList[cxpointer1:cxpointer2],
        ind1[0].andNodeInvertList[cxpointer1:cxpointer2],
    )
    # update the arrays seen by C code updateBool
    ind1[0].updateCpointers()
    ind2[0].updateCpointers()
    return ind1, ind2


# sets up GA toolbox from deap
def buildToolbox(individualLength, bitFlipProb, model, params):

    toolbox = base.Toolbox()  # build baseline toolbox
    weightTup = (-1.0,)  # specify weights of the errors
    for i in range(len(model.nodeList) - 1):
        weightTup += (-1.0,)
    creator.create(
        "FitnessMin", base.Fitness, weights=weightTup
    )  # make a fitness minimization function
    creator.create(
        "Individual", list, fitness=creator.FitnessMin
    )  # create a class of individuals that are lists

    # register our bitsring generator and how to create an individual, population
    toolbox.register("genRandomBitString", genBits, model=model)
    toolbox.register(
        "Individual", tools.initIterate, creator.Individual, toolbox.genRandomBitString
    )
    toolbox.register("population", tools.initRepeat, list, toolbox.Individual)
    # create statistics toolbox and give it functions
    stats = tools.Statistics(key=lambda ind: ind.fitness.values)
    stats.register("avg", numpy.mean)
    stats.register("std", numpy.std)
    stats.register("min", numpy.min)
    stats.register("max", numpy.max)

    # finish registering the toolbox functions
    toolbox.register("mate", tools.cxTwoPoint)
    toolbox.register("mutate", tools.mutFlipBit, indpb=bitFlipProb)
    toolbox.register("select", selNSGA2)
    toolbox.register("similar", numpy.array_equal)
    return toolbox, stats


# calculate  fitness for an individual
def evaluateByNode(individual, cells, model, sss, params, KOlist, KIlist, boolC):
    boolValues = [
        NP(
            list(individual),
            model,
            cells,
            model.initValueList[i],
            params,
            KOlist[i],
            KIlist[i],
            boolC,
        )
        for i in range(len(sss))
    ]
    return tuple(
        [
            numpy.sum(
                [
                    (boolValues[j][i] - sss[j][model.nodeList[i]]) ** 2
                    for j in range(0, len(sss))
                ]
            )
            for i in range(0, len(model.nodeList))
        ]
    )


# generates a random set of samples made up of cells by using parameteris from probInit seq
# to set up then iterating using strict Boolean modeling.
def runProbabilityBooleanSims(
    individual, model, sampleNum, cells, params, KOlist, KIlist, boolC
):
    seeds = []
    for i in range(0, sampleNum):
        seeds.append(random())
    samples = [
        sampler(
            individual, model, sampleNum, seeds[i], params, KOlist[i], KIlist[i], boolC
        )
        for i in range(sampleNum)
    ]
    return samples


# generates random seed samples... i.e. generates random starting states then runs EBN
def sampler(individual, model, cells, seeder, params, KOs, KIs, boolC):
    seed(seeder)
    cellArray = []
    sampleProbs = []
    # generate random proportions for each node to start
    for j in range(0, len(model.nodeList)):
        sampleProbs.append(random())
    return NP(individual, model, cells, sampleProbs, params, KOs, KIs, boolC)


# generates list of offspring to be compared... decides to do crossover or mutation
def varOrAdaptive(population, toolbox, model, lambda_, cxpb, mutpb, genfrac, mutModel):
    # algorithm for generating a list of offspring... copied and pasted from DEAP with modification for adaptive mutation
    assert (cxpb + mutpb) <= 1.0, (
        "The sum of the crossover and mutation "
        "probabilities must be smaller or equal to 1.0."
    )
    offspring = []
    for _ in range(lambda_):
        op_choice = random()
        if op_choice < cxpb:  # Apply crossover
            ind1, ind2 = list(map(toolbox.clone, sample(population, 2)))
            ind1, ind2 = cxTwoPointNode(ind1, ind2)
            del ind1.fitness.values
            offspring.append(ind1)
        elif op_choice < cxpb + mutpb:  # Apply mutation
            ind = toolbox.clone(choice(population))
            (ind,) = mutFlipBitAdapt(ind, genfrac, mutModel)
            del ind.fitness.values
            offspring.append(ind)
        else:  # shouldn't happen... clone existing individual
            offspring.append(choice(population))
    return offspring


# select node to mutate
def selectMutNode(errors):
    normerrors = [
        1.0 * error / numpy.sum(errors) for error in errors
    ]  # normalize errors to get a probability that the node  is modified
    probs = numpy.cumsum(normerrors)
    randy = random()  # randomly select a node to mutate
    return next(i for i in range(len(probs)) if probs[i] > randy)


# mutation algorithm
def mutFlipBitAdapt(indyIn, genfrac, mutModel):
    errors = list(indyIn.fitness.values)  # get errors
    individual = indyIn[1]
    model = indyIn[0]
    # get rid of errors in nodes that can't be changed
    errorNodes = 0
    for j in range(len(errors)):
        if model.andLenList[j] < 2:
            errors[j] = 0
        else:
            errorNodes = errorNodes + 1

    if numpy.sum(errors) < 0.05 * errorNodes or errorNodes == 0:
        # condition selection on number of incoming edges + downstream edges
        pseudoerrors = [
            len(model.possibilityList[i])
            if model.successorNums[i] == 0
            else len(model.possibilityList[i]) * model.successorNums[i]
            for i in range(len(model.nodeList))
        ]
        # zero out nodes that can't be changed
        for j in range(len(pseudoerrors)):
            if model.andLenList[j] < 2:
                pseudoerrors[j] = 0
        focusNode = selectMutNode(pseudoerrors)
    else:
        # if errors are relatively high, focus on nodes that fit the worst and have highest in-degree
        # calculate probabilities for mutating each node
        for i in range(len(errors)):
            temper = model.successorNums[i]
            if temper == 0:
                errors[i] = errors[i] * len(model.possibilityList[i])
            else:
                errors[i] = errors[i] * len(model.possibilityList[i]) * temper
        focusNode = selectMutNode(errors)
    # perform mutation
    if model.andLenList[focusNode] > 1:
        # find ends of the node of interest in the individual
        start = model.individualParse[focusNode]
        end = findEnd(focusNode, model)
        # mutate the inputs some of the time
        if len(model.possibilityList[focusNode]) > 3 and random() < mutModel:
            temppermup = []
            upstreamAdders = list(model.possibilityList[focusNode])
            rvals = list(model.rvalues[focusNode])
            while len(temppermup) < 3:
                randy = random()  # randomly select a node to mutate
                tempsum = sum(rvals)
                if tempsum == 0:
                    addNoder = int(math.floor(random() * len(upstreamAdders)))
                else:
                    recalc = numpy.cumsum([1.0 * rval / tempsum for rval in rvals])
                    addNoder = next(i for i in range(len(recalc)) if recalc[i] > randy)
                temppermup.append(upstreamAdders.pop(addNoder))
                rvals.pop(addNoder)
            model.update_upstream(focusNode, temppermup)
            model.updateCpointers()
        for i in range(start, end):
            if random() < 2 / (end - start + 1):
                individual[i] = 1
            else:
                individual[i] = 0
        # ensure that there is at least one shadow and node turned on
        if numpy.sum(individual[start:end]) == 0:
            individual[start] = 1
        indyIn[0] = model
        indyIn[1] = individual
    else:
        print("did not actually check")
    return (indyIn,)


def selNSGA2(individuals, k):
    # NSGA2 selection taken from deap
    """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 = sortNondominatedAdapt(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


# taken from deap and modified slightly to make pareto sorting less strict
def sortNondominatedAdapt(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 = list(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 dominated(fit_i, fit_j):
                dominating_fits[fit_j] += 1
                dominated_fits[fit_i].append(fit_j)
            elif dominated(fit_j, 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


# taken from deap and modified slightly to make pareto sorting less strict
def dominated(ind1, ind2):
    """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
    mean1 = numpy.mean(ind1.wvalues)
    mean2 = numpy.mean(ind2.wvalues)
    std1 = numpy.std(ind1.wvalues)
    if mean1 > mean2:
        not_equal = True
    elif mean1 < mean2:
        return False
    return not_equal


# taken from deap
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 range(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]] += 1.0 * (next[0][i] - prev[0][i]) / norm

    for i, dist in enumerate(distances):
        individuals[i].fitness.crowding_dist = dist


# master GA algorithm
def eaMuPlusLambdaAdaptive(
    toolbox,
    model,
    mu,
    lambda_,
    cxpb,
    mutpb,
    ngen,
    namer,
    newSSS,
    KOlist,
    KIlist,
    params,
    boolC,
    stats=None,
    verbose=__debug__,
):
    # population=[[copy.deepcopy(model),genBits(model)]for i in range(params.popSize)]
    population = toolbox.population(n=params.popSize)
    mutModel = params.mutModel
    logbook = tools.Logbook()
    logbook.header = ["gen", "nevals"] + (stats.fields if stats else [])
    lastcheck = []
    modellist = []
    fitnesslist = []
    popList = []
    # Evaluate the individuals with an invalid fitness
    invalid_ind = [ind for ind in population if not ind.fitness.valid]
    fitnesses = [
        evaluateByNode(
            indy[1], params.cells, indy[0], newSSS, params, KOlist, KIlist, boolC
        )
        for indy in invalid_ind
    ]

    for ind, fit in zip(invalid_ind, fitnesses):
        ind.fitness.values = fit
    fitnesslist.append([list(ind.fitness.values) for ind in population])
    popList.append([list(inder[1]) for inder in population])
    modellist.append(
        [
            [
                (modeler[0].size),
                list(modeler[0].nodeList),
                list(modeler[0].individualParse),
                list(modeler[0].andNodeList),
                list(modeler[0].andNodeInvertList),
                list(modeler[0].andLenList),
                list(modeler[0].nodeList),
                dict(modeler[0].nodeDict),
                list(modeler[0].initValueList),
            ]
            for modeler in 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))

    breaker = True
    for j in range(len(model.andLenList)):
        if model.andLenList[j] > 1:
            breaker = False
    for ind in population:
        if numpy.sum(ind.fitness.values) < 0.01 * len(ind.fitness.values):
            breaker = True
    if breaker:
        # outputList=[fitnesslist, popList, modellist]
        # pickle.dump( outputList, open( namer+"_pops.pickle", "wb" ) )
        return population, logbook

    # Begin the generational process
    for gen in range(1, ngen + 1):
        offspring = varOrAdaptive(
            population,
            toolbox,
            model,
            lambda_,
            0.5 + 0.5 * (1.0 - 1.0 * gen / ngen),
            (0.5 * gen / ngen),
            (1.0 * gen / ngen),
            mutModel,
        )
        # Evaluate the individuals with an invalid fitness
        invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
        fitnesses = [
            evaluateByNode(
                indy[1], params.cells, indy[0], newSSS, params, KOlist, KIlist, boolC
            )
            for indy in invalid_ind
        ]
        for ind, fit in zip(invalid_ind, fitnesses):
            ind.fitness.values = fit
        # Select the next generation population
        population[:] = toolbox.select(population + offspring, mu)
        fitnesslist.append([list(ind.fitness.values) for ind in population])
        popList.append([list(inder[1]) for inder in population])
        modellist.append(
            [
                [
                    (modeler[0].size),
                    list(modeler[0].nodeList),
                    list(modeler[0].individualParse),
                    list(modeler[0].andNodeList),
                    list(modeler[0].andNodeInvertList),
                    list(modeler[0].andLenList),
                    list(modeler[0].nodeList),
                    dict(modeler[0].nodeDict),
                    list(modeler[0].initValueList),
                ]
                for modeler in population
            ]
        )

        # 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))
        breaker = False
        for ind in population:
            if numpy.sum(ind.fitness.values) < 0.01 * len(ind.fitness.values):
                breaker = True
                saveInd = ind
        if breaker:
            errorTemp = saveInd.fitness.values
            for value in errorTemp:
                if value > 0.1:
                    breaker = False
        if breaker:
            # outputList=[fitnesslist, popList, modellist]
            # pickle.dump( outputList, open( namer+"_pops.pickle", "wb" ) )
            return population, logbook

    # outputList=[fitnesslist, popList, modellist]
    # pickle.dump( outputList, open( namer+"_pops.pickle", "wb" ) )
    return population, logbook


# wrapper for GA. eaMuPlusLambdaAdaptive does heavy lifting
def GAsearchModel(model, sampleList, params, KOlist, KIlist, namer, boolC):
    newInitValueList = genInitValueList(sampleList, model)  # set up initial value list
    model.initValueList = newInitValueList  # append initial value list to model
    toolbox, stats = buildToolbox(
        model.size, params.bitFlipProb, model, params
    )  # set up toolbox
    # run GA, find best in population, return
    population, logbook = eaMuPlusLambdaAdaptive(
        toolbox,
        model,
        mu=params.mu,
        lambda_=params.lambd,
        stats=stats,
        cxpb=params.crossoverProb,
        mutpb=params.mutationProb,
        ngen=params.generations,
        namer=namer,
        newSSS=sampleList,
        KOlist=KOlist,
        KIlist=KIlist,
        params=params,
        verbose=params.verbose,
        boolC=boolC,
    )
    out1, out2, model = findPopBest(population)
    return model, out1, out2


# wrapper to do local search in parrallell manner
def localSearch(model, indy, newSSS, params, KOlist, KIlist, boolC):
    outputs = [
        checkNodePossibilities(
            node, indy, newSSS, params.cells, model, params, KOlist, KIlist, boolC
        )
        for node in range(len(model.nodeList))
    ]
    equivs = []
    individual = []
    devs = []
    for output in outputs:
        individual.extend(output[0])
        equivs.append(output[1])
        devs.append(output[2])
    return individual, equivs, devs


# local search function
def checkNodePossibilities(
    node, indy, newSSS, cellNum, model, params, KOlist, KIlist, boolC
):
    tol = 0.01 * len(newSSS)  # set tolerance for equivalence
    end = findEnd(node, model)  # find end of model for this node
    start = model.individualParse[node]  # find start of model for this node
    truth = list(indy[start:end])
    equivs = [truth]  # add if
    if (end - start) == 0:
        return truth, equivs, 0.0
    indOptions = []
    indErrors = []
    # iterate over possibilities for this node
    for i in range(1, 2 ** (end - start)):
        tempultimate = list(indy)
        tempInd = bitList(i, len(truth))
        tempultimate[start:end] = tempInd  # set rule to one being checked
        currentsumtemp = evaluateByNode(
            tempultimate, cellNum, model, newSSS, params, KOlist, KIlist, boolC
        )
        currentsum = currentsumtemp[node]  # save the error found
        indOptions.append(tempInd)
        indErrors.append(currentsum)
        gc.collect()
    minny = min(indErrors)
    equivs = []
    # find the minimum error individual
    for i in range(len(indOptions)):
        if indErrors[i] < minny + tol:
            equivs.append(indOptions[i])
    truth = equivs[0]
    return truth, equivs, minny