diff --git a/kernel_tuner/strategies/common.py b/kernel_tuner/strategies/common.py index d01eae937..76ad8a568 100644 --- a/kernel_tuner/strategies/common.py +++ b/kernel_tuner/strategies/common.py @@ -3,6 +3,7 @@ from time import perf_counter import numpy as np +from scipy.spatial import distance from kernel_tuner import util from kernel_tuner.searchspace import Searchspace @@ -88,8 +89,17 @@ def __call__(self, x, check_restrictions=True): # else check if this is a legal (non-restricted) configuration if check_restrictions and self.searchspace.restrictions: + legal = self.searchspace.is_param_config_valid(tuple(params)) params_dict = dict(zip(self.searchspace.tune_params.keys(), params)) - legal = util.check_restrictions(self.searchspace.restrictions, params_dict, self.tuning_options.verbose) + + if "constraint_aware" in self.tuning_options.strategy_options and self.tuning_options.strategy_options["constraint_aware"]: + # attempt to repair + new_params = unscale_and_snap_to_nearest_valid(x, params, self.searchspace, self.tuning_options.eps) + if new_params: + params = new_params + legal = True + x_int = ",".join([str(i) for i in params]) + if not legal: result = params_dict result[self.tuning_options.objective] = util.InvalidConfig() @@ -243,3 +253,28 @@ def scale_from_params(params, tune_params, eps): for i, v in enumerate(tune_params.values()): x[i] = 0.5 * eps + v.index(params[i])*eps return x + + + +def unscale_and_snap_to_nearest_valid(x, params, searchspace, eps): + """Helper func to snap to the nearest valid configuration""" + + # params is nearest unscaled point, but is not valid + neighbors = get_neighbors(params, searchspace) + + if neighbors: + # sort on distance to x + neighbors.sort(key=lambda y: distance.euclidean(x,scale_from_params(y, searchspace.tune_params, eps))) + + # return closest valid neighbor + return neighbors[0] + + return [] + + +def get_neighbors(params, searchspace): + for neighbor_method in ["strictly-adjacent", "adjacent", "Hamming"]: + neighbors = searchspace.get_neighbors_no_cache(tuple(params), neighbor_method=neighbor_method) + if len(neighbors) > 0: + return neighbors + return [] diff --git a/kernel_tuner/strategies/firefly_algorithm.py b/kernel_tuner/strategies/firefly_algorithm.py index dc43aae6f..9971df047 100644 --- a/kernel_tuner/strategies/firefly_algorithm.py +++ b/kernel_tuner/strategies/firefly_algorithm.py @@ -13,7 +13,8 @@ maxiter=("Maximum number of iterations", 100), B0=("Maximum attractiveness", 1.0), gamma=("Light absorption coefficient", 1.0), - alpha=("Randomization parameter", 0.2)) + alpha=("Randomization parameter", 0.2), + constraint_aware=("constraint-aware optimization (True/False)", True)) def tune(searchspace: Searchspace, runner, tuning_options): @@ -23,7 +24,7 @@ def tune(searchspace: Searchspace, runner, tuning_options): # using this instead of get_bounds because scaling is used bounds, _, eps = cost_func.get_bounds_x0_eps() - num_particles, maxiter, B0, gamma, alpha = common.get_options(tuning_options.strategy_options, _options) + num_particles, maxiter, B0, gamma, alpha, constraint_aware = common.get_options(tuning_options.strategy_options, _options) best_score_global = sys.float_info.max best_position_global = [] @@ -34,9 +35,10 @@ def tune(searchspace: Searchspace, runner, tuning_options): swarm.append(Firefly(bounds)) # ensure particles start from legal points - population = list(list(p) for p in searchspace.get_random_sample(num_particles)) - for i, particle in enumerate(swarm): - particle.position = scale_from_params(population[i], searchspace.tune_params, eps) + if constraint_aware: + population = list(list(p) for p in searchspace.get_random_sample(num_particles)) + for i, particle in enumerate(swarm): + particle.position = scale_from_params(population[i], searchspace.tune_params, eps) # compute initial intensities for j in range(num_particles): diff --git a/kernel_tuner/strategies/genetic_algorithm.py b/kernel_tuner/strategies/genetic_algorithm.py index c29c150b5..3932baaa1 100644 --- a/kernel_tuner/strategies/genetic_algorithm.py +++ b/kernel_tuner/strategies/genetic_algorithm.py @@ -11,6 +11,7 @@ _options = dict( popsize=("population size", 20), maxiter=("maximum number of generations", 100), + constraint_aware=("constraint-aware optimization (True/False)", True), method=("crossover method to use, choose any from single_point, two_point, uniform, disruptive_uniform", "uniform"), mutation_chance=("chance to mutate is 1 in mutation_chance", 10), ) @@ -19,13 +20,14 @@ def tune(searchspace: Searchspace, runner, tuning_options): options = tuning_options.strategy_options - pop_size, generations, method, mutation_chance = common.get_options(options, _options) - crossover = supported_methods[method] + pop_size, generations, constraint_aware, method, mutation_chance = common.get_options(options, _options) + + GA = GeneticAlgorithm(pop_size, searchspace, constraint_aware, method, mutation_chance) best_score = 1e20 cost_func = CostFunc(searchspace, tuning_options, runner) - population = list(list(p) for p in searchspace.get_random_sample(pop_size)) + population = GA.generate_population() for generation in range(generations): @@ -33,7 +35,8 @@ def tune(searchspace: Searchspace, runner, tuning_options): weighted_population = [] for dna in population: try: - time = cost_func(dna, check_restrictions=False) + # if we are not constraint-aware we should check restrictions upon evaluation + time = cost_func(dna, check_restrictions=not constraint_aware) except util.StopCriterionReached as e: if tuning_options.verbose: print(e) @@ -51,18 +54,19 @@ def tune(searchspace: Searchspace, runner, tuning_options): if tuning_options.verbose: print("Generation %d, best_score %f" % (generation, best_score)) + # build new population for next generation population = [] # crossover and mutate while len(population) < pop_size: - dna1, dna2 = weighted_choice(weighted_population, 2) + dna1, dna2 = GA.weighted_choice(weighted_population, 2) - children = crossover(dna1, dna2) + children = GA.crossover(dna1, dna2) for child in children: - child = mutate(child, mutation_chance, searchspace) + child = GA.mutate(child) - if child not in population and searchspace.is_param_config_valid(tuple(child)): + if child not in population: population.append(child) if len(population) >= pop_size: @@ -75,57 +79,117 @@ def tune(searchspace: Searchspace, runner, tuning_options): tune.__doc__ = common.get_strategy_docstring("Genetic Algorithm", _options) +class GeneticAlgorithm: + + def __init__(self, pop_size, searchspace, constraint_aware=False, method="uniform", mutation_chance=10): + self.pop_size = pop_size + self.searchspace = searchspace + self.tune_params = searchspace.tune_params.copy() + self.constraint_aware = constraint_aware + self.crossover_method = supported_methods[method] + self.mutation_chance = mutation_chance -def weighted_choice(population, n): - """Randomly select n unique individuals from a weighted population, fitness determines probability of being selected.""" - - def random_index_betavariate(pop_size): - # has a higher probability of returning index of item at the head of the list - alpha = 1 - beta = 2.5 - return int(random.betavariate(alpha, beta) * pop_size) - - def random_index_weighted(pop_size): - """Use weights to increase probability of selection.""" - weights = [w for _, w in population] - # invert because lower is better - inverted_weights = [1.0 / w for w in weights] - prefix_sum = np.cumsum(inverted_weights) - total_weight = sum(inverted_weights) - randf = random.random() * total_weight - # return first index of prefix_sum larger than random number - return next(i for i, v in enumerate(prefix_sum) if v > randf) - - random_index = random_index_betavariate - - indices = [random_index(len(population)) for _ in range(n)] - chosen = [] - for ind in indices: - while ind in chosen: - ind = random_index(len(population)) - chosen.append(ind) - - return [population[ind][0] for ind in chosen] - - -def mutate(dna, mutation_chance, searchspace: Searchspace, cache=True): - """Mutate DNA with 1/mutation_chance chance.""" - # this is actually a neighbors problem with Hamming distance, choose randomly from returned searchspace list - if int(random.random() * mutation_chance) == 0: - if cache: - neighbors = searchspace.get_neighbors(tuple(dna), neighbor_method="Hamming") + def generate_population(self): + """ Constraint-aware population creation method """ + if self.constraint_aware: + pop = list(list(p) for p in self.searchspace.get_random_sample(self.pop_size)) else: - neighbors = searchspace.get_neighbors_no_cache(tuple(dna), neighbor_method="Hamming") - if len(neighbors) > 0: - return list(random.choice(neighbors)) - return dna + pop = [] + dna_size = len(self.tune_params) + for _ in range(self.pop_size): + dna = [] + for key in self.tune_params: + dna.append(random.choice(self.tune_params[key])) + pop.append(dna) + return pop + + def crossover(self, dna1, dna2): + """ Apply selected crossover method, repair dna if constraint-aware """ + dna1, dna2 = self.crossover_method(dna1, dna2) + if self.constraint_aware: + return self.repair(dna1), self.repair(dna2) + return dna1, dna2 + + def weighted_choice(self, population, n): + """Randomly select n unique individuals from a weighted population, fitness determines probability of being selected.""" + + def random_index_betavariate(pop_size): + # has a higher probability of returning index of item at the head of the list + alpha = 1 + beta = 2.5 + return int(random.betavariate(alpha, beta) * pop_size) + + def random_index_weighted(pop_size): + """Use weights to increase probability of selection.""" + weights = [w for _, w in population] + # invert because lower is better + inverted_weights = [1.0 / w for w in weights] + prefix_sum = np.cumsum(inverted_weights) + total_weight = sum(inverted_weights) + randf = random.random() * total_weight + # return first index of prefix_sum larger than random number + return next(i for i, v in enumerate(prefix_sum) if v > randf) + + random_index = random_index_betavariate + + indices = [random_index(len(population)) for _ in range(n)] + chosen = [] + for ind in indices: + while ind in chosen: + ind = random_index(len(population)) + chosen.append(ind) + + return [population[ind][0] for ind in chosen] + + + def mutate(self, dna, cache=False): + """Mutate DNA with 1/mutation_chance chance.""" + # this is actually a neighbors problem with Hamming distance, choose randomly from returned searchspace list + if int(random.random() * self.mutation_chance) == 0: + if self.constraint_aware: + if cache: + neighbors = self.searchspace.get_neighbors(tuple(dna), neighbor_method="Hamming") + else: + neighbors = self.searchspace.get_neighbors_no_cache(tuple(dna), neighbor_method="Hamming") + if len(neighbors) > 0: + return list(random.choice(neighbors)) + else: + # select a tunable parameter at random + mutate_index = random.randint(0, len(self.tune_params)-1) + mutate_key = list(self.tune_params.keys())[mutate_index] + # get all possible values for this parameter and remove current value + new_val_options = self.tune_params[mutate_key].copy() + new_val_options.remove(dna[mutate_index]) + # pick new value at random + if len(new_val_options) > 0: + new_val = random.choice(new_val_options) + dna[mutate_index] = new_val + return dna + + + def repair(self, dna): + """ It is possible that crossover methods yield a configuration that is not valid. """ + if not self.searchspace.is_param_config_valid(tuple(dna)): + # dna is not valid, try to repair it + # search for valid configurations neighboring this config + # start from strictly-adjacent to increasingly allowing more neighbors + for neighbor_method in ["strictly-adjacent", "adjacent", "Hamming"]: + neighbors = self.searchspace.get_neighbors_no_cache(tuple(dna), neighbor_method=neighbor_method) + + # if we have found valid neighboring configurations, select one at random + if len(neighbors) > 0: + new_dna = list(random.choice(neighbors)) + print(f"GA crossover resulted in invalid config {dna=}, repaired dna to {new_dna=}") + return new_dna + + return dna def single_point_crossover(dna1, dna2): """Crossover dna1 and dna2 at a random index.""" # check if you can do the crossovers using the neighbor index: check which valid parameter configuration is closest to the crossover, probably best to use "adjacent" as it is least strict? pos = int(random.random() * (len(dna1))) - return (dna1[:pos] + dna2[pos:], dna2[:pos] + dna1[pos:]) + return dna1[:pos] + dna2[pos:], dna2[:pos] + dna1[pos:] def two_point_crossover(dna1, dna2): @@ -137,7 +201,7 @@ def two_point_crossover(dna1, dna2): pos1, pos2 = sorted(random.sample(list(range(start, end)), 2)) child1 = dna1[:pos1] + dna2[pos1:pos2] + dna1[pos2:] child2 = dna2[:pos1] + dna1[pos1:pos2] + dna2[pos2:] - return (child1, child2) + return child1, child2 def uniform_crossover(dna1, dna2): @@ -168,7 +232,7 @@ def disruptive_uniform_crossover(dna1, dna2): child1[ind] = dna2[ind] child2[ind] = dna1[ind] swaps += 1 - return (child1, child2) + return child1, child2 supported_methods = { @@ -177,3 +241,4 @@ def disruptive_uniform_crossover(dna1, dna2): "uniform": uniform_crossover, "disruptive_uniform": disruptive_uniform_crossover, } + diff --git a/kernel_tuner/strategies/greedy_ils.py b/kernel_tuner/strategies/greedy_ils.py index a4c521746..c620ab925 100644 --- a/kernel_tuner/strategies/greedy_ils.py +++ b/kernel_tuner/strategies/greedy_ils.py @@ -1,9 +1,9 @@ """A simple greedy iterative local search algorithm for parameter search.""" +import random from kernel_tuner import util from kernel_tuner.searchspace import Searchspace from kernel_tuner.strategies import common from kernel_tuner.strategies.common import CostFunc -from kernel_tuner.strategies.genetic_algorithm import mutate from kernel_tuner.strategies.hillclimbers import base_hillclimb _options = dict(neighbor=("Method for selecting neighboring nodes, choose from Hamming or adjacent", "Hamming"), @@ -58,9 +58,14 @@ def tune(searchspace: Searchspace, runner, tuning_options): tune.__doc__ = common.get_strategy_docstring("Greedy Iterative Local Search (ILS)", _options) +def mutate(indiv, searchspace: Searchspace): + neighbors = searchspace.get_neighbors_no_cache(tuple(indiv), neighbor_method="Hamming") + return list(random.choice(neighbors)) + + def random_walk(indiv, permutation_size, no_improve, last_improve, searchspace: Searchspace): if last_improve >= no_improve: return searchspace.get_random_sample(1)[0] for _ in range(permutation_size): - indiv = mutate(indiv, 0, searchspace, cache=False) + indiv = mutate(indiv, searchspace) return indiv diff --git a/kernel_tuner/strategies/pso.py b/kernel_tuner/strategies/pso.py index 5b0df1429..efcd63815 100644 --- a/kernel_tuner/strategies/pso.py +++ b/kernel_tuner/strategies/pso.py @@ -13,7 +13,8 @@ maxiter=("Maximum number of iterations", 100), w=("Inertia weight constant", 0.5), c1=("Cognitive constant", 2.0), - c2=("Social constant", 1.0)) + c2=("Social constant", 1.0), + constraint_aware=("constraint-aware optimization (True/False)", False)) def tune(searchspace: Searchspace, runner, tuning_options): @@ -24,7 +25,7 @@ def tune(searchspace: Searchspace, runner, tuning_options): bounds, _, eps = cost_func.get_bounds_x0_eps() - num_particles, maxiter, w, c1, c2 = common.get_options(tuning_options.strategy_options, _options) + num_particles, maxiter, w, c1, c2, constraint_aware = common.get_options(tuning_options.strategy_options, _options) best_score_global = sys.float_info.max best_position_global = [] @@ -35,9 +36,10 @@ def tune(searchspace: Searchspace, runner, tuning_options): swarm.append(Particle(bounds)) # ensure particles start from legal points - population = list(list(p) for p in searchspace.get_random_sample(num_particles)) - for i, particle in enumerate(swarm): - particle.position = scale_from_params(population[i], searchspace.tune_params, eps) + if constraint_aware: + population = list(list(p) for p in searchspace.get_random_sample(num_particles)) + for i, particle in enumerate(swarm): + particle.position = scale_from_params(population[i], searchspace.tune_params, eps) # start optimization for i in range(maxiter): diff --git a/kernel_tuner/strategies/simulated_annealing.py b/kernel_tuner/strategies/simulated_annealing.py index dce929b7b..b380e5efb 100644 --- a/kernel_tuner/strategies/simulated_annealing.py +++ b/kernel_tuner/strategies/simulated_annealing.py @@ -10,16 +10,17 @@ from kernel_tuner.strategies.common import CostFunc _options = dict(T=("Starting temperature", 1.0), - T_min=("End temperature", 0.001), - alpha=("Alpha parameter", 0.995), - maxiter=("Number of iterations within each annealing step", 1)) + T_min=("End temperature", 0.001), + alpha=("Alpha parameter", 0.995), + maxiter=("Number of iterations within each annealing step", 1), + constraint_aware=("constraint-aware optimization (True/False)", True)) def tune(searchspace: Searchspace, runner, tuning_options): # SA works with real parameter values and does not need scaling cost_func = CostFunc(searchspace, tuning_options, runner) # optimization parameters - T, T_min, alpha, niter = common.get_options(tuning_options.strategy_options, _options) + T, T_min, alpha, niter, constraint_aware = common.get_options(tuning_options.strategy_options, _options) T_start = T # compute how many iterations would be needed to complete the annealing schedule @@ -30,7 +31,7 @@ def tune(searchspace: Searchspace, runner, tuning_options): max_feval = tuning_options.strategy_options.get("max_fevals", max_iter) # get random starting point and evaluate cost - pos = list(searchspace.get_random_sample(1)[0]) + pos = generate_starting_point(searchspace, constraint_aware) old_cost = cost_func(pos, check_restrictions=False) # main optimization loop @@ -46,9 +47,9 @@ def tune(searchspace: Searchspace, runner, tuning_options): for _ in range(niter): - new_pos = neighbor(pos, searchspace) + new_pos = neighbor(pos, searchspace, constraint_aware) try: - new_cost = cost_func(new_pos, check_restrictions=False) + new_cost = cost_func(new_pos, check_restrictions=not constraint_aware) except util.StopCriterionReached as e: if tuning_options.verbose: print(e) @@ -73,7 +74,7 @@ def tune(searchspace: Searchspace, runner, tuning_options): stuck = 0 c_old = c if stuck > 100: - pos = list(searchspace.get_random_sample(1)[0]) + pos = generate_starting_point(searchspace, constraint_aware) stuck = 0 # safeguard @@ -103,11 +104,49 @@ def acceptance_prob(old_cost, new_cost, T, tuning_options): return np.exp(((old_cost-new_cost)/old_cost)/T) -def neighbor(pos, searchspace: Searchspace): +def neighbor(pos, searchspace: Searchspace, constraint_aware=True): """Return a random neighbor of pos.""" - # Note: this is not the same as the previous implementation, because it is possible that non-edge parameters remain the same, but suggested configurations will all be within restrictions - neighbors = searchspace.get_neighbors(tuple(pos), neighbor_method='Hamming') if random.random() < 0.2 else searchspace.get_neighbors(tuple(pos), neighbor_method='strictly-adjacent') - if len(neighbors) > 0: - return list(random.choice(neighbors)) - # if there are no neighbors, return a random configuration - return list(searchspace.get_random_sample(1)[0]) + + if constraint_aware: + # Note: this is not the same as the previous implementation, because it is possible that non-edge parameters remain the same, but suggested configurations will all be within restrictions + neighbors = searchspace.get_neighbors(tuple(pos), neighbor_method='Hamming') if random.random() < 0.2 else searchspace.get_neighbors(tuple(pos), neighbor_method='strictly-adjacent') + if len(neighbors) > 0: + return list(random.choice(neighbors)) + # if there are no neighbors, return a random configuration + return list(searchspace.get_random_sample(1)[0]) + + else: + tune_params = searchspace.tune_params + size = len(pos) + pos_out = [] + # random mutation + # expected value is set that values all dimensions attempt to get mutated + for i in range(size): + key = list(tune_params.keys())[i] + values = tune_params[key] + + if random.random() < 0.2: #replace with random value + new_value = random_val(i, tune_params) + else: #adjacent value + ind = values.index(pos[i]) + if random.random() > 0.5: + ind += 1 + else: + ind -= 1 + ind = min(max(ind, 0), len(values)-1) + new_value = values[ind] + + pos_out.append(new_value) + return pos_out + +def random_val(index, tune_params): + """return a random value for a parameter""" + key = list(tune_params.keys())[index] + return random.choice(tune_params[key]) + +def generate_starting_point(searchspace: Searchspace, constraint_aware=True): + if constraint_aware: + return list(searchspace.get_random_sample(1)[0]) + else: + tune_params = searchspace.tune_params + return [random_val(i, tune_params) for i in range(len(tune_params))] diff --git a/test/strategies/test_genetic_algorithm.py b/test/strategies/test_genetic_algorithm.py index cb07f8d7f..d16ad11ce 100644 --- a/test/strategies/test_genetic_algorithm.py +++ b/test/strategies/test_genetic_algorithm.py @@ -14,10 +14,12 @@ def test_weighted_choice(): pop = searchspace.get_random_sample(pop_size) weighted_pop = [[p, i] for i, p in enumerate(pop)] - result = ga.weighted_choice(weighted_pop, 1) + GA = ga.GeneticAlgorithm(pop_size, searchspace) + + result = GA.weighted_choice(weighted_pop, 1) assert result[0] in pop - result = ga.weighted_choice(weighted_pop, 2) + result = GA.weighted_choice(weighted_pop, 2) print(result) assert result[0] in pop assert result[1] in pop @@ -41,9 +43,12 @@ def test_random_population(): def test_mutate(): - pop = searchspace.get_random_sample(1) - mutant = ga.mutate(pop[0], 10, searchspace) + GA = ga.GeneticAlgorithm(1, searchspace) + + pop = GA.generate_population() + + mutant = GA.mutate(pop[0]) assert len(pop[0]) == len(mutant) assert mutant[0] in tune_params["x"] assert mutant[1] in tune_params["y"] diff --git a/test/test_runners.py b/test/test_runners.py index 527c1d252..dd4a7f52b 100644 --- a/test/test_runners.py +++ b/test/test_runners.py @@ -140,6 +140,22 @@ def test_diff_evo(env): assert len(result) > 0 +def test_constraint_aware_GA(env): + options = dict(method="uniform", + constraint_aware=True, + popsize=5, + maxiter=2, + mutation_chance=10, + max_fevals=10) + result, _ = tune_kernel(*env, + strategy="genetic_algorithm", + strategy_options=options, + verbose=True, + cache=cache_filename, + simulation_mode=True) + assert len(result) > 0 + + @skip_if_no_pycuda def test_time_keeping(env): kernel_name, kernel_string, size, args, tune_params = env