<jupyter_start><jupyter_text>Optimizing a scikit model with strange waysWe'll now be trying to maximize some aspect of a specific machine learning model. I've chosen for that to be accuracy.Requires the Diabetes dataset file from [Kaggle](https://www.kaggle.com/uciml/pima-indians-diabetes-database).Author: <NAME><jupyter_code>import random
from collections import defaultdict
import numpy as np
import pandas as pd

import matplotlib.pyplot as plt
import seaborn as sns
sns.set_theme(style="darkgrid")

from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score
data = pd.read_csv('diabetes.csv')  # Just going to use these like globals. YOLO
X = data.drop('Outcome', axis=1).values
y = data[['Outcome']].values.ravel()<jupyter_output><empty_output><jupyter_text>Common stuff and commentsWe've got a few more challenges here, such as discrete variables in the search space, limited range variables in the search space (errors if we'd somehow end up testing zero predictors), and a way more expensive fitness calculation than we're used to (since it is essentially necessary to train and test a model to be able to get it).I'll try to do some weird kind of discretization in the case of continuous variables to reduce search space a bit - and to possibly make memoization of results feasible. So each numeric feature will have their own min, max, and step sizes. Might actually even work okay.The parameters of the RF that we'll use are max_depth, ccp_alpha, min_impurity_decrease, so three dimensions in the search space.PS: Absolutely not as flexible as it might first look - it's still written around just this one model.<jupyter_code># Common to all...
# Ranges takes care of the proper input ranges. If we leave these bounds, we project back.
# Minimum, maximum, 'step size' in decimal points.
ranges = [[2, 20, 0], [0.0, 1.0, 2], [0.0, 0.5, 2]]

# Creates n parameter sets with given ranges for the discussed numbers.
def create_population(n, ranges):
    models = []
    for _ in range(n):
        models.append([gen_param(*dist) for dist in ranges])
    return models

# Generates a random parameter
def gen_param(start, end, decimals):
    rnd = round(random.random() * (end - start) + start, decimals)
    return int(rnd) if decimals == 0 else rnd

# Creates a model based on generated parameters
def create_model(max_depth, ccp_alpha, min_impurity_decrease):
    return RandomForestClassifier(max_depth=max_depth, ccp_alpha=ccp_alpha, min_impurity_decrease=min_impurity_decrease,
                                 random_state=42)

# Updates memory, returns the best new one it found (each method will use same kind of memory)
# This is, of course, an awful name for the function, for its return values aren't really going to be expected
def update_memory(memory, pop):
    best, bestacc = None, 0.0
    for params in pop:
        # If model hasn't been evaluated, evaluate it.
        a,b,c = params
        if memory[a][b][c] == None:
            model = create_model(a,b,c)
            score = cross_val_score(model, X, y, cv=4).mean()
            memory[a][b][c] = score
            if score > bestacc:
                best, bestacc = params, score
    return best, bestacc

# Given some next generation, fix it to be within the limited search space
def fix_generation(pop, ranges):
    for params in pop:
        for i in range(len(params)):
            # Rounds each set of parameters to be within the wanted steps
            params[i] = round(params[i], ranges[i][2])
            # Projects them inside if they are outside the search space
            params[i] = max(params[i], ranges[i][0])  # Raises it to the minimum at least
            params[i] = min(params[i], ranges[i][1])  # Lowers it to the maximum at most

def plot(data, title):
    sns.lineplot(data = data)
    plt.xlabel('Generation')
    plt.ylabel('Best value')
    plt.title(title)
    plt.show()
# Example population of model parameters
create_population(4, ranges)
# Any of these can be fed into create_model as create_model(*elem).
create_model(*create_population(1, ranges)[0])<jupyter_output><empty_output><jupyter_text>Differential Evolution Implementation<jupyter_code># ranges  - ranges for the model parameters
# n       - population size
# scaling - scaling during generating next gen
# loops   - how many loops with no improvement to try
def de(ranges, n=30, scaling=0.5, loops=10):
    # Since some things are discretized now, we may be able to avoid
    # recalculating accuracy by remembering it.
    memory = defaultdict(lambda: defaultdict(lambda: defaultdict(lambda: None)))

    # Initial model set
    pop = create_population(n, ranges)
    best, bestacc = update_memory(memory, pop)
    history = [bestacc]
    
    # The main loop...
    loops_since_improvement = 0
    while loops_since_improvement < loops:
        loops_since_improvement += 1
        nextpop = de_next(pop, ranges, n, scaling)
        nextbest, nextbestacc = update_memory(memory, nextpop)
        if nextbestacc > bestacc:
            loops_since_improvement = 0
            best, bestacc = nextbest, nextbestacc
        history.append(bestacc)
        # overwrite pop with better ones from nextpop
        de_updatepop(pop, nextpop, memory)

    return best, bestacc, history

# ...
def de_next(pop, ranges, n, scaling):
    nextpop = []
    for i in range(n):
        a, b = random.choice(pop), random.choice(pop)  # Picking two random elements
        nextpop.append([pop[i][j]+scaling*(b[j]-a[j]) for j in range(len(a))])  # Add the scaled ab vector to pop[i]...
        
    fix_generation(nextpop, ranges)  # Solve some problems
    return nextpop

# ...
def de_updatepop(pop, nextpop, memory):
    # Now that everything in pop and nextpop are guaranteed to be in memory,
    # this couldn't be easier.
    for i in range(len(pop)):
        a1,b1,c1 = pop[i]
        a2,b2,c2 = nextpop[i]
        if memory[a2][b2][c2] > memory[a1][b1][c1]:
            pop[i] = nextpop[i]<jupyter_output><empty_output><jupyter_text>Sanity checkSeems to do okay and get a considered-good accuracy in the end, however the scaling has to be quite high (even above one for best outcome speed wise) for it to really go to a good place (and to be truly performant - by smashing into the corners of some variable range it's a lot more likely to try to compute something we've already computed). Feels strange to do it like this, I'll test it with several configs later.<jupyter_code>best, bestval, his = de(ranges, 15, 1.0, 10)
(best, bestval)<jupyter_output><empty_output><jupyter_text>Accuracy graph over timeWe find we get the best model VERY early though - in one or two generations it won't move any more. May be because of the chosen parameters, may be the data set.<jupyter_code>plot(his, "Differential Evolution")<jupyter_output><empty_output><jupyter_text>Measuring timeSo we know we spend most of the time for no reason. So let's measure and contrast the time this method takes with waiting ~3 loops for improvement and ~10 like above to see if we check many models we haven't made yet. (Also freeze the seed because otherwise things may get weird)<jupyter_code>random.seed(42)
%timeit -r 1 -n 1 de(ranges, 15, 1.0, 3)
random.seed(42)
%timeit -r 1 -n 1 de(ranges, 15, 1.0, 10)<jupyter_output>1min 21s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)<jupyter_text>Doesn't seem like adding seven extra loops changed runtime so much, suggesting we do get a lot of collisions.<jupyter_code>random.seed(42)
%timeit -r 1 -n 1 de(ranges, 15, 3.0, 3)  # Very high scaling - makes big jumps<jupyter_output>42.6 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)<jupyter_text>Bigger jumps seem to have cut down on processing time some.  Particle Swarm Optimization Implementation<jupyter_code># PSO: still broken as ever.
# ranges  - ranges for the model parameters
# n       - population size
# loops   - how many loops with no improvement to try
# lr      - how much velocity affects future velocity (between 0 and 1)
# c1 and c2   - weighting for personal and overall best when moving  (between 0 and 1)
# res_speed   - the speed we may reset to randomly
# res_speed_p - probability of reset
def pso(ranges, n=15, loops=5, lr=1.0, c1=0.5, c2=0.5, res_speed=5, res_speed_p=0.1):
    # Creating memory in hopes that it helps us catch potential collisions in the discretized search space
    memory = defaultdict(lambda: defaultdict(lambda: defaultdict(lambda: None)))

    # Initial model set
    pop     = create_population(n, ranges)
    popbest = pop[::]
    velocity      = pso_init_velocities(n)
    best, bestacc = update_memory(memory, pop)
    history = [bestacc]
    
    # The main loop...
    loops_since_improvement = 0
    while loops_since_improvement < loops:
        loops_since_improvement += 1
        
        pso_iterate(ranges, velocity, pop, popbest, best, lr, c1, c2, res_speed, res_speed_p)
        nextbest, nextbestacc = update_memory(memory, pop)
        pso_update_personal_best(memory, pop, popbest)
        if nextbestacc > bestacc:
            loops_since_improvement = 0
            best, bestacc = nextbest, nextbestacc
        history.append(bestacc)

    return best, bestacc, history

# ...
def pso_init_velocities(n):
    choices = list(range(-2,3))
    velocities = []
    for _ in range(n):
        velocities.append([random.choice(choices) for _ in range(3)])
    return velocities

# ...
def normalise(vector, length=1):
    l = np.dot(vector, vector)**0.5
    return [v/l for v in vector]

# ...
# Mutates pop and velocity
def pso_iterate(ranges, velocity, pop, popbest, best, lr, c1, c2, res_speed, res_speed_p):
    for i in range(len(pop)):
        z1, z2 = random.random(), random.random()
        velocity[i] = list(np.add(lr * np.array(velocity[i]),
                                  np.add(c1*z1*np.subtract(popbest[i], pop[i]), c2*z2*np.subtract(best, pop[i]))))
        pop[i] = list(np.add(pop[i], velocity[i]))
        # Reset speed for each particle individually this time around.
        if random.random() < res_speed_p:
            velocity[i] = normalise(velocity[i], res_speed)
    # Now we have to fix up any particles that have left the search space
    fix_generation(pop, ranges)

# ...
# Mutates popbest
def pso_update_personal_best(memory, pop, popbest):
    for i in range(len(pop)):
        a1,b1,c1 = pop[i]
        a2,b2,c2 = popbest[i]
        if memory[a1][b1][c1] > memory[a2][b2][c2]:
            popbest[i] = pop[i]<jupyter_output><empty_output><jupyter_text>Sanity checkFinds a pretty good solution as well. What really matters now is the speed.<jupyter_code>best, bestval, his = pso(ranges)
(best, bestval)<jupyter_output><empty_output><jupyter_text>Accuracy graph over timeSame 'issue' of converging very quickly.<jupyter_code>plot(his, "Particle Swarm Optimization")<jupyter_output><empty_output><jupyter_text>Measuring timeDoesn't seem to be too impressive, even given the kinda low pop we start with. If we cut pop, will it converge quicker?<jupyter_code>random.seed(42)
%timeit -r 1 -n 1 pso(ranges)
random.seed(42)
%timeit -r 1 -n 1 print(max(pso(ranges, n=10)[2]))<jupyter_output>0.7695312500000001
42.8 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)<jupyter_text>And we beat DE time with similar end result... What if we go a lot further again?<jupyter_code>random.seed(42)
%timeit -r 1 -n 1 print(max(pso(ranges, n=5)[2]))<jupyter_output>0.76171875
25.4 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)<jupyter_text>It's starting to hurt a tiny bit, but the speed is real. Let's try DE once more, trying to do as little work as possible...<jupyter_code>random.seed(42)
%timeit -r 1 -n 1 print(max(de(ranges, 12, 6.0, 3)[2]))  # Super big jumps and smaller pop<jupyter_output>0.7604166666666667
24.2 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)<jupyter_text>Same time for same result. This definitely makes life hard for those of us that just want to call something the best.  Genetic AlgorithmFinally, it is time for something simple. Implementation<jupyter_code># ranges  - ranges for the model parameters
# n       - population size
# bign    - candidates population size
# loops   - how many loops with no improvement to try
def ga(ranges, n=20, bign=60, loops=5):
    # Creating memory in hopes that it helps us catch potential collisions in the discretized search space
    memory = defaultdict(lambda: defaultdict(lambda: defaultdict(lambda: None)))

    # Initial model set
    pop     = create_population(n, ranges)
    best, bestacc = update_memory(memory, pop)
    history = [bestacc]
    
    # The main loop...
    loops_since_improvement = 0
    while loops_since_improvement < loops:
        loops_since_improvement += 1
        
        bigpop = ga_nextpop(ranges, pop, bign)
        nextbest, nextbestacc = update_memory(memory, bigpop)
        if nextbestacc > bestacc:
            loops_since_improvement = 0
            best, bestacc = nextbest, nextbestacc
        
        # Sort the models based on their accuracy and take only the best ones.
        pop = sorted(pop + bigpop, key=lambda x: memory[x[0]][x[1]][x[2]], reverse=True)[:n]
        
        history.append(bestacc)

    return best, bestacc, history

# ...
def ga_nextpop(ranges, pop, bign):
    nextpop = []
    for _ in range(bign):
        x1, x2 = random.choice(pop), random.choice(pop)
        nextpop.append(mutate(crossover(x1, x2)))
    fix_generation(nextpop, ranges)  # Solve out of range/legal search space issues
    return nextpop

# ...
def mutate(vector):
    if random.random() < 0.25:
        return vector  # just don't bother with some
    newvector = [elem * (random.random() * 1.5 + 0.5) for elem in vector]
    return newvector

# ...
def crossover(x1, x2):
    return [x1[i] if random.random() < 0.5 else x2[i] for i in range(len(x1))]<jupyter_output><empty_output><jupyter_text>Sanity checkIt kind of gets close but still I notice a difference of nearly a percentage point. DE and PSO found the same results, but GA finds something different...<jupyter_code>best, bestval, his = ga(ranges, n=10, bign=15)
(best, bestval)<jupyter_output><empty_output><jupyter_text>Accuracy graph over timeBut yeah, still does the same thing. Gets results fast and stays there.<jupyter_code>plot(his, "Genetic Algorithm")<jupyter_output><empty_output><jupyter_text>Measuring timeLong story short, not better for more time.<jupyter_code>random.seed(42)
%timeit -r 1 -n 1 ga(ranges, n=10, bign=15)<jupyter_output>45.6 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)
{'fraction_non_alphanumeric': 0.09663951817124085, 'fraction_numerical': 0.019368968585312434, 'mean_word_length': 3.6080100914538002, 'pattern_counts': {'":': 0, '<': 70, '<?xml version=': 0, '>': 70, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 2, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '72815', 'license': 'MIT', 'n_tokens_mistral': 4808, 'n_tokens_neox': 4561, 'n_words': 1672, 'path': 'executed/ScikitModels.ipynb', 'repo_name': 'Slopcorn/xeric', 'size': '14618'}