diff --git a/examples/example_5.py b/examples/example_5.py new file mode 100644 index 00000000..ed517835 --- /dev/null +++ b/examples/example_5.py @@ -0,0 +1,141 @@ +import numpy as np +import sklearn +import scipy.stats as stats + +from sklearn.preprocessing import StandardScaler, MinMaxScaler +from sklearn.model_selection import cross_validate +from sklearn.datasets import fetch_openml + +from matplotlib import pyplot as plt +from sklearn.utils import Bunch + +from suprb import SupRB +from suprb.optimizer.rule.es import ES1xLambda +from suprb.optimizer.solution.nsga2 import NonDominatedSortingGeneticAlgorithm2 +from suprb.optimizer.solution.sampler import BetaSolutionSampler, DiversitySolutionSampler +from suprb.optimizer.solution.spea2 import StrengthParetoEvolutionaryAlgorithm2 +from suprb.optimizer.solution.nsga3 import NonDominatedSortingGeneticAlgorithm3 +from suprb.optimizer.solution.ga import GeneticAlgorithm +from suprb.optimizer.solution.ts import TwoStageSolutionComposition +from suprb.logging.multi_objective import MOLogger + +from utils import log_scores + +import time + +algo_names = { + StrengthParetoEvolutionaryAlgorithm2: "SPEA2", + NonDominatedSortingGeneticAlgorithm2: "NSGA-II", + NonDominatedSortingGeneticAlgorithm3: "NSGA-III", + TwoStageSolutionComposition: "TS", +} + + +def plot_pareto_front(pareto_front: np.ndarray, title: str): + x = pareto_front[:, 0] + y = pareto_front[:, 1] + plt.step(x, y, linestyle="-", marker="x") + plt.title(title) + plt.xlabel("Complexity") + plt.ylabel("Error") + plt.xlim(0, 1) + plt.ylim(0, 1) + + +if __name__ == "__main__": + random_state = 42 + suprb_iter = 2 + sc_iter = 1 + + spea2 = StrengthParetoEvolutionaryAlgorithm2( + n_iter=sc_iter, + population_size=32, + sampler=BetaSolutionSampler(), + early_stopping_delta=0, + early_stopping_patience=10, + ) + nsga2 = NonDominatedSortingGeneticAlgorithm2( + n_iter=sc_iter, + population_size=32, + sampler=BetaSolutionSampler(), + early_stopping_delta=0, + early_stopping_patience=10, + ) + nsga3 = NonDominatedSortingGeneticAlgorithm3( + n_iter=sc_iter, + population_size=32, + sampler=BetaSolutionSampler(), + early_stopping_delta=0, + early_stopping_patience=-1, + ) + nsga3_es = NonDominatedSortingGeneticAlgorithm3( + n_iter=sc_iter, + population_size=32, + sampler=BetaSolutionSampler(), + early_stopping_delta=0, + early_stopping_patience=-1, + ) + ga1 = GeneticAlgorithm(n_iter=sc_iter) + ga2 = GeneticAlgorithm(n_iter=sc_iter * 2) + ts = TwoStageSolutionComposition(algorithm_1=ga1, algorithm_2=nsga3, switch_iteration=suprb_iter, warm_start=False) + sc_algos = (nsga3, nsga2, spea2, ts, nsga3_es) + logger_list = [] + time_list = [] + + plt.rcParams.update( + { + "text.usetex": True, + } + ) + + for i, sc in enumerate(sc_algos): + + data, _ = fetch_openml(name="Concrete_Data", version=1, return_X_y=True) + data = data.to_numpy() + + X, y = data[:, :8], data[:, 8] + X, y = sklearn.utils.shuffle(X, y, random_state=random_state) + warm_start: bool = (True,) + + X = MinMaxScaler(feature_range=(-1, 1)).fit_transform(X) + y = StandardScaler().fit_transform(y.reshape((-1, 1))).reshape((-1,)) + + model = SupRB( + n_iter=suprb_iter, + rule_discovery=ES1xLambda(n_iter=1000, delay=10), + solution_composition=sc, + logger=MOLogger(), + random_state=random_state, + ) + start_time = time.time() + scores = cross_validate( + model, + X, + y, + cv=2, + n_jobs=8, + verbose=10, + scoring=["r2", "neg_mean_squared_error"], + return_estimator=True, + ) + end_time = time.time() + time_list.append(end_time - start_time) + log_scores(scores) + logger_list.append(scores["estimator"][0].logger_) + print("Finished!") + + for t in time_list: + print(f"Time taken: {t}") + axes, plots = plt.subplots() + for l in logger_list: + ##### Plot Pareto Fronts ##### + pareto_front = l.pareto_fronts_ + pareto_front = np.array(pareto_front[suprb_iter - 1]) + hvs = l.metrics_["hypervolume"] + hv = hvs[suprb_iter - 1] + spreads = l.metrics_["spread"] + spread = spreads[suprb_iter - 1] + plot_pareto_front(pareto_front, f"$HV = {hv:.2f}, \Delta = {spread:.2f}$") + plt.show() + + ##### Plot Hypervolume ##### diff --git a/suprb/logging/default.py b/suprb/logging/default.py index 7a30761d..762f642c 100644 --- a/suprb/logging/default.py +++ b/suprb/logging/default.py @@ -36,8 +36,12 @@ def log_iteration(self, X: np.ndarray, y: np.ndarray, estimator: BaseRegressor, def log_metric(key, value): self.log_metric(key=key, value=value, step=estimator.step_) - def log_metric_stats(metric_name: str, attribute_name: str, lst: list): - comprehension = [getattr(e, attribute_name) for e in lst] + def log_metric_stats(metric_name: str, attribute_name: str, lst: list, index=None): + if index is None: + comprehension = [getattr(e, attribute_name) for e in lst] + else: + comprehension = [getattr(e, attribute_name)[index] for e in lst] + log_metric(metric_name + "_min", min(comprehension)) log_metric(metric_name + "_mean", sum(comprehension) / len(comprehension)) log_metric(metric_name + "_max", max(comprehension)) @@ -73,13 +77,26 @@ def log_metric_stats(metric_name: str, attribute_name: str, lst: list): population = estimator.solution_composition_.population_ log_metric("population_size", len(population)) # log_metric("population_diversity", genome_diversity(population)) # When using default logging, not all approaches are compatible with this - log_metric_stats("population_fitness", "fitness_", population) + + # This is a bit of a hack to support multidimensional objective functions ~Felix + if population[0].fitness_.__class__.__name__ == "list": + for i in range(len(population[0].fitness_)): + log_metric_stats(f"population_fitness_o_{i}", "fitness_", population, i) + else: + log_metric_stats("population_fitness", "fitness_", population) + log_metric_stats("population_error", "error_", population) log_metric_stats("population_complexity", "complexity_", population) # Log elitist elitist = estimator.solution_composition_.elitist() - log_metric("elitist_fitness", elitist.fitness_) + + if elitist.fitness_.__class__.__name__ == "list": + for i in range(len(elitist.fitness_)): + log_metric(f"elitist_fitness_o_{i}", elitist.fitness_[i]) + else: + log_metric("elitist_fitness", elitist.fitness_) + log_metric("elitist_error", elitist.error_) log_metric("elitist_complexity", elitist.complexity_) # log_metric("elitist_matched", matched_training_samples(elitist.subpopulation)) # When using default logging, not all approaches are compatible with this diff --git a/suprb/logging/metrics.py b/suprb/logging/metrics.py index 275e1529..74dc3ff4 100644 --- a/suprb/logging/metrics.py +++ b/suprb/logging/metrics.py @@ -27,3 +27,23 @@ def matched_training_samples(pool: list[Rule]): matched = np.stack([rule.match_set_ for rule in pool]).any(axis=0).nonzero()[0].shape[0] total = pool[0].match_set_.shape[0] return matched / total + + +def hypervolume(pareto_front: np.ndarray, reference_point: np.ndarray = None): + """Assumes that the pareto front is sorted by complexity in ascending order.""" + last_x = reference_point[0] + volume = 0 + for fitness in pareto_front: + volume += (last_x - fitness[0]) * np.prod(reference_point[1:] - fitness[1:]) + last_x = fitness[0] + return volume + + +def spread(pareto_front: np.ndarray): + if len(pareto_front) <= 1: + return 0 + """Assumes that the pareto front is sorted by complexity in ascending order.""" + distances = np.linalg.norm(pareto_front[:-1] - pareto_front[1:], axis=1) + avg_distance = np.mean(distances) + delta = np.sum(np.abs(distances - avg_distance)) / ((len(distances)) * avg_distance) + return delta diff --git a/suprb/logging/multi_objective.py b/suprb/logging/multi_objective.py new file mode 100644 index 00000000..6bac5daf --- /dev/null +++ b/suprb/logging/multi_objective.py @@ -0,0 +1,37 @@ +import numpy as np + +from . import DefaultLogger +from ..base import BaseRegressor + +from suprb.solution.fitness import c_norm, pseudo_accuracy + +from .metrics import hypervolume, spread +from ..optimizer.solution.base import MOSolutionComposition, SolutionComposition +from ..optimizer.solution.ts import TwoStageSolutionComposition + + +class MOLogger(DefaultLogger): + + pareto_fronts_: dict + + def log_init(self, X: np.ndarray, y: np.ndarray, estimator: BaseRegressor): + super().log_init(X, y, estimator) + self.pareto_fronts_ = {} + + def log_iteration(self, X: np.ndarray, y: np.ndarray, estimator: BaseRegressor, iteration: int): + super().log_iteration(X, y, estimator, estimator.step_) + sc = estimator.solution_composition_ + current_front = sc.pareto_front() + n = current_front[0].fitness.max_genome_length_ + if hasattr(current_front[0].fitness, "hv_reference"): + reference_points = current_front[0].fitness.hv_reference + else: + reference_points = np.array([1.0, 1.0]) + current_front = [ + [1 - c_norm(solution.complexity_, n), 1 - pseudo_accuracy(solution.error_)] for solution in current_front + ] + self.pareto_fronts_[iteration] = current_front + current_front = np.array(current_front) + self.log_metric("hypervolume", hypervolume(current_front, reference_points), estimator.step_) + self.log_metric("spread", spread(current_front), estimator.step_) + self.log_metric("sc_iterations", sc.step_, estimator.step_) diff --git a/suprb/optimizer/solution/__init__.py b/suprb/optimizer/solution/__init__.py index 8da2bc8e..8f89e07f 100644 --- a/suprb/optimizer/solution/__init__.py +++ b/suprb/optimizer/solution/__init__.py @@ -1,2 +1,3 @@ from .archive import SolutionArchive from .base import SolutionComposition +from .sampler import SolutionSampler diff --git a/suprb/optimizer/solution/base.py b/suprb/optimizer/solution/base.py index 5c9ebc5a..8b5bca81 100644 --- a/suprb/optimizer/solution/base.py +++ b/suprb/optimizer/solution/base.py @@ -7,7 +7,8 @@ from suprb.optimizer import BaseOptimizer from suprb.rule import Rule from suprb.utils import check_random_state -from . import SolutionArchive +from .archive import SolutionArchive +from .sampler import SolutionSampler class SolutionComposition(BaseOptimizer, metaclass=ABCMeta): @@ -145,3 +146,87 @@ def _reset(self): super()._reset() if hasattr(self, "population_"): del self.population_ + + +def hypervolume(pareto_front: list[Solution]): + pareto_front = sorted(pareto_front, key=lambda solution: solution.fitness_[0], reverse=True) + fitness_values = np.array([solution.fitness_ for solution in pareto_front]) + # Needs a MultiObjectiveSolutionFitness + reference_point = pareto_front[0].fitness.hv_reference_ + last_x = reference_point[0] + volume = 0 + for fitness in fitness_values: + volume += (last_x - fitness[0]) * np.prod(reference_point[1:] - fitness[1:]) + last_x = fitness[0] + return volume + + +class MOSolutionComposition(PopulationBasedSolutionComposition, metaclass=ABCMeta): + def __init__( + self, + population_size: int, + n_iter: int, + init: SolutionInit, + archive: SolutionArchive, + sampler: SolutionSampler, + random_state: int, + n_jobs: int, + warm_start: bool, + early_stopping_patience: int = -1, + early_stopping_delta: float = 0, + ): + super().__init__( + population_size=population_size, + n_iter=n_iter, + init=init, + archive=archive, + random_state=random_state, + n_jobs=n_jobs, + warm_start=warm_start, + ) + self.sampler = sampler + self.early_stopping_patience = early_stopping_patience + self.early_stopping_delta = early_stopping_delta + self._best_hypervolume = 0 + self._best_pareto_front = None + self._early_stopping_counter = 0 + self.step_ = 0 + + def check_early_stopping(self): + self.step_ += 1 + if self.early_stopping_patience < 0: + return False + hv = self.hypervolume() + hv_diff = hv - self._best_hypervolume + if self._best_hypervolume < hv: + self._best_hypervolume = hv + self._best_pareto_front = self.pareto_front() + + if hv_diff > self.early_stopping_delta: + self._early_stopping_counter = 0 + else: + self._early_stopping_counter += 1 + if self.early_stopping_patience <= self._early_stopping_counter: + return True + return False + + @abstractmethod + def pareto_front(self) -> list[Solution]: + pass + + def hypervolume(self) -> float: + return hypervolume(self.pareto_front()) + + def elitist(self) -> Optional[Solution]: + """Sample an elitist from the Pareto front""" + pf = self.pareto_front() + if len(pf) == 0: + return None + return self.sampler(pf, random_state=self.random_state_) + + def optimize(self, X: np.ndarray, y: np.ndarray, **kwargs) -> Union[Solution, list[Solution], None]: + self._best_hypervolume = 0 + self._best_pareto_front = None + self._early_stopping_counter = 0 + self.step_ = 0 + super().optimize(X, y, **kwargs) diff --git a/suprb/optimizer/solution/nsga2/__init__.py b/suprb/optimizer/solution/nsga2/__init__.py new file mode 100644 index 00000000..6afc7d01 --- /dev/null +++ b/suprb/optimizer/solution/nsga2/__init__.py @@ -0,0 +1,4 @@ +from .base import NonDominatedSortingGeneticAlgorithm2 +from .crossover import SolutionCrossover +from .mutation import SolutionMutation +from .selection import SolutionSelection diff --git a/suprb/optimizer/solution/nsga2/base.py b/suprb/optimizer/solution/nsga2/base.py new file mode 100644 index 00000000..8aad99c6 --- /dev/null +++ b/suprb/optimizer/solution/nsga2/base.py @@ -0,0 +1,142 @@ +import numpy as np +import scipy.stats as stats + +from suprb import Solution +from suprb.solution.initialization import SolutionInit, RandomInit +from ..base import MOSolutionComposition +from suprb.solution.fitness import NormalizedMOSolutionFitness +from suprb.utils import flatten + + +from .mutation import SolutionMutation, BitFlips +from .selection import SolutionSelection, BinaryTournament +from .crossover import SolutionCrossover, NPoint +from .sorting import fast_non_dominated_sort, calculate_crowding_distances +from ..sampler import SolutionSampler, BetaSolutionSampler + + +class NonDominatedSortingGeneticAlgorithm2(MOSolutionComposition): + """A fast and elitist multiobjective genetic algorithm. + + Implemented as described in 10.1109/4235.996017 + + Parameters + ---------- + n_iter: int + Iterations the metaheuristic will perform. + population_size: int + Number of solutions in the population. + mutation: SolutionMutation + crossover: SolutionCrossover + selection: SolutionSelection + init: SolutionInit + random_state : int, RandomState instance or None, default=None + Pass an int for reproducible results across multiple function calls. + warm_start: bool + If False, solutions are generated new for every `optimize()` call. + If True, solutions are used from previous runs. + n_jobs: int + The number of threads / processes the optimization uses. + """ + + def __init__( + self, + n_iter: int = 32, + population_size: int = 32, + mutation: SolutionMutation = BitFlips(), + crossover: SolutionCrossover = NPoint(n=3), + selection: SolutionSelection = BinaryTournament(), + sampler: SolutionSampler = BetaSolutionSampler(1.5, 1.5), + mutation_rate: float = 0.025, + crossover_rate: float = 0.75, + init: SolutionInit = RandomInit(fitness=NormalizedMOSolutionFitness()), + random_state: int = None, + n_jobs: int = 1, + warm_start: bool = True, + early_stopping_patience: int = -1, + early_stopping_delta: float = 0, + ): + super().__init__( + n_iter=n_iter, + population_size=population_size, + init=init, + archive=None, + random_state=random_state, + sampler=sampler, + n_jobs=n_jobs, + warm_start=warm_start, + early_stopping_patience=early_stopping_patience, + early_stopping_delta=early_stopping_delta, + ) + + self.mutation = mutation + self.crossover = crossover + self.selection = selection + + self.mutation_rate = mutation_rate + self.crossover_rate = crossover_rate + + def _optimize(self, X: np.ndarray, y: np.ndarray): + self.fit_population(X, y) + + for _ in range(self.n_iter): + fitness_values = np.array([solution.fitness_ for solution in self.population_]) + pareto_ranks = fast_non_dominated_sort(fitness_values) + crowding_distances = calculate_crowding_distances(fitness_values, pareto_ranks) + + parents = self.selection( + population=self.population_, + n=self.population_size, + random_state=self.random_state_, + pareto_ranks=pareto_ranks, + crowding_distances=crowding_distances, + ) + + # Note that this expression swallows the last element, if `population_size` is odd + parent_pairs = map(lambda *x: x, *([iter(parents)] * 2)) + + # Crossover + children = list( + flatten( + [ + ( + self.crossover(A, B, random_state=self.random_state_), + self.crossover(B, A, random_state=self.random_state_), + ) + for A, B in parent_pairs + ] + ) + ) + # If `population_size` is odd, we add the solution not selected for reproduction directly + if self.population_size % 2 != 0: + children.append(parents[-1]) + + # Mutation + mutated_children = [self.mutation(child, random_state=self.random_state_).fit(X, y) for child in children] + intermediate_pop = self.population_ + mutated_children + + # Selecting the next generation by the elitist of pareto rank and crowding score + intermediate_fitness_values = np.array([solution.fitness_ for solution in intermediate_pop]) + intermediate_pareto_ranks = fast_non_dominated_sort(intermediate_fitness_values) + intermediate_cds = calculate_crowding_distances(intermediate_fitness_values, intermediate_pareto_ranks) + # Sorting -intermediate_cds to get the sorting permutation that sorts in ascending order + sorting_permutation = np.argsort(-intermediate_cds) + intermediate_pop = [intermediate_pop[index] for index in sorting_permutation] + + # At this point the intermediate population is sorted by cds and can now be sorted with a + # stable sorting algorithm by rank, to get the result we want. + intermediate_pareto_ranks = intermediate_pareto_ranks[sorting_permutation] + sorting_permutation = np.argsort(intermediate_pareto_ranks, kind="stable") + intermediate_pop = [intermediate_pop[index] for index in sorting_permutation] + self.population_ = intermediate_pop[: self.population_size] + + if self.check_early_stopping(): + break + + def pareto_front(self) -> list[Solution]: + if not hasattr(self, "population_") or not self.population_: + return [] + fitness_values = np.array([solution.fitness_ for solution in self.population_]) + pareto_ranks = fast_non_dominated_sort(fitness_values) + pareto_front = np.array(self.population_)[pareto_ranks == 0] + return sorted(pareto_front, key=lambda x: x.fitness_[0], reverse=True) diff --git a/suprb/optimizer/solution/nsga2/crossover.py b/suprb/optimizer/solution/nsga2/crossover.py new file mode 100644 index 00000000..f9ca683c --- /dev/null +++ b/suprb/optimizer/solution/nsga2/crossover.py @@ -0,0 +1,55 @@ +from abc import ABCMeta, abstractmethod + +import numpy as np + +from suprb.base import BaseComponent +from suprb.solution import Solution +from suprb.utils import RandomState + + +class SolutionCrossover(BaseComponent, metaclass=ABCMeta): + + def __init__(self, crossover_rate: float = 0.9): + self.crossover_rate = crossover_rate + + def __call__(self, A: Solution, B: Solution, random_state: RandomState) -> Solution: + if random_state.random() < self.crossover_rate: + return self._crossover(A=A, B=B, random_state=random_state) + else: + # Just return the primary parent + return A + + @abstractmethod + def _crossover(self, A: Solution, B: Solution, random_state: RandomState) -> Solution: + pass + + +class NPoint(SolutionCrossover): + """Cut the genome at N points and alternate the pieces from solution A and B.""" + + def __init__(self, crossover_rate: float = 0.9, n: int = 2): + super().__init__(crossover_rate=crossover_rate) + self.n = n + + @staticmethod + def _single_point(A: Solution, B: Solution, index: int) -> Solution: + return A.clone(genome=np.append(A.genome[:index], B.genome[index:])) + + def _crossover(self, A: Solution, B: Solution, random_state: RandomState) -> Solution: + indices = random_state.choice(np.arange(len(A.genome)), size=min(self.n, len(A.genome)), replace=False) + for index in indices: + A = self._single_point(A, B, index) + B = self._single_point(B, A, index) + return A + + +class Uniform(SolutionCrossover): + """Decide for every bit with uniform probability if the bit in genome A or B is used.""" + + def _crossover(self, A: Solution, B: Solution, random_state: RandomState) -> Solution: + indices = random_state.random(size=len(A.genome)) <= 0.5 + genome = np.empty(A.genome.shape) + genome[indices] = A.genome[indices] + genome[~indices] = B.genome[~indices] + + return A.clone(genome=genome) diff --git a/suprb/optimizer/solution/nsga2/mutation.py b/suprb/optimizer/solution/nsga2/mutation.py new file mode 100644 index 00000000..5f1e1e61 --- /dev/null +++ b/suprb/optimizer/solution/nsga2/mutation.py @@ -0,0 +1,25 @@ +from abc import ABCMeta + +import numpy as np + +from suprb.base import BaseComponent +from suprb.solution import Solution +from suprb.utils import RandomState + + +class SolutionMutation(BaseComponent, metaclass=ABCMeta): + + def __init__(self, mutation_rate: float = 0.05): + self.mutation_rate = mutation_rate + + def __call__(self, solution: Solution, random_state: RandomState) -> Solution: + pass + + +class BitFlips(SolutionMutation): + """Flips every bit in the genome with probability `mutation_rate`.""" + + def __call__(self, solution: Solution, random_state: RandomState) -> Solution: + bit_flips = random_state.random(solution.genome.shape) < self.mutation_rate + genome = np.logical_xor(solution.genome, bit_flips) + return solution.clone(genome=genome) diff --git a/suprb/optimizer/solution/nsga2/selection.py b/suprb/optimizer/solution/nsga2/selection.py new file mode 100644 index 00000000..cf78e1fd --- /dev/null +++ b/suprb/optimizer/solution/nsga2/selection.py @@ -0,0 +1,50 @@ +from abc import ABCMeta + +import numpy as np + +from suprb.base import BaseComponent +from suprb.solution import Solution +from suprb.utils import RandomState + + +class SolutionSelection(BaseComponent, metaclass=ABCMeta): + + def __call__( + self, + population: list[Solution], + n: int, + random_state: RandomState, + pareto_ranks: np.ndarray, + crowding_distances: np.ndarray, + ) -> list[Solution]: + pass + + +class BinaryTournament(SolutionSelection): + """Draw 2 solutions n_parents times and select the best solution from each pair.""" + + def __init__(self): + pass + + def __call__( + self, + population: list[Solution], + n: int, + random_state: RandomState, + pareto_ranks: np.ndarray, + crowding_distances: np.ndarray, + ) -> list[Solution]: + + selection = [] + pairs = random_state.integers(low=0, high=len(population), size=(n, 2)) + for a, b in pairs: + if pareto_ranks[a] == pareto_ranks[b]: + if crowding_distances[a] >= crowding_distances[b]: + selection.append(population[a]) + else: + selection.append(population[b]) + elif pareto_ranks[a] > pareto_ranks[b]: + selection.append(population[a]) + else: + selection.append(population[b]) + return selection diff --git a/suprb/optimizer/solution/nsga2/sorting.py b/suprb/optimizer/solution/nsga2/sorting.py new file mode 100644 index 00000000..3d82ad5b --- /dev/null +++ b/suprb/optimizer/solution/nsga2/sorting.py @@ -0,0 +1,89 @@ +from typing import Tuple + +import numpy as np + + +def fast_non_dominated_sort(fitness_values: np.ndarray) -> np.ndarray: + """Sorts the fitness values into multiple levels of non domination. + + Parameters + ---------- + fitness_values: np.ndarray + numpy of shape (solution_count, objective_count) with all fitness values for every solution + + Returns + ------- + pareto ranks: level of non-dominated front for each solution (non sorted) + + """ + solution_count, objective_count = fitness_values.shape + pareto_ranks = np.ones(solution_count, dtype=np.int32) * -1 + dominated_count = np.zeros(solution_count, dtype=np.int32) + dominates_indices = [[] for _ in range(solution_count)] + + # Calculating the amount of solutions, that dominate i and + # the solutions i is dominating + for i in range(solution_count): + fitness_i = fitness_values[i][None, :] + + less_equal_count = np.sum(fitness_values <= fitness_i, axis=1) + less_than_count = np.sum(fitness_values < fitness_i, axis=1) + dominates_i_mask = (less_equal_count == objective_count) & (less_than_count >= 1) + dominated_count[i] = np.sum(dominates_i_mask) + + greater_equal_count = np.sum(fitness_values >= fitness_i, axis=1) + greater_than_count = np.sum(fitness_values > fitness_i, axis=1) + dominated_by_i_mask = (greater_equal_count == objective_count) & (greater_than_count >= 1) + dominates_indices[i] = np.argwhere(dominated_by_i_mask).flatten().tolist() + + pareto_ranks[dominated_count == 0] = 0 + front_rank = 0 + + # Actual non dominated sort + while -1 in pareto_ranks: + current_front = np.argwhere(pareto_ranks == front_rank).flatten() + front_rank += 1 + for solution in current_front: + solution_mask = np.zeros(solution_count, dtype=bool) + solution_mask[dominates_indices[solution]] = 1 + solution_mask = solution_mask & (dominated_count == 1) + dominated_count[dominates_indices[solution]] -= 1 + pareto_ranks[solution_mask] = front_rank + + return pareto_ranks + + +def calculate_crowding_distances(fitness_values: np.ndarray, pareto_ranks: np.ndarray) -> np.ndarray: + solution_count, objective_count = fitness_values.shape + crowding_distances = np.zeros(solution_count) + + for rank in np.unique(pareto_ranks): + front_indices = np.where(pareto_ranks == rank)[0] + front_size = len(front_indices) + + if front_size <= 2: + crowding_distances[front_indices] = np.inf + continue + + front_fitness = fitness_values[front_indices] + + for m in range(objective_count): + sorting_permutation = np.argsort(front_fitness[:, m]) + sorted_front = front_indices[sorting_permutation] + + crowding_distances[sorted_front[0]] = np.inf + crowding_distances[sorted_front[-1]] = np.inf + + min_f = front_fitness[sorting_permutation[0], m] + max_f = front_fitness[sorting_permutation[-1], m] + + # if max_f == min_f the crowding distance parts that result from objective_m are all 0 as all + # solution share the same coordinate in this dimension of the fitness function + if max_f > min_f: + normalized_range = max_f - min_f + for i in range(1, front_size - 1): + crowding_distances[sorted_front[i]] += ( + front_fitness[sorting_permutation[i + 1], m] - front_fitness[sorting_permutation[i - 1], m] + ) / normalized_range + + return crowding_distances diff --git a/suprb/optimizer/solution/nsga3/__init__.py b/suprb/optimizer/solution/nsga3/__init__.py new file mode 100644 index 00000000..39164975 --- /dev/null +++ b/suprb/optimizer/solution/nsga3/__init__.py @@ -0,0 +1,4 @@ +from .base import NonDominatedSortingGeneticAlgorithm3 +from .crossover import SolutionCrossover +from .mutation import SolutionMutation +from .selection import SolutionSelection diff --git a/suprb/optimizer/solution/nsga3/base.py b/suprb/optimizer/solution/nsga3/base.py new file mode 100644 index 00000000..2f05d5ee --- /dev/null +++ b/suprb/optimizer/solution/nsga3/base.py @@ -0,0 +1,202 @@ +import numpy as np +import scipy.stats as stats + +from suprb import Solution +from suprb.solution.initialization import SolutionInit, RandomInit +from ..base import MOSolutionComposition +from suprb.solution.fitness import NormalizedMOSolutionFitness +from suprb.utils import flatten + + +from .mutation import SolutionMutation, BitFlips +from .selection import SolutionSelection, ReferenceBasedBinaryTournament +from .crossover import SolutionCrossover, NPoint +from .sorting import fast_non_dominated_sort +from .reference import das_dennis_points, calc_ref_direction_distances +from .normalise import NSGAIIINormaliser, HyperPlaneNormaliser +from ..sampler import SolutionSampler, BetaSolutionSampler + + +class NonDominatedSortingGeneticAlgorithm3(MOSolutionComposition): + """A fast and elitist multiobjective genetic algorithm. + + Implemented as described in 10.1007/978-3-319-15892-1_3 + + Parameters + ---------- + n_iter: int + Iterations the metaheuristic will perform. + population_size: int + Number of solutions in the population. + mutation: SolutionMutation + crossover: SolutionCrossover + selection: SolutionSelection + init: SolutionInit + random_state : int, RandomState instance or None, default=None + Pass an int for reproducible results across multiple function calls. + warm_start: bool + If False, solutions are generated new for every `optimize()` call. + If True, solutions are used from previous runs. + n_jobs: int + The number of threads / processes the optimization uses. + """ + + def __init__( + self, + n_iter: int = 32, + n_reference_points: int = 15, + n_dimensions: int = 2, + population_size: int = 32, + normaliser: NSGAIIINormaliser = None, + mutation: SolutionMutation = BitFlips(), + crossover: SolutionCrossover = NPoint(n=3), + selection: SolutionSelection = ReferenceBasedBinaryTournament(), + sampler: SolutionSampler = BetaSolutionSampler(1.5, 1.5), + mutation_rate: float = 0.025, + crossover_rate: float = 0.75, + init: SolutionInit = RandomInit(fitness=NormalizedMOSolutionFitness()), + random_state: int = None, + n_jobs: int = 1, + warm_start: bool = True, + early_stopping_patience: int = -1, + early_stopping_delta: float = 0, + ): + super().__init__( + n_iter=n_iter, + population_size=population_size, + init=init, + archive=None, + sampler=sampler, + random_state=random_state, + n_jobs=n_jobs, + warm_start=warm_start, + early_stopping_patience=early_stopping_patience, + early_stopping_delta=early_stopping_delta, + ) + + self.n_dimensions = n_dimensions + self.n_reference_points = n_reference_points + self.mutation = mutation + self.crossover = crossover + self.selection = selection + self.sampler = sampler + fitness = self.init.fitness + if normaliser is None: + normaliser = HyperPlaneNormaliser(len(fitness.objective_func_)) + self.normaliser = normaliser + self.mutation_rate = mutation_rate + self.crossover_rate = crossover_rate + + def _optimize(self, X: np.ndarray, y: np.ndarray): + self.normaliser.reset() + self.fit_population(X, y) + fitness_values = np.array([solution.fitness_ for solution in self.population_]) + pareto_ranks = fast_non_dominated_sort(fitness_values) + normalised_fitness = self.normaliser(fitness_values, pareto_ranks) + + for i in range(len(self.population_)): + self.population_[i].internal_fitness_ = normalised_fitness[i] + reference_points = das_dennis_points(self.n_reference_points, 2) + + for _ in range(self.n_iter): + fitness_values = np.array([solution.internal_fitness_ for solution in self.population_]) + pareto_ranks = fast_non_dominated_sort(fitness_values) + ref_direction_distance, closest_ref_direction = calc_ref_direction_distances( + fitness_values, reference_points + ) + + parents = self.selection( + population=self.population_, + n=self.population_size, + random_state=self.random_state_, + pareto_ranks=pareto_ranks, + closest_ref_direction=closest_ref_direction, + ref_direction_distance=ref_direction_distance, + ) + + # Note that this expression swallows the last element, if `population_size` is odd + parent_pairs = map(lambda *x: x, *([iter(parents)] * 2)) + + # Crossover + children = list( + flatten( + [ + ( + self.crossover(A, B, random_state=self.random_state_), + self.crossover(B, A, random_state=self.random_state_), + ) + for A, B in parent_pairs + ] + ) + ) + # If `population_size` is odd, we add the solution not selected for reproduction directly + if self.population_size % 2 != 0: + children.append(parents[-1]) + + # Mutation + mutated_children = [self.mutation(child, random_state=self.random_state_).fit(X, y) for child in children] + union_pop = self.population_ + mutated_children + union_fitness_values = np.array([solution.fitness_ for solution in union_pop]) + union_pareto_ranks = fast_non_dominated_sort(union_fitness_values) + # I think we need to normalise the population in every iteration outside of the following if clause + # As the parent selection depends on objective values in U-NSGA-III + union_fitness_values = self.normaliser(union_fitness_values, union_pareto_ranks) + for l in range(len(union_pop)): + union_pop[l].internal_fitness_ = union_fitness_values[l] + union_pop = np.array(union_pop) + tmp_pop = [] + l = 0 + while len(tmp_pop) + np.sum(union_pareto_ranks == l) < self.population_size: + tmp_pop += union_pop[union_pareto_ranks == l].tolist() + l += 1 + next_pop = tmp_pop.copy() + tmp_pop += union_pop[union_pareto_ranks == l].tolist() + + if len(tmp_pop) == self.population_size: + self.population_ = tmp_pop + # Here we have to check for early stopping as the check at the end of the loop might not be triggered + if self.check_early_stopping(): + break + continue + + solutions_left_count = self.population_size - len(next_pop) + + union_ref_dist, union_closest_ref = calc_ref_direction_distances(union_fitness_values, reference_points) + niche_count = np.zeros(reference_points.shape[0]) + values, counts = np.unique(union_closest_ref[union_pareto_ranks < l], return_counts=True) + niche_count[values] = counts + + front_l = union_pop[union_pareto_ranks == l] + front_l_closest_ref = union_closest_ref[union_pareto_ranks == l] + front_l_ref_dist = union_ref_dist[union_pareto_ranks == l] + # Niching + k = 1 + while k < solutions_left_count: + min_index = self.random_state_.choice(np.nonzero(niche_count == niche_count.min())[0]) + candidates = front_l[front_l_closest_ref == min_index] + if len(candidates) == 0: + niche_count[min_index] = np.inf + continue + if niche_count[min_index] == 0: + candidate_distances = front_l_ref_dist[front_l_closest_ref == min_index] + next_pop.append(candidates[np.argmin(candidate_distances)]) + else: + next_pop.append(self.random_state_.choice(candidates)) + niche_count[min_index] = niche_count[min_index] + 1 + front_l_closest_ref = front_l_closest_ref[front_l != next_pop[-1]] + front_l_ref_dist = front_l_ref_dist[front_l != next_pop[-1]] + front_l = front_l[front_l != next_pop[-1]] + k += 1 + + self.population_ = next_pop + + if self.check_early_stopping(): + break + + def pareto_front(self) -> list[Solution]: + if not hasattr(self, "population_") or not self.population_: + return [] + fitness_values = np.array([solution.fitness_ for solution in self.population_]) + pareto_ranks = fast_non_dominated_sort(fitness_values) + pareto_front = np.array(self.population_)[pareto_ranks == 0] + return sorted(pareto_front, key=lambda x: x.fitness_[0], reverse=True) diff --git a/suprb/optimizer/solution/nsga3/crossover.py b/suprb/optimizer/solution/nsga3/crossover.py new file mode 100644 index 00000000..f9ca683c --- /dev/null +++ b/suprb/optimizer/solution/nsga3/crossover.py @@ -0,0 +1,55 @@ +from abc import ABCMeta, abstractmethod + +import numpy as np + +from suprb.base import BaseComponent +from suprb.solution import Solution +from suprb.utils import RandomState + + +class SolutionCrossover(BaseComponent, metaclass=ABCMeta): + + def __init__(self, crossover_rate: float = 0.9): + self.crossover_rate = crossover_rate + + def __call__(self, A: Solution, B: Solution, random_state: RandomState) -> Solution: + if random_state.random() < self.crossover_rate: + return self._crossover(A=A, B=B, random_state=random_state) + else: + # Just return the primary parent + return A + + @abstractmethod + def _crossover(self, A: Solution, B: Solution, random_state: RandomState) -> Solution: + pass + + +class NPoint(SolutionCrossover): + """Cut the genome at N points and alternate the pieces from solution A and B.""" + + def __init__(self, crossover_rate: float = 0.9, n: int = 2): + super().__init__(crossover_rate=crossover_rate) + self.n = n + + @staticmethod + def _single_point(A: Solution, B: Solution, index: int) -> Solution: + return A.clone(genome=np.append(A.genome[:index], B.genome[index:])) + + def _crossover(self, A: Solution, B: Solution, random_state: RandomState) -> Solution: + indices = random_state.choice(np.arange(len(A.genome)), size=min(self.n, len(A.genome)), replace=False) + for index in indices: + A = self._single_point(A, B, index) + B = self._single_point(B, A, index) + return A + + +class Uniform(SolutionCrossover): + """Decide for every bit with uniform probability if the bit in genome A or B is used.""" + + def _crossover(self, A: Solution, B: Solution, random_state: RandomState) -> Solution: + indices = random_state.random(size=len(A.genome)) <= 0.5 + genome = np.empty(A.genome.shape) + genome[indices] = A.genome[indices] + genome[~indices] = B.genome[~indices] + + return A.clone(genome=genome) diff --git a/suprb/optimizer/solution/nsga3/mutation.py b/suprb/optimizer/solution/nsga3/mutation.py new file mode 100644 index 00000000..5f1e1e61 --- /dev/null +++ b/suprb/optimizer/solution/nsga3/mutation.py @@ -0,0 +1,25 @@ +from abc import ABCMeta + +import numpy as np + +from suprb.base import BaseComponent +from suprb.solution import Solution +from suprb.utils import RandomState + + +class SolutionMutation(BaseComponent, metaclass=ABCMeta): + + def __init__(self, mutation_rate: float = 0.05): + self.mutation_rate = mutation_rate + + def __call__(self, solution: Solution, random_state: RandomState) -> Solution: + pass + + +class BitFlips(SolutionMutation): + """Flips every bit in the genome with probability `mutation_rate`.""" + + def __call__(self, solution: Solution, random_state: RandomState) -> Solution: + bit_flips = random_state.random(solution.genome.shape) < self.mutation_rate + genome = np.logical_xor(solution.genome, bit_flips) + return solution.clone(genome=genome) diff --git a/suprb/optimizer/solution/nsga3/normalise.py b/suprb/optimizer/solution/nsga3/normalise.py new file mode 100644 index 00000000..d1a4e2c7 --- /dev/null +++ b/suprb/optimizer/solution/nsga3/normalise.py @@ -0,0 +1,165 @@ +import traceback +from abc import ABCMeta, abstractmethod + +from suprb.solution import Solution + +import numpy as np + +from suprb.base import BaseComponent + + +class NSGAIIINormaliser(BaseComponent, metaclass=ABCMeta): + + def __init__(self) -> None: + self._ideal_point = None + self._nadir_point = None + + def reset(self): + self._ideal_point = None + self._nadir_point = None + + def _update_ideal_point(self, fitness_values: np.ndarray) -> np.ndarray: + if self._ideal_point is not None: + fitness_values = np.concatenate([fitness_values, self._ideal_point[np.newaxis, :]]) + + self._ideal_point = np.min(fitness_values, axis=0) + return self._ideal_point + + @abstractmethod + def _update_nadir_point(self, fitness_values: np.ndarray, pareto_ranks: np.ndarray) -> np.ndarray: + pass + + def __call__(self, fitness_values: np.ndarray, pareto_ranks: np.ndarray) -> np.ndarray: + """ + Parameters + ---------- + fitness_values : np.ndarray + Fitness values from which to update ideal and nadir points. + pareto_ranks: np.ndarray + Associated levels of pareto domination used to update nadir points. + + Returns + ------- + np.ndarray + Normalised fitness values. + """ + self._update_ideal_point(fitness_values) + self._update_nadir_point(fitness_values, pareto_ranks) + return (fitness_values - self._ideal_point) / (self._nadir_point - self._ideal_point) + + +class HyperPlaneNormaliser(NSGAIIINormaliser): + """ + Implemented as described in 10.1007/978-3-030-12598-1_19 + """ + + def __init__(self, objective_count, epsilon_asf: float = 10e-6, epsilon_nad: float = 10e-6) -> None: + super().__init__() + self.objective_count = objective_count + self.epsilon_asf = epsilon_asf + self._extreme_points: np.ndarray | None = None + self._worst_points_estimate = np.zeros(self.objective_count) + self.epsilon_nad = epsilon_nad + + def reset(self): + super().reset() + self._extreme_points = None + self._worst_points_estimate = np.zeros(self.objective_count) + + def _update_extreme_points(self, fitness_values: np.ndarray) -> np.ndarray: + if self._extreme_points is not None: + fitness_values = np.concatenate([fitness_values, self._extreme_points]) + else: + self._extreme_points = np.zeros((self.objective_count, self.objective_count)) + + for obj_idx in range(self.objective_count): + weights = np.ones(self.objective_count) * self.epsilon_asf + weights[obj_idx] = 1 + scalarised = asf(fitness_values, self._ideal_point, weights) + k = np.argmin(scalarised) + self._extreme_points[obj_idx] = fitness_values[k] + return self._extreme_points + + def _update_nadir_point(self, fitness_values: np.ndarray, pareto_ranks: np.ndarray) -> np.ndarray: + self._worst_points_estimate = np.max( + np.concatenate([fitness_values, self._worst_points_estimate[None, :]]), axis=0 + ) + self._update_extreme_points(fitness_values) + try: + normal, distance = find_plane_from_points(self._extreme_points - self._ideal_point) + intercepts = find_hyperplane_axis_intercepts(normal, distance) + nadir_point = np.zeros(self.objective_count) + for k in range(self.objective_count): + nadir_k = intercepts[k] + self._ideal_point[k] + if intercepts[k] < self.epsilon_nad: + raise ValueError("Negative intercept") + if nadir_k > self._worst_points_estimate[k]: + raise ValueError("Nadir point outside worst point") + nadir_point[k] = intercepts[k] + self._ideal_point[k] + + except Exception as e: + nadir_point = np.max(fitness_values, axis=0) + + fallback_nadir = np.max(fitness_values[pareto_ranks == 0], axis=0) + mask = nadir_point - self._ideal_point < self.epsilon_nad + nadir_point[mask] = fallback_nadir[mask] + self._nadir_point = nadir_point + return nadir_point + + +def find_hyperplane_axis_intercepts(normal_vector, distance) -> np.ndarray: + """ + Compute intercepts of a hyperplane with coordinate axes. + + Parameters + ---------- + normal_vector: np.ndarray + Normal vector of the hyperplane [n1, n2, ..., nn]. + distance: float + Distance of the hyperplane from the origin (d in equation). + + Returns + ------- + intercepts: np.ndarray + Intercepts with each axis as an n dimensional vector. + """ + normal_vector = np.array(normal_vector) + objective_count = normal_vector.shape[0] + intercepts = np.zeros(objective_count) + + for k in range(objective_count): + if normal_vector[k] != 0: + # Set all other coordinates to zero and solve for this axis + point = np.zeros_like(normal_vector) + intercepts[k] = -distance / normal_vector[k] + else: + # No intersection if normal component is zero + raise ValueError("Normal vector must not contain zero components.") + return intercepts + + +def find_plane_from_points(points: np.ndarray) -> tuple[np.ndarray, float]: + """ + Compute a n-dimensional hyperplane from a set of n points. + + Parameters + ---------- + points: np.ndarray + Set of points [n1, n2, ..., nn]. from which the Hyperplane is computed. + + Returns + ------- + normal, distance: tuple[np.ndarray, float] + The normal vector of the hyperplane as an n dimensional numpy array and the orthogonal distance. + """ + center = np.mean(points, axis=0) + centered_points = points - center + _, _, vh = np.linalg.svd(centered_points) + normal = vh[-1] + normal = normal / np.linalg.norm(normal) + distance = -np.dot(normal, center) + return normal, distance + + +def asf(fitness_values: np.ndarray, ideal_point: np.ndarray, weights: np.ndarray) -> float: + return np.max((fitness_values - ideal_point) / weights, axis=1) diff --git a/suprb/optimizer/solution/nsga3/reference.py b/suprb/optimizer/solution/nsga3/reference.py new file mode 100644 index 00000000..09676be7 --- /dev/null +++ b/suprb/optimizer/solution/nsga3/reference.py @@ -0,0 +1,65 @@ +from itertools import combinations_with_replacement +import numpy as np + + +def das_dennis_points(num_partitions: int, num_dimensions: int = 2) -> np.ndarray: + """ + Generate evenly distributed points within a simplex using the Das-Dennis methodology. + + Parameters + ---------- + num_partitions: int + The number of partitions (grids) to divide the simplex. + num_dimensions: + int The number of dimensions for the simplex. Default is 2. + + Returns + ------- + np.ndarray + A 2D NumPy array where each row represents a point within the simplex. + """ + + partitions = combinations_with_replacement(range(num_dimensions), num_partitions) + points = [] + for part in partitions: + coord = np.zeros(num_dimensions) + for idx in part: + coord[idx] += 1 + points.append(coord / num_partitions) + return np.array(points) + + +def calc_ref_direction_distances( + fitness_values: np.ndarray, reference_points: np.ndarray +) -> tuple[np.ndarray, np.ndarray]: + """ + Calculate distances between fitness values and reference points, and find the closest reference points. + + Parameters + ---------- + fitness_values : np.ndarray + A 2D array where each row represents a solution's fitness values. + reference_points : np.ndarray + A 2D array where each row represents a predefined reference direction. + + Returns + ------- + tuple[np.ndarray, np.ndarray] + A tuple containing: + - min_distances: A 1D array representing the minimum distance from each solution to the closest reference point. + - min_ref_points: A 1D array containing the index of the closest reference point for each solution. + """ + solution_count, _ = fitness_values.shape + reference_count, _ = reference_points.shape + min_distances = np.zeros(solution_count) + min_ref_points = np.zeros(solution_count, dtype=np.int32) + + for fit_idx in range(solution_count): + distances = np.zeros(reference_count) + for ref_idx in range(reference_points.shape[0]): + numerator = np.dot(np.dot(reference_points[ref_idx].T, fitness_values[fit_idx]), reference_points[ref_idx]) + denominator = np.linalg.norm(reference_points[ref_idx]) ** 2 + distances[ref_idx] = np.linalg.norm(fitness_values[fit_idx] - (numerator / denominator)) + min_distances[fit_idx] = np.min(distances) + min_ref_points[fit_idx] = np.argmin(distances) + return min_distances, min_ref_points diff --git a/suprb/optimizer/solution/nsga3/selection.py b/suprb/optimizer/solution/nsga3/selection.py new file mode 100644 index 00000000..40a5291e --- /dev/null +++ b/suprb/optimizer/solution/nsga3/selection.py @@ -0,0 +1,63 @@ +from abc import ABCMeta + +import numpy as np + +from suprb.base import BaseComponent +from suprb.solution import Solution +from suprb.utils import RandomState + + +class SolutionSelection(BaseComponent, metaclass=ABCMeta): + + def __call__( + self, + population: list[Solution], + n: int, + random_state: RandomState, + pareto_ranks: np.ndarray, + closest_ref_direction: np.ndarray, + ref_direction_distance: np.ndarray, + ) -> list[Solution]: + pass + + +class ReferenceBasedBinaryTournament(SolutionSelection): + """Draw 2 solutions n_parents times and select the best solution from each pair.""" + + def __init__(self): + pass + + def __call__( + self, + population: list[Solution], + n: int, + random_state: RandomState, + pareto_ranks: np.ndarray, + closest_ref_direction: np.ndarray, + ref_direction_distance: np.ndarray, + ) -> list[Solution]: + + selection = [] + pairs = random_state.integers(low=0, high=len(population), size=(n, 2)) + for a, b in pairs: + # select one solution randomly if they are not associated to the same reference direction + if closest_ref_direction[a] != closest_ref_direction[b]: + selection.append(population[a]) + continue + # select the solution that is closer to the reference line if pareto ranks are equal + if pareto_ranks[a] == pareto_ranks[b]: + if ref_direction_distance[a] < ref_direction_distance[b]: + selection.append(population[a]) + continue + else: + selection.append(population[b]) + continue + # If pareto ranks are not equal select the less dominated solution + if pareto_ranks[a] < pareto_ranks[b]: + selection.append(population[a]) + continue + else: + selection.append(population[b]) + continue + + return selection diff --git a/suprb/optimizer/solution/nsga3/sorting.py b/suprb/optimizer/solution/nsga3/sorting.py new file mode 100644 index 00000000..3d82ad5b --- /dev/null +++ b/suprb/optimizer/solution/nsga3/sorting.py @@ -0,0 +1,89 @@ +from typing import Tuple + +import numpy as np + + +def fast_non_dominated_sort(fitness_values: np.ndarray) -> np.ndarray: + """Sorts the fitness values into multiple levels of non domination. + + Parameters + ---------- + fitness_values: np.ndarray + numpy of shape (solution_count, objective_count) with all fitness values for every solution + + Returns + ------- + pareto ranks: level of non-dominated front for each solution (non sorted) + + """ + solution_count, objective_count = fitness_values.shape + pareto_ranks = np.ones(solution_count, dtype=np.int32) * -1 + dominated_count = np.zeros(solution_count, dtype=np.int32) + dominates_indices = [[] for _ in range(solution_count)] + + # Calculating the amount of solutions, that dominate i and + # the solutions i is dominating + for i in range(solution_count): + fitness_i = fitness_values[i][None, :] + + less_equal_count = np.sum(fitness_values <= fitness_i, axis=1) + less_than_count = np.sum(fitness_values < fitness_i, axis=1) + dominates_i_mask = (less_equal_count == objective_count) & (less_than_count >= 1) + dominated_count[i] = np.sum(dominates_i_mask) + + greater_equal_count = np.sum(fitness_values >= fitness_i, axis=1) + greater_than_count = np.sum(fitness_values > fitness_i, axis=1) + dominated_by_i_mask = (greater_equal_count == objective_count) & (greater_than_count >= 1) + dominates_indices[i] = np.argwhere(dominated_by_i_mask).flatten().tolist() + + pareto_ranks[dominated_count == 0] = 0 + front_rank = 0 + + # Actual non dominated sort + while -1 in pareto_ranks: + current_front = np.argwhere(pareto_ranks == front_rank).flatten() + front_rank += 1 + for solution in current_front: + solution_mask = np.zeros(solution_count, dtype=bool) + solution_mask[dominates_indices[solution]] = 1 + solution_mask = solution_mask & (dominated_count == 1) + dominated_count[dominates_indices[solution]] -= 1 + pareto_ranks[solution_mask] = front_rank + + return pareto_ranks + + +def calculate_crowding_distances(fitness_values: np.ndarray, pareto_ranks: np.ndarray) -> np.ndarray: + solution_count, objective_count = fitness_values.shape + crowding_distances = np.zeros(solution_count) + + for rank in np.unique(pareto_ranks): + front_indices = np.where(pareto_ranks == rank)[0] + front_size = len(front_indices) + + if front_size <= 2: + crowding_distances[front_indices] = np.inf + continue + + front_fitness = fitness_values[front_indices] + + for m in range(objective_count): + sorting_permutation = np.argsort(front_fitness[:, m]) + sorted_front = front_indices[sorting_permutation] + + crowding_distances[sorted_front[0]] = np.inf + crowding_distances[sorted_front[-1]] = np.inf + + min_f = front_fitness[sorting_permutation[0], m] + max_f = front_fitness[sorting_permutation[-1], m] + + # if max_f == min_f the crowding distance parts that result from objective_m are all 0 as all + # solution share the same coordinate in this dimension of the fitness function + if max_f > min_f: + normalized_range = max_f - min_f + for i in range(1, front_size - 1): + crowding_distances[sorted_front[i]] += ( + front_fitness[sorting_permutation[i + 1], m] - front_fitness[sorting_permutation[i - 1], m] + ) / normalized_range + + return crowding_distances diff --git a/suprb/optimizer/solution/sampler.py b/suprb/optimizer/solution/sampler.py new file mode 100644 index 00000000..bebdbd31 --- /dev/null +++ b/suprb/optimizer/solution/sampler.py @@ -0,0 +1,84 @@ +from abc import ABCMeta, abstractmethod +import scipy.stats as stats +from typing import Optional, Callable +from sklearn.utils import Bunch + +from suprb.base import BaseComponent +from suprb.solution import Solution +from suprb.utils import RandomState + +import numpy as np + + +def calculate_crowding_distances(fitness_values: np.ndarray, pareto_ranks: np.ndarray) -> np.ndarray: + solution_count, objective_count = fitness_values.shape + crowding_distances = np.zeros(solution_count) + + for rank in np.unique(pareto_ranks): + front_indices = np.where(pareto_ranks == rank)[0] + front_size = len(front_indices) + + if front_size <= 2: + crowding_distances[front_indices] = np.inf + continue + + front_fitness = fitness_values[front_indices] + + for m in range(objective_count): + sorting_permutation = np.argsort(front_fitness[:, m]) + sorted_front = front_indices[sorting_permutation] + + crowding_distances[sorted_front[0]] = np.inf + crowding_distances[sorted_front[-1]] = np.inf + + min_f = front_fitness[sorting_permutation[0], m] + max_f = front_fitness[sorting_permutation[-1], m] + + # if max_f == min_f the crowding distance parts that result from objective_m are all 0 as all + # solution share the same coordinate in this dimension of the fitness function + if max_f > min_f: + normalized_range = max_f - min_f + for i in range(1, front_size - 1): + crowding_distances[sorted_front[i]] += ( + front_fitness[sorting_permutation[i + 1], m] - front_fitness[sorting_permutation[i - 1], m] + ) / normalized_range + + return crowding_distances + + +class SolutionSampler(BaseComponent, metaclass=ABCMeta): + + @abstractmethod + def __call__(self, pareto_front: list[Solution], random_state: RandomState) -> Solution: + pass + + +class BetaSolutionSampler(SolutionSampler): + + def __init__(self, a: float = 1.5, b: float = 1.5, projected: bool = True): + self.a = a + self.b = b + self.projected = projected + + def __call__(self, pareto_front: list[Solution], random_state: RandomState) -> Solution: + if self.projected: + points = np.array([solution.fitness_ for solution in pareto_front]) + points = points / np.sum(points, axis=1, keepdims=True) + points = points[:, 0] + else: + points = np.linspace(0.0001, 1 - 0.0001, len(pareto_front)) + + weights = stats.beta.pdf(points, a=self.a, b=self.b) + weights[weights == np.inf] = 0 + weights = (weights / np.sum(weights)) if np.sum(weights) > 0 else (np.ones(len(weights)) / len(weights)) + return random_state.choice(pareto_front, p=weights) + + +class DiversitySolutionSampler(SolutionSampler): + + def __call__(self, pareto_front: list[Solution], random_state: RandomState) -> Solution: + fitness_values = np.array([solution.fitness_ for solution in pareto_front]) + weights = calculate_crowding_distances(fitness_values, np.zeros(len(pareto_front))) + weights[weights == np.inf] = 0 + weights = (weights / np.sum(weights)) if np.sum(weights) > 0 else (np.ones(len(weights)) / len(weights)) + return random_state.choice(pareto_front, p=weights) diff --git a/suprb/optimizer/solution/spea2/__init__.py b/suprb/optimizer/solution/spea2/__init__.py new file mode 100644 index 00000000..f3ca7b29 --- /dev/null +++ b/suprb/optimizer/solution/spea2/__init__.py @@ -0,0 +1,5 @@ +from .base import StrengthParetoEvolutionaryAlgorithm2 +from .crossover import SolutionCrossover +from .mutation import SolutionMutation +from .selection import SolutionSelection +from .archive import EnvironmentalArchive diff --git a/suprb/optimizer/solution/spea2/archive.py b/suprb/optimizer/solution/spea2/archive.py new file mode 100644 index 00000000..6ec4284b --- /dev/null +++ b/suprb/optimizer/solution/spea2/archive.py @@ -0,0 +1,50 @@ +from suprb.optimizer.solution.archive import SolutionArchive +from suprb.solution import Solution +from .internal_fitness import calculate_raw_internal_fitness, calculate_density, distance_to_kth +import numpy as np + + +class EnvironmentalArchive(SolutionArchive): + + def __init__(self, max_population_size, kth_nearest): + super().__init__() + self.max_population_size = max_population_size + self.kth_nearest = kth_nearest + + def __call__(self, new_population: list[Solution]): + pop_size = len(new_population) + fitness_values = np.array([solution.fitness_ for solution in new_population + self.population_]) + pop_and_arch = np.array([solution for solution in new_population + self.population_]) + raw_internal_fitness_values = calculate_raw_internal_fitness(fitness_values) + density_values = calculate_density(fitness_values, self.kth_nearest) + internal_fitness_values = raw_internal_fitness_values + density_values + + # Pass the internal fitness values to the optimizer to avoid recalculation ...? + for i, solution in enumerate(pop_and_arch): + solution.internal_fitness_ = internal_fitness_values[i] + + sorting_permutation = np.argsort(internal_fitness_values) + internal_fitness_values = internal_fitness_values[sorting_permutation] + pop_and_arch = pop_and_arch[sorting_permutation].tolist() + # We can index the 0 directly as there is always a pareto dominant element in a finite set. + non_dominated_count = np.unique(internal_fitness_values, return_counts=True)[1][0] + + if non_dominated_count < self.max_population_size: + self.population_ = pop_and_arch[: self.max_population_size] + else: + # In this case the population needs to be truncated! + self.population_ = pop_and_arch[:non_dominated_count] + # TODO: Vectorize this to achieve better performance + while len(self.population_) > self.max_population_size: + archive_fitness_values = np.array([solution.fitness_ for solution in self.population_]) + distances = np.zeros(len(self.population_)) + candidates_mask = np.ones(len(distances), dtype=np.bool) + for k in range(1, len(self.population_)): + for i in range(len(self.population_)): + distances[i] = distance_to_kth(archive_fitness_values[i], archive_fitness_values, k) + unique = np.unique(distances[candidates_mask])[0] + candidates_mask = (distances == unique) & candidates_mask + if np.sum(candidates_mask) == 1 or k == len(self.population_) - 1: + index = np.argmax(candidates_mask) + self.population_.pop(index) + break diff --git a/suprb/optimizer/solution/spea2/base.py b/suprb/optimizer/solution/spea2/base.py new file mode 100644 index 00000000..ae84656c --- /dev/null +++ b/suprb/optimizer/solution/spea2/base.py @@ -0,0 +1,130 @@ +import numpy as np +import scipy.stats as stats + +from suprb import Solution +from suprb.solution.fitness import NormalizedMOSolutionFitness +from suprb.solution.initialization import SolutionInit, RandomInit +from suprb.utils import flatten +from .sorting import fast_non_dominated_sort + +from ..base import MOSolutionComposition +from .mutation import SolutionMutation, BitFlips +from .selection import SolutionSelection, BinaryTournament +from .crossover import SolutionCrossover, NPoint +from .archive import EnvironmentalArchive +from ..sampler import SolutionSampler, BetaSolutionSampler +from .internal_fitness import calculate_raw_internal_fitness + + +class StrengthParetoEvolutionaryAlgorithm2(MOSolutionComposition): + """Stringth Pareto Evolutionary Algorithm 2. + + Implemented as described in 10.3929/ethz-a-004284029 + + Parameters + ---------- + n_iter: int + Iterations the metaheuristic will perform. + population_size: int + Number of solutions in the population. + mutation: SolutionMutation + crossover: SolutionCrossover + selection: SolutionSelection + init: SolutionInit + random_state : int, RandomState instance or None, default=None + Pass an int for reproducible results across multiple function calls. + warm_start: bool + If False, solutions are generated new for every `optimize()` call. + If True, solutions are used from previous runs. + n_jobs: int + The number of threads / processes the optimization uses. + """ + + def __init__( + self, + n_iter: int = 32, + population_size: int = 32, + archive_size: int = 32, + mutation: SolutionMutation = BitFlips(), + crossover: SolutionCrossover = NPoint(n=3), + selection: SolutionSelection = BinaryTournament(), + sampler: SolutionSampler = BetaSolutionSampler(1.5, 1.5), + mutation_rate: float = 0.025, + crossover_rate: float = 0.75, + kth_nearest: int = -1, + init: SolutionInit = RandomInit(fitness=NormalizedMOSolutionFitness()), + random_state: int = None, + n_jobs: int = 1, + warm_start: bool = True, + early_stopping_patience: int = -1, + early_stopping_delta: float = 0, + ): + super().__init__( + n_iter=n_iter, + population_size=population_size, + init=init, + archive=EnvironmentalArchive( + archive_size, kth_nearest if kth_nearest != -1 else int((population_size + archive_size) ** 0.5) + ), + sampler=sampler, + random_state=random_state, + n_jobs=n_jobs, + warm_start=warm_start, + early_stopping_patience=early_stopping_patience, + early_stopping_delta=early_stopping_delta, + ) + self.mutation = mutation + self.crossover = crossover + self.selection = selection + + self.mutation_rate = mutation_rate + self.crossover_rate = crossover_rate + + self.kth_nearest = kth_nearest + self.archive_size = archive_size + + def _optimize(self, X: np.ndarray, y: np.ndarray): + self.fit_population(X, y) + for _ in range(self.n_iter): + self.archive(self.population_) + internal_fitness_values = np.array([solution.internal_fitness_ for solution in self.archive.population_]) + parents = self.selection( + population=self.archive.population_, + n=self.population_size, + random_state=self.random_state_, + internal_fitness=internal_fitness_values, + ) + + # Note that this expression swallows the last element, if `population_size` is odd + parent_pairs = map(lambda *x: x, *([iter(parents)] * 2)) + + # Crossover + children = list( + flatten( + [ + ( + self.crossover(A, B, random_state=self.random_state_), + self.crossover(B, A, random_state=self.random_state_), + ) + for A, B in parent_pairs + ] + ) + ) + # If `population_size` is odd, we add the solution not selected for reproduction directly + if self.population_size % 2 != 0: + children.append(parents[-1]) + + # Mutation + mutated_children = [self.mutation(child, random_state=self.random_state_).fit(X, y) for child in children] + self.population_ = mutated_children + + if self.check_early_stopping(): + break + + def pareto_front(self) -> list[Solution]: + if not hasattr(self, "population_") or not self.population_: + return [] + fitness_values = np.array([solution.fitness_ for solution in self.archive.population_]) + pareto_ranks = fast_non_dominated_sort(fitness_values) + pareto_front = np.array(self.archive.population_)[pareto_ranks == 0] + return sorted(pareto_front, key=lambda x: x.fitness_[0], reverse=True) diff --git a/suprb/optimizer/solution/spea2/crossover.py b/suprb/optimizer/solution/spea2/crossover.py new file mode 100644 index 00000000..f9ca683c --- /dev/null +++ b/suprb/optimizer/solution/spea2/crossover.py @@ -0,0 +1,55 @@ +from abc import ABCMeta, abstractmethod + +import numpy as np + +from suprb.base import BaseComponent +from suprb.solution import Solution +from suprb.utils import RandomState + + +class SolutionCrossover(BaseComponent, metaclass=ABCMeta): + + def __init__(self, crossover_rate: float = 0.9): + self.crossover_rate = crossover_rate + + def __call__(self, A: Solution, B: Solution, random_state: RandomState) -> Solution: + if random_state.random() < self.crossover_rate: + return self._crossover(A=A, B=B, random_state=random_state) + else: + # Just return the primary parent + return A + + @abstractmethod + def _crossover(self, A: Solution, B: Solution, random_state: RandomState) -> Solution: + pass + + +class NPoint(SolutionCrossover): + """Cut the genome at N points and alternate the pieces from solution A and B.""" + + def __init__(self, crossover_rate: float = 0.9, n: int = 2): + super().__init__(crossover_rate=crossover_rate) + self.n = n + + @staticmethod + def _single_point(A: Solution, B: Solution, index: int) -> Solution: + return A.clone(genome=np.append(A.genome[:index], B.genome[index:])) + + def _crossover(self, A: Solution, B: Solution, random_state: RandomState) -> Solution: + indices = random_state.choice(np.arange(len(A.genome)), size=min(self.n, len(A.genome)), replace=False) + for index in indices: + A = self._single_point(A, B, index) + B = self._single_point(B, A, index) + return A + + +class Uniform(SolutionCrossover): + """Decide for every bit with uniform probability if the bit in genome A or B is used.""" + + def _crossover(self, A: Solution, B: Solution, random_state: RandomState) -> Solution: + indices = random_state.random(size=len(A.genome)) <= 0.5 + genome = np.empty(A.genome.shape) + genome[indices] = A.genome[indices] + genome[~indices] = B.genome[~indices] + + return A.clone(genome=genome) diff --git a/suprb/optimizer/solution/spea2/internal_fitness.py b/suprb/optimizer/solution/spea2/internal_fitness.py new file mode 100644 index 00000000..1f6eed0b --- /dev/null +++ b/suprb/optimizer/solution/spea2/internal_fitness.py @@ -0,0 +1,73 @@ +import numpy as np + + +def calculate_raw_internal_fitness(fitness_values: np.ndarray) -> np.ndarray: + """Calculates the raw internal fitness values R(i) for each solution. + These are the added strength values of all dominators of i. + + Parameters + ---------- + fitness_values: np.ndarray + numpy array of shape (solution_count, objective_count) with all fitness values for every solution + + Returns + ------- + raw_internal_fitness_values: np.ndarray + 1D numpy array of length solution_count with the raw internal fitness values for every solution + """ + solution_count, objective_count = fitness_values.shape + strength_values = np.zeros(solution_count, dtype=np.int32) + dominator_indices = [[]] * solution_count + raw_internal_fitness_values = np.zeros(solution_count) + + # First we calculate the strength values S(i) and the indices of Solutions that dominate i. + for i in range(solution_count): + fitness_i = fitness_values[i][None, :] + greater_equal_count = np.sum(fitness_values >= fitness_i, axis=1) + greater_than_count = np.sum(fitness_values > fitness_i, axis=1) + dominated_by_i_mask = (greater_equal_count == objective_count) & (greater_than_count >= 1) + strength_values[i] = np.sum(dominated_by_i_mask) + + less_equal_count = np.sum(fitness_values <= fitness_i, axis=1) + less_than_count = np.sum(fitness_values < fitness_i, axis=1) + dominates_i_mask = (less_equal_count == objective_count) & (less_than_count >= 1) + dominator_indices[i] = np.argwhere(dominates_i_mask).flatten().tolist() + + # The raw internal fitness values are calculated by summing up the strength values of + # all dominators of one solution. + + for i in range(solution_count): + if len(dominator_indices[i]) > 0: + raw_internal_fitness_values[i] = np.sum(strength_values[np.array(dominator_indices[i])]) + + return raw_internal_fitness_values + + +def calculate_density(fitness_values: np.ndarray, k: int) -> np.ndarray: + """Calculates the density D(i) for each solution. + + Parameters + ---------- + fitness_values: np.ndarray + numpy array of shape (solution_count, objective_count) with all fitness values for every solution + k: int + The distance to the k-th nearest neighbour is selected + + Returns + ------- + density_values: np.ndarray + 1D numpy array of length solution_count with the density values for every solution + """ + solution_count, objective_count = fitness_values.shape + density_values = np.zeros(solution_count, dtype=np.int32) + + for i in range(solution_count): + density_values[i] = distance_to_kth(fitness_values[i], fitness_values, k) + + return density_values + + +def distance_to_kth(fitness_i: np.ndarray, fitness_values: np.ndarray, k: int) -> np.ndarray: + distances = np.linalg.norm(fitness_i[None, :] - fitness_values, axis=-1) + distances = np.sort(distances) + return distances[k] diff --git a/suprb/optimizer/solution/spea2/mutation.py b/suprb/optimizer/solution/spea2/mutation.py new file mode 100644 index 00000000..5f1e1e61 --- /dev/null +++ b/suprb/optimizer/solution/spea2/mutation.py @@ -0,0 +1,25 @@ +from abc import ABCMeta + +import numpy as np + +from suprb.base import BaseComponent +from suprb.solution import Solution +from suprb.utils import RandomState + + +class SolutionMutation(BaseComponent, metaclass=ABCMeta): + + def __init__(self, mutation_rate: float = 0.05): + self.mutation_rate = mutation_rate + + def __call__(self, solution: Solution, random_state: RandomState) -> Solution: + pass + + +class BitFlips(SolutionMutation): + """Flips every bit in the genome with probability `mutation_rate`.""" + + def __call__(self, solution: Solution, random_state: RandomState) -> Solution: + bit_flips = random_state.random(solution.genome.shape) < self.mutation_rate + genome = np.logical_xor(solution.genome, bit_flips) + return solution.clone(genome=genome) diff --git a/suprb/optimizer/solution/spea2/selection.py b/suprb/optimizer/solution/spea2/selection.py new file mode 100644 index 00000000..0a7a92f6 --- /dev/null +++ b/suprb/optimizer/solution/spea2/selection.py @@ -0,0 +1,39 @@ +from abc import ABCMeta + +import numpy as np + +from suprb.base import BaseComponent +from suprb.solution import Solution +from suprb.utils import RandomState + + +class SolutionSelection(BaseComponent, metaclass=ABCMeta): + + def __call__( + self, population: list[Solution], n: int, random_state: RandomState, internal_fitness: np.ndarray + ) -> list[Solution]: + pass + + +class BinaryTournament(SolutionSelection): + """Draw 2 solutions n_parents times and select the best solution from each pair.""" + + def __init__(self): + pass + + def __call__( + self, + population: list[Solution], + n: int, + random_state: RandomState, + internal_fitness: np.ndarray, + ) -> list[Solution]: + + selection = [] + pairs = random_state.integers(low=0, high=len(population), size=(n, 2)) + for a, b in pairs: + if internal_fitness[a] < internal_fitness[b]: + selection.append(population[a]) + else: + selection.append(population[b]) + return selection diff --git a/suprb/optimizer/solution/spea2/sorting.py b/suprb/optimizer/solution/spea2/sorting.py new file mode 100644 index 00000000..84f91262 --- /dev/null +++ b/suprb/optimizer/solution/spea2/sorting.py @@ -0,0 +1,53 @@ +from typing import Tuple + +import numpy as np + + +def fast_non_dominated_sort(fitness_values: np.ndarray) -> np.ndarray: + """Sorts the fitness values into multiple levels of non domination. + + Parameters + ---------- + fitness_values: np.ndarray + numpy of shape (solution_count, objective_count) with all fitness values for every solution + + Returns + ------- + pareto ranks: level of non-dominated front for each solution (non sorted) + + """ + solution_count, objective_count = fitness_values.shape + pareto_ranks = np.ones(solution_count, dtype=np.int32) * -1 + dominated_count = np.zeros(solution_count, dtype=np.int32) + dominates_indices = [[] for _ in range(solution_count)] + + # Calculating the amount of solutions, that dominate i and + # the solutions i is dominating + for i in range(solution_count): + fitness_i = fitness_values[i][None, :] + + less_equal_count = np.sum(fitness_values <= fitness_i, axis=1) + less_than_count = np.sum(fitness_values < fitness_i, axis=1) + dominates_i_mask = (less_equal_count == objective_count) & (less_than_count >= 1) + dominated_count[i] = np.sum(dominates_i_mask) + + greater_equal_count = np.sum(fitness_values >= fitness_i, axis=1) + greater_than_count = np.sum(fitness_values > fitness_i, axis=1) + dominated_by_i_mask = (greater_equal_count == objective_count) & (greater_than_count >= 1) + dominates_indices[i] = np.argwhere(dominated_by_i_mask).flatten().tolist() + + pareto_ranks[dominated_count == 0] = 0 + front_rank = 0 + + # Actual non dominated sort + while -1 in pareto_ranks: + current_front = np.argwhere(pareto_ranks == front_rank).flatten() + front_rank += 1 + for solution in current_front: + solution_mask = np.zeros(solution_count, dtype=bool) + solution_mask[dominates_indices[solution]] = 1 + solution_mask = solution_mask & (dominated_count == 1) + dominated_count[dominates_indices[solution]] -= 1 + pareto_ranks[solution_mask] = front_rank + + return pareto_ranks diff --git a/suprb/optimizer/solution/ts/__init__.py b/suprb/optimizer/solution/ts/__init__.py new file mode 100644 index 00000000..15366697 --- /dev/null +++ b/suprb/optimizer/solution/ts/__init__.py @@ -0,0 +1 @@ +from .base import TwoStageSolutionComposition diff --git a/suprb/optimizer/solution/ts/base.py b/suprb/optimizer/solution/ts/base.py new file mode 100644 index 00000000..7f30ca00 --- /dev/null +++ b/suprb/optimizer/solution/ts/base.py @@ -0,0 +1,96 @@ +from typing import Optional, Union +import copy + +import numpy as np + +from suprb import Solution +from suprb.solution import SolutionInit, SolutionFitness +from suprb.solution.fitness import NormalizedMOSolutionFitness +from suprb.solution.initialization import RandomInit +from ..base import SolutionComposition, SolutionArchive, MOSolutionComposition, PopulationBasedSolutionComposition +from suprb.utils import check_random_state + + +class TwoStageSolutionComposition(SolutionComposition): + """ + Wrapper class for switching between different solution composition algorithms + + Parameters + ---------- + algorithm_1: PopulationBasedSolutionComposition + Algorithm to be used in the first stage. + algorithm_2: PopulationBasedSolutionComposition + Algorithm to be used in the second stage. + switch_iteration: int + Iteration at which to switch from algorithm_1 to algorithm_2. + + """ + + def __init__( + self, + algorithm_1: PopulationBasedSolutionComposition, + algorithm_2: Optional[PopulationBasedSolutionComposition] = None, + switch_iteration: Optional[int] = None, + n_iter: int = None, + init: SolutionInit = RandomInit(), + archive: SolutionArchive = None, + random_state: int = None, + n_jobs: int = None, + warm_start: bool = True, + ): + super().__init__( + n_iter=n_iter, init=init, archive=archive, random_state=random_state, n_jobs=n_jobs, warm_start=warm_start + ) + self.algorithm_1 = algorithm_1 + self.algorithm_2 = algorithm_2 + self.switch_iteration = switch_iteration + + self._suprb_step = 0 + self.current_algo_ = algorithm_1 + + def optimize(self, X: np.ndarray, y: np.ndarray, **kwargs) -> Union[Solution, list[Solution], None]: + # This is mega hacky but necessary if the initialization is to be kept in the suprb fit method + if self._suprb_step == 0: + self.algorithm_1.pool_ = self.pool_ + self.algorithm_1.init.fitness.max_genome_length_ = self.init.fitness.max_genome_length_ + if self.algorithm_2 is not None: + self.algorithm_2.pool_ = self.pool_ + self.algorithm_2.init.fitness.max_genome_length_ = self.init.fitness.max_genome_length_ + + self.algorithm_1.random_state = self.random_state + if self.algorithm_2 is not None: + self.algorithm_2.random_state = self.random_state + self.random_state_ = check_random_state(self.random_state) + + if self.algorithm_2 is not None and self._suprb_step == self.switch_iteration - 1: + self.current_algo_ = self.algorithm_2 + if self.warm_start: + pop = [] + algo_2_fitness = self.algorithm_2.init.fitness + for solution in self.algorithm_1.population_: + solution.fitness = algo_2_fitness + solution = self.init.pad(solution, self.random_state_) + solution.fit(X, y) + pop.append(solution) + self.algorithm_2.population_ = pop + + self.current_algo_.random_state = self.random_state + self.current_algo_.optimize(X, y) + + self.population_ = self.current_algo_.population_ + + # If Early Stopping is supported, get the actual iteration count + if hasattr(self.current_algo_, "step_"): + self.step_ = self.current_algo_.step_ + else: + self.step_ = self.current_algo_.n_iter + self._suprb_step += 1 + + def elitist(self) -> Optional[Solution]: + return self.current_algo_.elitist() + + def pareto_front(self) -> list[Solution]: + if isinstance(self.current_algo_, MOSolutionComposition): + return self.current_algo_.pareto_front() + else: + return [self.current_algo_.elitist()] diff --git a/suprb/solution/__init__.py b/suprb/solution/__init__.py index 495bdeb5..f05f0428 100644 --- a/suprb/solution/__init__.py +++ b/suprb/solution/__init__.py @@ -1,3 +1,3 @@ from .base import Solution, MixingModel, SolutionFitness -from .fitness import ComplexitySolutionFitness +from .fitness import ComplexitySolutionFitness, MultiObjectiveSolutionFitness from .initialization import SolutionInit diff --git a/suprb/solution/base.py b/suprb/solution/base.py index 5cbe2931..fb036fb7 100644 --- a/suprb/solution/base.py +++ b/suprb/solution/base.py @@ -2,6 +2,7 @@ import itertools from abc import ABCMeta, abstractmethod +from typing import Tuple, Union import numpy as np from sklearn.base import RegressorMixin @@ -29,7 +30,7 @@ class SolutionFitness(BaseFitness, metaclass=ABCMeta): max_genome_length_: int @abstractmethod - def __call__(self, solution: Solution) -> float: + def __call__(self, solution: Solution) -> Union[float, Tuple[float, float]]: pass @@ -51,8 +52,8 @@ def __init__( self.mixing = mixing self.fitness = fitness - def fit(self, X: np.ndarray, y: np.ndarray) -> Solution: - pred = self.predict(X, cache=True) + def fit(self, X: np.ndarray, y: np.ndarray, cache=True) -> Solution: + pred = self.predict(X, cache=cache) self.error_ = max(mean_squared_error(y, pred), 1e-4) self.input_size_ = self.genome.shape[0] self.complexity_ = np.sum(self.genome).item() # equivalent to np.count_nonzero, but possibly faster @@ -67,7 +68,7 @@ def predict(self, X: np.ndarray, cache=False) -> np.ndarray: The attribute `cache` decides if the prediction should use the cached data from the rules, which includes rule error, fitness, predictions and the binary match string. It is True while fitting, because the data is identical there and caching saves a good amount of time. - For predictions after fitting, it is false, because all data needs to be recalculated from scratch. + For predictions after fitting, it is false because all data needs to be recalculated from scratch. """ return self.mixing(X=X, subpopulation=self.subpopulation, cache=cache) diff --git a/suprb/solution/fitness.py b/suprb/solution/fitness.py index 6fb25c95..c2529df9 100644 --- a/suprb/solution/fitness.py +++ b/suprb/solution/fitness.py @@ -1,5 +1,5 @@ from abc import ABCMeta -from typing import Callable +from typing import Callable, Tuple import numpy as np @@ -38,6 +38,41 @@ def __call__(self, solution: Solution) -> float: ) +class MultiObjectiveSolutionFitness(SolutionFitness, metaclass=ABCMeta): + """Passes on fitness and complexity measures for MOOAs to minimize complexity along with error.""" + + objective_func_: list[Callable] + hv_reference_: np.ndarray + worst_point_estimate_: np.ndarray + + def __init__(self): + pass + + def __call__(self, solution: Solution) -> list: + return [self.objective_func_[i](solution) for i in range(len(self.objective_func_))] + + +class NormalizedMOSolutionFitness(MultiObjectiveSolutionFitness): + objective_func_: list[Callable] = [c_norm, pseudo_accuracy] + hv_reference_: np.ndarray = np.array([1.0, 1.0]) + worst_point_estimate_ = np.array([1.0, 1.0]) + + def __call__(self, solution: Solution) -> list: + return [ + 1 - self.objective_func_[0](solution.complexity_, self.max_genome_length_), + 1 - self.objective_func_[1](solution.error_), + ] + + def __init__(self): + super().__init__() + + +class SimpleMOSolutionFitness(MultiObjectiveSolutionFitness): + + def __call__(self, solution: Solution): + return [solution.complexity_, solution.error_] + + class ComplexityEmary(ComplexitySolutionFitness): """ Mixes the pseudo-accuracy as primary objective x1 diff --git a/suprb/suprb.py b/suprb/suprb.py index 7b9869a7..e8860599 100644 --- a/suprb/suprb.py +++ b/suprb/suprb.py @@ -146,7 +146,6 @@ def fit(self, X: np.ndarray, y: np.ndarray, cleanup=False): cleanup : bool Optional cleanup of unused rules and components after fitting. Can be used to reduce size if only the final model is relevant. Note that all information about the fitting process itself is removed. - Returns ------- self : BaseEstimator @@ -281,7 +280,6 @@ def _compose_solution(self, X: np.ndarray, y: np.ndarray): # Update the random state self.solution_composition_.random_state = self.solution_composition_seeds_[self.step_] - # Optimize self.solution_composition_.optimize(X, y) diff --git a/tests/test_solution.py b/tests/test_solution.py index c3983199..4ae41dd3 100644 --- a/tests/test_solution.py +++ b/tests/test_solution.py @@ -5,6 +5,8 @@ import suprb import suprb.logging.stdout from suprb.optimizer.solution.ga import GeneticAlgorithm +from suprb.optimizer.solution.nsga2 import NonDominatedSortingGeneticAlgorithm2 +from suprb.optimizer.solution.nsga3 import NonDominatedSortingGeneticAlgorithm3 from suprb.optimizer.solution.saga import ( SelfAdaptingGeneticAlgorithm1, SelfAdaptingGeneticAlgorithm2, @@ -12,6 +14,8 @@ SasGeneticAlgorithm, ) from suprb.optimizer.rule.es import ES1xLambda +from suprb.optimizer.solution.spea2 import StrengthParetoEvolutionaryAlgorithm2 +from suprb.optimizer.solution.ts import TwoStageSolutionComposition class TestSolution(unittest.TestCase): @@ -84,3 +88,64 @@ def test_check_sas(self): estimator.fit(X, y) check_estimator(estimator) + + def test_check_nsga2(self): + estimator = suprb.SupRB( + n_iter=4, + rule_discovery=ES1xLambda(n_iter=4, lmbda=1, delay=2), + solution_composition=NonDominatedSortingGeneticAlgorithm2(n_iter=2, population_size=2), + logger=suprb.logging.stdout.StdoutLogger(), + verbose=10, + ) + + X, y = _regression_dataset() + estimator.fit(X, y) + + # check_estimator(estimator) + + def test_check_nsga3(self): + estimator = suprb.SupRB( + n_iter=4, + rule_discovery=ES1xLambda(n_iter=4, lmbda=1, delay=2), + solution_composition=NonDominatedSortingGeneticAlgorithm3(n_iter=2, population_size=2), + logger=suprb.logging.stdout.StdoutLogger(), + verbose=10, + ) + + X, y = _regression_dataset() + estimator.fit(X, y) + + # check_estimator(estimator) + + def test_check_spea2(self): + estimator = suprb.SupRB( + n_iter=4, + rule_discovery=ES1xLambda(n_iter=4, lmbda=1, delay=2), + solution_composition=StrengthParetoEvolutionaryAlgorithm2(n_iter=2, population_size=2), + logger=suprb.logging.stdout.StdoutLogger(), + verbose=10, + ) + + X, y = _regression_dataset() + estimator.fit(X, y) + + # check_estimator(estimator) + + def test_check_ts(self): + estimator = suprb.SupRB( + n_iter=4, + rule_discovery=ES1xLambda(n_iter=4, lmbda=1, delay=2), + solution_composition=TwoStageSolutionComposition( + algorithm_1=GeneticAlgorithm(n_iter=2, population_size=2), + algorithm_2=NonDominatedSortingGeneticAlgorithm3(n_iter=2, population_size=2), + switch_iteration=3, + warm_start=False, + ), + logger=suprb.logging.stdout.StdoutLogger(), + verbose=10, + ) + + X, y = _regression_dataset() + estimator.fit(X, y) + + # check_estimator(estimator)