上下滑動查看更多源碼

# to maximize the objective
# fitness_weights = 1.0
# to minimize the objective
fitness_weights = -1.0

# copy the original dataframes into local copies, once
X_ga = X.copy()
y_ga = y.copy()

# 'const' (the first column) is not an actual feature, do not include it
X_features = X_ga.columns.to_list()[1:]

try:
    del creator.FitnessMax
    del creator.Individual
except Exception as e:
    pass

creator.create("FitnessMax", base.Fitness, weights=(fitness_weights,))
creator.create(
    "Individual", array.array, typecode='b', fitness=creator.FitnessMax
)

try:
    del toolbox
except Exception as e:
    pass

toolbox = base.Toolbox()
# Attribute generator
toolbox.register("attr_bool", random.randint, 0, 1)
# Structure initializers
toolbox.register(
    "individual",
    tools.initRepeat,
    creator.Individual,
    toolbox.attr_bool,
    len(X_features),
)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)


def evalOneMax(individual):
    # objective function
    # create True/False selector list for features
    # and add True at the start for 'const'
    cols_select = [True] + [i == 1 for i in list(individual)]
    # fit model using the features selected from the individual
    lin_mod = sm.OLS(y_ga, X_ga.loc[:, cols_select], hascnotallow=True).fit()
    return (lin_mod.bic,)


toolbox.register("evaluate", evalOneMax)
toolbox.register("mate", tools.cxTwoPoint)
toolbox.register("mutate", tools.mutFlipBit, indpb=0.05)
toolbox.register("select", tools.selTournament, tournsize=3)

random.seed(0)
pop = toolbox.population(n=300)
hof = tools.HallOfFame(1)
pop, log = algorithms.eaSimple(
    pop, toolbox, cxpb=0.5, mutpb=0.2, ngen=10, halloffame=hof, verbose=True
)

best_individual_ga_small = list(hof[0])
best_features_ga_small = [
    X_features[i] for i, val in enumerate(best_individual_ga_small) if val == 1
]
best_objective_ga_small = (
    sm.OLS(y_ga, X_ga[['const'] + best_features_ga_small], hascnotallow=True)
    .fit()
    .bic
)
print(f'best objective: {best_objective_ga_small}')
print(f'best features:  {best_features_ga_small}')

# to maximize the objective
# fitness_weights = 1.0
# to minimize the objective
fitness_weights = -1.0

# copy the original dataframes into local copies, once
X_ga = X.astype(float).copy()
y_ga = y.astype(float).copy()

# 'const' (the first column) is not an actual feature, do not include it
X_features = X_ga.columns.to_list()[1:]

try:
    del creator.FitnessMax
    del creator.Individual
except Exception as e:
    pass

creator.create("FitnessMax", base.Fitness, weights=(fitness_weights,))
creator.create(
    "Individual", array.array, typecode='b', fitness=creator.FitnessMax
)

try:
    del toolbox
except Exception as e:
    pass

toolbox = base.Toolbox()
# Attribute generator
toolbox.register("attr_bool", random.randint, 0, 1)
# Structure initializers
toolbox.register(
    "individual",
    tools.initRepeat,
    creator.Individual,
    toolbox.attr_bool,
    len(X_features),
)
toolbox.register("population", tools.initRepeat, list, toolbox.individual)


def evalMany(individuals):
    # individuals is an array of shape (n_individuals, n_features)
    # transform it into a list of lists:
    # - a list of individuals
    # - each individual is a list of feature selectors
    ind_list = [list(i) for i in list(individuals)]
    ret = []
    for ind in ind_list:
        # create a list of True/False feature selectors from each individual
        # pre-pend True to always select the 'const' feature
        cols_select = [True] + [i == 1 for i in list(ind)]
        # fit model using the features selected from the individual
        lin_mod = sm.OLS(y_ga, X_ga.loc[:, cols_select], hascnotallow=True).fit()
        ret.append((lin_mod.bic,))
    return ret


# multiprocess pool to evaluate individuals
def joblib_map(f, njobs, *iters):
    return Parallel(n_jobs=njobs)(delayed(f)(*args) for args in zip(*iters))


def selElitistAndTournament(
    individuals, k_tournament, k_elitist=0, tournsize=3
):
    # elitist tournament
    # in addition to the regular tournament, ensure the top #k_elistist individuals are preserved
    return tools.selBest(individuals, k_elitist) + tools.selTournament(
        individuals, k_tournament, tournsize
    )


# Hyperparameters
population_size = 1000
crossover_probability = 0.5
individual_mutation_probability = 0.2
gene_mutation_probability = 0.1
tournament_size = 3
elite_size = 0

n_jobs = os.cpu_count()

# register the pool as the mapper
# and the custom function as the evaluator
toolbox.register("map_multi", joblib_map)
toolbox.register("evaluate", evalMany)

toolbox.register("mate", tools.cxTwoPoint)
toolbox.register("mutate", tools.mutFlipBit, indpb=gene_mutation_probability)
# toolbox.register("select", tools.selTournament, tournsize=tournament_size)
# selection with tournament and optional elitism
toolbox.register(
    "select",
    selElitistAndTournament,
    k_tournament=population_size - elite_size,
    k_elitist=elite_size,
    tournsize=tournament_size,
)

random.seed(0)
population = toolbox.population(n=population_size)
hall_of_fame = tools.HallOfFame(1)
objective_runs_ga = 0

# Evaluate the individuals with an invalid fitness
invalid_ind = [ind for ind in population if not ind.fitness.valid]
# split the population in a list with n_jobs elements
# each element is an array containing multiple individuals
# send individuals to the evaluator
fitnesses_nested = toolbox.map_multi(
    toolbox.evaluate, n_jobs, np.array_split(invalid_ind, n_jobs)
)
objective_runs_ga += len(invalid_ind)
fitnesses = []
for l in fitnesses_nested:
    fitnesses += l
for ind, fit in zip(invalid_ind, fitnesses):
    ind.fitness.values = fit

hall_of_fame.update(population)

n_gen = 1000
best_objective_per_gen_ga = np.full(n_gen, np.nan)
best_objective_ga = np.nan
best_generation_ga = 0
gene_values_mean = np.zeros((n_gen, len(X_features)))
gene_maes = np.full(n_gen, np.nan)
time_to_best_ga = np.inf
# Begin the generational process
iterator = tqdm(range(1, n_gen + 1), desc='generation')
t_start = time.time()
for gen in iterator:
    t_start_loop = time.time()

    # Select the next generation of individuals
    offspring = toolbox.select(population)
    # Vary the pool of individuals via cross-over and mutation
    offspring = algorithms.varAnd(
        offspring,
        toolbox,
        crossover_probability,
        individual_mutation_probability,
    )
    # Evaluate the individuals with an invalid fitness
    invalid_ind = [ind for ind in offspring if not ind.fitness.valid]
    # split the population in a list with n_jobs elements
    # each element is an array containing multiple individuals
    # send list to the evaluator pool
    fitnesses_nested = toolbox.map_multi(
        toolbox.evaluate, n_jobs, np.array_split(invalid_ind, n_jobs)
    )
    objective_runs_ga += len(invalid_ind)
    fitnesses = []
    for l in fitnesses_nested:
        fitnesses += l
    for ind, fit in zip(invalid_ind, fitnesses):
        ind.fitness.values = fit
    # Update the hall of fame with the generated individuals
    hall_of_fame.update(offspring)
    # Replace the current population by the offspring
    population[:] = offspring

    t_end_loop = time.time()

    # record the mean gene values across the population
    gene_values_mean[gen - 1, :] = np.array(population).mean(axis=0)
    if gen >= 2:
        gene_maes[gen - 1] = mean_absolute_error(
            gene_values_mean[gen - 2, :], gene_values_mean[gen - 1, :]
        )
    # pick best individual for stats recording
    best_individual_ga = tools.selBest(population, 1)[0]
    best_objective_per_gen_ga[gen - 1] = best_individual_ga.fitness.values[0]
    if (
        best_objective_ga is np.nan
        or fitness_weights * best_individual_ga.fitness.values[0]
        > fitness_weights * best_objective_ga
    ):
        best_objective_ga = best_individual_ga.fitness.values[0]
        best_generation_ga = gen
        time_to_best_ga = t_end_loop - t_start
        print(
            f'gen: {gen:4n}, curr/prev gene MAE: {gene_maes[gen - 1]:.4f}, new best objective: {best_objective_ga:.4f}, time to best: {time_to_best_ga:.4f}'
        )

    if os.path.isfile('break'):
        # to gracefully break the loop, manually create a file called 'break'
        print(f'Found break file, stopping now.')
        iterator.close()
        break
g_completed_ga = gen

best_individual_ga = list(hall_of_fame[0])
best_features_ga = [
    X_features[i] for i, val in enumerate(best_individual_ga) if val == 1
]
best_objective_ga = (
    sm.OLS(y_ga, X_ga[['const'] + best_features_ga], hascnotallow=True).fit().bic
)  

## 強(qiáng)烈建議關(guān)注#公眾號：數(shù)據(jù)STUDIO 更多優(yōu)質(zhì)內(nèi)容每日更新！

best objective:  33705.569572544795
best generation: 787
objective runs:  600525
time to best:    158.027 sec

上下滑動查看更多源碼

print(f'best objective:  {best_objective_ga}')
print(f'best generation: {best_generation_ga}')
print(f'objective runs:  {objective_runs_ga}')
print(f'time to best:    {time_to_best_ga:.3f} sec')
print(f'best features:   {best_features_ga}')

gvm_df = (
    pd.DataFrame(
        gene_values_mean,
        columns=X_features,
        index=range(1, gene_values_mean.shape[0] + 1),
    )
    .sort_index(axis=1)
    .iloc[:g_completed_ga]
)
if g_completed_ga > 20:
    x_ticks = list(range(1, g_completed_ga + 1, round(g_completed_ga / 20)))
else:
    x_ticks = list(range(1, g_completed_ga + 1))
if x_ticks[-1] != g_completed_ga:
    x_ticks.append(g_completed_ga)

fig, ax = plt.subplots(
    3,
    1,
    sharex=True,
    height_ratios=[24, 1, 1],
    figsize=(12, 48),
    layout='constrained',
)

sns.heatmap(
    gvm_df.sort_index(axis=1).T,
    vmin=0.0,
    vmax=1.0,
    cmap='viridis',
    cbar=True,
    cbar_kws={'fraction': 0.01, 'anchor': (0.0, 1.0)},
    ax=ax[0],
)
ax[0].set_title('Population average of gene values after each generation')
ax[0].axvline(x=best_generation_ga, color='C7')
ax[0].set_xlabel('generation')
ax[0].tick_params(axis='both', reset=True)
ax[0].set_xticks(x_ticks)
ax[0].set_xticklabels(x_ticks)

ax[1].set_xlabel('generation')
ax[1].plot(list(range(2, g_completed_ga + 1)), gene_maes[1:g_completed_ga])
ax[1].vlines(
    x=best_generation_ga,
    ymin=gene_maes[1:g_completed_ga].min(),
    ymax=gene_maes[1:g_completed_ga].max(),
    colors='C7',
)
ax[1].set_xlabel('generation')
ax[1].set_ylabel('MAE')
ax[1].set_title(
    'Mean absolute error of average gene values between current generation and previous generation'
)
ax[1].tick_params(axis='both', reset=True)

ax[2].set_xlabel('generation')
ax[2].plot(
    list(range(1, g_completed_ga + 1)),
    best_objective_per_gen_ga[:g_completed_ga],
)
ax[2].vlines(
    x=best_generation_ga,
    ymin=min(best_objective_per_gen_ga[:g_completed_ga]),
    ymax=max(best_objective_per_gen_ga[:g_completed_ga]),
    colors='C7',
)
ax[2].set_xlabel('generation')
ax[2].set_ylabel('best objective')
ax[2].set_title('Best objective')

fig.suptitle('Genetic algorithm')
fig.savefig('ga-performance.png')
fig.show()

best generation: 787
objective runs:  600525
time to best:    158.027 sec
best features:   ['BedroomAbvGr', 'BldgType_1Fam', 'BsmtFinSF1', 'BsmtFinSF2', 'BsmtUnfSF', 'Condition1_Norm', 'Condition2_PosN', 'ExterQual_Ex', 'Exterior1st_BrkFace', 'Functional_Sev', 'Functional_Typ', 'GarageArea', 'GrLivArea', 'HeatingQC_Ex', 'KitchenQual_Ex', 'LandContour_HLS', 'LotArea', 'LotConfig_CulDSac', 'LowQualFinSF', 'Neighborhood_BrkSide', 'Neighborhood_Crawfor', 'Neighborhood_NWAmes', 'Neighborhood_NoRidge', 'Neighborhood_NridgHt', 'Neighborhood_Somerst', 'Neighborhood_StoneBr', 'OverallCond', 'OverallQual', 'PoolArea', 'RoofMatl_WdShngl', 'SaleCondition_Normal', 'SaleType_New', 'TotalBsmtSF', 'WoodDeckSF', 'YearBuilt']

上下滑動查看更多

SFS:    42.448 sec
GA:     158.027 sec
CMA-ES: 48.326 sec

SFS:    22791
GA:     600525
CMA-ES: 20000

baseline BIC: 34570.1662
SFS:          33708.9860
GA:           33705.5696
CMA-ES:       33703.0705