diff --git a/selection.py b/selection.py
new file mode 100644
index 0000000000..99c6614ed9
--- /dev/null
+++ b/selection.py
@@ -0,0 +1,211 @@
+import numpy as np
+
+# # dummy values for testing
+# lowerbound = -1
+# upperbound = 1
+# population_size = 200
+# n_vars = 265
+# generation_num=100
+
+# # population is (100, 265) with real values ranging from -1 to 1
+# population = np.random.uniform(lowerbound, upperbound, (population_size, n_vars))
+# fitness_population = np.random.uniform(50, 100, population_size)
+# # end dummy values for testing
+
+
+
+def roulette_wheel_selection(population, fitnesses):
+    selected_individuals = []
+    selected_fitness = []
+    n = population.shape[0]
+
+    total_fitness = np.sum(fitnesses)
+
+    for _ in range(n): 
+        alpha = np.random.uniform(0, total_fitness)
+        cumulative_sum = 0
+        j = 0
+        
+        # selection
+        while cumulative_sum < alpha and j < n:
+            cumulative_sum += fitnesses[j]
+            j += 1
+            
+        selected_individuals.append(population[j-1])
+        selected_fitness.append(fitnesses[j-1])
+    
+    return np.array(selected_individuals), np.array(selected_fitness)
+
+
+# Does not really work
+# def stochastic_universal_sampling(population, fitnesses):
+#     selected_individuals = []
+#     selected_fitness = []
+#     n = population.shape[0]
+
+#     mean_fitness = np.mean(fitnesses)
+#     alpha = np.random.rand()
+
+#     cumulative_sum = fitnesses[0]
+#     delta = alpha * mean_fitness
+
+#     j = 0
+#     while j < n - 1:
+#         if delta < cumulative_sum:
+#             selected_individuals.append(population[j])
+#             selected_fitness.append(fitnesses[j])
+#             delta += cumulative_sum
+#             break
+#         else:
+#             j += 1
+#             cumulative_sum += fitnesses[j] 
+
+#     return np.array(selected_individuals), np.array(selected_fitness)
+
+
+
+
+def linear_rank_selection(population, fitnesses):
+    selected_individuals = []
+    selected_fitness = []
+    n = population.shape[0]
+
+    sorted_indices = np.argsort(fitnesses)[::-1]
+    ranks = np.zeros_like(sorted_indices)
+
+    for rank, index in enumerate(sorted_indices, start=1):
+        ranks[index] = rank
+
+    probs = ranks / (n * (n - 1))
+
+    value = 1 / (n - 2.001)
+
+    while len(selected_individuals) < n:
+        for i in range(n):
+            alpha = np.random.uniform(0, value)
+            for j in range(n):
+                if probs[j] <= alpha:
+                    if len(selected_individuals) < n:
+                        selected_individuals.append(population[j])
+                        selected_fitness.append(fitnesses[j])
+                    break
+    
+    return np.array(selected_individuals), np.array(selected_fitness)
+
+
+
+def exponential_rank_selection(population, fitnesses):
+    selected_individuals = []
+    selected_fitness = []
+    n = population.shape[0]
+
+    sorted_indices = np.argsort(fitnesses)[::-1]
+    ranks = np.zeros_like(sorted_indices)
+
+    for rank, index in enumerate(sorted_indices, start=1):
+        ranks[index] = rank
+
+    probs = np.zeros(n)
+    c = (n * 2 * (n - 1)) / (6 * (n - 1) + n)
+    for i in range(n):
+        probs[i] = 1.0 * np.exp( - ranks[i] / c)
+
+
+    for i in range(n):
+        alpha = np.random.uniform(1 / 9 * c, 2 / c)
+        for j in range(n):
+            if probs[j] <= alpha:
+                selected_individuals.append(population[j])
+                selected_fitness.append(fitnesses[j])
+                break
+
+    return np.array(selected_individuals), np.array(selected_fitness)
+
+
+
+def tournament_selection(population, fitnesses):
+    selected_individuals = []
+    selected_fitness = []
+    n = population.shape[0]
+
+    k = 20
+    for _ in range(n):
+        temp = list(zip(population, fitnesses))
+        np.random.shuffle(temp)
+        res1, res2 = zip(*temp)
+        shuffled_population, shuffled_fitnesses = np.array(res1), np.array(res2)
+
+        # compare k individuals
+        best_out_of_k = np.argmax(shuffled_fitnesses[0:k])
+        selected_individuals.append(shuffled_population[best_out_of_k])
+        selected_fitness.append(shuffled_fitnesses[best_out_of_k])
+
+    return np.array(selected_individuals), np.array(selected_fitness)
+
+
+
+def selection_score(population, fitness_population, generation):
+    ''' 
+    Selection based on the dynamic approach from
+    'Parent Selection Operators for Genetic Algorithms'
+    Input: current population, current generation number, fitness of current population
+    Output: selection criterion
+    '''
+    best_idx = np.argmax(fitness_population)
+    best = population[best_idx]
+    criteria1 = 0
+    pop_size = population.shape[0]
+
+    # Hamming distance is binary so we use Euclidean distance instead
+    for individual in population:
+        criteria1 += np.linalg.norm(best - individual)
+    criteria1 /=  pop_size # normalize
+    criteria1 = np.exp(- criteria1 / generation) # decrease over generations
+
+    max_fitness = np.max(fitness_population)
+    min_fitness = np.min(fitness_population)
+    criteria2 = max_fitness / (max_fitness **2 + min_fitness**2) # maximisation problem
+    
+    criterion = 1/generation * criteria1 + ((generation-1)/generation) * criteria2
+
+    return criterion
+
+
+
+def dynamic_selection(population, fitnesses, generation):
+
+    rws_population, rws_fitness = roulette_wheel_selection(population, fitnesses)
+    rws_score = selection_score(rws_population, rws_fitness, generation)
+
+    # sus_population, sus_fitness = stochastic_universal_sampling(population, fitnesses)
+    # sus_score = selection_score(sus_population, sus_fitness, generation)
+
+    lrs_population, lrs_fitness = linear_rank_selection(population, fitnesses)
+    lrs_score = selection_score(lrs_population, lrs_fitness, generation)
+
+    ers_population, ers_fitness = exponential_rank_selection(population, fitnesses)
+    ers_score = selection_score(ers_population, ers_fitness, generation)
+
+    tos_population, tos_fitness = tournament_selection(population, fitnesses)
+    tos_score = selection_score(tos_population, tos_fitness, generation)
+
+    scores = [rws_score, lrs_score, ers_score, tos_score]
+    new_populations = [rws_population, lrs_population, ers_population, tos_population]
+    new_fitnesses = [rws_fitness, lrs_fitness, ers_fitness, tos_fitness]
+    best = np.argmax(scores)
+    # print(len(new_fitnesses[0]), len(new_fitnesses[1]), len(new_fitnesses[2]), len(new_fitnesses[3]))
+
+    return new_populations[best], new_fitnesses[best]
+
+
+# leave out truncation selection because it is not often used in practice and only for 
+# very large populations
+
+# for generation in range(1,100):
+#     pop, pop_fit = dynamic_selection(population, fitness_population, generation)
+#     print(generation, np.mean(pop_fit), np.std(pop_fit))
+# print(parent_selection(population, generation_num, fitness_population))
+# print(roulette_wheel_selection(population, fitness_population).shape)
+# dynamic_selection(population, fitness_population, generation_num)
+
+# print(population[0], population[1])
\ No newline at end of file
diff --git a/specialist_arco.py b/specialist_arco.py
new file mode 100644
index 0000000000..abf898faba
--- /dev/null
+++ b/specialist_arco.py
@@ -0,0 +1,226 @@
+import sys
+
+from evoman.environment import Environment
+from demo_controller import player_controller
+
+# imports other libs
+import numpy as np
+import os
+import tqdm
+import multiprocessing as mp
+import dill
+from multiprocessing import Pool
+
+# Set multiprocessing to use dill
+import multiprocessing.reduction
+multiprocessing.reduction.ForkingPickler = dill
+
+# runs simulation
+def simulation(env,x):
+    f,p,e,t = env.play(pcont=x)
+    return f
+
+# # evaluation
+# def evaluate(env, x):
+#     return np.array(list(map(lambda y: simulation(env,y), x)))
+
+class Specialist():
+    def __init__(self, n_hidden_neurons=10,
+                 experiment_name = 'optimization_test_arco',
+                 upperbound = 1,
+                 lowerbound = -1,
+                 population_size = 100,
+                 ) -> None:
+        self.headless = True
+        if self.headless:
+            os.environ["SDL_VIDEODRIVER"] = "dummy"
+
+        self.experiment_name = experiment_name
+        if not os.path.exists(self.experiment_name):
+            os.makedirs(self.experiment_name)
+
+        self.n_hidden_neurons = n_hidden_neurons
+
+        self.env = Environment(experiment_name=experiment_name,
+                    enemies=[2],
+                    playermode="ai",
+                    player_controller=player_controller(n_hidden_neurons), # you  can insert your own controller here
+                    enemymode="static",
+                    level=2,
+                    speed="fastest",
+                    visuals=False)
+        self.n_vars = (self.env.get_num_sensors()+1)*n_hidden_neurons + (n_hidden_neurons+1)*5
+        self.upperbound = upperbound
+        self.lowerbound = lowerbound
+        self.population_size = population_size
+
+    def simulation(self, neuron_values):
+        f, p, e, t = self.env.play(pcont=neuron_values)
+        return f
+
+    def fitness_eval(self, population):
+        return np.array([simulation(self.env, individual) for individual in population])
+    
+    
+    # def fitness_eval(self, population):
+    #     with Pool(mp.cpu_count()) as pool:
+    #         fitness_population = pool.starmap(self.simulation, [(self.env, individual) for individual in population])
+    #     return np.array(fitness_population)
+
+    def initialize(self):
+        return np.random.uniform(self.lowerbound,
+                                self.upperbound,
+                                (self.population_size, self.n_vars))
+
+    def mutation(self, child, p_mutation=0.2):
+
+        if np.random.uniform() > p_mutation:
+            #no mutation
+            return child
+        else:
+            child = np.array(child)
+            swap1, swap2 = np.random.choice(np.arange(1,10), size=2)
+            child[[swap1, swap2]] = child[[swap2, swap1]]
+            child_mutated = list(child)
+        return child_mutated
+
+    def limits(self, x):
+
+        if x>self.upperbound:
+            return self.upperbound
+        elif x<self.lowerbound:
+            return self.lowerbound
+        else:
+            return x
+
+    def tournament(self, pop):
+        c1 =  np.random.randint(0,pop.shape[0], 1)
+        c2 =  np.random.randint(0,pop.shape[0], 1)
+
+        c1_fit = self.fitness_eval(pop[c1])
+        c2_fit = self.fitness_eval(pop[c2])
+
+        if c1_fit > c2_fit:
+            return pop[c1][0]
+        else:
+            return pop[c2][0]
+
+    def crossover(self, pop, p_mutation):
+
+        total_offspring = np.zeros((0,self.n_vars))
+
+        for p in range(0, pop.shape[0], 2):
+            p1 = self.tournament(pop)
+            p2 = self.tournament(pop)
+
+            n_offspring = np.random.randint(1,3+1, 1)[0]
+            offspring = np.zeros( (n_offspring, self.n_vars) )
+
+            for f in range(0,n_offspring):
+
+                cross_prop = np.random.uniform(0,1)
+                offspring[f] = p1*cross_prop+p2*(1-cross_prop)
+
+                # mutation
+                for i in range(0,len(offspring[f])):
+                    if np.random.uniform(0 ,1)<=p_mutation:
+                        offspring[f][i] =   offspring[f][i]+np.random.normal(0, 1)
+
+                offspring[f] = np.array(list(map(lambda y: self.limits(y), offspring[f])))
+
+                total_offspring = np.vstack((total_offspring, offspring[f]))
+
+        return total_offspring
+
+    def normalize(self, pop, fit):
+
+        if (max(fit) - min(fit) ) > 0:
+            x_norm = (pop - min(fit))/(max(fit) - min(fit))
+        else:
+            x_norm = 0
+        if x_norm <= 0:
+            x_norm = 0.0000000001
+        return x_norm
+
+    def selection(self, new_population, new_fitness_population):
+        # TODO: find a better metric
+        p = np.array([self.normalize(pop, new_fitness_population) for pop in range(len(new_population))])
+        p = p / p.sum()
+        # print(sorted(p))
+        # chosen_idx = np.random.choice(np.arange(new_population.shape[0]),
+        #                         size=self.population_size,
+        #                         replace=False,
+        #                         p=p)
+
+        # chosen_idx = np.random.choice(np.arange(new_population.shape[0]),
+        #                         size=self.population_size,
+        #                         replace=True)
+        # return new_population[chosen_idx, :], new_fitness_population[chosen_idx]
+
+        fit_pop_cp = new_fitness_population
+        fit_pop_norm =  np.array(list(map(lambda y: self.normalize(y,fit_pop_cp), new_fitness_population))) # avoiding negative probabilities, as fitness is ranges from negative numbers
+        probs = (fit_pop_norm)/(fit_pop_norm).sum()
+        chosen = np.random.choice(new_population.shape[0], self.population_size , p=probs, replace=False)
+        print(chosen)
+        chosen = np.append(chosen[1:],np.argmax(new_fitness_population))
+        print(chosen)
+        print(chosen.shape)
+        pop = new_population[chosen]
+        fit_pop = new_fitness_population[chosen]
+
+        return pop, fit_pop
+
+
+    def train(self, total_generations=100):
+        # if no earlier training is done:
+        if not os.path.exists(self.experiment_name+'/evoman_solstate'):
+            population = self.initialize()
+            generation_number = 0
+        else:
+            print("Found earlier state")
+            self.env.load_state()
+            population = self.env.solutions[0]
+
+            # find generation we left off at:
+            with open(self.experiment_name + '/results.txt', 'r') as f:
+                for line in f:
+                    l = line
+                generation_number = int(l.strip().split()[1][:-1])
+
+        # Log mean, best, std
+        fitness_population = self.fitness_eval(population)
+        # mean = np.mean(fitness_population)
+        # best = np.argmax(fitness_population)
+        # std = np.std(fitness_population)
+
+        # Evolution loop
+        for gen_idx in tqdm.tqdm(range(generation_number, total_generations)):
+            # create offspring
+            offspring = self.crossover(population, p_mutation=0.2)# PLACEHOLDER
+            # mutated_offspring = [self.mutation(springie) for springie in offspring]
+
+            new_population = np.vstack((population, offspring))
+
+            # evaluate new population
+            new_fitness_population = self.fitness_eval(new_population)
+
+            # select population to continue to next generation
+            population, fitness_population = self.selection(new_population, new_fitness_population)
+
+            # save metrics for post-hoc evaluation
+            best = np.argmax(fitness_population)
+            mean = np.mean(fitness_population)
+            std  = np.std(fitness_population)
+
+            with open(self.experiment_name + '/results.txt', 'a') as f:
+                # save as best, mean, std
+                print(f"Generation {gen_idx}: {fitness_population[best]:.5f} {mean:.5f} {std:.5f}" )
+                f.write(f"Generation {gen_idx}: {fitness_population[best]:.5f} {mean:.5f} {std:.5f}\n")
+
+            np.savetxt(self.experiment_name + '/best.txt', population[best])
+            solutions = [population, fitness_population]
+            self.env.update_solutions(solutions)
+            self.env.save_state()
+
+if __name__ == '__main__':
+    Specialist(experiment_name = 'optimization_test_arco2', population_size=50).train(total_generations=100)
diff --git a/specialist_hannah.py b/specialist_hannah.py
new file mode 100644
index 0000000000..e20cc4e21e
--- /dev/null
+++ b/specialist_hannah.py
@@ -0,0 +1,103 @@
+# imports framework
+import sys, os
+
+from evoman.environment import Environment
+
+import numpy as np
+from demo_controller import player_controller
+
+# runs simulation
+def simulation(env,x):
+    # fitness, playerlife, enemylife, time
+    f,p,e,t = env.play(pcont=x)
+    return f
+
+# evaluation
+def evaluate(env, x):
+    return np.array(list(map(lambda y: simulation(env,y), x)))
+
+
+def main():
+    # choose this for not using visuals and thus making experiments faster
+    headless = True
+    if headless:
+        os.environ["SDL_VIDEODRIVER"] = "dummy"
+
+
+    experiment_name = 'optimization_test_hannah'
+    if not os.path.exists(experiment_name):
+        os.makedirs(experiment_name)
+
+    n_hidden_neurons = 10
+
+    # initializes simulation in individual evolution mode, for single static enemy.
+    env = Environment(experiment_name=experiment_name,
+                    enemies=[2],
+                    playermode="ai",
+                    player_controller=player_controller(n_hidden_neurons), # you  can insert your own controller here
+                    enemymode="static",
+                    level=2,
+                    speed="fastest",
+                    visuals=False)
+
+
+    # number of weights for multilayer with 10 hidden neurons
+    n_vars = (env.get_num_sensors()+1)*n_hidden_neurons + (n_hidden_neurons+1)*5
+    
+
+    # start writing your own code from here
+
+    # Create population
+        # Bounds for weights
+        # Store everything in two dimensional array n_population x n_weights
+        # Evaluate population (multiple of times) -> use network here
+
+    # Parent Selection
+        # Pairs of two so you can use for crossover
+        # Rank from 0 but some stochasticity
+
+    # Variation
+        # Crossover
+        # Mutation
+
+    # Survivor Selection
+        # Deterministically
+
+
+    # Full algorithm
+        # Parameters
+        # population_size = 50
+        # n_evaluations = 3
+        # n_offspring = 50
+        # weight_upper_bound = 2
+        # weight_lower_bound = -2
+        # mutation_sigma = .1
+        # generations = 10
+
+        # # Initialize environment, network and population. Perform an initial evaluation
+        # env = gym.make("MountainCar-v0")
+        # net = Perceptron (2, 3)
+        # pop = initialize_population(population_size, weight_lower_bound, weight_upper_bound)
+        # pop_fit = evaluate_population(pop, n_evaluations, net, env)
+
+        # for i in range (generations):
+        #     parents = parent_selection(pop, pop_fit, n_offspring)
+        #     offspring = crossover (parents)
+        #     offspring = mutate (offspring, weight_lower_bound, weight_upper_bound, mutation_sigma)
+
+        #     offspring_fit = evaluate_population(offspring, n_evaluations, net, env)
+
+        #     # concatenating to form a new population
+        #     pop = np.vstack((pop,offspring))
+        #     pop_fit = np.concatenate([pop_fit,offspring_fit])
+
+        #     pop, pop_fit = survivor_selection(pop, pop_fit, population_size)
+
+        #     print (f"Gen {i} - Best: {np.max (pop_fit)} - Mean: {np.mean(pop_fit)}")
+        #     clear_output(wait=True)
+        # env.close()
+
+
+
+if __name__ == '__main__':
+    main()