Repository: GheodeAI/metaheuristic-designer (language: python)

Issue:
Seeding solutions should accept a list of vectors instead of Individuals



Discussion / hints:

Fix patch:
diff --git a/docs/conf.py b/docs/conf.py
index a447185..4f17938 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -24,7 +24,7 @@
 author = "Eugenio Lorente-Ramos"
 
 # The full version, including alpha/beta/rc tags
-release = "0.1.3"
+release = "0.1.5"
 
 
 # -- General configuration ---------------------------------------------------
diff --git a/docs/metaheuristic_designer.algorithms.rst b/docs/metaheuristic_designer.algorithms.rst
index 8404a6a..3ac0738 100644
--- a/docs/metaheuristic_designer.algorithms.rst
+++ b/docs/metaheuristic_designer.algorithms.rst
@@ -16,3 +16,19 @@ Memetic Algorithm
    :members:
    :undoc-members:
    :show-inheritance:
+
+Algorithm Selection
+---------------------------------------------
+
+.. autoclass:: metaheuristic_designer.algorithms.AlgorithmSelection
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+Strategy Selection
+---------------------------------------------
+
+.. autoclass:: metaheuristic_designer.algorithms.StrategySelection
+   :members:
+   :undoc-members:
+   :show-inheritance:
diff --git a/docs/metaheuristic_designer.encodings.rst b/docs/metaheuristic_designer.encodings.rst
index b524f45..89256ab 100644
--- a/docs/metaheuristic_designer.encodings.rst
+++ b/docs/metaheuristic_designer.encodings.rst
@@ -33,6 +33,14 @@ ImageEncoding
    :undoc-members:
    :show-inheritance:
 
+AdaptionEncoding
+----------------
+
+.. autoclass:: metaheuristic_designer.encodings.AdaptionEncoding
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 CMAEncoding
 -----------
 
diff --git a/docs/metaheuristic_designer.operators.rst b/docs/metaheuristic_designer.operators.rst
index 0a89ade..7802f66 100644
--- a/docs/metaheuristic_designer.operators.rst
+++ b/docs/metaheuristic_designer.operators.rst
@@ -49,6 +49,23 @@ OperatorMeta
    :undoc-members:
    :show-inheritance:
 
+OperatorAdaptative
+----------------------------------------
+
+.. automodule:: metaheuristic_designer.operators.OperatorAdaptative
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+
+OperatorNull
+----------------------------------------
+
+.. automodule:: metaheuristic_designer.operators.OperatorNull
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 OperatorLambda
 ----------------------------------------
 
diff --git a/examples/algorithm_seleciton.py b/examples/algorithm_seleciton.py
new file mode 100644
index 0000000..8747fd7
--- /dev/null
+++ b/examples/algorithm_seleciton.py
@@ -0,0 +1,68 @@
+from metaheuristic_designer import ObjectiveFunc, ParamScheduler, simple
+from metaheuristic_designer.algorithms import GeneralAlgorithm, MemeticAlgorithm, AlgorithmSelection
+from metaheuristic_designer.operators import OperatorReal, OperatorInt, OperatorBinary
+from metaheuristic_designer.initializers import UniformVectorInitializer
+from metaheuristic_designer.selectionMethods import ParentSelection, SurvivorSelection
+from metaheuristic_designer.strategies import *
+
+from metaheuristic_designer.benchmarks import *
+
+import argparse
+
+from copy import copy
+import scipy as sp
+import numpy as np
+
+
+def run_algorithm(save_report):
+    params = {
+        "stop_cond": "neval",
+        "neval": 1e4,
+        "encoding": "real",
+        "verbose": False,
+    }
+
+    objfunc = HappyCat(3, "min")
+
+    # Define algorithms to be tested
+    algorithms = [
+        simple.hill_climb(objfunc, params),
+        simple.simulated_annealing(objfunc, params),
+        simple.evolution_strategy(objfunc, params),
+        simple.differential_evolution(objfunc, params),
+        simple.genetic_algorithm(objfunc, params),
+        simple.particle_swarm(objfunc, params),
+        simple.random_search(objfunc, params),
+    ]
+
+    algorithm_search = AlgorithmSelection(algorithms)
+
+    solution, best_fitness, report = algorithm_search.optimize()
+    print(f"solution: {solution}")
+    print(f"with fitness: {best_fitness}")
+    print(report)
+    if save_report:
+        report.to_csv("./examples/results/algorithm_selection_report.csv")
+
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "-s",
+        "--save-report",
+        dest="save_report",
+        action="store_true",
+        help="Saves the state of the search strategy",
+    )
+    args = parser.parse_args()
+
+    save_report = False
+
+    if args.save_report:
+        save_report = True
+
+    run_algorithm(save_report=save_report)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/examples/exec_basic.py b/examples/exec_basic.py
index d130d7b..a363dca 100644
--- a/examples/exec_basic.py
+++ b/examples/exec_basic.py
@@ -34,9 +34,7 @@ def run_algorithm(alg_name, memetic, save_state):
     # objfunc = Rastrigin(30, "min")
     objfunc = Rosenbrock(30, "min")
     # objfunc = Weierstrass(30, "min")
-    pop_initializer = UniformVectorInitializer(
-        objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=100
-    )
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=100)
 
     parent_params = ParamScheduler("Linear", {"amount": 20})
     # select_params = ParamScheduler("Linear")
@@ -54,10 +52,7 @@ def run_algorithm(alg_name, memetic, save_state):
         OperatorReal("DE/current-to-rand/1", DEparams),
     ]
 
-    neighborhood_structures = [
-        OperatorReal("Gauss", {"F": f}, name=f"Gauss(s={f:0.5e})")
-        for f in np.logspace(-6, 0, base=10, num=80)
-    ]
+    neighborhood_structures = [OperatorReal("Gauss", {"F": f}, name=f"Gauss(s={f:0.5e})") for f in np.logspace(-6, 0, base=10, num=80)]
 
     parent_sel_op = ParentSelection("Best", parent_params)
     selection_op = SurvivorSelection("(m+n)")
@@ -74,9 +69,7 @@ def run_algorithm(alg_name, memetic, save_state):
         search_strat = LocalSearch(pop_initializer, mutation_op, params={"iters": 20})
     elif alg_name == "SA":
         pop_initializer.pop_size = 1
-        search_strat = SA(
-            pop_initializer, mutation_op, {"iter": 100, "temp_init": 1, "alpha": 0.997}
-        )
+        search_strat = SA(pop_initializer, mutation_op, {"iter": 100, "temp_init": 1, "alpha": 0.997})
     elif alg_name == "ES":
         search_strat = ES(
             pop_initializer,
@@ -98,9 +91,7 @@ def run_algorithm(alg_name, memetic, save_state):
     elif alg_name == "HS":
         search_strat = HS(pop_initializer, {"HMCR": 0.8, "BW": 0.5, "PAR": 0.2})
     elif alg_name == "DE":
-        search_strat = DE(
-            pop_initializer, OperatorReal("DE/best/1", {"F": 0.8, "Cr": 0.8})
-        )
+        search_strat = DE(pop_initializer, OperatorReal("DE/best/1", {"F": 0.8, "Cr": 0.8}))
     elif alg_name == "PSO":
         search_strat = PSO(pop_initializer, {"w": 0.7, "c1": 1.5, "c2": 1.5})
     elif alg_name == "CRO":
@@ -151,9 +142,7 @@ def run_algorithm(alg_name, memetic, save_state):
         # local_search = HillClimb(pop_initializer)
         # search_strat = VNS(pop_initializer, neighborhood_structures, local_search, params={"iters": 500})
     elif alg_name == "GVNS":
-        local_search = VND(
-            pop_initializer, neighborhood_structures, params={"nchange": "cyclic"}
-        )
+        local_search = VND(pop_initializer, neighborhood_structures, params={"nchange": "cyclic"})
         search_strat = VNS(
             pop_initializer,
             neighborhood_structures,
@@ -171,9 +160,7 @@ def run_algorithm(alg_name, memetic, save_state):
         exit()
 
     if memetic:
-        alg = MemeticAlgorithm(
-            objfunc, search_strat, local_search, mem_select, params=params
-        )
+        alg = MemeticAlgorithm(objfunc, search_strat, local_search, mem_select, params=params)
     else:
         alg = GeneralAlgorithm(objfunc, search_strat, params=params)
 
diff --git a/examples/image_deblur.py b/examples/image_deblur.py
index 2ca3618..9a5916d 100644
--- a/examples/image_deblur.py
+++ b/examples/image_deblur.py
@@ -102,9 +102,7 @@ def run_algorithm(alg_name, img_file_name, memetic):
         )
         for i in range(100)
     ]
-    pop_initializer = DirectInitializer(
-        pop_initializer, init_population, encoding=encoding
-    )
+    pop_initializer = DirectInitializer(pop_initializer, init_population, encoding=encoding)
 
     mutation_op = OperatorReal("MutRand", {"method": "Cauchy", "F": 4, "N": 2})
     cross_op = OperatorReal("Multicross", {"Nindiv": 4})
@@ -123,9 +121,7 @@ def run_algorithm(alg_name, img_file_name, memetic):
         search_strat = LocalSearch(pop_initializer, mutation_op, {"iters": 20})
     elif alg_name == "SA":
         pop_initializer.pop_size = 1
-        search_strat = SA(
-            pop_initializer, mutation_op, {"iter": 100, "temp_init": 2, "alpha": 0.998}
-        )
+        search_strat = SA(pop_initializer, mutation_op, {"iter": 100, "temp_init": 2, "alpha": 0.998})
     elif alg_name == "ES":
         search_strat = ES(
             pop_initializer,
@@ -147,9 +143,7 @@ def run_algorithm(alg_name, img_file_name, memetic):
     elif alg_name == "HS":
         search_strat = HS(pop_initializer, {"HMCR": 0.8, "BW": 0.5, "PAR": 0.2})
     elif alg_name == "DE":
-        search_strat = DE(
-            pop_initializer, OperatorReal("DE/best/1", {"F": 0.8, "Cr": 0.8})
-        )
+        search_strat = DE(pop_initializer, OperatorReal("DE/best/1", {"F": 0.8, "Cr": 0.8}))
     elif alg_name == "PSO":
         search_strat = PSO(pop_initializer, {"w": 0.7, "c1": 1.5, "c2": 1.5})
     elif alg_name == "RandomSearch":
@@ -233,20 +227,14 @@ def run_algorithm(alg_name, img_file_name, memetic):
         # full_texture = image_blur
         render(full_texture, display_dim, src)
     alg.display_report(show_plots=True)
-    save_to_image(
-        image, f"{img_name}_{image_shape[0]}x{image_shape[1]}_deblured_{alg_name}.png"
-    )
+    save_to_image(image, f"{img_name}_{image_shape[0]}x{image_shape[1]}_deblured_{alg_name}.png")
 
 
 def main():
     parser = argparse.ArgumentParser()
     parser.add_argument("-a", "--algorithm", dest="alg", help="Specify an algorithm")
-    parser.add_argument(
-        "-i", "--image", dest="img", help="Specify an image as reference"
-    )
-    parser.add_argument(
-        "-m", "--memetic", dest="mem", action="store_true", help="Specify an algorithm"
-    )
+    parser.add_argument("-i", "--image", dest="img", help="Specify an image as reference")
+    parser.add_argument("-m", "--memetic", dest="mem", action="store_true", help="Specify an algorithm")
     args = parser.parse_args()
 
     algorithm_name = "SA"
diff --git a/examples/image_evolution.py b/examples/image_evolution.py
index 6c4566b..3eec555 100644
--- a/examples/image_evolution.py
+++ b/examples/image_evolution.py
@@ -23,9 +23,7 @@
 
 
 def render(image, display_dim, src):
-    texture = cv2.resize(
-        image.transpose([1, 0, 2]), display_dim, interpolation=cv2.INTER_NEAREST
-    )
+    texture = cv2.resize(image.transpose([1, 0, 2]), display_dim, interpolation=cv2.INTER_NEAREST)
     pygame.surfarray.blit_array(src, texture)
     pygame.display.flip()
 
@@ -50,6 +48,9 @@ def run_algorithm(alg_name, img_file_name, memetic):
         "patience": 200,
         "verbose": True,
         "v_timer": 0.5,
+        # Parallel
+        "parallel": False,
+        "threads": 8,
     }
 
     display = True
@@ -65,9 +66,9 @@ def run_algorithm(alg_name, img_file_name, memetic):
     img_name = img_file_name.split("/")[-1]
     img_name = img_name.split(".")[0]
 
-    # objfunc = ImgApprox(image_shape, reference_img, img_name=img_name)
+    objfunc = ImgApprox(image_shape, reference_img, img_name=img_name)
     # objfunc = ImgEntropy(image_shape, 256)
-    objfunc = ImgExperimental(image_shape, reference_img, img_name=img_name)
+    # objfunc = ImgExperimental(image_shape, reference_img, img_name=img_name)
 
     encoding = ImageEncoding(image_shape, color=True)
     pop_initializer = UniformVectorInitializer(
@@ -78,16 +79,12 @@ def run_algorithm(alg_name, img_file_name, memetic):
         encoding=encoding,
     )
 
-    mutation_op = OperatorReal(
-        "MutRand", {"method": "Uniform", "Low": -20, "Up": 20, "N": 15}
-    )
+    mutation_op = OperatorReal("MutRand", {"method": "Uniform", "Low": -20, "Up": 20, "N": 15})
     cross_op = OperatorReal("Multipoint")
 
     op_list = [
         OperatorReal("Multipoint"),
-        OperatorReal(
-            "MutRand", {"method": "Cauchy", "F": 5, "N": 10}, name="MutCauchy"
-        ),
+        OperatorReal("MutRand", {"method": "Cauchy", "F": 5, "N": 10}, name="MutCauchy"),
         OperatorReal("MutRand", {"method": "Gauss", "F": 5, "N": 10}, name="MutGauss"),
         OperatorReal(
             "MutSample",
@@ -113,7 +110,7 @@ def run_algorithm(alg_name, img_file_name, memetic):
     selection_op = SurvivorSelection("Elitism", {"amount": 10})
 
     mem_select = ParentSelection("Best", {"amount": 5})
-    neihbourhood_op = OperatorInt("MutRand", {"method": "Uniform", "Low": -3, "Up": -3, "N": 3})
+    neihbourhood_op = OperatorInt("MutRand", {"method": "Uniform", "Low": -10, "Up": -10, "N": 3})
     local_search = LocalSearch(pop_initializer, neihbourhood_op, params={"iters": 10})
 
     if alg_name == "HillClimb":
@@ -121,12 +118,10 @@ def run_algorithm(alg_name, img_file_name, memetic):
         search_strat = HillClimb(pop_initializer, mutation_op)
     elif alg_name == "LocalSearch":
         pop_initializer.pop_size = 1
-        search_strat = LocalSearch(pop_initializer, mutation_op, params={"iters": 20})
+        search_strat = LocalSearch(pop_initializer, mutation_op, params={"iters": 40})
     elif alg_name == "SA":
         pop_initializer.pop_size = 1
-        search_strat = SA(
-            pop_initializer, mutation_op, {"iter": 100, "temp_init": 1, "alpha": 0.997}
-        )
+        search_strat = SA(pop_initializer, mutation_op, {"iter": 100, "temp_init": 1, "alpha": 0.997})
     elif alg_name == "ES":
         search_strat = ES(
             pop_initializer,
@@ -148,9 +143,7 @@ def run_algorithm(alg_name, img_file_name, memetic):
     elif alg_name == "HS":
         search_strat = HS(pop_initializer, {"HMCR": 0.8, "BW": 0.5, "PAR": 0.2})
     elif alg_name == "DE":
-        search_strat = DE(
-            pop_initializer, OperatorReal("DE/best/1", {"F": 0.8, "Cr": 0.8})
-        )
+        search_strat = DE(pop_initializer, OperatorReal("DE/best/1", {"F": 0.8, "Cr": 0.8}))
     elif alg_name == "PSO":
         search_strat = PSO(pop_initializer, {"w": 0.7, "c1": 1.5, "c2": 1.5})
     elif alg_name == "CRO":
@@ -211,9 +204,7 @@ def run_algorithm(alg_name, img_file_name, memetic):
         exit()
 
     if memetic:
-        alg = MemeticAlgorithm(
-            objfunc, search_strat, local_search, mem_select, params=params
-        )
+        alg = MemeticAlgorithm(objfunc, search_strat, local_search, mem_select, params=params)
     else:
         alg = GeneralAlgorithm(objfunc, search_strat, params)
 
@@ -258,12 +249,8 @@ def run_algorithm(alg_name, img_file_name, memetic):
 def main():
     parser = argparse.ArgumentParser()
     parser.add_argument("-a", "--algorithm", dest="alg", help="Specify an algorithm")
-    parser.add_argument(
-        "-i", "--image", dest="img", help="Specify an image as reference"
-    )
-    parser.add_argument(
-        "-m", "--memetic", dest="mem", action="store_true", help="Specify an algorithm"
-    )
+    parser.add_argument("-i", "--image", dest="img", help="Specify an image as reference")
+    parser.add_argument("-m", "--memetic", dest="mem", action="store_true", help="Specify an algorithm")
     args = parser.parse_args()
 
     algorithm_name = "SA"
diff --git a/examples/np_problem_example.py b/examples/np_problem_example.py
index f42aaf8..2125a01 100644
--- a/examples/np_problem_example.py
+++ b/examples/np_problem_example.py
@@ -151,9 +151,7 @@ def run_algorithm(alg_name, problem_name, memetic, save_state):
         objfunc = BinKnapsack(weights, values, capacity)
 
         encoding = TypeCastEncoding(int, bool)
-        pop_initializer = UniformVectorInitializer(
-            objfunc.vecsize, 0, 1, pop_size=100, dtype=int, encoding=encoding
-        )
+        pop_initializer = UniformVectorInitializer(objfunc.vecsize, 0, 1, pop_size=100, dtype=int, encoding=encoding)
 
         parent_params = ParamScheduler("Linear", {"amount": 20})
         # select_params = ParamScheduler("Linear")
@@ -170,18 +168,13 @@ def run_algorithm(alg_name, problem_name, memetic, save_state):
             OperatorBinary("MutSample", {"method": "Bernoulli", "p": 0.2, "N": 2}),
         ]
 
-        neighborhood_structures = [
-            OperatorBinary("Flip", {"N": n}, name=f"Flip(n={n})")
-            for n in range(objfunc.vecsize)
-        ]
+        neighborhood_structures = [OperatorBinary("Flip", {"N": n}, name=f"Flip(n={n})") for n in range(objfunc.vecsize)]
 
     elif problem_name == "3SAT":
         objfunc = ThreeSAT.from_cnf_file("./data/sat_examples/uf50-03.cnf")
 
         encoding = TypeCastEncoding(int, bool)
-        pop_initializer = UniformVectorInitializer(
-            objfunc.vecsize, 0, 1, pop_size=100, dtype=int, encoding=encoding
-        )
+        pop_initializer = UniformVectorInitializer(objfunc.vecsize, 0, 1, pop_size=100, dtype=int, encoding=encoding)
 
         parent_params = ParamScheduler("Linear", {"amount": 20})
         # select_params = ParamScheduler("Linear")
@@ -198,10 +191,7 @@ def run_algorithm(alg_name, problem_name, memetic, save_state):
             OperatorBinary("MutSample", {"method": "Bernoulli", "p": 0.2, "N": 2}),
         ]
 
-        neighborhood_structures = [
-            OperatorBinary("Flip", {"N": n}, name=f"Flip(n={n})")
-            for n in range(objfunc.vecsize)
-        ]
+        neighborhood_structures = [OperatorBinary("Flip", {"N": n}, name=f"Flip(n={n})") for n in range(objfunc.vecsize)]
     elif problem_name == "MaxClique":
         g = nx.gnp_random_graph(100, 0.8)
         adj_mat = nx.adjacency_matrix(g).todense()
@@ -229,10 +219,7 @@ def run_algorithm(alg_name, problem_name, memetic, save_state):
             OperatorPerm("PMX"),
         ]
 
-        neighborhood_structures = [
-            OperatorPerm("Flip", {"N": n}, name=f"Flip(n={n})")
-            for n in range(objfunc.vecsize)
-        ]
+        neighborhood_structures = [OperatorPerm("Flip", {"N": n}, name=f"Flip(n={n})") for n in range(objfunc.vecsize)]
     else:
         raise ValueError(f"The problem '{problem_name}' does not exist.")
 
@@ -277,9 +264,7 @@ def run_algorithm(alg_name, problem_name, memetic, save_state):
     elif alg_name == "HS":
         search_strat = HS(pop_initializer, params={"HMCR": 0.8, "BW": 0.5, "PAR": 0.2})
     elif alg_name == "DE":
-        search_strat = DE(
-            pop_initializer, OperatorReal("DE/best/1", params={"F": 0.8, "Cr": 0.8})
-        )
+        search_strat = DE(pop_initializer, OperatorReal("DE/best/1", params={"F": 0.8, "Cr": 0.8}))
     elif alg_name == "PSO":
         search_strat = PSO(pop_initializer, params={"w": 0.7, "c1": 1.5, "c2": 1.5})
     elif alg_name == "CRO":
@@ -330,9 +315,7 @@ def run_algorithm(alg_name, problem_name, memetic, save_state):
         # local_search = HillClimb(pop_initializer)
         # search_strat = VNS(pop_initializer, neighborhood_structures, local_search, params={"iters": 500})
     elif alg_name == "GVNS":
-        local_search = VND(
-            pop_initializer, neighborhood_structures, params={"nchange": "cyclic"}
-        )
+        local_search = VND(pop_initializer, neighborhood_structures, params={"nchange": "cyclic"})
         search_strat = VNS(
             pop_initializer,
             neighborhood_structures,
@@ -348,9 +331,7 @@ def run_algorithm(alg_name, problem_name, memetic, save_state):
         exit()
 
     if memetic:
-        alg = MemeticAlgorithm(
-            objfunc, search_strat, local_search, mem_select, params=params
-        )
+        alg = MemeticAlgorithm(objfunc, search_strat, local_search, mem_select, params=params)
     else:
         alg = GeneralAlgorithm(objfunc, search_strat, params=params)
 
@@ -365,9 +346,7 @@ def run_algorithm(alg_name, problem_name, memetic, save_state):
 def main():
     parser = argparse.ArgumentParser()
     parser.add_argument("-a", "--algorithm", dest="alg", help="Specify an algorithm")
-    parser.add_argument(
-        "-p", "--problem", dest="prob", help="Specify an problem to solve"
-    )
+    parser.add_argument("-p", "--problem", dest="prob", help="Specify an problem to solve")
     parser.add_argument(
         "-m",
         "--memetic",
diff --git a/examples/obj_from_lambda.py b/examples/obj_from_lambda.py
index 460868a..9167946 100644
--- a/examples/obj_from_lambda.py
+++ b/examples/obj_from_lambda.py
@@ -25,9 +25,7 @@ def run_algorithm():
 
     objfunc = ObjectiveFromLambda(func, 10, "min")
 
-    pop_init = UniformVectorInitializer(
-        objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=1
-    )
+    pop_init = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=1)
 
     mutation_op = OperatorReal("RandNoise", {"method": "Cauchy", "F": 0.0001})
 
diff --git a/examples/strategy_seleciton.py b/examples/strategy_seleciton.py
new file mode 100644
index 0000000..4dcc9c6
--- /dev/null
+++ b/examples/strategy_seleciton.py
@@ -0,0 +1,119 @@
+from metaheuristic_designer import ObjectiveFunc, ParamScheduler, simple
+from metaheuristic_designer.algorithms import GeneralAlgorithm, MemeticAlgorithm, StrategySelection
+from metaheuristic_designer.operators import OperatorReal, OperatorInt, OperatorBinary
+from metaheuristic_designer.initializers import UniformVectorInitializer
+from metaheuristic_designer.selectionMethods import ParentSelection, SurvivorSelection
+from metaheuristic_designer.strategies import *
+
+from metaheuristic_designer.benchmarks import *
+
+import argparse
+
+from copy import copy
+import scipy as sp
+import numpy as np
+
+
+def run_algorithm(save_report):
+    objfunc = HappyCat(3, "min")
+    single_initializer = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=1)
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=100)
+
+    # Define algorithms to be tested
+    strategies = [
+        HillClimb(single_initializer, OperatorReal("RandNoise", {"method": "Gauss", "F": 1e-4}), name="HillClimb-Gauss"),
+        HillClimb(single_initializer, OperatorReal("RandNoise", {"method": "Cauchy", "F": 1e-4}), name="HillClimb-Cauchy"),
+        LocalSearch(single_initializer, OperatorReal("RandNoise", {"method": "Cauchy", "F": 1e-4}), params={"iters": 20}, name="LocalSearch-Cauchy"),
+        LocalSearch(single_initializer, OperatorReal("RandNoise", {"method": "Gauss", "F": 1e-4}), params={"iters": 20}, name="LocalSearch-Gauss"),
+        SA(
+            single_initializer,
+            OperatorReal("RandNoise", {"method": "Gauss", "F": 1e-4}),
+            params={"iter": 100, "temp_init": 1, "alpha": 0.997},
+            name="SA-Gauss",
+        ),
+        SA(
+            single_initializer,
+            OperatorReal("RandNoise", {"method": "Cauchy", "F": 1e-4}),
+            params={"iter": 100, "temp_init": 1, "alpha": 0.997},
+            name="SA-Cauchy",
+        ),
+        SA(
+            pop_initializer,
+            OperatorReal("RandNoise", {"method": "Gauss", "F": 1e-4}),
+            params={"iter": 100, "temp_init": 1, "alpha": 0.997},
+            name="ParallelSA-Gauss",
+        ),
+        ES(
+            pop_initializer,
+            OperatorReal("RandNoise", {"method": "Gauss", "F": 1e-4}),
+            OperatorReal("Nothing"),
+            ParentSelection("Nothing"),
+            SurvivorSelection("(m+n)"),
+            params={"offspringSize": 150},
+            name="ES-(100+150)",
+        ),
+        ES(
+            pop_initializer,
+            OperatorReal("RandNoise", {"method": "Gauss", "F": 1e-4}),
+            OperatorReal("Nothing"),
+            ParentSelection("Nothing"),
+            SurvivorSelection("(m,n)"),
+            params={"offspringSize": 400},
+            name="ES-(100,400)",
+        ),
+        GA(
+            pop_initializer,
+            OperatorReal("RandNoise", {"method": "Gauss", "F": 1e-4}),
+            OperatorReal("Multipoint"),
+            ParentSelection("Tournament", {"amount": 60, "p": 0.1}),
+            SurvivorSelection("Elitism", {"amount": 10}),
+            params={"pcross": 0.8, "pmut": 0.1},
+            name="GA",
+        ),
+        PSO(pop_initializer, {"w": 0.7, "c1": 1.5, "c2": 1.5}, name="PSO"),
+        DE(pop_initializer, OperatorReal("DE/best/1", {"F": 0.8, "Cr": 0.8}), name="DE/best/1"),
+        DE(pop_initializer, OperatorReal("DE/rand/1", {"F": 0.8, "Cr": 0.8}), name="DE/rand/1"),
+        DE(pop_initializer, OperatorReal("DE/current-to-best/1", {"F": 0.8, "Cr": 0.8}), name="DE/current-to-best/1"),
+        RandomSearch(pop_initializer),
+    ]
+
+    algorithm_search = StrategySelection(
+        objfunc,
+        strategies,
+        algorithm_params={
+            "stop_cond": "neval",
+            "neval": 1e4,
+            "verbose": False,
+        },
+        params={"verbose": False, "repetitions": 10},
+    )
+
+    solution, best_fitness, report = algorithm_search.optimize()
+    print(f"solution: {solution}")
+    print(f"with fitness: {best_fitness}")
+    print(report)
+    if save_report:
+        report.to_csv("./examples/results/strategy_selection_report.csv")
+
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument(
+        "-s",
+        "--save-report",
+        dest="save_report",
+        action="store_true",
+        help="Saves the state of the search strategy",
+    )
+    args = parser.parse_args()
+
+    save_report = False
+
+    if args.save_report:
+        save_report = True
+
+    run_algorithm(save_report=save_report)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/noxfile.py b/noxfile.py
index 01b01b1..e2b92ae 100644
--- a/noxfile.py
+++ b/noxfile.py
@@ -8,6 +8,4 @@
 def test(session):
     session.install(".")
     session.install("pytest", "coverage")
-    session.run(
-        "coverage", "run", "--data-file", ".coverage.nox", "--parallel", "-m", "pytest"
-    )
+    session.run("coverage", "run", "--data-file", ".coverage.nox", "--parallel", "-m", "pytest")
diff --git a/pyproject.toml b/pyproject.toml
index f3e7313..61f71d1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "metaheuristic-designer"
-version = "0.1.4"
+version = "0.1.5"
 authors = [
     {name = "eugenioLR", email = "eugeniolorente0@gmail.com"}
 ]
@@ -31,9 +31,11 @@ dev = [
     "coverage>=7.3.1",
     "sphinx_material>=0.0.36",
     "sphinx>=7.1.1",
+    "black[jupyter]>=23.10.1",
     "pygame>=2.5.0",    
     "networkx>=3.1",
-    "opencv-python>=4.8"
+    "opencv-python>=4.8",
+    "pandas>=2.1.2"
 ]
 
 examples = [
diff --git a/src/metaheuristic_designer/Algorithm.py b/src/metaheuristic_designer/Algorithm.py
index 41fecc5..4820cc8 100644
--- a/src/metaheuristic_designer/Algorithm.py
+++ b/src/metaheuristic_designer/Algorithm.py
@@ -5,11 +5,14 @@
 import numpy as np
 import pyparsing as pp
 from .utils import NumpyEncoder
+import matplotlib.pyplot as plt
 
 
 class Algorithm(ABC):
     """
-    General framework for metaheuristic algorithms.
+    Abstract Algorithm class.
+
+    This class defines the structure of all optimization algorithms.
 
     Parameters
     ----------
@@ -22,12 +25,7 @@ class Algorithm(ABC):
         Dictionary of parameters to define the stopping condition and output of the algorithm.
     """
 
-    def __init__(
-        self,
-        objfunc: ObjectiveFunc,
-        search_strategy: SearchStrategy,
-        params: Union[ParamScheduler, dict] = None,
-    ):
+    def __init__(self, objfunc: ObjectiveFunc, search_strategy: SearchStrategy, params: Union[ParamScheduler, dict] = None, name: str = None):
         """
         Constructor of the Search class
         """
@@ -35,6 +33,7 @@ def __init__(
         self.params = params
         self.search_strategy = search_strategy
         self.objfunc = objfunc
+        self._name = name
 
         if params is None:
             params = {}
@@ -48,20 +47,22 @@ def __init__(
         self.stop_cond_parsed = parse_stopping_cond(self.stop_cond)
 
         self.progress_metric = params.get("progress_metric", self.stop_cond)
-        self.progress_metric_parsed = (
-            parse_stopping_cond(self.progress_metric)
-            if "progress_metric" in params
-            else self.stop_cond_parsed
-        )
+        self.progress_metric_parsed = parse_stopping_cond(self.progress_metric) if "progress_metric" in params else self.stop_cond_parsed
 
-        self.Ngen = params.get("ngen", 100)
-        self.Neval = params.get("neval", 1e5)
+        self.ngen = params.get("ngen", 100)
+        self.neval = params.get("neval", 1e5)
         self.time_limit = params.get("time_limit", 10.0)
         self.cpu_time_limit = params.get("cpu_time_limit", 10.0)
         self.fit_target = params.get("fit_target", 1e-10)
+        if self.fit_target == 0:
+            self.fit_target = 1e-10
         self.max_patience = params.get("patience", 1)
         self.patience_left = self.max_patience
 
+        # Parallel parameters
+        self.parallel = params.get("parallel", False)
+        self.threads = params.get("threads", 8)
+
         # Metrics
         self.fit_history = []
         self.best_history = []
@@ -73,6 +74,14 @@ def __init__(
         self.real_time_spent = 0
         self.converged_steps = 0
 
+    @property
+    def name(self):
+        return self._name if self._name else self.search_strategy.name
+
+    @name.setter
+    def name(self, new_name: str):
+        self._name = new_name
+
     def restart(self):
         """
         Resets the internal values of the algorithm and the number of evaluations of the fitness function.
@@ -113,9 +122,7 @@ def best_solution(self) -> Tuple[Individual, float]:
 
         return self.search_strategy.best_solution()
 
-    def stopping_condition(
-        self, gen: int, real_time_start: float, cpu_time_start: float
-    ) -> bool:
+    def stopping_condition(self, gen: int, real_time_start: float, cpu_time_start: float) -> bool:
         """
         Given the state of the algorithm, returns wether we have finished or not.
 
@@ -134,9 +141,9 @@ def stopping_condition(
             Whether the algorithm has reached its end
         """
 
-        neval_reached = self.objfunc.counter >= self.Neval
+        neval_reached = self.objfunc.counter >= self.neval
 
-        ngen_reached = gen >= self.Ngen
+        ngen_reached = gen >= self.ngen
 
         real_time_reached = time.time() - real_time_start >= self.time_limit
 
@@ -159,9 +166,7 @@ def stopping_condition(
             patience_reached,
         )
 
-    def get_progress(
-        self, gen: int, real_time_start: float, cpu_time_start: float
-    ) -> float:
+    def get_progress(self, gen: int, real_time_start: float, cpu_time_start: float) -> float:
         """
         Given the state of the algorithm, returns a number between 0 and 1 indicating
         how close to the end of the algorithm we are, 0 when starting and 1 when finished.
@@ -181,21 +186,28 @@ def get_progress(
             Indicator of how close it the algorithm to finishing, 1 means the algorithm should be stopped.
         """
 
-        neval_reached = self.objfunc.counter / self.Neval
+        neval_reached = self.objfunc.counter / self.neval
 
-        ngen_reached = gen / self.Ngen
+        ngen_reached = gen / self.ngen
 
         real_time_reached = (time.time() - real_time_start) / self.time_limit
 
         cpu_time_reached = (time.process_time() - cpu_time_start) / self.cpu_time_limit
 
         best_fitness = self.best_solution()[1]
+
+        # if self.objfunc.mode == "max":
+        #     if self.fit_target == 0:
+        #         self.fit_target = 1e-40
+        #     target_reached = (best_fitness / self.fit_target
+        # else:
+        #     if best_fitness == 0:
+        #         best_fitness = 1e-40
+        #     target_reached = self.fit_target / best_fitness
         if self.objfunc.mode == "max":
-            target_reached = best_fitness / self.fit_target
+            target_reached = 1 - (self.best_solution()[1] - self.fit_target) / self.fit_target
         else:
-            if best_fitness == 0:
-                best_fitness = 1e-40
-            target_reached = self.fit_target / best_fitness
+            target_reached = 1 - (self.fit_target - self.best_solution()[1]) / self.fit_target
 
         patience_prec = 1 - self.patience_left / self.max_patience
 
@@ -209,9 +221,7 @@ def get_progress(
             patience_prec,
         )
 
-    def update(
-        self, real_time_start: float, cpu_time_start: float, pass_step: bool = True
-    ):
+    def update(self, real_time_start: float, cpu_time_start: float, pass_step: bool = True):
         """
         Updates the attributes of the optimization algorithm.
         This function should be called once per iteration of the algorithm.
@@ -238,9 +248,7 @@ def update(
 
         self.progress = self.get_progress(self.steps, real_time_start, cpu_time_start)
 
-        self.ended = self.stopping_condition(
-            self.steps, real_time_start, cpu_time_start
-        )
+        self.ended = self.stopping_condition(self.steps, real_time_start, cpu_time_start)
 
     def initialize(self):
         """
@@ -248,12 +256,12 @@ def initialize(self):
         """
 
         self.restart()
-        self.search_strategy.initialize(self.objfunc)
+        initial_population = self.search_strategy.initialize(self.objfunc)
+        initial_population = self.search_strategy.evaluate_population(initial_population, self.objfunc, self.parallel, self.threads)
+        self.search_strategy.population = initial_population
 
     @abstractmethod
-    def step(
-        self, time_start: float = 0, verbose: bool = False
-    ) -> Tuple[Individual, float]:
+    def step(self, time_start: float = 0, verbose: bool = False) -> Tuple[Individual, float]:
         """
         Performs an iteration of the algorithm.
 
@@ -359,9 +367,7 @@ def get_state(
 
         if show_best_solution:
             data["best_fitness"] = self.best_solution()[1]
-            data["best_individual"] = self.search_strategy.best.get_state(
-                show_speed=False, show_best=False
-            )
+            data["best_individual"] = self.search_strategy.best.get_state(show_speed=False, show_best=False)
 
         if show_fit_history:
             data["fit_history"] = self.fit_history
@@ -369,9 +375,7 @@ def get_state(
         if show_gen_history:
             data["best_history"] = self.best_history
 
-        data["search_strat_state"] = self.search_strategy.get_state(
-            show_pop, show_pop_details
-        )
+        data["search_strat_state"] = self.search_strategy.get_state(show_pop, show_pop_details)
 
         return data
 
@@ -422,15 +426,10 @@ def store_state(
             fp.write(dumped)
 
     def init_info(self):
-        print(
-            f"Initializing optimization of {self.objfunc.name} using {self.search_strategy.name}"
-        )
-        print(
-            f"-----------------------------{'-'*len(self.objfunc.name)}-------{'-'*len(self.search_strategy.name)}"
-        )
+        print(f"Initializing optimization of {self.objfunc.name} using {self.search_strategy.name}")
+        print(f"-----------------------------{'-'*len(self.objfunc.name)}-------{'-'*len(self.search_strategy.name)}")
         print()
 
-    @abstractmethod
     def step_info(self, start_time: float = 0):
         """
         Displays information about the current state of the algotithm.
@@ -441,7 +440,17 @@ def step_info(self, start_time: float = 0):
             Indicates to the algorihm how much time has already passed.
         """
 
-    @abstractmethod
+        print(f"Optimizing {self.objfunc.name} using {self.name}:")
+        print(f"\tReal time Spent: {round(time.time() - start_time,2)} s")
+        print(f"\tCPU time Spent:  {round(time.time() - start_time,2)} s")
+        print(f"\tGeneration: {self.steps}")
+        best_fitness = self.best_solution()[1]
+        print(f"\tBest fitness: {best_fitness}")
+        print(f"\tEvaluations of fitness: {self.objfunc.counter}")
+        print()
+        self.search_strategy.extra_step_info()
+        print()
+
     def display_report(self, show_plots: bool = True):
         """
         Shows a summary of the execution of the algorithm.
@@ -452,6 +461,26 @@ def display_report(self, show_plots: bool = True):
             Whether to display plots about the algorithm or not.
         """
 
+        print("Number of generations:", len(self.fit_history))
+        print("Real time spent: ", round(self.real_time_spent, 5), "s", sep="")
+        print("CPU time spent: ", round(self.cpu_time_spent, 5), "s", sep="")
+        print("Number of fitness evaluations:", self.objfunc.counter)
+
+        best_fitness = self.best_solution()[1]
+        print("Best fitness:", best_fitness)
+
+        if show_plots:
+            # Plot fitness history
+            plt.axhline(y=0, color="black", alpha=0.9)
+            plt.axvline(x=0, color="black", alpha=0.9)
+            plt.plot(self.fit_history, "blue")
+            plt.xlabel("generations")
+            plt.ylabel("fitness")
+            plt.title(f"{self.search_strategy.name} fitness")
+            plt.show()
+
+        self.search_strategy.extra_report(show_plots)
+
 
 def parse_stopping_cond(condition_str: str) -> List[str | List]:
     """
@@ -471,13 +500,9 @@ def parse_stopping_cond(condition_str: str) -> List[str | List]:
 
     orop = pp.Literal("and")
     andop = pp.Literal("or")
-    condition = pp.oneOf(
-        ["neval", "ngen", "time_limit", "cpu_time_limit", "fit_target", "convergence"]
-    )
+    condition = pp.oneOf(["neval", "ngen", "time_limit", "cpu_time_limit", "fit_target", "convergence"])
 
-    expr = pp.infixNotation(
-        condition, [(orop, 2, pp.opAssoc.RIGHT), (andop, 2, pp.opAssoc.RIGHT)]
-    )
+    expr = pp.infixNotation(condition, [(orop, 2, pp.opAssoc.RIGHT), (andop, 2, pp.opAssoc.RIGHT)])
 
     return expr.parse_string(condition_str).as_list()
 
@@ -522,12 +547,8 @@ def process_condition(
 
     if isinstance(cond_parsed, list):
         if len(cond_parsed) == 3:
-            cond1 = process_condition(
-                cond_parsed[0], neval, ngen, real_time, cpu_time, target, patience
-            )
-            cond2 = process_condition(
-                cond_parsed[2], neval, ngen, real_time, cpu_time, target, patience
-            )
+            cond1 = process_condition(cond_parsed[0], neval, ngen, real_time, cpu_time, target, patience)
+            cond2 = process_condition(cond_parsed[2], neval, ngen, real_time, cpu_time, target, patience)
 
             if cond_parsed[1] == "or":
                 result = cond1 or cond2
@@ -535,9 +556,7 @@ def process_condition(
                 result = cond1 and cond2
 
         elif len(cond_parsed) == 1:
-            result = process_condition(
-                cond_parsed[0], neval, ngen, real_time, cpu_time, target, patience
-            )
+            result = process_condition(cond_parsed[0], neval, ngen, real_time, cpu_time, target, patience)
 
     else:
         if cond_parsed == "neval":
@@ -596,12 +615,8 @@ def process_progress(
 
     if isinstance(cond_parsed, list):
         if len(cond_parsed) == 3:
-            progress1 = process_progress(
-                cond_parsed[0], neval, ngen, real_time, cpu_time, target, patience
-            )
-            progress2 = process_progress(
-                cond_parsed[2], neval, ngen, real_time, cpu_time, target, patience
-            )
+            progress1 = process_progress(cond_parsed[0], neval, ngen, real_time, cpu_time, target, patience)
+            progress2 = process_progress(cond_parsed[2], neval, ngen, real_time, cpu_time, target, patience)
 
             if cond_parsed[1] == "or":
                 result = max(progress1, progress2)
@@ -609,9 +624,7 @@ def process_progress(
                 result = min(progress1, progress2)
 
         elif len(cond_parsed) == 1:
-            result = process_progress(
-                cond_parsed[0], neval, ngen, real_time, cpu_time, target, patience
-            )
+            result = process_progress(cond_parsed[0], neval, ngen, real_time, cpu_time, target, patience)
     else:
         if cond_parsed == "neval":
             result = neval
diff --git a/src/metaheuristic_designer/Encoding.py b/src/metaheuristic_designer/Encoding.py
index 4b27a58..1599134 100644
--- a/src/metaheuristic_designer/Encoding.py
+++ b/src/metaheuristic_designer/Encoding.py
@@ -4,10 +4,11 @@
 
 class Encoding(ABC):
     """
-    Base class for transforming between phenotype and genotype.
+    Abstract Encoding class
+
+    This class transforms between phenotype and genotype.
     """
 
-    @abstractmethod
     def encode(self, phenotype) -> Any:
         """
         Encodes a viable solution to our problem to the encoding used in each individual of the algorithm.
diff --git a/src/metaheuristic_designer/Individual.py b/src/metaheuristic_designer/Individual.py
index 4ea7729..b10b184 100644
--- a/src/metaheuristic_designer/Individual.py
+++ b/src/metaheuristic_designer/Individual.py
@@ -55,9 +55,7 @@ def __copy__(self) -> Individual:
         Returns a copy of the Individual.
         """
 
-        copied_ind = Individual(
-            self.objfunc, copy(self._genotype), copy(self.speed), self.encoding
-        )
+        copied_ind = Individual(self.objfunc, copy(self._genotype), copy(self.speed), self.encoding)
         copied_ind._fitness = self._fitness
         copied_ind.fitness_calculated = self.fitness_calculated
         copied_ind.best = copy(self.best)
@@ -90,20 +88,26 @@ def apply_speed(self) -> Individual:
             Individual with the speed applied.
         """
 
-        return Individual(
-            self.objfunc, self._genotype + self.speed, self.speed, self.encoding
-        )
+        return Individual(self.objfunc, self._genotype + self.speed, self.speed, self.encoding)
 
-    @property
-    def fitness(self) -> float:
+    def calculate_fitness(self) -> float:
         """
-        The fitness of the individual, optimized to be calculated only once per individual.
+        Calculates the fitness of the individual if it has not been calculated before
         """
 
         if not self.fitness_calculated:
+            # This calls the setter function
             self.fitness = self.objfunc(self)
         return self._fitness
 
+    @property
+    def fitness(self) -> float:
+        """
+        The fitness of the individual, optimized to be calculated only once per individual.
+        """
+
+        return self.calculate_fitness()
+
     @fitness.setter
     def fitness(self, fit: float):
         """
diff --git a/src/metaheuristic_designer/Initializer.py b/src/metaheuristic_designer/Initializer.py
index 26fb002..67f31ca 100644
--- a/src/metaheuristic_designer/Initializer.py
+++ b/src/metaheuristic_designer/Initializer.py
@@ -58,9 +58,7 @@ def generate_individual(self, objfunc: ObjectiveFunc) -> Individual:
 
         return self.generate_random(objfunc)
 
-    def generate_population(
-        self, objfunc: ObjectiveFunc, n_indiv: int = None
-    ) -> List[Individual]:
+    def generate_population(self, objfunc: ObjectiveFunc, n_indiv: int = None) -> List[Individual]:
         """
         Generate n_indiv Individuals using the generate_individual method.
 
diff --git a/src/metaheuristic_designer/Operator.py b/src/metaheuristic_designer/Operator.py
index 415175c..5f40e36 100644
--- a/src/metaheuristic_designer/Operator.py
+++ b/src/metaheuristic_designer/Operator.py
@@ -7,6 +7,8 @@ class Operator(ABC):
     """
     Abstract Operator class.
 
+    This class modifies the genotype of one individual in order to perform some optimization task.
+
     Parameters
     ----------
     params: ParamScheduler or dict, optional
diff --git a/src/metaheuristic_designer/ParamScheduler.py b/src/metaheuristic_designer/ParamScheduler.py
index 2f850b3..5ddc768 100644
--- a/src/metaheuristic_designer/ParamScheduler.py
+++ b/src/metaheuristic_designer/ParamScheduler.py
@@ -61,6 +61,9 @@ def get(key: str, def_value: Any = None) -> Any:
 
         return self.current_params.get(key, def_value)
 
+    def update(update_dict: dict):
+        self.current_params.update(update_dict)
+
     def reset(self):
         """
         Sets all the parameters to their initial values.
@@ -111,9 +114,7 @@ def step(self, progress: float):
                 if type(self.param_schedule[key]) in (list, tuple):
                     start_param = self.param_schedule[key][0]
                     end_param = self.param_schedule[key][1]
-                    self.current_params[key] = (
-                        1 - progress
-                    ) * start_param + progress * end_param
+                    self.current_params[key] = (1 - progress) * start_param + progress * end_param
 
         elif self.strategy == "exp":
             # with f(x) = k·e^{a·x}+b,  f(0) = p[0],  f(1) = p[1]
diff --git a/src/metaheuristic_designer/SearchStrategy.py b/src/metaheuristic_designer/SearchStrategy.py
index c2560be..af06811 100644
--- a/src/metaheuristic_designer/SearchStrategy.py
+++ b/src/metaheuristic_designer/SearchStrategy.py
@@ -4,10 +4,19 @@
 from .ParamScheduler import ParamScheduler
 from .selectionMethods import ParentSelection, SurvivorSelection
 from .Operator import Operator
+from multiprocessing import Pool
+
+
+def evaluate_indiv(indiv):
+    calculation_done = not indiv.fitness_calculated
+    indiv.calculate_fitness()
+    return indiv, calculation_done
 
 
 class SearchStrategy(ABC):
     """
+    Abstract Search Strategy class.
+    
     This is the class that defines how the optimization will be carried out.
 
     Parameters
@@ -116,11 +125,25 @@ def initialize(self, objfunc: ObjectiveFunc):
 
         self.population = self.pop_init.generate_population(objfunc)
 
-        self.best = max(self.population, key=lambda x: x.fitness)
+        return self.population
+
+    def evaluate_population(self, population, objfunc, parallel=False, threads=8):
+        if parallel:
+            with Pool(threads) as p:
+                result_pairs = p.map(evaluate_indiv, population)
+            population, calculated = map(list, zip(*result_pairs))
+            objfunc.counter += sum(calculated)
+        else:
+            [indiv.calculate_fitness() for indiv in population]
+
+        current_best = max(population, key=lambda x: x.fitness)
+
+        if not self.best or self.best.fitness < current_best.fitness:
+            self.best = current_best
+
+        return population
 
-    def select_parents(
-        self, population: List[Individual], **kwargs
-    ) -> Tuple[List[Individual], List[int]]:
+    def select_parents(self, population: List[Individual], **kwargs) -> Tuple[List[Individual], List[int]]:
         """
         Selects the individuals that will be perturbed in this generation to generate the offspring.
 
@@ -138,9 +161,7 @@ def select_parents(
         return population
 
     @abstractmethod
-    def perturb(
-        self, parent_list: List[Individual], objfunc: ObjectiveFunc, **kwargs
-    ) -> List[Individual]:
+    def perturb(self, parent_list: List[Individual], objfunc: ObjectiveFunc, **kwargs) -> List[Individual]:
         """
         Applies operators to the population to get the next generation of individuals.
 
@@ -157,9 +178,7 @@ def perturb(
             The list of individuals modified by the operators of the search strategy.
         """
 
-    def select_individuals(
-        self, population: List[Individual], offspring: List[Individual], **kwargs
-    ) -> List[Individual]:
+    def select_individuals(self, population: List[Individual], offspring: List[Individual], **kwargs) -> List[Individual]:
         """
         Selects the individuals that will pass to the next generation.
 
@@ -221,10 +240,7 @@ def get_state(self, show_pop: bool = False, show_pop_details: bool = False) -> d
             data["survivor_sel"] = [surv.get_state() for surv in self.surv_sel]
 
         if show_pop:
-            data["population"] = [
-                ind.get_state(show_speed=show_pop_details, show_best=show_pop_details)
-                for ind in self.population
-            ]
+            data["population"] = [ind.get_state(show_speed=show_pop_details, show_best=show_pop_details) for ind in self.population]
 
         return data
 
diff --git a/src/metaheuristic_designer/SelectionMethod.py b/src/metaheuristic_designer/SelectionMethod.py
index ad016f3..b1f3d9b 100644
--- a/src/metaheuristic_designer/SelectionMethod.py
+++ b/src/metaheuristic_designer/SelectionMethod.py
@@ -5,7 +5,9 @@
 
 class SelectionMethod(ABC):
     """
-    Survivor selection methods
+    Abstract Selection Method class.
+
+    This class defines the structure for individual selection methods.
 
     Parameters
     ----------
@@ -38,9 +40,7 @@ def __init__(
         else:
             self.params = params
 
-    def __call__(
-        self, popul: List[Individual], offspring: List[Individual] = None
-    ) -> List[Individual]:
+    def __call__(self, popul: List[Individual], offspring: List[Individual] = None) -> List[Individual]:
         """
         Shorthand for calling the 'select' method
         """
@@ -48,9 +48,7 @@ def __call__(
         return self.select(popul, offspring)
 
     @abstractmethod
-    def select(
-        self, population: List[Individual], offspring: List[Individual] = None
-    ) -> List[Individual]:
+    def select(self, population: List[Individual], offspring: List[Individual] = None) -> List[Individual]:
         """
         Takes a population with its offspring and returns the individuals that survive
         to produce the next generation.
@@ -83,6 +81,9 @@ def step(self, progress: float):
             if "maxPopSize" in self.params:
                 self.params["maxPopSize"] = round(self.params["maxPopSize"])
 
+    def set_param(self, **kwargs):
+        self.params.update(kwargs)
+
     def get_state(self):
         """
         Gets the current state of the algorithm as a dictionary.
diff --git a/src/metaheuristic_designer/algorithms/AlgorithmSelection.py b/src/metaheuristic_designer/algorithms/AlgorithmSelection.py
new file mode 100644
index 0000000..153d008
--- /dev/null
+++ b/src/metaheuristic_designer/algorithms/AlgorithmSelection.py
@@ -0,0 +1,131 @@
+from __future__ import annotations
+import time
+import pandas as pd
+from matplotlib import pyplot as plt
+from ..Algorithm import Algorithm
+from collections import Counter
+import enlighten
+
+
+class AlgorithmSelection:
+    """
+    Utility to evaluate and compare the performance of different algorithms.
+
+    Parameters
+    ----------
+
+    algorithm_list: Iterable[Algorithm]
+        List of algorithms to evaluate.
+    params: ParamScheduler or dict, optional
+        Indicates whether to show progress bars with 'verbose' and the number of times to repeat each algorithm with 'repetitions'
+    """
+
+    def __init__(
+        self,
+        algorithm_list: Iterable[Algorithm],
+        params: Union[ParamScheduler, dict] = None,
+    ):
+        if params is None:
+            params = {}
+
+        self.repetitions = params.get("repetitions", 10)
+
+        self.algorithm_list = algorithm_list
+
+        # Avoid repeating names
+        name_counter = Counter()
+        for alg in algorithm_list:
+            prev_name = alg.name
+            if alg.name in name_counter:
+                alg.name = alg.name + str(name_counter[alg.name] + 1)
+            name_counter.update([prev_name])
+
+        self.solutions = []
+        self.verbose = params.get("verbose", True)
+
+    def optimize(self) -> Tuple[Any, float, pd.DataFrame]:
+        """
+        Evaluates all the provided search strategies and returns the best overall solution
+        """
+
+        if self.verbose:
+            print(f"Running {len(self.algorithm_list)} algorithms {self.repetitions} times each.")
+
+        best_solution = None
+        best_fitness = 0
+        report_raw = pd.DataFrame(columns=["name", "realtime", "cputime", "fitness"])
+
+        # Create progress bar manager and global progress bar
+        if self.verbose:
+            bar_manager = enlighten.get_manager()
+            algorithm_bar = bar_manager.counter(total=len(self.algorithm_list), desc="Launching algorithms", color="red")
+
+        for idx, algorithm in enumerate(self.algorithm_list):
+            # Create new progress bar for the new algorithm
+            if self.verbose:
+                repetition_bar = bar_manager.counter(total=self.repetitions, desc=f"Evaluating {algorithm.name}", color="green", leave=False)
+
+            for rep in range(self.repetitions):
+                # Optimize using the algorithm
+                solution, fitness = algorithm.optimize()
+
+                # Get the dataframe row
+                report_raw.loc[len(report_raw.index)] = {
+                    "name": algorithm.name,
+                    "realtime": algorithm.real_time_spent,
+                    "cputime": algorithm.cpu_time_spent,
+                    "fitness": fitness,
+                }
+
+                # Save the solution if it improves the previous one
+                if (
+                    best_solution is None
+                    or (algorithm.objfunc.mode == "min" and best_fitness > fitness)
+                    or (algorithm.objfunc.mode == "max" and best_fitness < fitness)
+                ):
+                    best_solution = solution
+                    best_fitness = fitness
+
+                # Update progress bar
+                if self.verbose:
+                    repetition_bar.update()
+
+                # Reset the algorithm data
+                algorithm.restart()
+
+            # Update progress bar
+            if self.verbose:
+                algorithm_bar.update()
+
+        # Stop displaying progress bars
+        if self.verbose:
+            bar_manager.stop()
+
+        # Obtain statistics about the executions
+        report_gropued = report_raw.groupby("name", sort=False)
+        report = pd.DataFrame()
+        for group_name, group in report_gropued:
+            report = pd.concat(
+                [
+                    report,
+                    pd.DataFrame(
+                        {
+                            "name": [group_name],
+                            "realtime_min": [group["realtime"].min()],
+                            "realtime_avg": [group["realtime"].mean()],
+                            "realtime_max": [group["realtime"].max()],
+                            "realtime_std": [group["realtime"].std()],
+                            "cputime_min": [group["cputime"].min()],
+                            "cputime_avg": [group["cputime"].mean()],
+                            "cputime_max": [group["cputime"].max()],
+                            "cputime_std": [group["cputime"].std()],
+                            "fitness_min": [group["fitness"].min()],
+                            "fitness_avg": [group["fitness"].mean()],
+                            "fitness_max": [group["fitness"].max()],
+                            "fitness_std": [group["fitness"].std()],
+                        }
+                    ),
+                ]
+            )
+
+        return best_solution, best_fitness, report.reset_index(drop=True)
diff --git a/src/metaheuristic_designer/algorithms/GeneralAlgorithm.py b/src/metaheuristic_designer/algorithms/GeneralAlgorithm.py
index 69151a2..efcd37e 100644
--- a/src/metaheuristic_designer/algorithms/GeneralAlgorithm.py
+++ b/src/metaheuristic_designer/algorithms/GeneralAlgorithm.py
@@ -6,7 +6,9 @@
 
 class GeneralAlgorithm(Algorithm):
     """
-    General framework for metaheuristic algorithms
+    General framework for metaheuristic algorithms.
+
+    Performs a loop of parent selection, perturbation, and survivor selection until a stopping condition is reached. 
 
     Parameters
     ----------
@@ -19,36 +21,28 @@ class GeneralAlgorithm(Algorithm):
         Dictionary of parameters to define the stopping condition and output of the algorithm.
     """
 
-    def __init__(
-        self,
-        objfunc: ObjectiveFunc,
-        search_strategy: SearchStrategy,
-        params: Union[ParamScheduler, dict] = None,
-    ):
+    def __init__(self, objfunc: ObjectiveFunc, search_strategy: SearchStrategy, params: Union[ParamScheduler, dict] = None, name: str = None):
         """
         Constructor of the Metaheuristic class
         """
 
-        super().__init__(objfunc, search_strategy, params)
+        super().__init__(objfunc, search_strategy, params, name)
 
     def step(self, time_start=0, verbose=False):
         # Get the population of this generation
         population = self.search_strategy.population
 
         # Generate their parents
-        parents = self.search_strategy.select_parents(
-            population, progress=self.progress, history=self.best_history
-        )
+        parents = self.search_strategy.select_parents(population, progress=self.progress, history=self.best_history)
 
         # Evolve the selected parents
-        offspring = self.search_strategy.perturb(
-            parents, self.objfunc, progress=self.progress, history=self.best_history
-        )
+        offspring = self.search_strategy.perturb(parents, self.objfunc, progress=self.progress, history=self.best_history)
+
+        # Get the fitness of the individuals
+        offspring = self.search_strategy.evaluate_population(offspring, self.objfunc, self.parallel, self.threads)
 
         # Select the individuals that remain for the next generation
-        population = self.search_strategy.select_individuals(
-            population, offspring, progress=self.progress, history=self.best_history
-        )
+        population = self.search_strategy.select_individuals(population, offspring, progress=self.progress, history=self.best_history)
 
         # Assign the newly generate population
         self.search_strategy.population = population
@@ -62,35 +56,3 @@ def step(self, time_start=0, verbose=False):
         self.fit_history.append(best_fitness)
 
         return (best_individual, best_fitness)
-
-    def step_info(self, start_time=0):
-        print(f"Optimizing {self.objfunc.name} using {self.search_strategy.name}:")
-        print(f"\tReal time Spent: {round(time.time() - start_time,2)} s")
-        print(f"\tCPU time Spent:  {round(time.time() - start_time,2)} s")
-        print(f"\tGeneration: {self.steps}")
-        best_fitness = self.best_solution()[1]
-        print(f"\tBest fitness: {best_fitness}")
-        print(f"\tEvaluations of fitness: {self.objfunc.counter}")
-        self.search_strategy.extra_step_info()
-        print()
-
-    def display_report(self, show_plots=True):
-        print("Number of generations:", len(self.fit_history))
-        print("Real time spent: ", round(self.real_time_spent, 5), "s", sep="")
-        print("CPU time spent: ", round(self.cpu_time_spent, 5), "s", sep="")
-        print("Number of fitness evaluations:", self.objfunc.counter)
-
-        best_fitness = self.best_solution()[1]
-        print("Best fitness:", best_fitness)
-
-        if show_plots:
-            # Plot fitness history
-            plt.axhline(y=0, color="black", alpha=0.9)
-            plt.axvline(x=0, color="black", alpha=0.9)
-            plt.plot(self.fit_history, "blue")
-            plt.xlabel("generations")
-            plt.ylabel("fitness")
-            plt.title(f"{self.search_strategy.name} fitness")
-            plt.show()
-
-        self.search_strategy.extra_report(show_plots)
diff --git a/src/metaheuristic_designer/algorithms/MemeticAlgorithm.py b/src/metaheuristic_designer/algorithms/MemeticAlgorithm.py
index f6ebfea..250107e 100644
--- a/src/metaheuristic_designer/algorithms/MemeticAlgorithm.py
+++ b/src/metaheuristic_designer/algorithms/MemeticAlgorithm.py
@@ -6,7 +6,7 @@
 
 class MemeticAlgorithm(Algorithm):
     """
-    General framework for metaheuristic algorithms
+    Framework for algorithms that incorporate a combination of a main search strategy and a local search step.
 
     Parameters
     ----------
@@ -23,18 +23,25 @@ class MemeticAlgorithm(Algorithm):
         Dictionary of parameters to define the stopping condition and output of the algorithm.
     """
 
-    def __init__(
-        self, objfunc, search_strategy, local_search, improve_choice, params=None
-    ):
+    def __init__(self, objfunc, search_strategy, local_search, improve_choice, params=None, name=None):
         """
         Constructor of the Metaheuristic class
         """
 
-        super().__init__(objfunc, search_strategy, params)
+        super().__init__(objfunc, search_strategy, params, name)
 
         self.local_search = local_search
         self.improve_choice = improve_choice
 
+    @property
+    def name(self):
+        backup_name = f"Memetic {self.search_strategy.name}"
+        return self._name if self._name else backup_name
+
+    @name.setter
+    def name(self, new_name: str):
+        self._name = new_name
+
     def initialize(self):
         super().initialize()
         self.local_search.initialize(self.objfunc)
@@ -52,34 +59,28 @@ def _do_local_search(self, offspring):
         for idx, val in enumerate(off_idxs):
             offspring[val] = improved[idx]
 
-        current_best = max(improved, key=lambda x: x.fitness)
-        if self.search_strategy.best.fitness < current_best.fitness:
-            self.search_strategy.best = current_best
+        offspring = self.search_strategy.evaluate_population(offspring, self.objfunc, self.parallel, self.threads)
 
         return offspring
 
     def step(self, time_start=0, verbose=False):
         population = self.search_strategy.population
 
-        parents = self.search_strategy.select_parents(
-            population, progress=self.progress, history=self.best_history
-        )
+        parents = self.search_strategy.select_parents(population, progress=self.progress, history=self.best_history)
 
-        offspring = self.search_strategy.perturb(
-            parents, self.objfunc, progress=self.progress, history=self.best_history
-        )
+        offspring = self.search_strategy.perturb(parents, self.objfunc, progress=self.progress, history=self.best_history)
+
+        # Get the fitness of the individuals
+        offspring = self.search_strategy.evaluate_population(offspring, self.objfunc, self.parallel, self.threads)
 
         offspring = self._do_local_search(offspring)
 
-        population = self.search_strategy.select_individuals(
-            population, offspring, progress=self.progress, history=self.best_history
-        )
+        population = self.search_strategy.select_individuals(population, offspring, progress=self.progress, history=self.best_history)
 
         self.search_strategy.population = population
 
         best_individual, best_fitness = self.search_strategy.best_solution()
         self.search_strategy.update_params(progress=self.progress)
-        self.steps += 1
 
         # Store information
         self.best_history.append(best_individual)
@@ -126,36 +127,6 @@ def get_state(
         return data
 
     def step_info(self, start_time=0):
-        print(
-            f"Optimizing {self.objfunc.name} using {self.search_strategy.name}+{self.local_search.name}:"
-        )
-        print(f"\tReal time Spent: {round(time.time() - start_time,2)} s")
-        print(f"\tCPU time Spent:  {round(time.time() - start_time,2)} s")
-        print(f"\tGeneration: {self.steps}")
-        best_fitness = self.best_solution()[1]
-        print(f"\tBest fitness: {best_fitness}")
-        print(f"\tEvaluations of fitness: {self.objfunc.counter}")
-        self.search_strategy.extra_step_info()
+        super().step_info(start_time)
         self.local_search.extra_step_info()
         print()
-
-    def display_report(self, show_plots=True):
-        print("Number of generations:", len(self.fit_history))
-        print("Real time spent: ", round(self.real_time_spent, 5), "s", sep="")
-        print("CPU time spent: ", round(self.cpu_time_spent, 5), "s", sep="")
-        print("Number of fitness evaluations:", self.objfunc.counter)
-
-        best_fitness = self.best_solution()[1]
-        print("Best fitness:", best_fitness)
-
-        if show_plots:
-            # Plot fitness history
-            plt.axhline(y=0, color="black", alpha=0.9)
-            plt.axvline(x=0, color="black", alpha=0.9)
-            plt.plot(self.fit_history, "blue")
-            plt.xlabel("generations")
-            plt.ylabel("fitness")
-            plt.title(f"{self.search_strategy.name} fitness")
-            plt.show()
-
-        self.search_strategy.extra_report(show_plots)
diff --git a/src/metaheuristic_designer/algorithms/StrategySelection.py b/src/metaheuristic_designer/algorithms/StrategySelection.py
new file mode 100644
index 0000000..1fa95f2
--- /dev/null
+++ b/src/metaheuristic_designer/algorithms/StrategySelection.py
@@ -0,0 +1,45 @@
+from __future__ import annotations
+import time
+import pandas as pd
+from matplotlib import pyplot as plt
+from ..SearchStrategy import SearchStrategy
+from .GeneralAlgorithm import GeneralAlgorithm
+from .AlgorithmSelection import AlgorithmSelection
+from collections import Counter
+
+
+class StrategySelection:
+    """
+    Utility to evaluate and compare the performance of different search strategies.
+
+    Parameters
+    ----------
+
+    objfunc: ObjectiveFunc
+        Objective function to evaluate
+    strategy_list: Iterable[SearchStrategy]
+        List of algorithms to evaluate.
+    algorithm_params: ParamScheduler or dict, optional
+        Parameters shared by all the algorithms being run.
+    params: ParamScheduler or dict, optional
+        Indicates whether to show progress bars with 'verbose' and the number of times to repeat each algorithm with 'repetitions'
+    """
+
+    def __init__(
+        self,
+        objfunc: ObjectiveFunc,
+        strategy_list: Iterable[SearchStrategy],
+        algorithm_params: Union[ParamScheduler, dict] = None,
+        params: Union[ParamScheduler, dict] = None,
+    ):
+        self.strategy_list = strategy_list
+        algorithm_list = [GeneralAlgorithm(objfunc, strategy, algorithm_params) for strategy in strategy_list]
+
+        self.algorithm_selection = AlgorithmSelection(algorithm_list, params)
+
+    def optimize(self) -> Tuple[Any, float, pd.DataFrame]:
+        """
+        Evaluates all the provided search strategies and returns the best overall solution
+        """
+
+        return self.algorithm_selection.optimize()
diff --git a/src/metaheuristic_designer/algorithms/__init__.py b/src/metaheuristic_designer/algorithms/__init__.py
index 057144a..ecdfddb 100644
--- a/src/metaheuristic_designer/algorithms/__init__.py
+++ b/src/metaheuristic_designer/algorithms/__init__.py
@@ -1,2 +1,4 @@
 from .GeneralAlgorithm import GeneralAlgorithm
 from .MemeticAlgorithm import MemeticAlgorithm
+from .AlgorithmSelection import AlgorithmSelection
+from .StrategySelection import StrategySelection
diff --git a/src/metaheuristic_designer/benchmarks/benchmark_funcs.py b/src/metaheuristic_designer/benchmarks/benchmark_funcs.py
index e889b0e..6ab34f7 100644
--- a/src/metaheuristic_designer/benchmarks/benchmark_funcs.py
+++ b/src/metaheuristic_designer/benchmarks/benchmark_funcs.py
@@ -3,6 +3,7 @@
 # from numba import jit
 from ..ObjectiveFunc import ObjectiveVectorFunc
 from ..utils import RAND_GEN
+import time
 
 
 class MaxOnes(ObjectiveVectorFunc):
@@ -43,6 +44,18 @@ def repair_solution(self, solution):
         return np.clip(solution.copy(), 0, 1)
 
 
+class SleepTest(ObjectiveVectorFunc):
+    def __init__(self, size, sleep_time=2, opt="min"):
+        self.size = size
+        self.sleep_time = sleep_time
+        super().__init__(self.size, opt, -100, 100, name="Sphere function")
+
+    def objective(self, solution):
+        time.sleep(self.sleep_time)
+        print("hello")
+        return _sphere(solution)
+
+
 ### Benchmark functions
 class Sphere(ObjectiveVectorFunc):
     def __init__(self, size, opt="min"):
@@ -56,9 +69,7 @@ def objective(self, solution):
 class HighCondElliptic(ObjectiveVectorFunc):
     def __init__(self, size, opt="min"):
         self.size = size
-        super().__init__(
-            self.size, opt, -5.12, 5.12, name="High condition elliptic function"
-        )
+        super().__init__(self.size, opt, -5.12, 5.12, name="High condition elliptic function")
 
     def objective(self, solution):
         return _high_cond_elipt_f(solution)
@@ -166,9 +177,7 @@ def objective(self, solution):
 class ExpandedGriewankPlusRosenbrock(ObjectiveVectorFunc):
     def __init__(self, size, opt="min"):
         self.size = size
-        super().__init__(
-            self.size, opt, -100, 100, name="Expanded Griewank + Rosenbrock"
-        )
+        super().__init__(self.size, opt, -100, 100, name="Expanded Griewank + Rosenbrock")
 
     def objective(self, solution):
         return _exp_griewank_plus_rosenbrock(solution)
@@ -192,9 +201,7 @@ def __init__(self, size, opt="min", lim_min=-1, lim_max=1):
         self.size = size
         self.lim_min = lim_min
         self.lim_max = lim_max
-        super().__init__(
-            self.size, opt, low_lim=lim_min, up_lim=lim_max, name="Sum Powell"
-        )
+        super().__init__(self.size, opt, low_lim=lim_min, up_lim=lim_max, name="Sum Powell")
 
     def objective(self, solution):
         return _sum_powell(solution)
@@ -204,19 +211,12 @@ def repair_solution(self, solution, parent=None):
         if parent:
             mask_inf = solution < self.lim_min
             mask_sup = solution > self.lim_max
-            solution[mask_inf] = parent[mask_inf] + RAND_GEN.random() * (
-                self.lim_min - parent[mask_inf]
-            )
-            solution[mask_sup] = parent[mask_sup] + RAND_GEN.random() * (
-                parent[mask_sup] - self.lim_max
-            )
+            solution[mask_inf] = parent[mask_inf] + RAND_GEN.random() * (self.lim_min - parent[mask_inf])
+            solution[mask_sup] = parent[mask_sup] + RAND_GEN.random() * (parent[mask_sup] - self.lim_max)
         # random in range
         else:
             mask = (solution < self.lim_min) | (solution > self.lim_max)
-            solution[mask] = (
-                RAND_GEN.random(len(mask[mask == True])) * (self.lim_max - self.lim_min)
-                - self.lim_min
-            )
+            solution[mask] = RAND_GEN.random(len(mask[mask == True])) * (self.lim_max - self.lim_min) - self.lim_min
         return solution
 
 
@@ -229,9 +229,7 @@ def __init__(self, size, opt="min", lim_min=-10, lim_max=10):
         self.size = size
         self.lim_min = lim_min
         self.lim_max = lim_max
-        super().__init__(
-            self.size, opt, low_lim=lim_min, up_lim=lim_max, name="N4 Xin-She Yang"
-        )
+        super().__init__(self.size, opt, low_lim=lim_min, up_lim=lim_max, name="N4 Xin-She Yang")
 
     def objective(self, solution):
         return _n4xinshe_yang(solution)
@@ -241,19 +239,12 @@ def repair_solution(self, solution, parent=None):
         if parent:
             mask_inf = solution < self.lim_min
             mask_sup = solution > self.lim_max
-            solution[mask_inf] = parent[mask_inf] + RAND_GEN.random() * (
-                self.lim_min - parent[mask_inf]
-            )
-            solution[mask_sup] = parent[mask_sup] + RAND_GEN.random() * (
-                parent[mask_sup] - self.lim_max
-            )
+            solution[mask_inf] = parent[mask_inf] + RAND_GEN.random() * (self.lim_min - parent[mask_inf])
+            solution[mask_sup] = parent[mask_sup] + RAND_GEN.random() * (parent[mask_sup] - self.lim_max)
         # random in range
         else:
             mask = (solution < self.lim_min) | (solution > self.lim_max)
-            solution[mask] = (
-                RAND_GEN.random(len(mask[mask == True])) * (self.lim_max - self.lim_min)
-                - self.lim_min
-            )
+            solution[mask] = RAND_GEN.random(len(mask[mask == True])) * (self.lim_max - self.lim_min) - self.lim_min
         return solution
 
 
@@ -296,14 +287,7 @@ def _ackley(solution):
 
 # @jit(nopython=False)
 def _weierstrass(solution, iter=20):
-    return np.sum(
-        np.array(
-            [
-                0.5**k * np.cos(2 * np.pi * 3**k * (solution + 0.5))
-                for k in range(iter)
-            ]
-        )
-    )
+    return np.sum(np.array([0.5**k * np.cos(2 * np.pi * 3**k * (solution + 0.5)) for k in range(iter)]))
 
 
 # @jit(nopython=True)
@@ -345,10 +329,7 @@ def _katsuura(solution):
         + (i + 1)
         * np.sum(
             (
-                np.abs(
-                    2 ** (np.arange(1, 32 + 1)) * solution[i]
-                    - np.round(2 ** (np.arange(1, 32 + 1)) * solution[i])
-                )
+                np.abs(2 ** (np.arange(1, 32 + 1)) * solution[i] - np.round(2 ** (np.arange(1, 32 + 1)) * solution[i]))
                 * 2 ** (-np.arange(1, 32 + 1, dtype=float))
             )
             ** (10 / solution.size**1.2)
@@ -396,9 +377,7 @@ def _exp_shafferF6(solution):
     term2 = 1 + 0.001 * (solution[:-1] ** 2 + solution[1:] ** 2).sum()
     temp = 0.5 + term1 / term2
 
-    term1 = (
-        np.sin(np.sqrt(np.sum((solution.size - 1) ** 2 + solution[0] ** 2))) ** 2 - 0.5
-    )
+    term1 = np.sin(np.sqrt(np.sum((solution.size - 1) ** 2 + solution[0] ** 2))) ** 2 - 0.5
     term2 = 1 + 0.001 * ((solution.size - 1) ** 2 + solution[0] ** 2)
 
     return temp + 0.5 + term1 / term2
diff --git a/src/metaheuristic_designer/benchmarks/classic_problems.py b/src/metaheuristic_designer/benchmarks/classic_problems.py
index 1f3c93d..988c65f 100644
--- a/src/metaheuristic_designer/benchmarks/classic_problems.py
+++ b/src/metaheuristic_designer/benchmarks/classic_problems.py
@@ -23,9 +23,7 @@ class ThreeSAT(ObjectiveVectorFunc):
 
     def __init__(self, clauses):
         if not isinstance(clauses, np.ndarray) or clauses.shape[1] != 3:
-            raise ValueError(
-                "The caluses must be represented as an array of size (n_clauses, 3)."
-            )
+            raise ValueError("The caluses must be represented as an array of size (n_clauses, 3).")
 
         self.clauses = clauses
         self.n_vars = np.abs(clauses).max()
@@ -50,14 +48,10 @@ def from_cnf_file(path):
         clauses_arr = np.asarray(clauses)
 
         if len(clauses) != n_clauses:
-            warnings.warn(
-                "The number of clauses in the file was incorrect.", stacklevel=2
-            )
+            warnings.warn("The number of clauses in the file was incorrect.", stacklevel=2)
 
         if np.abs(clauses_arr).max() != n_vars:
-            warnings.warn(
-                "The number of variables in the file was incorrect.", stacklevel=2
-            )
+            warnings.warn("The number of variables in the file was incorrect.", stacklevel=2)
 
         return ThreeSAT(clauses_arr)
 
@@ -104,9 +98,7 @@ def __init__(self, cost, value, max_weight):
         cost = np.asarray(cost)
         value = np.asarray(value)
         if cost.size != value.size:
-            raise ValueError(
-                "The value vector must have the same dimension as the cost vector."
-            )
+            raise ValueError("The value vector must have the same dimension as the cost vector.")
 
         self.cost = cost
         self.value = value
diff --git a/src/metaheuristic_designer/benchmarks/img_funcs.py b/src/metaheuristic_designer/benchmarks/img_funcs.py
index dad5412..778ea72 100644
--- a/src/metaheuristic_designer/benchmarks/img_funcs.py
+++ b/src/metaheuristic_designer/benchmarks/img_funcs.py
@@ -7,7 +7,7 @@
 
 
 class ImgApprox(ObjectiveVectorFunc):
-    def __init__(self, img_dim, reference, opt="min", img_name=""):
+    def __init__(self, img_dim, reference, opt="min", img_name="", diff_func="sts"):
         self.img_dim = tuple(img_dim) + (3,)
         self.size = img_dim[0] * img_dim[1] * 3
         self.reference = reference.resize((img_dim[0], img_dim[1]))
@@ -18,31 +18,30 @@ def __init__(self, img_dim, reference, opt="min", img_name=""):
         else:
             name = f'Approximating "{img_name}"'
 
-        opt = "max"
+        self.diff_func = diff_func
+        if diff_func in ["mse"]:
+            opt = "min"
+        else:
+            opt = "max"
 
         super().__init__(self.size, opt, 0, 256, name=name)
 
     def objective(self, solution):
-        # return imgdistance(solution, self.reference)
-        return (
-            sum(
-                metrics.structural_similarity(
-                    solution[:, :, i], self.reference[:, :, i]
-                )
-                for i in range(3)
-            )
-            / 3
-        )
+        error = 0
+        if self.diff_func == "mse":
+            error = imgdistance_mse(solution, self.reference)
+        elif self.diff_func == "sts":
+            error = 0
+            for i in range(3):
+                error += metrics.structural_similarity(solution[:, :, i], self.reference[:, :, i])
+            error /= 3
+
+        return error
 
     def repair_solution(self, solution):
         return np.clip(solution, 0, 255)
 
 
-# @jit(nopython=True)
-def imgdistance(img, reference):
-    return np.sum((img - reference) ** 2)
-
-
 class ImgStd(ObjectiveVectorFunc):
     def __init__(self, img_dim, opt="max"):
         self.size = img_dim[0] * img_dim[1] * 3
@@ -69,10 +68,7 @@ def __init__(self, img_dim, nbins=10, opt="min"):
 
     def objective(self, solution):
         solution_channels = solution.reshape([3, -1])
-        img_hists = [
-            np.histogram(solution_channels[i], bins=np.linspace(0, 256, self.nbins))[0]
-            for i in range(3)
-        ]
+        img_hists = [np.histogram(solution_channels[i], bins=np.linspace(0, 256, self.nbins))[0] for i in range(3)]
         img_hists = np.array(img_hists) / solution_channels.shape[1]
         img_hists_no_zeros = img_hists
         img_hists_no_zeros[img_hists == 0] = 1
@@ -86,14 +82,12 @@ class ImgExperimental(ObjectiveVectorFunc):
     def __init__(self, img_dim, reference, img_name, opt="min"):
         self.img_dim = tuple(img_dim) + (3,)
         self.size = img_dim[0] * img_dim[1] * 3
-        self.reference = np.asarray(reference.resize([img_dim[0], img_dim[1]]))[
-            :, :, :3
-        ].astype(np.uint32)
+        self.reference = np.asarray(reference.resize([img_dim[0], img_dim[1]]))[:, :, :3].astype(np.uint32)
 
         super().__init__(self.size, opt, 0, 256, name="Image approx and std")
 
     def objective(self, solution):
-        dist = imgdistance(solution, self.reference)
+        dist = imgdistance_mse(solution, self.reference)
         dist_norm = dist / (np.sqrt(self.size) * 255)
 
         solution_rounded = solution // 75
@@ -112,3 +106,8 @@ def repair_solution(self, solution):
 
     def repair_speed(self, solution):
         return np.clip(solution, -255, 255)
+
+
+# @jit(nopython=True)
+def imgdistance_mse(img, reference):
+    return np.sum((img - reference) ** 2)
diff --git a/src/metaheuristic_designer/encodings/AdaptionEncoding.py b/src/metaheuristic_designer/encodings/AdaptionEncoding.py
new file mode 100644
index 0000000..c7f2d5c
--- /dev/null
+++ b/src/metaheuristic_designer/encodings/AdaptionEncoding.py
@@ -0,0 +1,33 @@
+from __future__ import annotations
+from ..Encoding import Encoding
+from abc import abstractmethod
+from .DefaultEncoding import DefaultEncoding
+
+
+class AdaptionEncoding(Encoding):
+    """
+    Abstract Adaption Encoding class.
+
+    This kind of encoding will represent solutions as a vector with the solution and some parameters for the operators of the algorithm.
+    """
+
+    def __init__(self, vecsize: int, nparams: int, base_encoding: Encoding = None):
+        self.vecsize = vecsize
+        self.nparams = nparams
+        if base_encoding is None:
+            base_encoding = DefaultEncoding()
+        self.base_encoding = base_encoding
+
+    def encode(self, phenotype: Any, param_vec: np.ndarray = None) -> np.ndarray:
+        phenotype_encoded = self.base_encoding.encode(phenotype)
+        return np.concatenate([phenotype_encoded, param_vec])
+
+    def decode(self, genotype: np.ndarray) -> np.ndarray:
+        return self.base_encoding.decode(genotype[:self.vecsize])
+
+    def decode_param_vec(self, genotype):
+        return genotype[self.vecsize:]
+
+    @abstractmethod
+    def decode_param(self, genotype) -> dict:
+        pass
diff --git a/src/metaheuristic_designer/encodings/CMAEncoding.py b/src/metaheuristic_designer/encodings/CMAEncoding.py
index 9570d3e..111e178 100644
--- a/src/metaheuristic_designer/encodings/CMAEncoding.py
+++ b/src/metaheuristic_designer/encodings/CMAEncoding.py
@@ -1,27 +1,15 @@
 from __future__ import annotations
 import numpy as np
-from ..Encoding import Encoding
+
+# from ..Encoding import Encoding
+from .AdaptionEncoding import AdaptionEncoding
 from .DefaultEncoding import DefaultEncoding
 
 
-class CMAEncoding(Encoding):
+class CMAEncoding(AdaptionEncoding):
     """
-    Decoder used to implement the CMA-ES algorithm
+    Decoder used to implement the CMA-ES algorithm.
     """
 
-    def __init__(self, nparams: int, pre_encoding: Encoding = None):
-        self.nparams = nparams
-
-        if pre_encoding is None:
-            pre_encoding = DefaultDecoder()
-        self.pre_encoding = pre_encoding
-
-    def encode(self, phenotype: np.ndarray, param_vec: np.ndarray = None) -> np.ndarray:
-        encoded = self.pre_encoding.encode(phenotype)
-        if param_vec is None:
-            param_vec = np.ones(self.nparams)
-        return np.concatenate(encoded, param_vec)
-
-    def decode(self, genotype: np.ndarray) -> np.ndarray:
-        decoded = genotype[: -self.nparams]
-        return self.pre_encoding.decode(decoded)
+    def decode_param(self, genotype: np.ndarray) -> np.ndarray:
+        return {"F": self.decode_param_vec(genotype)}
diff --git a/src/metaheuristic_designer/encodings/DefaultEncoding.py b/src/metaheuristic_designer/encodings/DefaultEncoding.py
index f669967..e94539f 100644
--- a/src/metaheuristic_designer/encodings/DefaultEncoding.py
+++ b/src/metaheuristic_designer/encodings/DefaultEncoding.py
@@ -5,11 +5,11 @@
 
 class DefaultEncoding(Encoding):
     """
-    Default encoder that uses the input vector as from the individual as the solution
+    Default encoder that uses the genotype directly as a solution.
     """
 
-    def encode(self, phenotype: np.ndarray) -> np.ndarray:
+    def encode(self, phenotype: Any) -> Any:
         return phenotype
 
-    def decode(self, genotype: np.ndarray) -> np.ndarray:
+    def decode(self, genotype: Any) -> Any:
         return genotype
diff --git a/src/metaheuristic_designer/encodings/ImageEncoding.py b/src/metaheuristic_designer/encodings/ImageEncoding.py
index 6276207..0ee1d5d 100644
--- a/src/metaheuristic_designer/encodings/ImageEncoding.py
+++ b/src/metaheuristic_designer/encodings/ImageEncoding.py
@@ -5,7 +5,7 @@
 
 class ImageEncoding(Encoding):
     """
-    Decoder used to evolve images
+    Decoder used to evolve images.
     """
 
     def __init__(self, shape, color=True):
diff --git a/src/metaheuristic_designer/encodings/LambdaEncoding.py b/src/metaheuristic_designer/encodings/LambdaEncoding.py
index b7ea5e5..62ee879 100644
--- a/src/metaheuristic_designer/encodings/LambdaEncoding.py
+++ b/src/metaheuristic_designer/encodings/LambdaEncoding.py
@@ -5,7 +5,7 @@
 
 class LambdaEncoding(Encoding):
     """
-    Decoder that uses user specified functions
+    Decoder that uses user specified functions.
     """
 
     def __init__(self, encode_fn: callable, decode_fn: callable):
diff --git a/src/metaheuristic_designer/encodings/MatrixEncoding.py b/src/metaheuristic_designer/encodings/MatrixEncoding.py
index 43c12ce..34be3ed 100644
--- a/src/metaheuristic_designer/encodings/MatrixEncoding.py
+++ b/src/metaheuristic_designer/encodings/MatrixEncoding.py
@@ -5,7 +5,7 @@
 
 class MatrixEncoding(Encoding):
     """
-    Decoder used to evolve matrices
+    Decoder used to evolve matrices.
     """
 
     def __init__(self, shape):
diff --git a/src/metaheuristic_designer/encodings/__init__.py b/src/metaheuristic_designer/encodings/__init__.py
index c343637..aeec22f 100644
--- a/src/metaheuristic_designer/encodings/__init__.py
+++ b/src/metaheuristic_designer/encodings/__init__.py
@@ -4,3 +4,4 @@
 from .MatrixEncoding import MatrixEncoding
 from .TypeCastEncoding import TypeCastEncoding
 from .LambdaEncoding import LambdaEncoding
+from .AdaptionEncoding import AdaptionEncoding
diff --git a/src/metaheuristic_designer/initializers/GaussianInitializer.py b/src/metaheuristic_designer/initializers/GaussianInitializer.py
index ab9f67e..dc79436 100644
--- a/src/metaheuristic_designer/initializers/GaussianInitializer.py
+++ b/src/metaheuristic_designer/initializers/GaussianInitializer.py
@@ -26,18 +26,14 @@ class GaussianInitializer(Initializer):
         Data type used in each of the components of the vector in the individual.
     """
 
-    def __init__(
-        self, genotype_size, g_mean, g_std, pop_size=1, encoding=None, dtype=float
-    ):
+    def __init__(self, genotype_size, g_mean, g_std, pop_size=1, encoding=None, dtype=float):
         super().__init__(pop_size, encoding)
 
         self.genotype_size = genotype_size
 
         if type(g_mean) in [list, tuple, np.ndarray]:
             if len(g_mean) != genotype_size:
-                raise ValueError(
-                    f"If g_mean is a sequence it must be of length {genotype_size}."
-                )
+                raise ValueError(f"If g_mean is a sequence it must be of length {genotype_size}.")
 
             self.g_mean = g_mean
         else:
@@ -45,9 +41,7 @@ def __init__(
 
         if type(g_std) in [list, tuple, np.ndarray]:
             if len(g_std) != genotype_size:
-                raise ValueError(
-                    f"If g_std is a sequence it must be of length {genotype_size}."
-                )
+                raise ValueError(f"If g_std is a sequence it must be of length {genotype_size}.")
 
             self.g_std = g_std
         else:
@@ -58,9 +52,7 @@ def __init__(
 
 class GaussianVectorInitializer(GaussianInitializer):
     def generate_random(self, objfunc):
-        new_vector_float = RAND_GEN.normal(
-            self.g_mean, self.g_std, size=self.genotype_size
-        )
+        new_vector_float = RAND_GEN.normal(self.g_mean, self.g_std, size=self.genotype_size)
         if self.dtype is int:
             new_vector = np.round(new_vector_float).astype(self.dtype)
         else:
diff --git a/src/metaheuristic_designer/initializers/LambdaInitializer.py b/src/metaheuristic_designer/initializers/LambdaInitializer.py
index 116b26e..4309547 100644
--- a/src/metaheuristic_designer/initializers/LambdaInitializer.py
+++ b/src/metaheuristic_designer/initializers/LambdaInitializer.py
@@ -20,9 +20,7 @@ class LambdaInitializer(Initializer):
         Encoding that will be passed to each individual.
     """
 
-    def __init__(
-        self, generator: callable, pop_size: int = 1, encoding: Encoding = None
-    ):
+    def __init__(self, generator: callable, pop_size: int = 1, encoding: Encoding = None):
         self.pop_size = pop_size
         self.generator = generator
 
diff --git a/src/metaheuristic_designer/initializers/SeedInitializer.py b/src/metaheuristic_designer/initializers/SeedInitializer.py
index 35cbc85..502f683 100644
--- a/src/metaheuristic_designer/initializers/SeedInitializer.py
+++ b/src/metaheuristic_designer/initializers/SeedInitializer.py
@@ -14,8 +14,8 @@ class SeedProbInitializer(Initializer):
     ----------
     default_init: Initializer
         Initializer used to generate random individuals.
-    solutions: List[Individual]
-        The predefined individuals that will be inserted into the population.
+    solutions: List[Any]
+        The predefined solutions that will be inserted into the population.
     insert_prob: float
         Probability of inserting one of the predefined solutions into the population.
     """
@@ -23,7 +23,7 @@ class SeedProbInitializer(Initializer):
     def __init__(
         self,
         default_init: Initializer,
-        solutions: List[Individual],
+        solutions: List[Any],
         insert_prob: float = 0.1,
     ):
         super().__init__(default_init.pop_size)
@@ -38,7 +38,8 @@ def generate_random(self, objfunc):
     def generate_individual(self, objfunc):
         new_indiv = None
         if random.random() < self.insert_prob:
-            new_indiv = random.choice(self.solutions)
+            new_solution = random.choice(self.solutions)
+            new_indiv = Individual(objfunc, new_solution, encoding=self.default_init.encoding)
         else:
             new_indiv = self.default_init.generate_individual(objfunc)
         return new_indiv
@@ -53,8 +54,8 @@ class SeedDetermInitializer(Initializer):
     ----------
     default_init: Initializer
         Initializer used to generate random individuals.
-    solutions: List[Individual]
-        The predefined individuals that will be inserted into the population.
+    solutions: List[Any]
+        The predefined solutions that will be inserted into the population.
     n_to_insert: float
         Amount of predefined individuals to insert in the population.
     """
@@ -62,7 +63,7 @@ class SeedDetermInitializer(Initializer):
     def __init__(
         self,
         default_init: Initializer,
-        solutions: List[Individual],
+        solutions: List[Any],
         n_to_insert: int = None,
     ):
         super().__init__(default_init.pop_size)
@@ -82,7 +83,8 @@ def generate_random(self, objfunc):
     def generate_individual(self, objfunc):
         new_indiv = None
         if self.inserted < self.number_to_insert:
-            new_indiv = self.solutions[self.inserted % len(self.solutions)]
+            new_solution = self.solutions[self.inserted % len(self.solutions)]
+            new_indiv = Individual(objfunc, new_solution, encoding=self.default_init.encoding)
         else:
             new_indiv = self.default_init.generate_individual(objfunc)
 
diff --git a/src/metaheuristic_designer/initializers/UniformInitializer.py b/src/metaheuristic_designer/initializers/UniformInitializer.py
index 6840285..45a863b 100644
--- a/src/metaheuristic_designer/initializers/UniformInitializer.py
+++ b/src/metaheuristic_designer/initializers/UniformInitializer.py
@@ -3,6 +3,7 @@
 import random
 from ..Initializer import Initializer
 from ..Individual import Individual
+from ..encodings import AdaptionEncoding
 from ..utils import RAND_GEN
 
 
@@ -26,18 +27,16 @@ class UniformInitializer(Initializer):
         Data type used in each of the components of the vector in the individual.
     """
 
-    def __init__(
-        self, genotype_size, low_lim, up_lim, pop_size=1, encoding=None, dtype=float
-    ):
+    def __init__(self, genotype_size, low_lim, up_lim, pop_size=1, encoding=None, dtype=float):
         super().__init__(pop_size, encoding)
 
         self.genotype_size = genotype_size
+        if isinstance(encoding, AdaptionEncoding):
+            self.genotype_size = encoding.vecsize + encoding.nparams
 
         if type(low_lim) in [list, tuple, np.ndarray]:
             if len(low_lim) != genotype_size:
-                raise ValueError(
-                    f"If low_lim is a sequence it must be of length {genotype_size}."
-                )
+                raise ValueError(f"If low_lim is a sequence it must be of length {genotype_size}.")
 
             self.low_lim = low_lim
         else:
@@ -45,9 +44,7 @@ def __init__(
 
         if type(up_lim) in [list, tuple, np.ndarray]:
             if len(up_lim) != genotype_size:
-                raise ValueError(
-                    f"If up_lim is a sequence it must be of length {genotype_size}."
-                )
+                raise ValueError(f"If up_lim is a sequence it must be of length {genotype_size}.")
 
             self.up_lim = up_lim
         else:
@@ -58,9 +55,7 @@ def __init__(
 
 class UniformVectorInitializer(UniformInitializer):
     def generate_random(self, objfunc):
-        new_vector_float = RAND_GEN.uniform(
-            self.low_lim, self.up_lim, size=self.genotype_size
-        )
+        new_vector_float = RAND_GEN.uniform(self.low_lim, self.up_lim, size=self.genotype_size)
         if self.dtype is int:
             new_vector = np.round(new_vector_float).astype(self.dtype)
         else:
@@ -74,9 +69,7 @@ def generate_individual(self, objfunc):
 
 class UniformListInitializer(UniformInitializer):
     def generate_random(self, objfunc):
-        new_list = [
-            RAND_GEN.uniform(low, up) for low, up in zip(self.low_lim, self.up_lim)
-        ]
+        new_list = [RAND_GEN.uniform(low, up) for low, up in zip(self.low_lim, self.up_lim)]
         return Individual(objfunc, new_list, encoding=self.encoding)
 
     def generate_individual(self, objfunc):
diff --git a/src/metaheuristic_designer/operators/OperatorAdaptative.py b/src/metaheuristic_designer/operators/OperatorAdaptative.py
new file mode 100644
index 0000000..35d0c36
--- /dev/null
+++ b/src/metaheuristic_designer/operators/OperatorAdaptative.py
@@ -0,0 +1,69 @@
+from __future__ import annotations
+import random
+from ..Operator import Operator
+from ..encodings import AdaptionEncoding
+from .OperatorMeta import OperatorMeta
+from .OperatorReal import OperatorReal
+from .OperatorNull import OperatorNull
+from copy import copy
+import numpy as np
+import enum
+from enum import Enum
+
+
+class OperatorAdaptative(Operator):
+    """
+    Operator class that allow algorithms to self-adapt by mutating the operator's parameters.
+
+    Parameters
+    ----------
+        base_operator: Operator
+            Operator that will be applied to the solution we are evaluating.
+        param_operator: Operator
+            Operator that will be applied to the parameters of the base operator.
+        param_encoding: AdaptionEncoding
+            Encoding that divides the genotype into the solution and the operator's parameters.
+        params: Union[ParamScheduler, dict]
+            Optional parameters that are used by the operator.
+        name: str
+            Name of the operator.
+    """
+
+    def __init__(
+        self,
+        base_operator: Operator,
+        param_operator: Operator,
+        param_encoding: AdaptionEncoding,
+        params: Union[ParamScheduler, dict] = None,
+        name: str = None,
+    ):
+        """
+        Constructor for the OperatorAdaptative class
+        """
+
+        super().__init__(params, base_operator.name + "-Adaptative")
+
+        vecmask = np.concatenate([np.zeros(param_encoding.vecsize), np.ones(param_encoding.nparams)])
+        null_op = OperatorNull()
+
+        self.base_operator = base_operator
+        self.param_operator = param_operator
+        self.base_operator_split = OperatorMeta("Split", [base_operator, null_op], {"mask": vecmask})
+        self.param_operator_split = OperatorMeta("Split", [null_op, param_operator], {"mask": vecmask})
+        self.param_encoding = param_encoding
+
+    def evolve(self, indiv, population, objfunc, global_best, initializer=None):
+        # Evolve only parameters
+        indiv_conf_param = self.param_operator_split.evolve(indiv, population, objfunc, global_best, initializer)
+
+        # Apply parameters to operator
+        indiv_params = self.param_encoding.decode_param(indiv_conf_param.genotype)
+        self.base_operator.params.update(indiv_params)
+
+        # Evolve solution
+        new_indiv = self.base_operator_split.evolve(indiv_conf_param, population, objfunc, global_best, initializer)
+        return new_indiv
+
+    def step(self, progress: float):
+        super().step(progress)
+        self.base_operator.step(progress)
diff --git a/src/metaheuristic_designer/operators/OperatorBinary.py b/src/metaheuristic_designer/operators/OperatorBinary.py
index 85f8447..a0db15e 100644
--- a/src/metaheuristic_designer/operators/OperatorBinary.py
+++ b/src/metaheuristic_designer/operators/OperatorBinary.py
@@ -69,11 +69,9 @@ class OperatorBinary(Operator):
         Name that is associated with the operator.
     """
 
-    def __init__(
-        self, method: str, params: Union[ParamScheduler, dict] = None, name: str = None
-    ):
+    def __init__(self, method: str, params: Union[ParamScheduler, dict] = None, name: str = None):
         """
-        Constructor for the Operator class
+        Constructor for the OperatorBinary class
         """
 
         if name is None:
@@ -103,9 +101,7 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
         params = copy(self.params)
 
         if "Cr" in params and "N" not in params:
-            params["N"] = np.count_nonzero(
-                RAND_GEN.random(indiv.genotype.size) < params["Cr"]
-            )
+            params["N"] = np.count_nonzero(RAND_GEN.random(indiv.genotype.size) < params["Cr"])
 
         if "N" in params:
             params["N"] = round(params["N"])
@@ -121,9 +117,7 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
             new_indiv.genotype = cross_mp(new_indiv.genotype, indiv2.genotype.copy())
 
         elif self.method == BinOpMethods.MULTICROSS:
-            new_indiv.genotype = multi_cross(
-                new_indiv.genotype, others, params["Nindiv"]
-            )
+            new_indiv.genotype = multi_cross(new_indiv.genotype, others, params["Nindiv"])
 
         elif self.method == BinOpMethods.PERM:
             new_indiv.genotype = permutation(new_indiv.genotype, params["N"])
@@ -144,14 +138,10 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
             new_indiv = initializer.generate_random(objfunc)
 
         elif self.method == BinOpMethods.RANDOM_MASK:
-            mask_pos = np.hstack(
-                [np.ones(params["N"]), np.zeros(new_indiv.genotype.size - params["N"])]
-            ).astype(bool)
+            mask_pos = np.hstack([np.ones(params["N"]), np.zeros(new_indiv.genotype.size - params["N"])]).astype(bool)
             RAND_GEN.shuffle(mask_pos)
 
-            new_indiv.genotype[mask_pos] = initializer.generate_random(
-                objfunc
-            ).genotype[mask_pos]
+            new_indiv.genotype[mask_pos] = initializer.generate_random(objfunc).genotype[mask_pos]
 
         elif self.method == BinOpMethods.DUMMY:
             new_indiv.genotype = dummy_op(new_indiv.genotype, params["F"])
@@ -160,5 +150,8 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
             fn = params["function"]
             new_indiv.genotype = fn(indiv, population, objfunc, params)
 
+        elif self.method == BinOpMethods.NOTHING:
+            new_indiv = indiv
+
         new_indiv.genotype = (new_indiv.genotype != 0).astype(np.int32)
         return new_indiv
diff --git a/src/metaheuristic_designer/operators/OperatorInt.py b/src/metaheuristic_designer/operators/OperatorInt.py
index 77509b7..a45f668 100644
--- a/src/metaheuristic_designer/operators/OperatorInt.py
+++ b/src/metaheuristic_designer/operators/OperatorInt.py
@@ -107,11 +107,9 @@ class OperatorInt(Operator):
         Name that is associated with the operator.
     """
 
-    def __init__(
-        self, method: str, params: Union[ParamScheduler, dict] = None, name: str = None
-    ):
+    def __init__(self, method: str, params: Union[ParamScheduler, dict] = None, name: str = None):
         """
-        Constructor for the Operator class
+        Constructor for the OperatorInt class
         """
 
         if name is None:
@@ -152,9 +150,7 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
         params = copy(self.params)
 
         if "Cr" in params and "N" not in params:
-            params["N"] = np.count_nonzero(
-                RAND_GEN.random(indiv.genotype.size) < params["Cr"]
-            )
+            params["N"] = np.count_nonzero(RAND_GEN.random(indiv.genotype.size) < params["Cr"])
 
         if "N" in params:
             params["N"] = round(params["N"])
@@ -170,19 +166,13 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
             new_indiv.genotype = cross_mp(new_indiv.genotype, indiv2.genotype.copy())
 
         elif self.method == IntOpMethods.WEIGHTED_AVG:
-            new_indiv.genotype = weighted_average(
-                new_indiv.genotype, indiv2.genotype.copy(), params["F"]
-            )
+            new_indiv.genotype = weighted_average(new_indiv.genotype, indiv2.genotype.copy(), params["F"])
 
         elif self.method == IntOpMethods.BLXALPHA:
-            new_indiv.genotype = blxalpha(
-                new_indiv.genotype, indiv2.genotype.copy(), params["Cr"]
-            )
+            new_indiv.genotype = blxalpha(new_indiv.genotype, indiv2.genotype.copy(), params["Cr"])
 
         elif self.method == IntOpMethods.MULTICROSS:
-            new_indiv.genotype = multi_cross(
-                new_indiv.genotype, others, params["Nindiv"]
-            )
+            new_indiv.genotype = multi_cross(new_indiv.genotype, others, params["Nindiv"])
 
         elif self.method == IntOpMethods.XOR:
             new_indiv.genotype = xor_mask(new_indiv.genotype, params["N"])
@@ -191,9 +181,7 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
             new_indiv.genotype = xor_cross(new_indiv.genotype, indiv2.genotype.copy())
 
         elif self.method == IntOpMethods.CROSSINTERAVG:
-            new_indiv.genotype = cross_inter_avg(
-                new_indiv.genotype, others, params["N"]
-            )
+            new_indiv.genotype = cross_inter_avg(new_indiv.genotype, others, params["N"])
 
         elif self.method == IntOpMethods.PERM:
             new_indiv.genotype = permutation(new_indiv.genotype, params["N"])
@@ -208,9 +196,7 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
             new_indiv.genotype = cauchy(new_indiv.genotype, params["F"])
 
         elif self.method == IntOpMethods.UNIFORM:
-            new_indiv.genotype = uniform(
-                new_indiv.genotype, params["Low"], params["Up"]
-            )
+            new_indiv.genotype = uniform(new_indiv.genotype, params["Low"], params["Up"])
 
         elif self.method == IntOpMethods.POISSON:
             new_indiv.genotype = poisson(new_indiv.genotype, params["F"])
@@ -231,44 +217,28 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
             new_indiv.genotype = mutate_sample(new_indiv.genotype, others, params)
 
         elif self.method == IntOpMethods.DE_RAND_1:
-            new_indiv.genotype = DE_rand1(
-                new_indiv.genotype, others, params["F"], params["Cr"]
-            )
+            new_indiv.genotype = DE_rand1(new_indiv.genotype, others, params["F"], params["Cr"])
 
         elif self.method == IntOpMethods.DE_BEST_1:
-            new_indiv.genotype = DE_best1(
-                new_indiv.genotype, others, params["F"], params["Cr"]
-            )
+            new_indiv.genotype = DE_best1(new_indiv.genotype, others, params["F"], params["Cr"])
 
         elif self.method == IntOpMethods.DE_RAND_2:
-            new_indiv.genotype = DE_rand2(
-                new_indiv.genotype, others, params["F"], params["Cr"]
-            )
+            new_indiv.genotype = DE_rand2(new_indiv.genotype, others, params["F"], params["Cr"])
 
         elif self.method == IntOpMethods.DE_BEST_2:
-            new_indiv.genotype = DE_best2(
-                new_indiv.genotype, others, params["F"], params["Cr"]
-            )
+            new_indiv.genotype = DE_best2(new_indiv.genotype, others, params["F"], params["Cr"])
 
         elif self.method == IntOpMethods.DE_CTRAND_1:
-            new_indiv.genotype = DE_current_to_rand1(
-                new_indiv.genotype, others, params["F"], params["Cr"]
-            )
+            new_indiv.genotype = DE_current_to_rand1(new_indiv.genotype, others, params["F"], params["Cr"])
 
         elif self.method == IntOpMethods.DE_CTBEST_1:
-            new_indiv.genotype = DE_current_to_best1(
-                new_indiv.genotype, others, params["F"], params["Cr"]
-            )
+            new_indiv.genotype = DE_current_to_best1(new_indiv.genotype, others, params["F"], params["Cr"])
 
         elif self.method == IntOpMethods.DE_CTPBEST_1:
-            new_indiv.genotype = DE_current_to_pbest1(
-                new_indiv.genotype, others, params["F"], params["Cr"], params["P"]
-            )
+            new_indiv.genotype = DE_current_to_pbest1(new_indiv.genotype, others, params["F"], params["Cr"], params["P"])
 
         elif self.method == IntOpMethods.PSO:
-            new_indiv = pso_operator(
-                indiv, others, global_best, params["w"], params["c1"], params["c2"]
-            )
+            new_indiv = pso_operator(indiv, others, global_best, params["w"], params["c1"], params["c2"])
 
         elif self.method == IntOpMethods.FIREFLY:
             new_indiv.genotype = firefly(
@@ -285,14 +255,10 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
             new_indiv = initializer.generate_random(objfunc)
 
         elif self.method == IntOpMethods.RANDOM_MASK:
-            mask_pos = np.hstack(
-                [np.ones(params["N"]), np.zeros(new_indiv.genotype.size - params["N"])]
-            ).astype(bool)
+            mask_pos = np.hstack([np.ones(params["N"]), np.zeros(new_indiv.genotype.size - params["N"])]).astype(bool)
             RAND_GEN.shuffle(mask_pos)
 
-            new_indiv.genotype[mask_pos] = initializer.generate_random(
-                objfunc
-            ).genotype[mask_pos]
+            new_indiv.genotype[mask_pos] = initializer.generate_random(objfunc).genotype[mask_pos]
 
         elif self.method == IntOpMethods.DUMMY:
             new_indiv.genotype = dummy_op(new_indiv.genotype, params["F"])
@@ -301,5 +267,8 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
             fn = params["function"]
             new_indiv.genotype = fn(indiv, population, objfunc, params)
 
+        elif self.method == IntOpMethods.NOTHING:
+            new_indiv = indiv
+
         new_indiv.genotype = np.round(new_indiv.genotype)
         return new_indiv
diff --git a/src/metaheuristic_designer/operators/OperatorLambda.py b/src/metaheuristic_designer/operators/OperatorLambda.py
index 38e3dcb..e2417eb 100644
--- a/src/metaheuristic_designer/operators/OperatorLambda.py
+++ b/src/metaheuristic_designer/operators/OperatorLambda.py
@@ -23,11 +23,9 @@ class OperatorLambda(Operator):
         Name that is associated with the operator.
     """
 
-    def __init__(
-        self, fn: callable, params: Union[ParamScheduler, dict] = None, name: str = None
-    ):
+    def __init__(self, fn: callable, params: Union[ParamScheduler, dict] = None, name: str = None):
         """
-        Constructor for the OperatorReal class
+        Constructor for the OperatorLambda class
         """
 
         self.fn = fn
diff --git a/src/metaheuristic_designer/operators/OperatorList.py b/src/metaheuristic_designer/operators/OperatorList.py
index 027f502..c29c330 100644
--- a/src/metaheuristic_designer/operators/OperatorList.py
+++ b/src/metaheuristic_designer/operators/OperatorList.py
@@ -44,11 +44,9 @@ class OperatorList(Operator):
         Name that is associated with the operator.
     """
 
-    def __init__(
-        self, method: str, params: Union[ParamScheduler, dict] = None, name: str = None
-    ):
+    def __init__(self, method: str, params: Union[ParamScheduler, dict] = None, name: str = None):
         """
-        Constructor for the OperatorReal class
+        Constructor for the OperatorList class
         """
 
         if name is None:
@@ -71,9 +69,11 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
                 params["maxlen"],
                 params["generator"],
             )
+
         elif self.method == ListOpMethods.SHRINK:
-            nex_indiv.genotype = shrink(
-                new_indiv.genotype, params["N"], params["method"]
-            )
+            nex_indiv.genotype = shrink(new_indiv.genotype, params["N"], params["method"])
+
+        elif self.method == ListOpMethods.NOTHING:
+            new_indiv = indiv
 
         return new_indiv
diff --git a/src/metaheuristic_designer/operators/OperatorMeta.py b/src/metaheuristic_designer/operators/OperatorMeta.py
index 15a5634..e3a721b 100644
--- a/src/metaheuristic_designer/operators/OperatorMeta.py
+++ b/src/metaheuristic_designer/operators/OperatorMeta.py
@@ -55,7 +55,7 @@ def __init__(
         name: str = None,
     ):
         """
-        Constructor for the Operator class
+        Constructor for the OperatorMeta class
         """
 
         self.op_list = op_list
@@ -89,19 +89,12 @@ def __init__(
         self.mask = params.get("mask", 0)
 
         # If we have a branch with 2 operators and "p" is given as an input
-        if (
-            self.method == MetaOpMethods.BRANCH
-            and "weights" not in params
-            and "p" in params
-            and len(op_list) == 2
-        ):
+        if self.method == MetaOpMethods.BRANCH and "weights" not in params and "p" in params and len(op_list) == 2:
             params["weights"] = [params["p"], 1 - params["p"]]
 
     def evolve(self, indiv, population, objfunc, global_best, initializer=None):
         if self.method == MetaOpMethods.BRANCH:
-            self.chosen_idx = random.choices(
-                range(len(self.op_list)), k=1, weights=self.params["weights"]
-            )[0]
+            self.chosen_idx = random.choices(range(len(self.op_list)), k=1, weights=self.params["weights"])[0]
             chosen_op = self.op_list[self.chosen_idx]
             result = chosen_op(indiv, population, objfunc, global_best, initializer)
 
@@ -124,16 +117,12 @@ def evolve(self, indiv, population, objfunc, global_best, initializer=None):
             for idx_op, op in enumerate(self.op_list):
                 if np.any(self.mask == idx_op):
                     indiv_copy.genotype = indiv.genotype[self.mask == idx_op]
-                    global_best_copy.genotype = global_best.genotype[
-                        self.mask == idx_op
-                    ]
+                    global_best_copy.genotype = global_best.genotype[self.mask == idx_op]
 
                     for idx_pop, val in enumerate(population_copy):
                         val.genotype = population[idx_pop].genotype[self.mask == idx_op]
 
-                    aux_indiv = op(
-                        indiv_copy, population_copy, objfunc, global_best, initializer
-                    )
+                    aux_indiv = op(indiv_copy, population_copy, objfunc, global_best, initializer)
                     result.genotype[self.mask == idx_op] = aux_indiv.genotype
 
         return result
diff --git a/src/metaheuristic_designer/operators/OperatorNull.py b/src/metaheuristic_designer/operators/OperatorNull.py
new file mode 100644
index 0000000..fcb951d
--- /dev/null
+++ b/src/metaheuristic_designer/operators/OperatorNull.py
@@ -0,0 +1,38 @@
+from __future__ import annotations
+import numpy as np
+from ..Operator import Operator
+from .vector_operator_functions import *
+from ..ParamScheduler import ParamScheduler
+from copy import copy
+import enum
+from enum import Enum
+from ..utils import RAND_GEN
+
+
+class OperatorNull(Operator):
+    """
+    Operator class that returns the individual without changes.
+    Surprisingly useful.
+
+    Parameters
+    ----------
+    fn: callable
+        Function that will be applied when operating on an individual.
+    params: ParamScheduler or dict, optional
+        Dictionary of parameters to define the operator.
+    name: str, optional
+        Name that is associated with the operator.
+    """
+
+    def __init__(self, name: str = None):
+        """
+        Constructor for the OperatorNull class
+        """
+
+        if name is None:
+            name = "Nothing"
+
+        super().__init__({}, name)
+
+    def evolve(self, indiv, *args):
+        return indiv
diff --git a/src/metaheuristic_designer/operators/OperatorPerm.py b/src/metaheuristic_designer/operators/OperatorPerm.py
index f493b11..edcdc75 100644
--- a/src/metaheuristic_designer/operators/OperatorPerm.py
+++ b/src/metaheuristic_designer/operators/OperatorPerm.py
@@ -60,11 +60,9 @@ class OperatorPerm(Operator):
         Name that is associated with the operator.
     """
 
-    def __init__(
-        self, method: str, params: Union[ParamScheduler, dict] = None, name: str = None
-    ):
+    def __init__(self, method: str, params: Union[ParamScheduler, dict] = None, name: str = None):
         """
-        Constructor for the Operator class
+        Constructor for the OperatorPerm class
         """
 
         if name is None:
@@ -91,9 +89,7 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
         params = copy(self.params)
 
         if "Cr" in params and "N" not in params:
-            params["N"] = np.count_nonzero(
-                RAND_GEN.random(indiv.genotype.size) < params["Cr"]
-            )
+            params["N"] = np.count_nonzero(RAND_GEN.random(indiv.genotype.size) < params["Cr"])
 
         if "N" in params:
             params["N"] = round(params["N"])
@@ -130,4 +126,7 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
             fn = params["function"]
             new_indiv.genotype = fn(indiv, population, objfunc, params)
 
+        elif self.method == PermOpMethods.NOTHING:
+            new_indiv = indiv
+
         return new_indiv
diff --git a/src/metaheuristic_designer/operators/OperatorReal.py b/src/metaheuristic_designer/operators/OperatorReal.py
index 9a85dac..442d089 100644
--- a/src/metaheuristic_designer/operators/OperatorReal.py
+++ b/src/metaheuristic_designer/operators/OperatorReal.py
@@ -108,11 +108,9 @@ class OperatorReal(Operator):
         Name that is associated with the operator.
     """
 
-    def __init__(
-        self, method: str, params: Union[ParamScheduler, dict] = None, name: str = None
-    ):
+    def __init__(self, method: str, params: Union[ParamScheduler, dict] = None, name: str = None):
         """
-        Constructor for the OperatorReal class
+        Constructor for the OperatorPerm class
         """
 
         if name is None:
@@ -147,9 +145,7 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
         params = copy(self.params)
 
         if "Cr" in params and "N" not in params:
-            params["N"] = np.count_nonzero(
-                RAND_GEN.random(indiv.genotype.size) < params["Cr"]
-            )
+            params["N"] = np.count_nonzero(RAND_GEN.random(indiv.genotype.size) < params["Cr"])
 
         if "N" in params:
             params["N"] = round(params["N"])
@@ -165,34 +161,22 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
             new_indiv.genotype = cross_mp(new_indiv.genotype, indiv2.genotype.copy())
 
         elif self.method == RealOpMethods.WEIGHTED_AVG:
-            new_indiv.genotype = weighted_average(
-                new_indiv.genotype, indiv2.genotype.copy(), params["F"]
-            )
+            new_indiv.genotype = weighted_average(new_indiv.genotype, indiv2.genotype.copy(), params["F"])
 
         elif self.method == RealOpMethods.BLXALPHA:
-            new_indiv.genotype = blxalpha(
-                new_indiv.genotype, indiv2.genotype.copy(), params["Cr"]
-            )
+            new_indiv.genotype = blxalpha(new_indiv.genotype, indiv2.genotype.copy(), params["Cr"])
 
         elif self.method == RealOpMethods.SBX:
-            new_indiv.genotype = sbx(
-                new_indiv.genotype, indiv2.genotype.copy(), params["Cr"]
-            )
+            new_indiv.genotype = sbx(new_indiv.genotype, indiv2.genotype.copy(), params["Cr"])
 
         elif self.method == RealOpMethods.MULTICROSS:
-            new_indiv.genotype = multi_cross(
-                new_indiv.genotype, others, params["Nindiv"]
-            )
+            new_indiv.genotype = multi_cross(new_indiv.genotype, others, params["Nindiv"])
 
         elif self.method == RealOpMethods.CROSSINTERAVG:
-            new_indiv.genotype = cross_inter_avg(
-                new_indiv.genotype, others, params["N"]
-            )
+            new_indiv.genotype = cross_inter_avg(new_indiv.genotype, others, params["N"])
 
         elif self.method == RealOpMethods.MUTATE1SIGMA:
-            new_indiv.genotype = mutate_1_sigma(
-                new_indiv.genotype[0], params["epsilon"], params["tau"]
-            )
+            new_indiv.genotype = mutate_1_sigma(new_indiv.genotype[0], params["epsilon"], params["tau"])
 
         elif self.method == RealOpMethods.MUTATENSIGMAS:
             new_indiv.genotype = mutate_n_sigmas(
@@ -203,9 +187,7 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
             )
 
         elif self.method == RealOpMethods.SAMPLESIGMA:
-            new_indiv.genotype = sample_1_sigma(
-                new_indiv.genotype, params["N"], params["epsilon"], params["tau"]
-            )
+            new_indiv.genotype = sample_1_sigma(new_indiv.genotype, params["N"], params["epsilon"], params["tau"])
 
         elif self.method == RealOpMethods.PERM:
             new_indiv.genotype = permutation(new_indiv.genotype, params["N"])
@@ -220,9 +202,7 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
             new_indiv.genotype = cauchy(new_indiv.genotype, params["F"])
 
         elif self.method == RealOpMethods.UNIFORM:
-            new_indiv.genotype = uniform(
-                new_indiv.genotype, params["Low"], params["Up"]
-            )
+            new_indiv.genotype = uniform(new_indiv.genotype, params["Low"], params["Up"])
 
         elif self.method == RealOpMethods.MUTNOISE:
             new_indiv.genotype = mutate_rand(new_indiv.genotype, others, params)
@@ -237,44 +217,28 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
             new_indiv.genotype = rand_sample(new_indiv.genotype, others, params)
 
         elif self.method == RealOpMethods.DE_RAND_1:
-            new_indiv.genotype = DE_rand1(
-                new_indiv.genotype, others, params["F"], params["Cr"]
-            )
+            new_indiv.genotype = DE_rand1(new_indiv.genotype, others, params["F"], params["Cr"])
 
         elif self.method == RealOpMethods.DE_BEST_1:
-            new_indiv.genotype = DE_best1(
-                new_indiv.genotype, others, params["F"], params["Cr"]
-            )
+            new_indiv.genotype = DE_best1(new_indiv.genotype, others, params["F"], params["Cr"])
 
         elif self.method == RealOpMethods.DE_RAND_2:
-            new_indiv.genotype = DE_rand2(
-                new_indiv.genotype, others, params["F"], params["Cr"]
-            )
+            new_indiv.genotype = DE_rand2(new_indiv.genotype, others, params["F"], params["Cr"])
 
         elif self.method == RealOpMethods.DE_BEST_2:
-            new_indiv.genotype = DE_best2(
-                new_indiv.genotype, others, params["F"], params["Cr"]
-            )
+            new_indiv.genotype = DE_best2(new_indiv.genotype, others, params["F"], params["Cr"])
 
         elif self.method == RealOpMethods.DE_CTRAND_1:
-            new_indiv.genotype = DE_current_to_rand1(
-                new_indiv.genotype, others, params["F"], params["Cr"]
-            )
+            new_indiv.genotype = DE_current_to_rand1(new_indiv.genotype, others, params["F"], params["Cr"])
 
         elif self.method == RealOpMethods.DE_CTBEST_1:
-            new_indiv.genotype = DE_current_to_best1(
-                new_indiv.genotype, others, params["F"], params["Cr"]
-            )
+            new_indiv.genotype = DE_current_to_best1(new_indiv.genotype, others, params["F"], params["Cr"])
 
         elif self.method == RealOpMethods.DE_CTPBEST_1:
-            new_indiv.genotype = DE_current_to_pbest1(
-                new_indiv.genotype, others, params["F"], params["Cr"], params["P"]
-            )
+            new_indiv.genotype = DE_current_to_pbest1(new_indiv.genotype, others, params["F"], params["Cr"], params["P"])
 
         elif self.method == RealOpMethods.PSO:
-            new_indiv = pso_operator(
-                indiv, others, global_best, params["w"], params["c1"], params["c2"]
-            )
+            new_indiv = pso_operator(indiv, others, global_best, params["w"], params["c1"], params["c2"])
 
         elif self.method == RealOpMethods.FIREFLY:
             new_indiv.genotype = firefly(
@@ -291,14 +255,10 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
             new_indiv = initializer.generate_random(objfunc)
 
         elif self.method == RealOpMethods.RANDOM_MASK:
-            mask_pos = np.hstack(
-                [np.ones(params["N"]), np.zeros(new_indiv.genotype.size - params["N"])]
-            ).astype(bool)
+            mask_pos = np.hstack([np.ones(params["N"]), np.zeros(new_indiv.genotype.size - params["N"])]).astype(bool)
             RAND_GEN.shuffle(mask_pos)
 
-            new_indiv.genotype[mask_pos] = initializer.generate_random(
-                objfunc
-            ).genotype[mask_pos]
+            new_indiv.genotype[mask_pos] = initializer.generate_random(objfunc).genotype[mask_pos]
 
         elif self.method == RealOpMethods.DUMMY:
             new_indiv.genotype = dummy_op(new_indiv.genotype, params["F"])
@@ -307,4 +267,7 @@ def evolve(self, indiv, population, objfunc, global_best, initializer):
             fn = params["function"]
             new_indiv.genotype = fn(indiv, population, objfunc, params)
 
+        elif self.method == RealOpMethods.NOTHING:
+            new_indiv = indiv
+
         return new_indiv
diff --git a/src/metaheuristic_designer/operators/__init__.py b/src/metaheuristic_designer/operators/__init__.py
index 2183127..57e9ebc 100644
--- a/src/metaheuristic_designer/operators/__init__.py
+++ b/src/metaheuristic_designer/operators/__init__.py
@@ -4,4 +4,6 @@
 from .OperatorBinary import OperatorBinary, BinOpMethods, bin_ops_map
 from .OperatorMeta import OperatorMeta, MetaOpMethods, meta_ops_map
 from .OperatorLambda import OperatorLambda
+from .OperatorNull import OperatorNull
+from .OperatorAdaptative import OperatorAdaptative
 from .vector_operator_functions import *
diff --git a/src/metaheuristic_designer/operators/list_operator_functions.py b/src/metaheuristic_designer/operators/list_operator_functions.py
index f6dd0e5..ef9344c 100644
--- a/src/metaheuristic_designer/operators/list_operator_functions.py
+++ b/src/metaheuristic_designer/operators/list_operator_functions.py
@@ -24,9 +24,7 @@ def expand(input_list, n, method, maxlen, val_generator=None):
 def expand_rand(input_list, n, val_generator):
     for i in range(n):
         inject_idx = random.randrange(len(input_list))
-        input_list = (
-            input_list[:inject_idx] + [val_generator()] + input_list[inject_idx:]
-        )
+        input_list = input_list[:inject_idx] + [val_generator()] + input_list[inject_idx:]
 
     return input_list + new_values
 
diff --git a/src/metaheuristic_designer/operators/vector_operator_functions.py b/src/metaheuristic_designer/operators/vector_operator_functions.py
index 7e9662c..229f89a 100644
--- a/src/metaheuristic_designer/operators/vector_operator_functions.py
+++ b/src/metaheuristic_designer/operators/vector_operator_functions.py
@@ -74,9 +74,7 @@ def mutate_sample(vector, population, params):
     RAND_GEN.shuffle(mask_pos)
     popul_matrix = np.vstack([i.genotype for i in population])
     mean = popul_matrix.mean(axis=0)[mask_pos]
-    std = (
-        popul_matrix.std(axis=0)[mask_pos] + 1e-6
-    ) * strength  # ensure there will be some standard deviation
+    std = (popul_matrix.std(axis=0)[mask_pos] + 1e-6) * strength  # ensure there will be some standard deviation
 
     rand_vec = sample_distribution(method, n, mean, std, low, up)
 
@@ -97,9 +95,7 @@ def rand_sample(vector, population, params):
 
     popul_matrix = np.vstack([i.genotype for i in population])
     mean = popul_matrix.mean(axis=0)
-    std = (
-        popul_matrix.std(axis=0) + 1e-6
-    ) * strength  # ensure there will be some standard deviation
+    std = (popul_matrix.std(axis=0) + 1e-6) * strength  # ensure there will be some standard deviation
 
     rand_vec = sample_distribution(method, vector.shape, mean, std, low, up)
 
@@ -205,20 +201,15 @@ def mutate_1_sigma(sigma, epsilon, tau):
     Mutate a sigma value in base of tau param, where epsilon is de minimum value that a sigma can have.
     """
 
-    return max(epsilon, np.exp(tau * RAND_GEN.normal()))
+    return np.maximum(epsilon, np.exp(tau * RAND_GEN.normal()))
 
 
-def mutate_n_sigmas(list_sigmas, epsilon, tau, tau_multiple):
+def mutate_n_sigmas(sigmas, epsilon, tau, tau_multiple):
     """
     Mutate a list of sigmas values in base of tau and tau_multiple params, where epsilon is de minimum value that a sigma can have.
     """
 
-    base_tau = tau * RAND_GEN.normal()
-    new_sigmas = [
-        max(epsilon, sigma * np.exp(base_tau + tau_multiple * RAND_GEN.normal()))
-        for sigma in list_sigmas
-    ]
-    return new_sigmas
+    return np.maximum(epsilon, sigmas*np.exp(tau * RAND_GEN.normal() + tau_multiple*RAND_GEN.normal()))
 
 
 def xor_mask(vector, n, mode="byte"):
@@ -481,11 +472,7 @@ def DE_rand2(vector, population, F, CR):
     if len(population) > 5:
         r1, r2, r3, r4, r5 = random.sample(population, 5)
 
-        v = (
-            r1.genotype
-            + F * (r2.genotype - r3.genotype)
-            + F * (r4.genotype - r5.genotype)
-        )
+        v = r1.genotype + F * (r2.genotype - r3.genotype) + F * (r4.genotype - r5.genotype)
         mask_pos = RAND_GEN.random(vector.shape) <= CR
         vector[mask_pos] = v[mask_pos]
     return vector
@@ -501,11 +488,7 @@ def DE_best2(vector, population, F, CR):
         best = population[fitness.index(max(fitness))]
         r1, r2, r3, r4 = random.sample(population, 4)
 
-        v = (
-            best.genotype
-            + F * (r1.genotype - r2.genotype)
-            + F * (r3.genotype - r4.genotype)
-        )
+        v = best.genotype + F * (r1.genotype - r2.genotype) + F * (r3.genotype - r4.genotype)
         mask_pos = RAND_GEN.random(vector.shape) <= CR
         vector[mask_pos] = v[mask_pos]
     return vector
@@ -519,11 +502,7 @@ def DE_current_to_rand1(vector, population, F, CR):
     if len(population) > 3:
         r1, r2, r3 = random.sample(population, 3)
 
-        v = (
-            vector
-            + RAND_GEN.random() * (r1.genotype - vector)
-            + F * (r2.genotype - r3.genotype)
-        )
+        v = vector + RAND_GEN.random() * (r1.genotype - vector) + F * (r2.genotype - r3.genotype)
         mask_pos = RAND_GEN.random(vector.shape) <= CR
         vector[mask_pos] = v[mask_pos]
     return vector
@@ -571,11 +550,7 @@ def pso_operator(indiv, population, global_best, w, c1, c2):
     c1 = c1 * RAND_GEN.random(indiv.genotype.shape)
     c2 = c2 * RAND_GEN.random(indiv.genotype.shape)
 
-    indiv.speed = (
-        w * indiv.speed
-        + c1 * (indiv.best - indiv.genotype)
-        + c2 * (global_best.genotype - indiv.genotype)
-    )
+    indiv.speed = w * indiv.speed + c1 * (indiv.best - indiv.genotype) + c2 * (global_best.genotype - indiv.genotype)
     return indiv.apply_speed()
 
 
@@ -592,12 +567,7 @@ def firefly(solution, population, objfunc, alpha_0, beta_0, delta, gamma):
             r = np.linalg.norm(solution.genotype - ind.genotype)
             alpha = alpha_0 * delta**idx
             beta = beta_0 * np.exp(-gamma * (r / (sol_range * np.sqrt(n_dim))) ** 2)
-            new_vector = (
-                new_vector
-                + beta * (ind.genotype - new_vector)
-                + alpha * sol_range * random.random()
-                - 0.5
-            )
+            new_vector = new_vector + beta * (ind.genotype - new_vector) + alpha * sol_range * random.random() - 0.5
             new_vector = objfunc.repair_solution(new_vector)
 
     return new_vector
diff --git a/src/metaheuristic_designer/selectionMethods/ParentSelection.py b/src/metaheuristic_designer/selectionMethods/ParentSelection.py
index b97d0d6..c0e1582 100644
--- a/src/metaheuristic_designer/selectionMethods/ParentSelection.py
+++ b/src/metaheuristic_designer/selectionMethods/ParentSelection.py
@@ -38,6 +38,8 @@ class ParentSelection(SelectionMethod):
     """
     Parent selection methods.
 
+    Selects the individuals that will be perturbed in this generation.
+
     Parameters
     ----------
     method: str
@@ -73,16 +75,12 @@ def __init__(
             if "F" not in self.params:
                 self.params["F"] = None
 
-    def select(
-        self, population: List[Individual], offsping: List[Individual] = None
-    ) -> List[Individual]:
+    def select(self, population: List[Individual], offsping: List[Individual] = None) -> List[Individual]:
         population = population.copy()
         parents = []
 
         if self.method == ParentSelMethod.TOURNAMENT:
-            parents = prob_tournament(
-                population, self.params["amount"], self.params["p"]
-            )
+            parents = prob_tournament(population, self.params["amount"], self.params["p"])
 
         elif self.method == ParentSelMethod.BEST:
             parents = select_best(population, self.params["amount"])
diff --git a/src/metaheuristic_designer/selectionMethods/SurvivorSelection.py b/src/metaheuristic_designer/selectionMethods/SurvivorSelection.py
index 82b5843..f94253c 100644
--- a/src/metaheuristic_designer/selectionMethods/SurvivorSelection.py
+++ b/src/metaheuristic_designer/selectionMethods/SurvivorSelection.py
@@ -45,7 +45,9 @@ def from_str(str_input):
 
 class SurvivorSelection(SelectionMethod):
     """
-    Survivor selection methods
+    Survivor selection methods.
+
+    Selects the individuals that will remain for the next generation of our algorithm.
 
     Parameters
     ----------
@@ -77,9 +79,7 @@ def __init__(
 
         self.method = SurvSelMethod.from_str(method)
 
-    def select(
-        self, popul: List[Individual], offspring: List[Individual]
-    ) -> List[Individual]:
+    def select(self, popul: List[Individual], offspring: List[Individual]) -> List[Individual]:
         result = []
         if self.method == SurvSelMethod.ELITISM:
             result = elitism(popul, offspring, self.params["amount"])
diff --git a/src/metaheuristic_designer/selectionMethods/survivor_selection_functions.py b/src/metaheuristic_designer/selectionMethods/survivor_selection_functions.py
index fe25c7c..ae3b8dd 100644
--- a/src/metaheuristic_designer/selectionMethods/survivor_selection_functions.py
+++ b/src/metaheuristic_designer/selectionMethods/survivor_selection_functions.py
@@ -75,7 +75,7 @@ def prob_one_to_one(popul, offspring, p):
 
     new_population = []
     for parent, child in zip(popul, offspring):
-        if random.random() < p and child.fitness > parent.fitness:
+        if random.random() < p or child.fitness > parent.fitness:
             new_population.append(child)
         else:
             new_population.append(parent)
@@ -106,9 +106,7 @@ def elitism(popul, offspring, amount):
     survivors: List[Individual]
         The individuals selected for the next generation.
     """
-    selected_offspring = sorted(offspring, reverse=True, key=lambda x: x.fitness)[
-        : len(popul) - amount
-    ]
+    selected_offspring = sorted(offspring, reverse=True, key=lambda x: x.fitness)[: len(popul) - amount]
     best_parents = sorted(popul, reverse=True, key=lambda x: x.fitness)[:amount]
 
     return best_parents + selected_offspring
@@ -132,9 +130,7 @@ def cond_elitism(popul, offspring, amount):
     """
 
     best_parents = sorted(popul, reverse=True, key=lambda x: x.fitness)[:amount]
-    new_offspring = sorted(offspring, reverse=True, key=lambda x: x.fitness)[
-        : len(popul)
-    ]
+    new_offspring = sorted(offspring, reverse=True, key=lambda x: x.fitness)[: len(popul)]
     best_offspring = new_offspring[:amount]
 
     for idx, val in enumerate(best_parents):
diff --git a/src/metaheuristic_designer/simple/differential_evolution.py b/src/metaheuristic_designer/simple/differential_evolution.py
index 1f20838..2af1ff8 100644
--- a/src/metaheuristic_designer/simple/differential_evolution.py
+++ b/src/metaheuristic_designer/simple/differential_evolution.py
@@ -24,9 +24,7 @@ def differential_evolution(objfunc: ObjectiveVectorFunc, params: dict) -> Algori
     """
 
     if "encoding" not in params:
-        raise ValueError(
-            f'You must specify the encoding in the params structure, the options are "real", "int" and "bin"'
-        )
+        raise ValueError(f'You must specify the encoding in the params structure, the options are "real", "int" and "bin"')
 
     encoding_str = params["encoding"]
 
@@ -37,9 +35,7 @@ def differential_evolution(objfunc: ObjectiveVectorFunc, params: dict) -> Algori
     elif encoding_str.lower() == "bin":
         alg = _differential_evolution_bin_vec(objfunc, params)
     else:
-        raise ValueError(
-            f'The encoding "{encoding_str}" does not exist, try "real", "int" or "bin"'
-        )
+        raise ValueError(f'The encoding "{encoding_str}" does not exist, try "real", "int" or "bin"')
 
     return alg
 
@@ -66,9 +62,7 @@ def _differential_evolution_real_vec(objfunc, params):
     ]:
         raise ValueError(f'Differential evolution strategy "{de_type}" does not exist.')
 
-    pop_initializer = UniformVectorInitializer(
-        objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=float
-    )
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=float)
 
     de_op = OperatorReal(de_type, {"F": f, "Cr": cr})
 
diff --git a/src/metaheuristic_designer/simple/evolution_strategy.py b/src/metaheuristic_designer/simple/evolution_strategy.py
index 9ef426c..613754a 100644
--- a/src/metaheuristic_designer/simple/evolution_strategy.py
+++ b/src/metaheuristic_designer/simple/evolution_strategy.py
@@ -25,9 +25,7 @@ def evolution_strategy(objfunc: ObjectiveVectorFunc, params: dict) -> Algorithm:
     """
 
     if "encoding" not in params:
-        raise ValueError(
-            f'You must specify the encoding in the params structure, the options are "real", "int" and "bin"'
-        )
+        raise ValueError(f'You must specify the encoding in the params structure, the options are "real", "int" and "bin"')
 
     encoding_str = params["encoding"]
 
@@ -38,9 +36,7 @@ def evolution_strategy(objfunc: ObjectiveVectorFunc, params: dict) -> Algorithm:
     elif encoding_str.lower() == "real":
         alg = _evolution_strategy_real_vec(objfunc, params)
     else:
-        raise ValueError(
-            f'The encoding "{encoding_str}" does not exist, try "real", "int" or "bin"'
-        )
+        raise ValueError(f'The encoding "{encoding_str}" does not exist, try "real", "int" or "bin"')
 
     return alg
 
@@ -58,9 +54,7 @@ def _evolution_strategy_bin_vec(objfunc, params):
 
     encoding = TypeCastEncoding(int, bool)
 
-    pop_initializer = UniformVectorInitializer(
-        objfunc.vecsize, 0, 1, pop_size=pop_size, dtype=int, encoding=encoding
-    )
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, 0, 1, pop_size=pop_size, dtype=int, encoding=encoding)
 
     cross_op = OperatorBinary("Nothing")
     mutation_op = OperatorBinary("Flip", {"N": mutstr})
@@ -90,9 +84,7 @@ def _evolution_strategy_int_vec(objfunc, params):
     offspring_size = params.get("offspring_size", 150)
     mutstr = params.get("mutstr", 1)
 
-    pop_initializer = UniformVectorInitializer(
-        objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=int
-    )
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=int)
 
     cross_op = OperatorInt("Nothing")
     mutation_op = OperatorInt(
@@ -130,9 +122,7 @@ def _evolution_strategy_real_vec(objfunc, params):
     offspring_size = params.get("offspring_size", 150)
     mutstr = params.get("mutstr", 1e-5)
 
-    pop_initializer = UniformVectorInitializer(
-        objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=float
-    )
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=float)
 
     cross_op = OperatorReal("Nothing")
     mutation_op = OperatorReal("RandNoise", {"method": "Gauss", "F": mutstr})
diff --git a/src/metaheuristic_designer/simple/genetic_algorithm.py b/src/metaheuristic_designer/simple/genetic_algorithm.py
index 4165b5f..d964636 100644
--- a/src/metaheuristic_designer/simple/genetic_algorithm.py
+++ b/src/metaheuristic_designer/simple/genetic_algorithm.py
@@ -25,9 +25,7 @@ def genetic_algorithm(objfunc: ObjectiveVectorFunc, params: dict) -> Algorithm:
     """
 
     if "encoding" not in params:
-        raise ValueError(
-            f'You must specify the encoding in the params structure, the options are "real", "int" and "bin"'
-        )
+        raise ValueError(f'You must specify the encoding in the params structure, the options are "real", "int" and "bin"')
 
     encoding_str = params["encoding"]
 
@@ -38,9 +36,7 @@ def genetic_algorithm(objfunc: ObjectiveVectorFunc, params: dict) -> Algorithm:
     elif encoding_str.lower() == "real":
         alg = _genetic_algorithm_real_vec(objfunc, params)
     else:
-        raise ValueError(
-            f'The encoding "{encoding_str}" does not exist, try "real", "int" or "bin"'
-        )
+        raise ValueError(f'The encoding "{encoding_str}" does not exist, try "real", "int" or "bin"')
 
     return alg
 
@@ -60,9 +56,7 @@ def _genetic_algorithm_bin_vec(objfunc, params):
 
     encoding = TypeCastEncoding(int, bool)
 
-    pop_initializer = UniformVectorInitializer(
-        objfunc.vecsize, 0, 1, pop_size=pop_size, dtype=int, encoding=encoding
-    )
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, 0, 1, pop_size=pop_size, dtype=int, encoding=encoding)
 
     cross_op = OperatorBinary(cross_method)
     mutation_op = OperatorBinary("Flip", {"N": mutstr})
@@ -95,9 +89,7 @@ def _genetic_algorithm_int_vec(objfunc, params):
     pmut = params.get("pmut", 0.1)
     mutstr = params.get("mutstr", 1)
 
-    pop_initializer = UniformVectorInitializer(
-        objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=int
-    )
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=int)
 
     cross_op = OperatorInt(cross_method)
     mutation_op = OperatorInt(
@@ -138,9 +130,7 @@ def _genetic_algorithm_real_vec(objfunc, params):
     pmut = params.get("pmut", 0.1)
     mutstr = params.get("mutstr", 1e-5)
 
-    pop_initializer = UniformVectorInitializer(
-        objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=float
-    )
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=float)
 
     cross_op = OperatorReal(cross_method)
     mutation_op = OperatorReal("RandNoise", {"method": "Gauss", "F": mutstr, "N": 1})
diff --git a/src/metaheuristic_designer/simple/hill_climb.py b/src/metaheuristic_designer/simple/hill_climb.py
index 49e7b17..2a3a140 100644
--- a/src/metaheuristic_designer/simple/hill_climb.py
+++ b/src/metaheuristic_designer/simple/hill_climb.py
@@ -24,9 +24,7 @@ def hill_climb(objfunc: ObjectiveVectorFunc, params: dict) -> Algorithm:
     """
 
     if "encoding" not in params:
-        raise ValueError(
-            f'You must specify the encoding in the params structure, the options are "real", "int" and "bin"'
-        )
+        raise ValueError(f'You must specify the encoding in the params structure, the options are "real", "int" and "bin"')
 
     encoding_str = params["encoding"]
 
@@ -37,9 +35,7 @@ def hill_climb(objfunc: ObjectiveVectorFunc, params: dict) -> Algorithm:
     elif encoding_str.lower() == "real":
         alg = _hill_climb_real_vec(objfunc, params)
     else:
-        raise ValueError(
-            f'The encoding "{encoding_str}" does not exist, try "real", "int" or "bin"'
-        )
+        raise ValueError(f'The encoding "{encoding_str}" does not exist, try "real", "int" or "bin"')
 
     return alg
 
@@ -54,9 +50,7 @@ def _hill_climb_bin_vec(objfunc, params):
 
     encoding = TypeCastEncoding(int, bool)
 
-    pop_initializer = UniformVectorInitializer(
-        objfunc.vecsize, 0, 1, pop_size=1, dtype=int, encoding=encoding
-    )
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, 0, 1, pop_size=1, dtype=int, encoding=encoding)
 
     mutation_op = OperatorBinary("Flip", {"N": mutstr})
 
@@ -73,9 +67,7 @@ def _hill_climb_int_vec(objfunc, params):
 
     mutstr = params.get("mutstr", 1)
 
-    pop_initializer = UniformVectorInitializer(
-        objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=1, dtype=int
-    )
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=1, dtype=int)
 
     mutation_op = OperatorInt(
         "MutRand",
@@ -100,9 +92,7 @@ def _hill_climb_real_vec(objfunc, params):
 
     mutstr = params.get("mutstr", 1e-5)
 
-    pop_initializer = UniformVectorInitializer(
-        objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=1, dtype=float
-    )
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=1, dtype=float)
 
     mutation_op = OperatorReal("RandNoise", {"method": "Gauss", "F": mutstr})
 
diff --git a/src/metaheuristic_designer/simple/particle_swarm.py b/src/metaheuristic_designer/simple/particle_swarm.py
index c1f1f12..71ab54e 100644
--- a/src/metaheuristic_designer/simple/particle_swarm.py
+++ b/src/metaheuristic_designer/simple/particle_swarm.py
@@ -23,9 +23,7 @@ def particle_swarm(objfunc: ObjectiveVectorFunc, params: dict) -> Algorithm:
     """
 
     if "encoding" not in params:
-        raise ValueError(
-            f'You must specify the encoding in the params structure, the options are "real", "int" and "bin"'
-        )
+        raise ValueError(f'You must specify the encoding in the params structure, the options are "real", "int" and "bin"')
 
     encoding_str = params["encoding"]
 
@@ -36,9 +34,7 @@ def particle_swarm(objfunc: ObjectiveVectorFunc, params: dict) -> Algorithm:
     elif encoding_str.lower() == "bin":
         alg = _particle_swarm_bin_vec(objfunc, params)
     else:
-        raise ValueError(
-            f'The encoding "{encoding_str}" does not exist, try "real", "int" or "bin"'
-        )
+        raise ValueError(f'The encoding "{encoding_str}" does not exist, try "real", "int" or "bin"')
 
     return alg
 
@@ -54,9 +50,7 @@ def _particle_swarm_real_vec(objfunc, params):
     c1 = params.get("c1", 1.5)
     c2 = params.get("c2", 1.5)
 
-    pop_initializer = UniformVectorInitializer(
-        objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=float
-    )
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=float)
 
     search_strat = PSO(pop_initializer, {"w": w, "c1": c1, "c2": c2})
 
diff --git a/src/metaheuristic_designer/simple/random_search.py b/src/metaheuristic_designer/simple/random_search.py
index 422e026..c98fc61 100644
--- a/src/metaheuristic_designer/simple/random_search.py
+++ b/src/metaheuristic_designer/simple/random_search.py
@@ -23,9 +23,7 @@ def random_search(objfunc: ObjectiveVectorFunc, params: dict) -> Algorithm:
     """
 
     if "encoding" not in params:
-        raise ValueError(
-            f'You must specify the encoding in the params structure, the options are "real", "int" and "bin"'
-        )
+        raise ValueError(f'You must specify the encoding in the params structure, the options are "real", "int" and "bin"')
 
     encoding_str = params["encoding"]
 
@@ -36,9 +34,7 @@ def random_search(objfunc: ObjectiveVectorFunc, params: dict) -> Algorithm:
     elif encoding_str.lower() == "real":
         alg = _random_search_real_vec(objfunc, params)
     else:
-        raise ValueError(
-            f'The encoding "{encoding_str}" does not exist, try "real", "int" or "bin"'
-        )
+        raise ValueError(f'The encoding "{encoding_str}" does not exist, try "real", "int" or "bin"')
 
     return alg
 
@@ -51,9 +47,7 @@ def _random_search_bin_vec(objfunc, params):
 
     encoding = TypeCastEncoding(int, bool)
 
-    pop_initializer = UniformVectorInitializer(
-        objfunc.vecsize, 0, 1, pop_size=1, dtype=int, encoding=encoding
-    )
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, 0, 1, pop_size=1, dtype=int, encoding=encoding)
 
     search_strat = RandomSearch(pop_initializer)
 
@@ -66,9 +60,7 @@ def _random_search_int_vec(objfunc, params):
     This objective function should accept integer coded vectors.
     """
 
-    pop_initializer = UniformVectorInitializer(
-        objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=int
-    )
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=int)
 
     search_strat = RandomSearch(pop_initializer)
 
@@ -87,9 +79,7 @@ def _random_search_real_vec(objfunc, params):
     pmut = params.get("pmut", 0.1)
     mutstr = params.get("mutstr", 1e-5)
 
-    pop_initializer = UniformVectorInitializer(
-        objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=float
-    )
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=float)
 
     search_strat = RandomSearch(pop_initializer)
 
diff --git a/src/metaheuristic_designer/simple/simulated_annealing.py b/src/metaheuristic_designer/simple/simulated_annealing.py
index 10bc9f1..a7110fd 100644
--- a/src/metaheuristic_designer/simple/simulated_annealing.py
+++ b/src/metaheuristic_designer/simple/simulated_annealing.py
@@ -24,9 +24,7 @@ def simulated_annealing(objfunc: ObjectiveVectorFunc, params: dict) -> Algorithm
     """
 
     if "encoding" not in params:
-        raise ValueError(
-            f'You must specify the encoding in the params structure, the options are "real", "int" and "bin"'
-        )
+        raise ValueError(f'You must specify the encoding in the params structure, the options are "real", "int" and "bin"')
 
     encoding_str = params["encoding"]
 
@@ -37,9 +35,7 @@ def simulated_annealing(objfunc: ObjectiveVectorFunc, params: dict) -> Algorithm
     elif encoding_str.lower() == "real":
         alg = _simulated_annealing_real_vec(objfunc, params)
     else:
-        raise ValueError(
-            f'The encoding "{encoding_str}" does not exist, try "real", "int" or "bin"'
-        )
+        raise ValueError(f'The encoding "{encoding_str}" does not exist, try "real", "int" or "bin"')
 
     return alg
 
@@ -57,9 +53,7 @@ def _simulated_annealing_bin_vec(objfunc, params):
 
     encoding = TypeCastEncoding(int, bool)
 
-    pop_initializer = UniformVectorInitializer(
-        objfunc.vecsize, 0, 1, pop_size=1, dtype=int, encoding=encoding
-    )
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, 0, 1, pop_size=1, dtype=int, encoding=encoding)
 
     mutation_op = OperatorBinary("Flip", {"N": mutstr})
 
@@ -83,9 +77,7 @@ def _simulated_annealing_int_vec(objfunc, params):
     alpha = params.get("alpha", 0.99)
     mutstr = params.get("mutstr", 1)
 
-    pop_initializer = UniformVectorInitializer(
-        objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=1, dtype=int
-    )
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=1, dtype=int)
 
     mutation_op = OperatorInt(
         "MutRand",
@@ -117,9 +109,7 @@ def _simulated_annealing_real_vec(objfunc, params):
     alpha = params.get("alpha", 0.99)
     mutstr = params.get("mutstr", 1e-5)
 
-    pop_initializer = UniformVectorInitializer(
-        objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=1, dtype=float
-    )
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=1, dtype=float)
 
     mutation_op = OperatorReal("RandNoise", {"method": "Gauss", "F": mutstr})
 
diff --git a/src/metaheuristic_designer/strategies/CRO/CRO.py b/src/metaheuristic_designer/strategies/CRO/CRO.py
index 3b730ec..8964884 100644
--- a/src/metaheuristic_designer/strategies/CRO/CRO.py
+++ b/src/metaheuristic_designer/strategies/CRO/CRO.py
@@ -35,6 +35,4 @@ def __init__(
         pop_init = deepcopy(pop_init)
         pop_init.pop_size = round(pop_init.pop_size * params["rho"])
 
-        super().__init__(
-            pop_init, evolve_op, params=params, selection_op=selection_op, name=name
-        )
+        super().__init__(pop_init, evolve_op, params=params, selection_op=selection_op, name=name)
diff --git a/src/metaheuristic_designer/strategies/CRO/DPCRO_SL.py b/src/metaheuristic_designer/strategies/CRO/DPCRO_SL.py
index eda470d..5511308 100644
--- a/src/metaheuristic_designer/strategies/CRO/DPCRO_SL.py
+++ b/src/metaheuristic_designer/strategies/CRO/DPCRO_SL.py
@@ -31,9 +31,7 @@ def __init__(
         self.dyn_steps = params["dyn_steps"]
         self.prob_amp = params["prob_amp"]
 
-        self.operator_idx = random.choices(
-            range(len(self.operator_list)), k=self.maxpopsize
-        )
+        self.operator_idx = random.choices(range(len(self.operator_list)), k=self.maxpopsize)
         self.operator_weight = [1 / len(operator_list)] * len(operator_list)
 
         if self.dyn_method == "success":
@@ -50,72 +48,6 @@ def __init__(
         self.operator_metric = [0] * len(operator_list)
         self.operator_history = []
 
-    def perturb(self, parent_list, objfunc, **kwargs):
-        offspring = []
-
-        divided_population = [[] for i in self.operator_list]
-        for idx, indiv in enumerate(parent_list):
-            op_idx = self.operator_idx[idx]
-            divided_population[op_idx].append(indiv)
-
-        for idx, indiv in enumerate(parent_list):
-            # Select operator
-            op_idx = self.operator_idx[idx]
-
-            op = self.operator_list[op_idx]
-
-            new_parent_list = (
-                divided_population[op_idx] if self.group_subs else parent_list
-            )
-
-            # Apply operator
-            new_indiv = op(indiv, new_parent_list, objfunc, self.best, self.pop_init)
-            new_indiv.genotype = objfunc.repair_solution(new_indiv.genotype)
-            new_indiv.speed = objfunc.repair_speed(new_indiv.speed)
-
-            # Collect data about each operator
-            if self.dyn_method == "fitness" or self.dyn_method == "diff":
-                self.operator_data[op_idx].append(new_indiv.fitness)
-
-            # Add to offspring list
-            offspring.append(new_indiv)
-
-        # Update best solution
-        current_best = max(offspring, key=lambda x: x.fitness)
-        if self.best.fitness < current_best.fitness:
-            self.best = current_best
-
-        return offspring
-
-    def select_individuals(self, population, offspring, **kwargs):
-        offspring_ids = [i.id for i in offspring]
-        new_population = self.selection_op(population, offspring)
-        new_ids = [i.id for i in new_population]
-
-        for idx, off_id in enumerate(offspring_ids):
-            op_idx = self.operator_idx[idx]
-
-            self.larva_count[op_idx] += 1
-
-            if off_id in new_ids:
-                self.operator_data[op_idx][0] += 1
-
-        return new_population
-
-    def update_params(self, **kwargs):
-        progress = kwargs["progress"]
-
-        self._generate_substrates(progress)
-        super().update_params(progress=progress)
-
-    def extra_step_info(self):
-        print("\n\tSubstrate probability:")
-        op_names = [i.name for i in self.operator_list]
-        weights = [round(i, 6) for i in self.operator_weight]
-        adjust = max([len(i) for i in op_names])
-        for idx, val in enumerate(op_names):
-            print(f"\t\t{val}:".ljust(adjust + 3, " ") + f"{weights[idx]}")
-
     def _operator_metric(self, data):
         result = 0
 
@@ -152,9 +84,7 @@ def _operator_probability(self, values):
         # if there are numerical error default repeat with a default value
         if (amplified_vec == 0).any() or not np.isfinite(amplified_vec).all():
             if not self.prob_amp_warned:
-                print(
-                    "Warning: the probability amplification parameter is too small, defaulting to prob_amp = 1"
-                )
+                print("Warning: the probability amplification parameter is too small, defaulting to prob_amp = 1")
                 self.prob_amp_warned = True
             prob = exp_vec / exp_vec.sum()
         else:
@@ -172,9 +102,7 @@ def _evaluate_operators(self):
 
         # take reference data for the calculation of the difference of the next evaluation
         if self.dyn_method == "diff":
-            full_data = [
-                data_point for op_data in self.operator_data for data_point in op_data
-            ]
+            full_data = [data_point for op_data in self.operator_data for data_point in op_data]
             metric = self._operator_metric(full_data)
 
         # calculate the value of each operator with the data gathered
@@ -197,9 +125,7 @@ def _evaluate_operators(self):
                 # Calculate the difference of the fitness in this generation to the previous one and
                 # store the current value for the next evaluation
                 if self.dyn_method == "diff":
-                    self.operator_metric[idx] = (
-                        self.operator_metric[idx] - self.operator_metric_prev[idx]
-                    )
+                    self.operator_metric[idx] = self.operator_metric[idx] - self.operator_metric_prev[idx]
                     self.operator_metric_prev[idx] = metric
 
                 # Reset data for next iteration
@@ -217,9 +143,71 @@ def _generate_substrates(self, progress=0):
         self.operator_w_history.append(self.operator_weight)
 
         # Choose each operator with the weights chosen
-        self.operator_idx = random.choices(
-            range(n_operators), weights=self.operator_weight, k=self.maxpopsize
-        )
+        self.operator_idx = random.choices(range(n_operators), weights=self.operator_weight, k=self.maxpopsize)
 
         # save the evaluation of each operator
         self.operator_history.append(np.array(self.operator_metric))
+
+    def perturb(self, parent_list, objfunc, **kwargs):
+        offspring = []
+
+        divided_population = [[] for i in self.operator_list]
+        for idx, indiv in enumerate(parent_list):
+            op_idx = self.operator_idx[idx]
+            divided_population[op_idx].append(indiv)
+
+        for idx, indiv in enumerate(parent_list):
+            # Select operator
+            op_idx = self.operator_idx[idx]
+
+            op = self.operator_list[op_idx]
+
+            new_parent_list = divided_population[op_idx] if self.group_subs else parent_list
+
+            # Apply operator
+            new_indiv = op(indiv, new_parent_list, objfunc, self.best, self.pop_init)
+            new_indiv.genotype = objfunc.repair_solution(new_indiv.genotype)
+            new_indiv.speed = objfunc.repair_speed(new_indiv.speed)
+
+            # Collect data about each operator
+            if self.dyn_method == "fitness" or self.dyn_method == "diff":
+                self.operator_data[op_idx].append(new_indiv.fitness)
+
+            # Add to offspring list
+            offspring.append(new_indiv)
+
+        # Update best solution
+        current_best = max(offspring, key=lambda x: x.fitness)
+        if self.best.fitness < current_best.fitness:
+            self.best = current_best
+
+        return offspring
+
+    def select_individuals(self, population, offspring, **kwargs):
+        offspring_ids = [i.id for i in offspring]
+        new_population = self.selection_op(population, offspring)
+        new_ids = [i.id for i in new_population]
+
+        for idx, off_id in enumerate(offspring_ids):
+            op_idx = self.operator_idx[idx]
+
+            self.larva_count[op_idx] += 1
+
+            if off_id in new_ids:
+                self.operator_data[op_idx][0] += 1
+
+        return new_population
+
+    def update_params(self, **kwargs):
+        progress = kwargs["progress"]
+
+        self._generate_substrates(progress)
+        super().update_params(progress=progress)
+
+    def extra_step_info(self):
+        print("\n\tSubstrate probability:")
+        op_names = [i.name for i in self.operator_list]
+        weights = [round(i, 6) for i in self.operator_weight]
+        adjust = max([len(i) for i in op_names])
+        for idx, val in enumerate(op_names):
+            print(f"\t\t{val}:".ljust(adjust + 3, " ") + f"{weights[idx]}")
diff --git a/src/metaheuristic_designer/strategies/CRO/PCRO_SL.py b/src/metaheuristic_designer/strategies/CRO/PCRO_SL.py
index c7ffdab..8a4e161 100644
--- a/src/metaheuristic_designer/strategies/CRO/PCRO_SL.py
+++ b/src/metaheuristic_designer/strategies/CRO/PCRO_SL.py
@@ -23,15 +23,11 @@ def __init__(
         name: str = "PCRO-SL",
     ):
         super().__init__(pop_init, operator_list, params=params, name=name)
-        self.operator_idx = random.choices(
-            range(len(self.operator_list)), k=self.maxpopsize
-        )
+        self.operator_idx = random.choices(range(len(self.operator_list)), k=self.maxpopsize)
 
     def update_params(self, **kwargs):
         progress = kwargs["progress"]
 
-        self.operator_idx = random.choices(
-            range(len(self.operator_list)), k=self.maxpopsize
-        )
+        self.operator_idx = random.choices(range(len(self.operator_list)), k=self.maxpopsize)
 
         super().update_params(progress=progress)
diff --git a/src/metaheuristic_designer/strategies/Classic/CMA_ES.py b/src/metaheuristic_designer/strategies/Classic/CMA_ES.py
index be21196..848e95d 100644
--- a/src/metaheuristic_designer/strategies/Classic/CMA_ES.py
+++ b/src/metaheuristic_designer/strategies/Classic/CMA_ES.py
@@ -32,12 +32,8 @@ def __init__(
 
         mutate = OperatorMeta("Sequence", [mutate1, rand1])
 
-        super().__init__(
-            pop_init, mutate, cross, parent_select, selection, params, name
-        )
+        super().__init__(pop_init, mutate, cross, parent_select, selection, params, name)
 
     def initialize(self, objfunc):
-        objfunc.encoding = CMAEncoding(
-            self.params["nparams"], pre_encoding=objfunc.encoding
-        )
+        objfunc.encoding = CMAEncoding(self.params["nparams"], pre_encoding=objfunc.encoding)
         super().initialize(objfunc)
diff --git a/src/metaheuristic_designer/strategies/Classic/DE.py b/src/metaheuristic_designer/strategies/Classic/DE.py
index 09cac9e..3ec7eb7 100644
--- a/src/metaheuristic_designer/strategies/Classic/DE.py
+++ b/src/metaheuristic_designer/strategies/Classic/DE.py
@@ -20,6 +20,4 @@ def __init__(
         if selection_op is None:
             selection_op = SurvivorSelection("One-to-one")
 
-        super().__init__(
-            pop_init, de_op, selection_op=selection_op, params=params, name=name
-        )
+        super().__init__(pop_init, de_op, selection_op=selection_op, params=params, name=name)
diff --git a/src/metaheuristic_designer/strategies/Classic/ES.py b/src/metaheuristic_designer/strategies/Classic/ES.py
index 944aaf1..4bf2b79 100644
--- a/src/metaheuristic_designer/strategies/Classic/ES.py
+++ b/src/metaheuristic_designer/strategies/Classic/ES.py
@@ -1,6 +1,7 @@
 from __future__ import annotations
 from typing import Union
-from ...operators import OperatorMeta
+from ...operators import OperatorMeta, OperatorReal
+from ...selectionMethods import ParentSelection
 from ..VariablePopulation import VariablePopulation
 
 
@@ -13,14 +14,18 @@ def __init__(
         self,
         pop_init: Initializer,
         mutation_op: Operator,
-        cross_op: Operator,
-        parent_sel_op: ParentSelection,
-        selection_op: SurvivorSelection,
+        cross_op: Operator = None,
+        parent_sel_op: ParentSelection = None,
+        selection_op: SurvivorSelection = None,
         params: Union[ParamScheduler, dict] = {},
         name: str = "ES",
     ):
+        if cross_op is None:
+            evolve_op = mutation_op
+        else:
+            evolve_op = OperatorMeta("Sequence", [mutation_op, cross_op])
+
         offspring_size = params.get("offspringSize", pop_init.pop_size)
-        evolve_op = OperatorMeta("Sequence", [mutation_op, cross_op])
 
         super().__init__(
             pop_init,
diff --git a/src/metaheuristic_designer/strategies/Classic/GA.py b/src/metaheuristic_designer/strategies/Classic/GA.py
index 49dac4f..6e95b44 100644
--- a/src/metaheuristic_designer/strategies/Classic/GA.py
+++ b/src/metaheuristic_designer/strategies/Classic/GA.py
@@ -1,7 +1,7 @@
 from __future__ import annotations
 from typing import Union
 import numpy as np
-from ...operators import OperatorReal, OperatorMeta
+from ...operators import OperatorReal, OperatorMeta, OperatorNull
 from ..VariablePopulation import VariablePopulation
 
 
@@ -23,14 +23,10 @@ def __init__(
         self.pmut = params.get("pmut", 0.1)
         self.pcross = params.get("pcross", 0.9)
 
-        null_operator = OperatorReal("Nothing")
+        null_operator = OperatorNull()
 
-        prob_mut_op = OperatorMeta(
-            "Branch", [mutation_op, null_operator], {"p": self.pmut}
-        )
-        prob_cross_op = OperatorMeta(
-            "Branch", [cross_op, null_operator], {"p": self.pcross}
-        )
+        prob_mut_op = OperatorMeta("Branch", [mutation_op, null_operator], {"p": self.pmut})
+        prob_cross_op = OperatorMeta("Branch", [cross_op, null_operator], {"p": self.pcross})
 
         evolve_op = OperatorMeta("Sequence", [prob_mut_op, prob_cross_op])
 
diff --git a/src/metaheuristic_designer/strategies/Classic/PSO.py b/src/metaheuristic_designer/strategies/Classic/PSO.py
index a47674f..2a24ac4 100644
--- a/src/metaheuristic_designer/strategies/Classic/PSO.py
+++ b/src/metaheuristic_designer/strategies/Classic/PSO.py
@@ -34,9 +34,7 @@ def __init__(
 
         selection_op = SurvivorSelection("Generational")
 
-        super().__init__(
-            pop_init, pso_op, selection_op=selection_op, params=params, name=name
-        )
+        super().__init__(pop_init, pso_op, selection_op=selection_op, params=params, name=name)
 
     def extra_step_info(self):
         popul_matrix = np.array(list(map(lambda x: x.genotype, self.population)))
diff --git a/src/metaheuristic_designer/strategies/Classic/SA.py b/src/metaheuristic_designer/strategies/Classic/SA.py
index 73befb9..8917aa7 100644
--- a/src/metaheuristic_designer/strategies/Classic/SA.py
+++ b/src/metaheuristic_designer/strategies/Classic/SA.py
@@ -4,10 +4,12 @@
 import numpy as np
 from ...ParamScheduler import ParamScheduler
 from ...SearchStrategy import SearchStrategy
+from ...selectionMethods import SurvivorSelection
 from ...Operator import Operator
+from ..HillClimb import HillClimb
 
 
-class SA(SearchStrategy):
+class SA(HillClimb):
     """
     Simulated annealing
     """
@@ -20,31 +22,14 @@ def __init__(
         name: str = "SA",
     ):
         self.iter = params.get("iter", 100)
+        self.iter_count = 0
         self.temp_init = params.get("temp_init", 100)
         self.temp = self.temp_init
         self.alpha = params.get("alpha", 0.99)
 
-        self.perturb_op = perturb_op
+        selection_op = SurvivorSelection("ProbHillClimb", {"p": np.exp(-1 / self.temp_init)})
 
-        super().__init__(pop_init, params=params, name=name)
-
-    def perturb(self, indiv_list, objfunc, **kwargs):
-        indiv = indiv_list[0]
-        for _ in range(self.iter):
-            new_indiv = self.perturb_op(
-                indiv, indiv_list, objfunc, self.best, self.pop_init
-            )
-            new_indiv.genotype = objfunc.repair_solution(new_indiv.genotype)
-
-            # Accept the new solution even if it is worse with a probability
-            p = np.exp(-1 / self.temp)
-            if new_indiv.fitness > indiv.fitness or random.random() < p:
-                indiv = new_indiv
-
-            if new_indiv.fitness > self.best.fitness:
-                self.best = new_indiv
-
-        return [indiv]
+        super().__init__(pop_init, perturb_op, selection_op, params=params, name=name)
 
     def update_params(self, **kwargs):
         progress = kwargs["progress"]
@@ -58,9 +43,14 @@ def update_params(self, **kwargs):
             self.iter = round(self.params["iter"])
             self.alpha = self.params["alpha"]
 
-        self.temp = self.temp * self.alpha
+        self.iter_count += 1
+        if self.iter_count > self.iter:
+            self.temp = self.temp * self.alpha
+            self.selection_op.params["p"] = np.exp(-1 / self.temp)
+            self.iter_count = 0
 
     def extra_step_info(self):
         print()
+        print(f"\tTemp iters: {self.iter_count}/{self.iter}")
         print(f"\tTemperature: {float(self.temp):0.3}")
         print(f"\tAccept prob: {np.exp(-1 / self.temp):0.3}")
diff --git a/src/metaheuristic_designer/strategies/HillClimb.py b/src/metaheuristic_designer/strategies/HillClimb.py
index 908f07e..4e51254 100644
--- a/src/metaheuristic_designer/strategies/HillClimb.py
+++ b/src/metaheuristic_designer/strategies/HillClimb.py
@@ -4,7 +4,7 @@
 from ..ParamScheduler import ParamScheduler
 from ..SearchStrategy import SearchStrategy
 from ..Operator import Operator
-from ..operators import OperatorReal
+from ..operators import OperatorNull
 from ..selectionMethods import SurvivorSelection
 
 
@@ -24,36 +24,29 @@ def __init__(
         self.iterations = params.get("iters", 1)
 
         if perturb_op is None:
-            perturb_op = OperatorReal("Nothing")
+            perturb_op = OperatorNull()
         self.perturb_op = perturb_op
 
         if selection_op is None:
-            selection_op = SurvivorSelection("One-to-One")
+            selection_op = SurvivorSelection("HillClimb")
         self.selection_op = selection_op
 
         super().__init__(pop_init, params=params, name=name)
 
     def perturb(self, indiv_list, objfunc, **kwargs):
-        next_indiv_list = copy(indiv_list)
-        for _ in range(self.iterations):
-            offspring = []
-            for indiv in next_indiv_list:
-                # Perturb individual
-                new_indiv = self.perturb_op(
-                    indiv, next_indiv_list, objfunc, self.best, self.pop_init
-                )
-                new_indiv.genotype = objfunc.repair_solution(new_indiv.genotype)
-
-                offspring.append(new_indiv)
-
-            # Keep best individual regardless of selection method
-            current_best = max(offspring, key=lambda x: x.fitness)
-            if self.best.fitness <= current_best.fitness:
-                self.best = current_best
-
-            next_indiv_list = self.selection_op(next_indiv_list, offspring)
-
-        return next_indiv_list
+        offspring = []
+        for indiv in indiv_list:
+            # Perturb individual
+            new_indiv = self.perturb_op(indiv, indiv_list, objfunc, self.best, self.pop_init)
+            new_indiv.genotype = objfunc.repair_solution(new_indiv.genotype)
+            new_indiv.speed = objfunc.repair_speed(new_indiv.speed)
+
+            offspring.append(new_indiv)
+
+        return offspring
+
+    def select_individuals(self, population, offspring, **kwargs):
+        return self.selection_op(population, offspring)
 
     def update_params(self, **kwargs):
         progress = kwargs["progress"]
diff --git a/src/metaheuristic_designer/strategies/LocalSearch.py b/src/metaheuristic_designer/strategies/LocalSearch.py
index 4f0dd3a..a62e29f 100644
--- a/src/metaheuristic_designer/strategies/LocalSearch.py
+++ b/src/metaheuristic_designer/strategies/LocalSearch.py
@@ -4,7 +4,7 @@
 from ..ParamScheduler import ParamScheduler
 from ..SearchStrategy import SearchStrategy
 from ..Operator import Operator
-from ..operators import OperatorReal
+from ..operators import OperatorNull
 from ..selectionMethods import SurvivorSelection
 
 
@@ -24,36 +24,30 @@ def __init__(
         self.iterations = params.get("iters", 100)
 
         if perturb_op is None:
-            perturb_op = OperatorReal("Nothing")
+            perturb_op = OperatorNull()
         self.perturb_op = perturb_op
 
         if selection_op is None:
-            selection_op = SurvivorSelection("One-to-One")
+            selection_op = SurvivorSelection("KeepBest", {"amount": 1})
         self.selection_op = selection_op
 
         super().__init__(pop_init, params=params, name=name)
 
     def perturb(self, indiv_list, objfunc, **kwargs):
-        next_indiv_list = copy(indiv_list)
-        for _ in range(self.iterations):
-            offspring = []
-            for indiv in indiv_list:
-                # Perturb individual
-                new_indiv = self.perturb_op(
-                    indiv, indiv_list, objfunc, self.best, self.pop_init
-                )
-                new_indiv.genotype = objfunc.repair_solution(new_indiv.genotype)
-
-                offspring.append(new_indiv)
-
-            # Keep best individual regardless of selection method
-            current_best = max(offspring, key=lambda x: x.fitness)
-            if self.best.fitness < current_best.fitness:
-                self.best = current_best
-
-            next_indiv_list = self.selection_op(next_indiv_list, offspring)
-
-        return next_indiv_list
+        offspring = []
+        indiv = indiv_list[0]
+        for i in range(self.iterations):
+            # Perturb individual
+            new_indiv = self.perturb_op(indiv, indiv_list, objfunc, self.best, self.pop_init)
+            new_indiv.genotype = objfunc.repair_solution(new_indiv.genotype)
+            new_indiv.speed = objfunc.repair_speed(new_indiv.speed)
+
+            offspring.append(new_indiv)
+
+        return offspring
+
+    def select_individuals(self, population, offspring, **kwargs):
+        return self.selection_op(population, offspring)
 
     def update_params(self, **kwargs):
         progress = kwargs["progress"]
diff --git a/src/metaheuristic_designer/strategies/StaticPopulation.py b/src/metaheuristic_designer/strategies/StaticPopulation.py
index 632b0b7..416e65a 100644
--- a/src/metaheuristic_designer/strategies/StaticPopulation.py
+++ b/src/metaheuristic_designer/strategies/StaticPopulation.py
@@ -41,20 +41,13 @@ def perturb(self, parent_list, objfunc, **kwargs):
         offspring = []
         for indiv in parent_list:
             # Apply operator
-            new_indiv = self.operator(
-                indiv, parent_list, objfunc, self.best, self.pop_init
-            )
+            new_indiv = self.operator(indiv, parent_list, objfunc, self.best, self.pop_init)
             new_indiv.genotype = objfunc.repair_solution(new_indiv.genotype)
             new_indiv.speed = objfunc.repair_speed(new_indiv.speed)
 
             # Add to offspring list
             offspring.append(new_indiv)
 
-        # Update best solution
-        current_best = max(offspring, key=lambda x: x.fitness)
-        if self.best.fitness < current_best.fitness:
-            self.best = current_best
-
         return offspring
 
     def select_individuals(self, population, offspring, **kwargs):
diff --git a/src/metaheuristic_designer/strategies/VNS/RVNS.py b/src/metaheuristic_designer/strategies/VNS/RVNS.py
index 2472f61..7a629cf 100644
--- a/src/metaheuristic_designer/strategies/VNS/RVNS.py
+++ b/src/metaheuristic_designer/strategies/VNS/RVNS.py
@@ -46,18 +46,11 @@ def perturb(self, indiv_list, objfunc, **kwargs):
         offspring = []
         for indiv in indiv_list:
             # Perturb individual
-            new_indiv = self.perturb_op(
-                indiv, indiv_list, objfunc, self.best, self.pop_init
-            )
+            new_indiv = self.perturb_op(indiv, indiv_list, objfunc, self.best, self.pop_init)
             new_indiv.genotype = objfunc.repair_solution(new_indiv.genotype)
 
             offspring.append(new_indiv)
 
-        # Keep best individual regardless of selection method
-        current_best = max(offspring, key=lambda x: x.fitness)
-        if self.best.fitness < current_best.fitness:
-            self.best = current_best
-
         return offspring
 
     def select_individuals(self, population, offspring, **kwargs):
@@ -76,15 +69,10 @@ def update_params(self, **kwargs):
         if isinstance(self.perturb_op, Operator):
             self.perturb_op.step(progress)
 
-        if (
-            self.perturb_op.chosen_idx >= len(self.op_list)
-            or self.perturb_op.chosen_idx < 0
-        ):
+        if self.perturb_op.chosen_idx >= len(self.op_list) or self.perturb_op.chosen_idx < 0:
             self.perturb_op.chosen_idx = 0
 
     def extra_step_info(self):
         idx = self.perturb_op.chosen_idx
 
-        print(
-            f"\tCurrent Operator: {idx}/{len(self.op_list)}, {self.op_list[idx].name}"
-        )
+        print(f"\tCurrent Operator: {idx}/{len(self.op_list)}, {self.op_list[idx].name}")
diff --git a/src/metaheuristic_designer/strategies/VNS/VND.py b/src/metaheuristic_designer/strategies/VNS/VND.py
index 90884e9..a217e10 100644
--- a/src/metaheuristic_designer/strategies/VNS/VND.py
+++ b/src/metaheuristic_designer/strategies/VNS/VND.py
@@ -31,11 +31,7 @@ def __init__(
         self.op_list = op_list
         self.perturb_op = OperatorMeta("Pick", op_list, {"init_idx": 0})
 
-        self.nchange = (
-            NeighborhoodChange.from_str(params["nchange"])
-            if "nchange" in params
-            else NeighborhoodChange.SEQ
-        )
+        self.nchange = NeighborhoodChange.from_str(params["nchange"]) if "nchange" in params else NeighborhoodChange.SEQ
 
         self.current_op = 0
 
@@ -60,9 +56,7 @@ def perturb(self, indiv_list, objfunc, **kwargs):
             offspring = []
             for indiv in indiv_list:
                 # Perturb individual
-                new_indiv = self.perturb_op(
-                    indiv, indiv_list, objfunc, self.best, self.pop_init
-                )
+                new_indiv = self.perturb_op(indiv, indiv_list, objfunc, self.best, self.pop_init)
                 new_indiv.genotype = objfunc.repair_solution(new_indiv.genotype)
 
                 offspring.append(new_indiv)
@@ -79,9 +73,7 @@ def perturb(self, indiv_list, objfunc, **kwargs):
     def select_individuals(self, population, offspring, **kwargs):
         new_population = self.selection_op(population, offspring)
 
-        self.perturb_op.chosen_idx = next_neighborhood(
-            offspring[0], population[0], self.perturb_op.chosen_idx, self.nchange
-        )
+        self.perturb_op.chosen_idx = next_neighborhood(offspring[0], population[0], self.perturb_op.chosen_idx, self.nchange)
 
         return new_population
 
@@ -91,16 +83,11 @@ def update_params(self, **kwargs):
         if isinstance(self.perturb_op, Operator):
             self.perturb_op.step(progress)
 
-        if (
-            self.perturb_op.chosen_idx >= len(self.op_list)
-            or self.perturb_op.chosen_idx < 0
-        ):
+        if self.perturb_op.chosen_idx >= len(self.op_list) or self.perturb_op.chosen_idx < 0:
             self.perturb_op.chosen_idx = 0
 
     def extra_step_info(self):
         idx = self.perturb_op.chosen_idx
 
-        print(
-            f"\tCurrent Operator: {idx}/{len(self.op_list)}, {self.op_list[idx].name}"
-        )
+        print(f"\tCurrent Operator: {idx}/{len(self.op_list)}, {self.op_list[idx].name}")
         # time.sleep(0.25)
diff --git a/src/metaheuristic_designer/strategies/VNS/VNS.py b/src/metaheuristic_designer/strategies/VNS/VNS.py
index 33833c7..9d563f0 100644
--- a/src/metaheuristic_designer/strategies/VNS/VNS.py
+++ b/src/metaheuristic_designer/strategies/VNS/VNS.py
@@ -28,11 +28,7 @@ def __init__(
         self.op_list = op_list
         self.perturb_op = OperatorMeta("Pick", op_list, {"init_idx": 0})
 
-        self.nchange = (
-            NeighborhoodChange.from_str(params["nchange"])
-            if "nchange" in params
-            else NeighborhoodChange.SEQ
-        )
+        self.nchange = NeighborhoodChange.from_str(params["nchange"]) if "nchange" in params else NeighborhoodChange.SEQ
 
         self.local_search = local_search
 
@@ -57,9 +53,7 @@ def perturb(self, indiv_list, objfunc, **kwargs):
         offspring = []
         for indiv in indiv_list:
             # Perturb individual
-            new_indiv = self.perturb_op(
-                indiv, indiv_list, objfunc, self.best, self.pop_init
-            )
+            new_indiv = self.perturb_op(indiv, indiv_list, objfunc, self.best, self.pop_init)
             new_indiv.genotype = objfunc.repair_solution(new_indiv.genotype)
 
             # Local search
@@ -70,9 +64,7 @@ def perturb(self, indiv_list, objfunc, **kwargs):
 
                 offspring = self.local_search.perturb(parents, objfunc, kwargs)
 
-                population = self.local_search.select_individuals(
-                    population, offspring, kwargs
-                )
+                population = self.local_search.select_individuals(population, offspring, kwargs)
 
                 self.local_search.update_params()
 
@@ -90,9 +82,7 @@ def perturb(self, indiv_list, objfunc, **kwargs):
     def select_individuals(self, population, offspring, **kwargs):
         new_population = self.selection_op(population, offspring)
 
-        self.perturb_op.chosen_idx = next_neighborhood(
-            offspring[0], population[0], self.perturb_op.chosen_idx, self.nchange
-        )
+        self.perturb_op.chosen_idx = next_neighborhood(offspring[0], population[0], self.perturb_op.chosen_idx, self.nchange)
 
         return new_population
 
@@ -102,15 +92,10 @@ def update_params(self, **kwargs):
         if isinstance(self.perturb_op, Operator):
             self.perturb_op.step(progress)
 
-        if (
-            self.perturb_op.chosen_idx >= len(self.op_list)
-            or self.perturb_op.chosen_idx < 0
-        ):
+        if self.perturb_op.chosen_idx >= len(self.op_list) or self.perturb_op.chosen_idx < 0:
             self.perturb_op.chosen_idx = 0
 
     def extra_step_info(self):
         idx = self.perturb_op.chosen_idx
 
-        print(
-            f"\tCurrent Operator: {idx}/{len(self.op_list)}, {self.op_list[idx].name}"
-        )
+        print(f"\tCurrent Operator: {idx}/{len(self.op_list)}, {self.op_list[idx].name}")
diff --git a/src/metaheuristic_designer/strategies/VariablePopulation.py b/src/metaheuristic_designer/strategies/VariablePopulation.py
index 8da8d23..b92997a 100644
--- a/src/metaheuristic_designer/strategies/VariablePopulation.py
+++ b/src/metaheuristic_designer/strategies/VariablePopulation.py
@@ -49,20 +49,13 @@ def perturb(self, parent_list, objfunc, **kwargs):
         while len(offspring) < self.n_offspring:
             # Apply operator
             indiv = random.choice(parent_list)
-            new_indiv = self.operator(
-                indiv, parent_list, objfunc, self.best, self.pop_init
-            )
+            new_indiv = self.operator(indiv, parent_list, objfunc, self.best, self.pop_init)
             new_indiv.genotype = objfunc.repair_solution(new_indiv.genotype)
             new_indiv.speed = objfunc.repair_speed(new_indiv.speed)
 
             # Add to offspring list
             offspring.append(new_indiv)
 
-        # Update best solution
-        current_best = max(offspring, key=lambda x: x.fitness)
-        if self.best.fitness < current_best.fitness:
-            self.best = current_best
-
         return offspring
 
     def select_individuals(self, population, offspring, **kwargs):
diff --git a/tutorials/algorithm_selection.ipynb b/tutorials/algorithm_selection.ipynb
new file mode 100644
index 0000000..1b385b3
--- /dev/null
+++ b/tutorials/algorithm_selection.ipynb
@@ -0,0 +1,828 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Automatic algorithm/strategy selection\n",
+    "Here we will show you how you can perform an algorithm selection, helping you choose the best optimization algorithm for your specific problem."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 1,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from metaheuristic_designer import simple\n",
+    "from metaheuristic_designer.operators import OperatorReal\n",
+    "from metaheuristic_designer.selectionMethods import ParentSelection, SurvivorSelection\n",
+    "from metaheuristic_designer.algorithms import AlgorithmSelection, StrategySelection\n",
+    "from metaheuristic_designer.initializers import UniformVectorInitializer\n",
+    "from metaheuristic_designer.strategies import HillClimb, SA, ES, GA, DE, PSO, RandomSearch\n",
+    "from metaheuristic_designer.benchmarks import *\n",
+    "import numpy as np"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For the following examples, we will use the Rastrigin function for which we need to find the minimum value. This is a very hard optimization problem and will allow us to compare the performance of different algorithms."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 2,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "rastrigin_func = Rastrigin(size=2)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Algorithm selection\n",
+    "We will first explore the AlgorithmSelection class, which launches each of the specified algorithms multiple times and gets some statistics about the evaluation.\n",
+    "\n",
+    "We begin by defining the algorithm we are going to use. The ```simple``` package gives us premade implementations of well known optimization algorithms with a single function call, although customization is limited."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 3,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Running 7 algorithms 10 times each.\n"
+     ]
+    },
+    {
+     "data": {
+      "text/html": [
+       "<style>\n",
+       ".enlighten-fg-red {\n",
+       "  color: #cd0000;\n",
+       "}\n",
+       "</style>\n",
+       "<div class=\"enlighten\">\n",
+       "  <div class=\"enlighten-bar\">\n",
+       "    <pre>Launching algorithms   0%|<span class=\"enlighten-fg-red\">                                                   </span>| 0/7 [00:00&lt;&#63;, 0.00/s]</pre>\n",
+       "  </div>\n",
+       "  <div class=\"enlighten-bar\">\n",
+       "    <pre></pre>\n",
+       "  </div>\n",
+       "</div>\n"
+      ],
+      "text/plain": [
+       "<IPython.core.display.HTML object>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "# The parameters common to all algorithms\n",
+    "params = {\n",
+    "    # Stopping condition\n",
+    "    \"stop_cond\": \"neval\",\n",
+    "    \"neval\": 5e4,\n",
+    "    # Encoding\n",
+    "    \"encoding\": \"real\",\n",
+    "    # Disable verbose, the output gets very cluttered otherwise\n",
+    "    \"verbose\": False,\n",
+    "    # Number of times to launch each algorithm\n",
+    "    \"repetitions\": 10,\n",
+    "}\n",
+    "\n",
+    "# Define the algorithms we are going to use\n",
+    "algorithms = [\n",
+    "    simple.random_search(rastrigin_func, params),\n",
+    "    simple.hill_climb(rastrigin_func, params),\n",
+    "    simple.simulated_annealing(rastrigin_func, params),\n",
+    "    simple.evolution_strategy(rastrigin_func, params),\n",
+    "    simple.genetic_algorithm(rastrigin_func, params),\n",
+    "    simple.differential_evolution(rastrigin_func, params),\n",
+    "    simple.particle_swarm(rastrigin_func, params),\n",
+    "]\n",
+    "\n",
+    "# Instanciate the AlgorithmSelection class\n",
+    "algorithm_search = AlgorithmSelection(algorithms)\n",
+    "\n",
+    "# Launch all the algorithms and get the best overall solution\n",
+    "solution, best_fitness, report = algorithm_search.optimize()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The optimize function should have returned the best overall solution and a report with statistics about the execution"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 4,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "best solution: [-2.62463187e-09 -2.31890834e-09]\n",
+      "fitness: 0.0\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(f\"best solution: {solution}\")\n",
+    "print(f\"fitness: {best_fitness}\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 5,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>name</th>\n",
+       "      <th>realtime_min</th>\n",
+       "      <th>realtime_avg</th>\n",
+       "      <th>realtime_max</th>\n",
+       "      <th>realtime_std</th>\n",
+       "      <th>cputime_min</th>\n",
+       "      <th>cputime_avg</th>\n",
+       "      <th>cputime_max</th>\n",
+       "      <th>cputime_std</th>\n",
+       "      <th>fitness_min</th>\n",
+       "      <th>fitness_avg</th>\n",
+       "      <th>fitness_max</th>\n",
+       "      <th>fitness_std</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>5</th>\n",
+       "      <td>DE</td>\n",
+       "      <td>1.520876</td>\n",
+       "      <td>1.538080</td>\n",
+       "      <td>1.558256</td>\n",
+       "      <td>0.010276</td>\n",
+       "      <td>1.520177</td>\n",
+       "      <td>1.541681</td>\n",
+       "      <td>1.556947</td>\n",
+       "      <td>0.010985</td>\n",
+       "      <td>0.000000</td>\n",
+       "      <td>0.000000</td>\n",
+       "      <td>0.000000</td>\n",
+       "      <td>0.000000</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>6</th>\n",
+       "      <td>PSO</td>\n",
+       "      <td>1.113223</td>\n",
+       "      <td>1.127316</td>\n",
+       "      <td>1.141176</td>\n",
+       "      <td>0.008725</td>\n",
+       "      <td>1.112528</td>\n",
+       "      <td>1.128177</td>\n",
+       "      <td>1.145475</td>\n",
+       "      <td>0.009978</td>\n",
+       "      <td>0.000000</td>\n",
+       "      <td>0.000000</td>\n",
+       "      <td>0.000000</td>\n",
+       "      <td>0.000000</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>0</th>\n",
+       "      <td>RandomSearch</td>\n",
+       "      <td>1.251807</td>\n",
+       "      <td>1.279428</td>\n",
+       "      <td>1.313645</td>\n",
+       "      <td>0.020328</td>\n",
+       "      <td>1.252516</td>\n",
+       "      <td>1.280534</td>\n",
+       "      <td>1.315798</td>\n",
+       "      <td>0.020718</td>\n",
+       "      <td>0.004720</td>\n",
+       "      <td>0.031401</td>\n",
+       "      <td>0.119305</td>\n",
+       "      <td>0.046518</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>4</th>\n",
+       "      <td>GA</td>\n",
+       "      <td>1.118457</td>\n",
+       "      <td>1.142087</td>\n",
+       "      <td>1.209193</td>\n",
+       "      <td>0.028553</td>\n",
+       "      <td>1.117847</td>\n",
+       "      <td>1.141340</td>\n",
+       "      <td>1.208274</td>\n",
+       "      <td>0.028467</td>\n",
+       "      <td>0.800926</td>\n",
+       "      <td>1.389587</td>\n",
+       "      <td>2.562032</td>\n",
+       "      <td>0.497320</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>3</th>\n",
+       "      <td>ES</td>\n",
+       "      <td>1.199183</td>\n",
+       "      <td>1.214232</td>\n",
+       "      <td>1.239551</td>\n",
+       "      <td>0.011732</td>\n",
+       "      <td>1.198628</td>\n",
+       "      <td>1.213559</td>\n",
+       "      <td>1.238718</td>\n",
+       "      <td>0.011679</td>\n",
+       "      <td>1.412214</td>\n",
+       "      <td>2.467525</td>\n",
+       "      <td>5.448058</td>\n",
+       "      <td>1.712328</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1</th>\n",
+       "      <td>HillClimb</td>\n",
+       "      <td>1.127098</td>\n",
+       "      <td>1.147712</td>\n",
+       "      <td>1.191707</td>\n",
+       "      <td>0.020343</td>\n",
+       "      <td>1.127515</td>\n",
+       "      <td>1.148799</td>\n",
+       "      <td>1.193783</td>\n",
+       "      <td>0.020693</td>\n",
+       "      <td>1.989918</td>\n",
+       "      <td>11.466948</td>\n",
+       "      <td>33.794655</td>\n",
+       "      <td>11.957213</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>2</th>\n",
+       "      <td>SA</td>\n",
+       "      <td>1.149471</td>\n",
+       "      <td>1.170844</td>\n",
+       "      <td>1.203308</td>\n",
+       "      <td>0.017003</td>\n",
+       "      <td>1.149264</td>\n",
+       "      <td>1.171789</td>\n",
+       "      <td>1.202562</td>\n",
+       "      <td>0.016068</td>\n",
+       "      <td>12.178331</td>\n",
+       "      <td>19.627940</td>\n",
+       "      <td>54.316949</td>\n",
+       "      <td>12.808269</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "           name  realtime_min  realtime_avg  realtime_max  realtime_std  \\\n",
+       "5            DE      1.520876      1.538080      1.558256      0.010276   \n",
+       "6           PSO      1.113223      1.127316      1.141176      0.008725   \n",
+       "0  RandomSearch      1.251807      1.279428      1.313645      0.020328   \n",
+       "4            GA      1.118457      1.142087      1.209193      0.028553   \n",
+       "3            ES      1.199183      1.214232      1.239551      0.011732   \n",
+       "1     HillClimb      1.127098      1.147712      1.191707      0.020343   \n",
+       "2            SA      1.149471      1.170844      1.203308      0.017003   \n",
+       "\n",
+       "   cputime_min  cputime_avg  cputime_max  cputime_std  fitness_min  \\\n",
+       "5     1.520177     1.541681     1.556947     0.010985     0.000000   \n",
+       "6     1.112528     1.128177     1.145475     0.009978     0.000000   \n",
+       "0     1.252516     1.280534     1.315798     0.020718     0.004720   \n",
+       "4     1.117847     1.141340     1.208274     0.028467     0.800926   \n",
+       "3     1.198628     1.213559     1.238718     0.011679     1.412214   \n",
+       "1     1.127515     1.148799     1.193783     0.020693     1.989918   \n",
+       "2     1.149264     1.171789     1.202562     0.016068    12.178331   \n",
+       "\n",
+       "   fitness_avg  fitness_max  fitness_std  \n",
+       "5     0.000000     0.000000     0.000000  \n",
+       "6     0.000000     0.000000     0.000000  \n",
+       "0     0.031401     0.119305     0.046518  \n",
+       "4     1.389587     2.562032     0.497320  \n",
+       "3     2.467525     5.448058     1.712328  \n",
+       "1    11.466948    33.794655    11.957213  \n",
+       "2    19.627940    54.316949    12.808269  "
+      ]
+     },
+     "execution_count": 5,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "report.sort_values(\"fitness_avg\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "With this, we can clearly see that the best algorithm for this task seems to be the Differential Evolution algorithm, and surprisingly, RandomSearch performs better than HillClimb, SA, ES and even GA. This could be caused by some misconfiguration of parameters, but since we just used the default implementation, we have no way of changing this."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Strategy Selection\n",
+    "We could have specified the algorithms directly constructing them piece by piece, but unless you need to try some specific implementation of an algorithm, it is much simpler to search for the best Search Strategy. This is where the StrategySearch class is used. It is very similar to the AlgorithmSelection class, but using search strategies.\n",
+    "\n",
+    "This also has the advantage that you can choose a custom name for each strategy."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "Running 11 algorithms 10 times each.\n"
+     ]
+    },
+    {
+     "data": {
+      "text/html": [
+       "<style>\n",
+       ".enlighten-fg-red {\n",
+       "  color: #cd0000;\n",
+       "}\n",
+       ".enlighten-fg-green {\n",
+       "  color: #00cd00;\n",
+       "}\n",
+       "</style>\n",
+       "<div class=\"enlighten\">\n",
+       "  <div class=\"enlighten-bar\">\n",
+       "    <pre>Launching algorithms 100%|<span class=\"enlighten-fg-red\">█████████████████████████████████████████████</span>| 11/11 [02:17&lt;00:00, 0.08/s]</pre>\n",
+       "  </div>\n",
+       "  <div class=\"enlighten-bar\">\n",
+       "    <pre>Evaluating RandomSearch 100%|<span class=\"enlighten-fg-green\">██████████████████████████████████████████</span>| 10/10 [00:13&lt;00:00, 0.77/s]</pre>\n",
+       "  </div>\n",
+       "  <div class=\"enlighten-bar\">\n",
+       "    <pre>Evaluating HillClimb 100%|<span class=\"enlighten-fg-green\">█████████████████████████████████████████████</span>| 10/10 [00:11&lt;00:00, 0.87/s]</pre>\n",
+       "  </div>\n",
+       "  <div class=\"enlighten-bar\">\n",
+       "    <pre>Evaluating SA 100%|<span class=\"enlighten-fg-green\">████████████████████████████████████████████████████</span>| 10/10 [00:12&lt;00:00, 0.86/s]</pre>\n",
+       "  </div>\n",
+       "  <div class=\"enlighten-bar\">\n",
+       "    <pre>Evaluating ParallelSA 100%|<span class=\"enlighten-fg-green\">████████████████████████████████████████████</span>| 10/10 [00:09&lt;00:00, 1.09/s]</pre>\n",
+       "  </div>\n",
+       "  <div class=\"enlighten-bar\">\n",
+       "    <pre>Evaluating ES-(μ+λ) 100%|<span class=\"enlighten-fg-green\">██████████████████████████████████████████████</span>| 10/10 [00:09&lt;00:00, 1.07/s]</pre>\n",
+       "  </div>\n",
+       "  <div class=\"enlighten-bar\">\n",
+       "    <pre>Evaluating ES-(μ,λ) 100%|<span class=\"enlighten-fg-green\">██████████████████████████████████████████████</span>| 10/10 [00:09&lt;00:00, 1.06/s]</pre>\n",
+       "  </div>\n",
+       "  <div class=\"enlighten-bar\">\n",
+       "    <pre>Evaluating GA 100%|<span class=\"enlighten-fg-green\">████████████████████████████████████████████████████</span>| 10/10 [00:18&lt;00:00, 0.55/s]</pre>\n",
+       "  </div>\n",
+       "  <div class=\"enlighten-bar\">\n",
+       "    <pre>Evaluating DE/best/1 100%|<span class=\"enlighten-fg-green\">█████████████████████████████████████████████</span>| 10/10 [00:16&lt;00:00, 0.64/s]</pre>\n",
+       "  </div>\n",
+       "  <div class=\"enlighten-bar\">\n",
+       "    <pre>Evaluating DE/rand/1 100%|<span class=\"enlighten-fg-green\">█████████████████████████████████████████████</span>| 10/10 [00:12&lt;00:00, 0.86/s]</pre>\n",
+       "  </div>\n",
+       "  <div class=\"enlighten-bar\">\n",
+       "    <pre>Evaluating DE/current-to-best/1 100%|<span class=\"enlighten-fg-green\">██████████████████████████████████</span>| 10/10 [00:16&lt;00:00, 0.62/s]</pre>\n",
+       "  </div>\n",
+       "  <div class=\"enlighten-bar\">\n",
+       "    <pre>Evaluating PSO 100%|<span class=\"enlighten-fg-green\">███████████████████████████████████████████████████</span>| 10/10 [00:12&lt;00:00, 0.86/s]</pre>\n",
+       "  </div>\n",
+       "</div>\n"
+      ],
+      "text/plain": [
+       "<IPython.core.display.HTML object>"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "# Initializer for sinlge-solution algorithms\n",
+    "single_initializer = UniformVectorInitializer(rastrigin_func.vecsize, rastrigin_func.low_lim, rastrigin_func.up_lim, pop_size=1)\n",
+    "\n",
+    "# Initializer for population-based algorithms\n",
+    "pop_initializer = UniformVectorInitializer(rastrigin_func.vecsize, rastrigin_func.low_lim, rastrigin_func.up_lim, pop_size=100)\n",
+    "\n",
+    "# Define strategies to be tested\n",
+    "strategies = [\n",
+    "    RandomSearch(pop_initializer),\n",
+    "    HillClimb(\n",
+    "        single_initializer,\n",
+    "        perturb_op=OperatorReal(\"RandNoise\", {\"method\": \"Gauss\", \"F\": 1e-3}),\n",
+    "        name=\"HillClimb\",\n",
+    "    ),\n",
+    "    SA(\n",
+    "        single_initializer,\n",
+    "        perturb_op=OperatorReal(\"RandNoise\", {\"method\": \"Gauss\", \"F\": 1e-3}),\n",
+    "        params={\"iter\": 100, \"temp_init\": 1, \"alpha\": 0.995},\n",
+    "        name=\"SA\",\n",
+    "    ),\n",
+    "    SA(\n",
+    "        pop_initializer,\n",
+    "        perturb_op=OperatorReal(\"RandNoise\", {\"method\": \"Gauss\", \"F\": 1e-3}),\n",
+    "        params={\"iter\": 100, \"temp_init\": 1, \"alpha\": 0.995},\n",
+    "        name=\"ParallelSA\",\n",
+    "    ),\n",
+    "    ES(\n",
+    "        pop_initializer,\n",
+    "        mutation_op=OperatorReal(\"RandNoise\", {\"method\": \"Gauss\", \"F\": 5e-3}),\n",
+    "        selection_op=SurvivorSelection(\"(m+n)\"),\n",
+    "        params={\"offspringSize\": 150},\n",
+    "        name=\"ES-(μ+λ)\",\n",
+    "    ),\n",
+    "    ES(\n",
+    "        pop_initializer,\n",
+    "        mutation_op=OperatorReal(\"RandNoise\", {\"method\": \"Gauss\", \"F\": 5e-2}),\n",
+    "        selection_op=SurvivorSelection(\"(m,n)\"),\n",
+    "        params={\"offspringSize\": 700},\n",
+    "        name=\"ES-(μ,λ)\",\n",
+    "    ),\n",
+    "    GA(\n",
+    "        pop_initializer,\n",
+    "        mutation_op=OperatorReal(\"RandNoise\", {\"method\": \"Gauss\", \"F\": 1e-3}),\n",
+    "        cross_op=OperatorReal(\"Multipoint\"),\n",
+    "        parent_sel_op=ParentSelection(\"Tournament\", {\"amount\": 60, \"p\": 0.1}),\n",
+    "        selection_op=SurvivorSelection(\"Elitism\", {\"amount\": 10}),\n",
+    "        params={\"pcross\": 0.9, \"pmut\": 0.1},\n",
+    "        name=\"GA\",\n",
+    "    ),\n",
+    "    DE(\n",
+    "        pop_initializer,\n",
+    "        OperatorReal(\"DE/best/1\", {\"F\": 0.8, \"Cr\": 0.8}),\n",
+    "        name=\"DE/best/1\",\n",
+    "    ),\n",
+    "    DE(\n",
+    "        pop_initializer,\n",
+    "        OperatorReal(\"DE/rand/1\", {\"F\": 0.8, \"Cr\": 0.8}),\n",
+    "        name=\"DE/rand/1\",\n",
+    "    ),\n",
+    "    DE(\n",
+    "        pop_initializer,\n",
+    "        OperatorReal(\"DE/current-to-best/1\", {\"F\": 0.8, \"Cr\": 0.8}),\n",
+    "        name=\"DE/current-to-best/1\",\n",
+    "    ),\n",
+    "    PSO(pop_initializer, {\"w\": 0.7, \"c1\": 1.5, \"c2\": 1.5}, name=\"PSO\"),\n",
+    "]\n",
+    "\n",
+    "# Instanciate the StrategySelection class\n",
+    "algorithm_search = StrategySelection(\n",
+    "    rastrigin_func,\n",
+    "    strategies,\n",
+    "    algorithm_params={\n",
+    "        \"stop_cond\": \"neval\",\n",
+    "        \"neval\": 5e4,\n",
+    "        \"verbose\": False,\n",
+    "    },\n",
+    "    params={\"verbose\": True, \"repetitions\": 10},\n",
+    ")\n",
+    "\n",
+    "# Launch all the algorithms and get the best overall solution\n",
+    "solution, best_fitness, report = algorithm_search.optimize()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 10,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "best solution: [-1.26199997e-09 -4.53316307e-10]\n",
+      "fitness: 0.0\n"
+     ]
+    }
+   ],
+   "source": [
+    "print(f\"best solution: {solution}\")\n",
+    "print(f\"fitness: {best_fitness}\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 11,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/html": [
+       "<div>\n",
+       "<style scoped>\n",
+       "    .dataframe tbody tr th:only-of-type {\n",
+       "        vertical-align: middle;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe tbody tr th {\n",
+       "        vertical-align: top;\n",
+       "    }\n",
+       "\n",
+       "    .dataframe thead th {\n",
+       "        text-align: right;\n",
+       "    }\n",
+       "</style>\n",
+       "<table border=\"1\" class=\"dataframe\">\n",
+       "  <thead>\n",
+       "    <tr style=\"text-align: right;\">\n",
+       "      <th></th>\n",
+       "      <th>name</th>\n",
+       "      <th>realtime_min</th>\n",
+       "      <th>realtime_avg</th>\n",
+       "      <th>realtime_max</th>\n",
+       "      <th>realtime_std</th>\n",
+       "      <th>cputime_min</th>\n",
+       "      <th>cputime_avg</th>\n",
+       "      <th>cputime_max</th>\n",
+       "      <th>cputime_std</th>\n",
+       "      <th>fitness_min</th>\n",
+       "      <th>fitness_avg</th>\n",
+       "      <th>fitness_max</th>\n",
+       "      <th>fitness_std</th>\n",
+       "    </tr>\n",
+       "  </thead>\n",
+       "  <tbody>\n",
+       "    <tr>\n",
+       "      <th>7</th>\n",
+       "      <td>DE/best/1</td>\n",
+       "      <td>1.519729</td>\n",
+       "      <td>1.548633</td>\n",
+       "      <td>1.608609</td>\n",
+       "      <td>0.025094</td>\n",
+       "      <td>1.519781</td>\n",
+       "      <td>1.565137</td>\n",
+       "      <td>1.628799</td>\n",
+       "      <td>0.035645</td>\n",
+       "      <td>0.000000e+00</td>\n",
+       "      <td>0.000000e+00</td>\n",
+       "      <td>0.000000e+00</td>\n",
+       "      <td>0.000000</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>8</th>\n",
+       "      <td>DE/rand/1</td>\n",
+       "      <td>1.123156</td>\n",
+       "      <td>1.167058</td>\n",
+       "      <td>1.297242</td>\n",
+       "      <td>0.058852</td>\n",
+       "      <td>1.128656</td>\n",
+       "      <td>1.174887</td>\n",
+       "      <td>1.296135</td>\n",
+       "      <td>0.055647</td>\n",
+       "      <td>0.000000e+00</td>\n",
+       "      <td>0.000000e+00</td>\n",
+       "      <td>0.000000e+00</td>\n",
+       "      <td>0.000000</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>9</th>\n",
+       "      <td>DE/current-to-best/1</td>\n",
+       "      <td>1.584857</td>\n",
+       "      <td>1.613552</td>\n",
+       "      <td>1.653009</td>\n",
+       "      <td>0.020102</td>\n",
+       "      <td>1.597510</td>\n",
+       "      <td>1.635562</td>\n",
+       "      <td>1.727576</td>\n",
+       "      <td>0.037934</td>\n",
+       "      <td>0.000000e+00</td>\n",
+       "      <td>0.000000e+00</td>\n",
+       "      <td>0.000000e+00</td>\n",
+       "      <td>0.000000</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>10</th>\n",
+       "      <td>PSO</td>\n",
+       "      <td>1.115229</td>\n",
+       "      <td>1.156385</td>\n",
+       "      <td>1.321878</td>\n",
+       "      <td>0.061156</td>\n",
+       "      <td>1.115460</td>\n",
+       "      <td>1.169345</td>\n",
+       "      <td>1.320756</td>\n",
+       "      <td>0.061212</td>\n",
+       "      <td>0.000000e+00</td>\n",
+       "      <td>0.000000e+00</td>\n",
+       "      <td>0.000000e+00</td>\n",
+       "      <td>0.000000</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>6</th>\n",
+       "      <td>GA</td>\n",
+       "      <td>1.769814</td>\n",
+       "      <td>1.803421</td>\n",
+       "      <td>1.826834</td>\n",
+       "      <td>0.016155</td>\n",
+       "      <td>1.768550</td>\n",
+       "      <td>1.802076</td>\n",
+       "      <td>1.825370</td>\n",
+       "      <td>0.016126</td>\n",
+       "      <td>1.526814e-10</td>\n",
+       "      <td>1.526814e-10</td>\n",
+       "      <td>1.526814e-10</td>\n",
+       "      <td>0.000000</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>5</th>\n",
+       "      <td>ES-(μ,λ)</td>\n",
+       "      <td>0.930633</td>\n",
+       "      <td>0.943650</td>\n",
+       "      <td>0.991546</td>\n",
+       "      <td>0.017881</td>\n",
+       "      <td>0.936875</td>\n",
+       "      <td>0.948075</td>\n",
+       "      <td>0.992188</td>\n",
+       "      <td>0.016175</td>\n",
+       "      <td>1.081313e-06</td>\n",
+       "      <td>1.081313e-06</td>\n",
+       "      <td>1.081313e-06</td>\n",
+       "      <td>0.000000</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>0</th>\n",
+       "      <td>RandomSearch</td>\n",
+       "      <td>1.252175</td>\n",
+       "      <td>1.293047</td>\n",
+       "      <td>1.346313</td>\n",
+       "      <td>0.033340</td>\n",
+       "      <td>1.253476</td>\n",
+       "      <td>1.291402</td>\n",
+       "      <td>1.343438</td>\n",
+       "      <td>0.032443</td>\n",
+       "      <td>3.427405e-03</td>\n",
+       "      <td>8.270857e-03</td>\n",
+       "      <td>1.034662e-02</td>\n",
+       "      <td>0.003342</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>4</th>\n",
+       "      <td>ES-(μ+λ)</td>\n",
+       "      <td>0.922964</td>\n",
+       "      <td>0.929995</td>\n",
+       "      <td>0.944922</td>\n",
+       "      <td>0.005929</td>\n",
+       "      <td>0.924595</td>\n",
+       "      <td>0.945595</td>\n",
+       "      <td>1.008896</td>\n",
+       "      <td>0.027014</td>\n",
+       "      <td>8.391705e-07</td>\n",
+       "      <td>3.979843e-01</td>\n",
+       "      <td>1.989918e+00</td>\n",
+       "      <td>0.839023</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>3</th>\n",
+       "      <td>ParallelSA</td>\n",
+       "      <td>0.906290</td>\n",
+       "      <td>0.918373</td>\n",
+       "      <td>0.938245</td>\n",
+       "      <td>0.010232</td>\n",
+       "      <td>0.907828</td>\n",
+       "      <td>0.938450</td>\n",
+       "      <td>1.008967</td>\n",
+       "      <td>0.034142</td>\n",
+       "      <td>2.701125e-04</td>\n",
+       "      <td>5.977198e-01</td>\n",
+       "      <td>9.949638e-01</td>\n",
+       "      <td>0.426150</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>1</th>\n",
+       "      <td>HillClimb</td>\n",
+       "      <td>1.124578</td>\n",
+       "      <td>1.142825</td>\n",
+       "      <td>1.157456</td>\n",
+       "      <td>0.010962</td>\n",
+       "      <td>1.124174</td>\n",
+       "      <td>1.143407</td>\n",
+       "      <td>1.158496</td>\n",
+       "      <td>0.011480</td>\n",
+       "      <td>9.949591e-01</td>\n",
+       "      <td>4.676294e+00</td>\n",
+       "      <td>1.691420e+01</td>\n",
+       "      <td>4.475443</td>\n",
+       "    </tr>\n",
+       "    <tr>\n",
+       "      <th>2</th>\n",
+       "      <td>SA</td>\n",
+       "      <td>1.134191</td>\n",
+       "      <td>1.156348</td>\n",
+       "      <td>1.186226</td>\n",
+       "      <td>0.016692</td>\n",
+       "      <td>1.136885</td>\n",
+       "      <td>1.157306</td>\n",
+       "      <td>1.188024</td>\n",
+       "      <td>0.016494</td>\n",
+       "      <td>7.959662e+00</td>\n",
+       "      <td>1.174046e+01</td>\n",
+       "      <td>2.487372e+01</td>\n",
+       "      <td>6.375946</td>\n",
+       "    </tr>\n",
+       "  </tbody>\n",
+       "</table>\n",
+       "</div>"
+      ],
+      "text/plain": [
+       "                    name  realtime_min  realtime_avg  realtime_max  \\\n",
+       "7              DE/best/1      1.519729      1.548633      1.608609   \n",
+       "8              DE/rand/1      1.123156      1.167058      1.297242   \n",
+       "9   DE/current-to-best/1      1.584857      1.613552      1.653009   \n",
+       "10                   PSO      1.115229      1.156385      1.321878   \n",
+       "6                     GA      1.769814      1.803421      1.826834   \n",
+       "5               ES-(μ,λ)      0.930633      0.943650      0.991546   \n",
+       "0           RandomSearch      1.252175      1.293047      1.346313   \n",
+       "4               ES-(μ+λ)      0.922964      0.929995      0.944922   \n",
+       "3             ParallelSA      0.906290      0.918373      0.938245   \n",
+       "1              HillClimb      1.124578      1.142825      1.157456   \n",
+       "2                     SA      1.134191      1.156348      1.186226   \n",
+       "\n",
+       "    realtime_std  cputime_min  cputime_avg  cputime_max  cputime_std  \\\n",
+       "7       0.025094     1.519781     1.565137     1.628799     0.035645   \n",
+       "8       0.058852     1.128656     1.174887     1.296135     0.055647   \n",
+       "9       0.020102     1.597510     1.635562     1.727576     0.037934   \n",
+       "10      0.061156     1.115460     1.169345     1.320756     0.061212   \n",
+       "6       0.016155     1.768550     1.802076     1.825370     0.016126   \n",
+       "5       0.017881     0.936875     0.948075     0.992188     0.016175   \n",
+       "0       0.033340     1.253476     1.291402     1.343438     0.032443   \n",
+       "4       0.005929     0.924595     0.945595     1.008896     0.027014   \n",
+       "3       0.010232     0.907828     0.938450     1.008967     0.034142   \n",
+       "1       0.010962     1.124174     1.143407     1.158496     0.011480   \n",
+       "2       0.016692     1.136885     1.157306     1.188024     0.016494   \n",
+       "\n",
+       "     fitness_min   fitness_avg   fitness_max  fitness_std  \n",
+       "7   0.000000e+00  0.000000e+00  0.000000e+00     0.000000  \n",
+       "8   0.000000e+00  0.000000e+00  0.000000e+00     0.000000  \n",
+       "9   0.000000e+00  0.000000e+00  0.000000e+00     0.000000  \n",
+       "10  0.000000e+00  0.000000e+00  0.000000e+00     0.000000  \n",
+       "6   1.526814e-10  1.526814e-10  1.526814e-10     0.000000  \n",
+       "5   1.081313e-06  1.081313e-06  1.081313e-06     0.000000  \n",
+       "0   3.427405e-03  8.270857e-03  1.034662e-02     0.003342  \n",
+       "4   8.391705e-07  3.979843e-01  1.989918e+00     0.839023  \n",
+       "3   2.701125e-04  5.977198e-01  9.949638e-01     0.426150  \n",
+       "1   9.949591e-01  4.676294e+00  1.691420e+01     4.475443  \n",
+       "2   7.959662e+00  1.174046e+01  2.487372e+01     6.375946  "
+      ]
+     },
+     "execution_count": 11,
+     "metadata": {},
+     "output_type": "execute_result"
+    }
+   ],
+   "source": [
+    "report.sort_values(\"fitness_avg\")"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "evolcomp",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.4"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/tutorials/complete_tutorial.ipynb b/tutorials/complete_tutorial.ipynb
index ae9fa54..6eb4c92 100644
--- a/tutorials/complete_tutorial.ipynb
+++ b/tutorials/complete_tutorial.ipynb
@@ -48,9 +48,7 @@
     "    Takes an integer-valued vector and outputs the number of treasures collected.\n",
     "    \"\"\"\n",
     "\n",
-    "    def __init__(\n",
-    "        self, n_moves, n_tresures, grid_size=(10, 10), name: str = \"Tresure Hunt\"\n",
-    "    ):\n",
+    "    def __init__(self, n_moves, n_tresures, grid_size=(10, 10), name: str = \"Tresure Hunt\"):\n",
     "        self.n_moves = n_moves\n",
     "        self.n_tresures = n_tresures\n",
     "        self.grid_size = np.array(grid_size)\n",
@@ -91,12 +89,7 @@
     "            new_pos = pos + move_set[move]\n",
     "\n",
     "            # Check if it is on bounds\n",
-    "            if (\n",
-    "                new_pos[0] >= 0\n",
-    "                and new_pos[0] < self.grid_size[0]\n",
-    "                and new_pos[1] >= 0\n",
-    "                and new_pos[1] < self.grid_size[1]\n",
-    "            ):\n",
+    "            if new_pos[0] >= 0 and new_pos[0] < self.grid_size[0] and new_pos[1] >= 0 and new_pos[1] < self.grid_size[1]:\n",
     "                pos = new_pos\n",
     "\n",
     "                # Collect the treasure if there is one\n",
@@ -135,12 +128,7 @@
     "\n",
     "        for move in move_vector:\n",
     "            new_pos = pos + move_set[move]\n",
-    "            if (\n",
-    "                new_pos[0] >= 0\n",
-    "                and new_pos[0] < self.grid_size[0]\n",
-    "                and new_pos[1] >= 0\n",
-    "                and new_pos[1] < self.grid_size[1]\n",
-    "            ):\n",
+    "            if new_pos[0] >= 0 and new_pos[0] < self.grid_size[0] and new_pos[1] >= 0 and new_pos[1] < self.grid_size[1]:\n",
     "                pos = new_pos\n",
     "                if grid_copy[*pos] == 0:\n",
     "                    grid_copy[*pos] = 2\n",
@@ -300,9 +288,7 @@
     "params = {\"stop_cond\": \"neval\", \"neval\": 50000, \"verbose\": True, \"v_timer\": 2}\n",
     "\n",
     "# Define a population initialization module\n",
-    "pop_init = UniformVectorInitializer(\n",
-    "    objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=1, dtype=int\n",
-    ")\n",
+    "pop_init = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=1, dtype=int)\n",
     "\n",
     "# Define the operators to use\n",
     "mutate_op = OperatorInt(\"MutSample\", {\"method\": \"Uniform\", \"Low\": 0, \"Up\": 3, \"N\": 8})\n",
@@ -465,9 +451,7 @@
     "pop_size = 100\n",
     "\n",
     "# Define a population initialization module\n",
-    "pop_init = UniformVectorInitializer(\n",
-    "    objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=int\n",
-    ")\n",
+    "pop_init = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=int)\n",
     "\n",
     "# Define the operators to use\n",
     "mutate_op = OperatorInt(\"MutSample\", {\"method\": \"Uniform\", \"Low\": 0, \"Up\": 3, \"N\": 8})\n",
@@ -480,9 +464,7 @@
     "surv_sel = SurvivorSelection(\"KeepBest\")\n",
     "\n",
     "# Instanciate the algorithm\n",
-    "algorithm = GA(\n",
-    "    pop_init, mutate_op, cross_op, parent_sel, surv_sel, {\"pmut\": 0.1, \"pcross\": 0.8}\n",
-    ")\n",
+    "algorithm = GA(pop_init, mutate_op, cross_op, parent_sel, surv_sel, {\"pmut\": 0.1, \"pcross\": 0.8})\n",
     "\n",
     "# Define a Search method\n",
     "search = GeneralSearch(objfunc, algorithm, params)\n",
@@ -628,9 +610,7 @@
     "pop_size = 100\n",
     "\n",
     "# Define a population initialization module\n",
-    "pop_init = UniformVectorInitializer(\n",
-    "    objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=int\n",
-    ")\n",
+    "pop_init = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=pop_size, dtype=int)\n",
     "\n",
     "# Define the operators to use\n",
     "mutate_op = OperatorInt(\"MutSample\", {\"method\": \"Uniform\", \"Low\": 0, \"Up\": 3, \"N\": 8})\n",
diff --git a/tutorials/create_encoding.ipynb b/tutorials/create_encoding.ipynb
index 8e55ff3..aefb005 100644
--- a/tutorials/create_encoding.ipynb
+++ b/tutorials/create_encoding.ipynb
@@ -59,6 +59,7 @@
     "def encode_round(vector):\n",
     "    return np.round(vector).astype(int)\n",
     "\n",
+    "\n",
     "def decode_round(vector):\n",
     "    return vector.astype(float)"
    ]
@@ -137,9 +138,7 @@
     "objfunc = Sphere(3)\n",
     "\n",
     "# Define a population initialization module\n",
-    "pop_init = UniformVectorInitializer(\n",
-    "    objfunc.vecsize, objfunc.low_lim*100, objfunc.up_lim*100, pop_size=100, encoding=round_encoding\n",
-    ")\n",
+    "pop_init = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim * 100, objfunc.up_lim * 100, pop_size=100, encoding=round_encoding)\n",
     "\n",
     "# Define the operator to be used\n",
     "mutation = OperatorReal(\"MutNoise\", {\"method\": \"Normal\", \"F\": 1, \"N\": 1})\n",
@@ -219,6 +218,7 @@
     "    def decode(self, genotype):\n",
     "        return np.round(genotype).astype(int)\n",
     "\n",
+    "\n",
     "custom_encoding = RoundEncoding()"
    ]
   },
@@ -280,9 +280,7 @@
     "objfunc = Sphere(3)\n",
     "\n",
     "# Define a population initialization module\n",
-    "pop_init = UniformVectorInitializer(\n",
-    "    objfunc.vecsize, objfunc.low_lim*100, objfunc.up_lim*100, pop_size=100, encoding=custom_encoding\n",
-    ")\n",
+    "pop_init = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim * 100, objfunc.up_lim * 100, pop_size=100, encoding=custom_encoding)\n",
     "\n",
     "# Define the operator to be used\n",
     "mutation = OperatorReal(\"MutNoise\", {\"method\": \"Normal\", \"F\": 1, \"N\": 1})\n",
diff --git a/tutorials/create_operators.ipynb b/tutorials/create_operators.ipynb
index 86083f9..6c8757c 100644
--- a/tutorials/create_operators.ipynb
+++ b/tutorials/create_operators.ipynb
@@ -77,9 +77,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "operator = OperatorLambda(\n",
-    "    mutate_custom, name=\"GaussianNoise_lambda\", params={\"F\": 1e-3}\n",
-    ")"
+    "operator = OperatorLambda(mutate_custom, name=\"GaussianNoise_lambda\", params={\"F\": 1e-3})"
    ]
   },
   {
@@ -153,9 +151,7 @@
     "objfunc = Sphere(3)\n",
     "\n",
     "# Define a population initialization module\n",
-    "pop_init = UniformVectorInitializer(\n",
-    "    objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=1\n",
-    ")\n",
+    "pop_init = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=1)\n",
     "\n",
     "# Instanciate the newly created algorithm\n",
     "algorithm = HillClimb(pop_init, operator)\n",
@@ -322,9 +318,7 @@
     "objfunc = Sphere(3)\n",
     "\n",
     "# Define a population initialization module\n",
-    "pop_init = UniformVectorInitializer(\n",
-    "    objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=1\n",
-    ")\n",
+    "pop_init = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=1)\n",
     "\n",
     "# Instanciate the newly created algorithm\n",
     "algorithm = HillClimb(pop_init, operator)\n",
diff --git a/tutorials/create_strategy.ipynb b/tutorials/create_strategy.ipynb
index 0cc6a8f..56fd641 100644
--- a/tutorials/create_strategy.ipynb
+++ b/tutorials/create_strategy.ipynb
@@ -83,9 +83,7 @@
     "        self.offspring_size = offspring_size\n",
     "\n",
     "        # Create the mutation operator for use later\n",
-    "        self.mutate_op = OperatorReal(\n",
-    "            \"MutRand\", {\"method\": \"Gauss\", \"F\": sigma, \"N\": 1}\n",
-    "        )\n",
+    "        self.mutate_op = OperatorReal(\"MutRand\", {\"method\": \"Gauss\", \"F\": sigma, \"N\": 1})\n",
     "\n",
     "        # Force pop_init to agree with the population size, this is not mandatory.\n",
     "        pop_init.pop_size = parent_size\n",
@@ -117,9 +115,7 @@
     "            indiv = random.choice(parent_list)\n",
     "\n",
     "            # Apply the operator.\n",
-    "            new_indiv = self.mutate_op(\n",
-    "                indiv, parent_list, objfunc, self.best, self.pop_init\n",
-    "            )\n",
+    "            new_indiv = self.mutate_op(indiv, parent_list, objfunc, self.best, self.pop_init)\n",
     "\n",
     "            # Repair the solution to keep the vector within the bounds\n",
     "            # of the problem.\n",
@@ -528,9 +524,7 @@
     "objfunc = Sphere(3)\n",
     "\n",
     "# Define a population initialization module\n",
-    "pop_init = UniformVectorInitializer(\n",
-    "    objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=n_parents\n",
-    ")\n",
+    "pop_init = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=n_parents)\n",
     "\n",
     "# Instanciate the newly created algorithm\n",
     "algorithm = BasicES(pop_init, n_parents, n_offspring, 1e-3)\n",
@@ -577,13 +571,6 @@
     "print(f\"Solution fitness: {best_fitness}\\n\")\n",
     "search.display_report()"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {


Test patch:
diff --git a/examples/adaption_test.py b/examples/adaption_test.py
new file mode 100644
index 0000000..fc0b1a9
--- /dev/null
+++ b/examples/adaption_test.py
@@ -0,0 +1,120 @@
+from metaheuristic_designer import ObjectiveFunc, ParamScheduler
+from metaheuristic_designer.algorithms import GeneralAlgorithm, MemeticAlgorithm
+from metaheuristic_designer.operators import OperatorReal, OperatorAdaptative
+from metaheuristic_designer.initializers import UniformVectorInitializer
+from metaheuristic_designer.selectionMethods import ParentSelection, SurvivorSelection
+from metaheuristic_designer.encodings import AdaptionEncoding
+from metaheuristic_designer.strategies import *
+
+from metaheuristic_designer.benchmarks import *
+
+import argparse
+
+from copy import copy
+import scipy as sp
+import numpy as np
+
+
+class STDAdaptEncoding(AdaptionEncoding):
+    def decode_param(self, genotype):
+        # print(genotype)
+        return {"F": np.maximum(1e-7, self.decode_param_vec(genotype))}
+        
+
+def run_algorithm(alg_name, memetic, save_state):
+    params = {
+        "stop_cond": "convergence or time_limit",
+        "progress_metric": "time_limit",
+        "time_limit": 100.0,
+        "cpu_time_limit": 100.0,
+        "ngen": 1000,
+        "neval": 3e6,
+        "fit_target": 1e-10,
+        "patience": 10,
+        "verbose": True,
+        "v_timer": 0.5,
+    }
+
+    # objfunc = Sphere(30, "min")
+    objfunc = Rastrigin(2, "min")
+    # objfunc = Rosenbrock(2, "min")
+    # objfunc = Ackley(30, "min")
+    # objfunc = Weierstrass(30, "min")
+    # objfunc = HappyCat(3, "min")
+    
+    # mutation_op = OperatorReal("RandNoise", {"method": "Gauss"})
+    mutation_op = OperatorReal("MutNoise", {"method": "Gauss", "N": 1})
+
+    param_op = OperatorReal("Mutate1Sigma", {"tau": 1/np.sqrt(objfunc.vecsize), "epsilon": 1e-7})
+    adaption_encoding = STDAdaptEncoding(objfunc.vecsize, nparams=1)
+    # param_op = OperatorReal("MutateNSigmas", {"tau": 1/np.sqrt(2+objfunc.vecsize), "tau_multiple": 0.5/np.sqrt(objfunc.vecsize), "epsilon": 1e-7})
+    # adaption_encoding = STDAdaptEncoding(objfunc.vecsize, nparams=objfunc.vecsize)
+
+    ada_mutation_op = OperatorAdaptative(mutation_op, param_op, adaption_encoding)
+
+    pop_initializer = UniformVectorInitializer(objfunc.vecsize, objfunc.low_lim, objfunc.up_lim, pop_size=100, encoding=adaption_encoding)    
+
+    cross_op = OperatorReal("Multipoint")
+
+    parent_sel_op = ParentSelection("Nothing")
+    selection_op = SurvivorSelection("(m+n)")
+
+    search_strat = ES(
+        pop_initializer,
+        ada_mutation_op,
+        cross_op,
+        parent_sel_op,
+        selection_op,
+        {"offspringSize": 700},
+        name="Adaptative-ES"
+    )
+
+    alg = GeneralAlgorithm(objfunc, search_strat, params=params)
+
+    ind, fit = alg.optimize()
+    print(f"full solution vector: {ind}")
+    print(f"solution: {adaption_encoding.decode(ind)}")
+    print(f"params: {adaption_encoding.decode_param(ind)}")
+    alg.display_report(show_plots=True)
+
+    if save_state:
+        alg.store_state("./examples/results/test.json", readable=True, show_pop=True)
+
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument("-a", "--algorithm", dest="alg", help="Specify an algorithm")
+    parser.add_argument(
+        "-m",
+        "--memetic",
+        dest="mem",
+        action="store_true",
+        help="Does local search after mutation",
+    )
+    parser.add_argument(
+        "-s",
+        "--save-state",
+        dest="save_state",
+        action="store_true",
+        help="Saves the state of the search strategy",
+    )
+    args = parser.parse_args()
+
+    algorithm_name = "ES"
+    mem = False
+    save_state = False
+
+    if args.alg:
+        algorithm_name = args.alg
+
+    if args.mem:
+        mem = True
+
+    if args.save_state:
+        save_state = True
+
+    run_algorithm(alg_name=algorithm_name, memetic=mem, save_state=save_state)
+
+
+if __name__ == "__main__":
+    main()