diff --git a/experiments/brain-speciation/MorphologyCompatibility.py b/experiments/brain-speciation/MorphologyCompatibility.py
new file mode 100644
index 0000000000..a2fff685ee
--- /dev/null
+++ b/experiments/brain-speciation/MorphologyCompatibility.py
@@ -0,0 +1,138 @@
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+if TYPE_CHECKING:
+    from typing import Optional
+    from pyrevolve.evolution.individual import Individual
+
+
+class MorphologyCompatibility:
+    def __init__(self,
+                 total_threshold: float = 1.0,
+                 branching_modules_count: float = 0.0,
+                 branching: float = 0.0,
+                 extremities: float = 0.0,
+                 limbs: float = 0.0,
+                 extensiveness: float = 0.0,
+                 length_of_limbs: float = 0.0,
+                 coverage: float = 0.0,
+                 joints: float = 0.0,
+                 active_hinges_count: float = 0.0,
+                 proportion: float = 0.0,
+                 width: float = 0.0,
+                 height: float = 0.0,
+                 z_depth: float = 0.0,
+                 absolute_size: float = 0.0,
+                 size: float = 0.0,
+                 sensors: float = 0.0,
+                 symmetry: float = 0.0,
+                 hinge_count: float = 0.0,
+                 brick_count: float = 0.0,
+                 brick_sensor_count: float = 0.0,
+                 touch_sensor_count: float = 0.0,
+                 free_slots: float = 0.0,
+                 height_base_ratio: float = 0.0,
+                 max_permitted_modules: Optional[int] = None,
+                 symmetry_vertical: float = 0.0,
+                 base_density: float = 0.0,
+                 bottom_layer: float = 0.0,
+                 ):
+        # Total threshold
+        self.total_threshold: float = total_threshold
+
+        # Absolute branching
+        self.branching_modules_count: float = branching_modules_count
+        # Relative branching
+        self.branching: float = branching
+        # Absolute number of limbs
+        self.extremities: float = extremities
+        # Relative number of limbs
+        self.limbs: float = limbs
+        # Absolute length of limbs
+        self.extensiveness: float = extensiveness
+        # Relative length of limbs
+        self.length_of_limbs: float = length_of_limbs
+        # Coverage
+        self.coverage: float = coverage
+        # Relative number of effective active joints
+        self.joints: float = joints
+        # Absolute number of effective active joints
+        self.active_hinges_count: float = active_hinges_count
+        # Proportion
+        self.proportion: float = proportion
+        # Width
+        self.width: float = width
+        # Height
+        self.height: float = height
+        # Z depth
+        self.z_depth: float = z_depth
+        # Absolute size
+        self.absolute_size: float = absolute_size
+        # Relative size in respect of the max body size `max_permitted_modules`
+        self.size: float = size
+        # Proportion of sensor vs empty slots
+        self.sensors: float = sensors
+        # Body symmetry in the xy plane
+        self.symmetry: float = symmetry
+        # Number of active joints
+        self.hinge_count: float = hinge_count
+        # Number of bricks
+        self.brick_count: float = brick_count
+        # Number of brick sensors
+        self.brick_sensor_count: float = brick_sensor_count
+        # Number of touch sensors
+        self.touch_sensor_count: float = touch_sensor_count
+        # Number of free slots
+        self.free_slots: float = free_slots
+        # Ratio of the height over the root of the area of the base
+        self.height_base_ratio: float = height_base_ratio
+        # Maximum number of modules allowed (sensors excluded)
+        self.max_permitted_modules: Optional[int] = max_permitted_modules
+        # Vertical symmetry
+        self.symmetry_vertical: float = symmetry_vertical
+        # Base model density
+        self.base_density: float = base_density
+        # Bottom layer of the robot
+        self.bottom_layer: float = bottom_layer
+
+    def compatible_individuals(self,
+                               individual1: Individual,
+                               individual2: Individual) -> bool:
+        morph_measure_1 = individual1.phenotype.measure_body()
+        morph_measure_2 = individual2.phenotype.measure_body()
+        _1 = morph_measure_1
+        _2 = morph_measure_2
+
+        # TODO consider normalization of some of these values, some are already normalized by definition
+
+        total_distance: float = 0.0
+        total_distance += self.branching_modules_count * abs(_2.branching_modules_count - _1.branching_modules_count)
+        total_distance += self.branching * abs(_2.branching - _1.branching)
+        total_distance += self.extremities * abs(_2.extremities - _1.extremities)
+        total_distance += self.limbs * abs(_2.limbs - _1.limbs)
+        total_distance += self.extensiveness * abs(_2.extensiveness - _1.extensiveness)
+        total_distance += self.length_of_limbs * abs(_2.length_of_limbs - _1.length_of_limbs)
+        total_distance += self.coverage * abs(_2.coverage - _1.coverage)
+        total_distance += self.joints * abs(_2.joints - _1.joints)
+        total_distance += self.active_hinges_count * abs(_2.active_hinges_count - _1.active_hinges_count)
+        total_distance += self.proportion * abs(_2.proportion - _1.proportion)
+        total_distance += self.width * abs(_2.width - _1.width)
+        total_distance += self.height * abs(_2.height - _1.height)
+        total_distance += self.z_depth * abs(_2.z_depth - _1.z_depth)
+        total_distance += self.absolute_size * abs(_2.absolute_size - _1.absolute_size)
+        if self.max_permitted_modules is not None:
+            total_distance += self.size * \
+                              abs(_2.absolute_size - _1.absolute_size) / self.max_permitted_modules
+        total_distance += self.sensors * abs(_2.sensors - _1.sensors)
+        total_distance += self.symmetry * abs(_2.symmetry - _1.symmetry)
+        total_distance += self.hinge_count * abs(_2.hinge_count - _1.hinge_count)
+        total_distance += self.brick_count * abs(_2.brick_count - _1.brick_count)
+        total_distance += self.brick_sensor_count * abs(_2.brick_sensor_count - _1.brick_sensor_count)
+        total_distance += self.touch_sensor_count * abs(_2.touch_sensor_count - _1.touch_sensor_count)
+        total_distance += self.free_slots * abs(_2.free_slots - _1.free_slots)
+        total_distance += self.height_base_ratio * abs(_2.height_base_ratio - _1.height_base_ratio)
+        total_distance += self.symmetry_vertical * abs(_2.symmetry_vertical - _1.symmetry_vertical)
+        total_distance += self.base_density * abs(_2.base_density - _1.base_density)
+        total_distance += self.bottom_layer * abs(_2.bottom_layer - _1.bottom_layer)
+
+        return total_distance <= self.total_threshold
diff --git a/experiments/brain-speciation/manager.py b/experiments/brain-speciation/manager.py
index b9c58e8ae7..83aeeff1e5 100644
--- a/experiments/brain-speciation/manager.py
+++ b/experiments/brain-speciation/manager.py
@@ -1,7 +1,9 @@
 #!/usr/bin/env python3
+from __future__ import annotations
+
 from pyrevolve import parser
 from pyrevolve.evolution import fitness
-from pyrevolve.evolution.selection import multiple_selection, tournament_selection
+from pyrevolve.evolution.selection import multiple_selection_with_duplicates, tournament_selection
 from pyrevolve.evolution.speciation.population_speciated import PopulationSpeciated
 from pyrevolve.evolution.speciation.population_speciated_config import PopulationSpeciatedConfig
 from pyrevolve.evolution.speciation.population_speciated_management import steady_state_speciated_population_management
@@ -11,13 +13,17 @@
 from pyrevolve.genotype.lsystem_neat.mutation import LSystemNeatMutationConf as lMutationConfig
 from pyrevolve.genotype.plasticoding.mutation.mutation import MutationConfig as plasticMutationConfig
 from pyrevolve.genotype.lsystem_neat.mutation import standard_mutation as lmutation
-
 from pyrevolve.util.supervisor.analyzer_queue import AnalyzerQueue
 from pyrevolve.util.supervisor.simulator_queue import SimulatorQueue
 from pyrevolve.custom_logging.logger import logger
 from pyrevolve.genotype.plasticoding import PlasticodingConfig
 from pyrevolve.genotype.lsystem_neat.lsystem_neat_genotype import LSystemCPGHyperNEATGenotype, LSystemCPGHyperNEATGenotypeConfig
 from pyrevolve.genotype.neat_brain_genome.neat_brain_genome import NeatBrainGenomeConfig
+from .MorphologyCompatibility import MorphologyCompatibility
+
+from typing import TYPE_CHECKING
+if TYPE_CHECKING:
+    from pyrevolve.evolution.individual import Individual
 
 
 async def run():
@@ -31,7 +37,7 @@ async def run():
     offspring_size = 50
 
     body_conf = PlasticodingConfig(
-        max_structural_modules=20,  # TODO increase
+        max_structural_modules=20,
         allow_vertical_brick=False,
         use_movement_commands=True,
         use_rotation_commands=False,
@@ -67,17 +73,32 @@ async def run():
     crossover_conf = lCrossoverConfig(
         crossover_prob=0.8,
     )
+
+    compatibitity_tester = MorphologyCompatibility(
+        total_threshold=1.0,
+        size=1.0,
+        brick_count=1.0,
+        proportion=1.0,
+        coverage=1.0,
+        joints=1.5,
+        branching=1.0,
+        symmetry=0.0,
+        max_permitted_modules=body_conf.max_structural_modules,
+    )
+
     # experiment params #
 
     # Parse command line / file input arguments
     args = parser.parse_args()
     experiment_management = ExperimentManagement(args)
+    has_offspring = False
     do_recovery = args.recovery_enabled and not experiment_management.experiment_is_new()
 
     logger.info(f'Activated run {args.run} of experiment {args.experiment_name}')
 
     if do_recovery:
-        gen_num, has_offspring, next_robot_id, next_species_id = experiment_management.read_recovery_state(population_size, offspring_size)
+        gen_num, has_offspring, next_robot_id, next_species_id = \
+            experiment_management.read_recovery_state(population_size, offspring_size, species=True)
 
         if gen_num == num_generations-1:
             logger.info('Experiment is already complete.')
@@ -97,9 +118,15 @@ async def run():
     if next_species_id < 0:
         next_species_id = 1
 
-    def are_genomes_compatible_fn(genotype1: LSystemCPGHyperNEATGenotype,
-                                  genotype2: LSystemCPGHyperNEATGenotype) -> bool:
-        return genotype1.is_brain_compatible(genotype2, genotype_conf)
+    def are_individuals_brains_compatible_fn(individual1: Individual,
+                                             individual2: Individual) -> bool:
+        assert isinstance(individual1.genotype, LSystemCPGHyperNEATGenotype)
+        assert isinstance(individual2.genotype, LSystemCPGHyperNEATGenotype)
+        return individual1.genotype.is_brain_compatible(individual2.genotype, genotype_conf)
+
+    def are_individuals_morphologies_compatible_fn(individual1: Individual,
+                                                   individual2: Individual) -> bool:
+        return compatibitity_tester.compatible_individuals(individual1, individual2)
 
     population_conf = PopulationSpeciatedConfig(
         population_size=population_size,
@@ -111,7 +138,7 @@ def are_genomes_compatible_fn(genotype1: LSystemCPGHyperNEATGenotype,
         crossover_operator=lcrossover,
         crossover_conf=crossover_conf,
         selection=lambda individuals: tournament_selection(individuals, 2),
-        parent_selection=lambda individuals: multiple_selection(individuals, 2, tournament_selection),
+        parent_selection=lambda individuals: multiple_selection_with_duplicates(individuals, 2, tournament_selection),
         population_management=steady_state_speciated_population_management,
         population_management_selector=tournament_selection,
         evaluation_time=args.evaluation_time,
@@ -119,10 +146,11 @@ def are_genomes_compatible_fn(genotype1: LSystemCPGHyperNEATGenotype,
         experiment_name=args.experiment_name,
         experiment_management=experiment_management,
         # species stuff
-        are_genomes_compatible_fn=are_genomes_compatible_fn,
+        # are_individuals_compatible_fn=are_individuals_brains_compatible_fn,
+        are_individuals_compatible_fn=are_individuals_morphologies_compatible_fn,
         young_age_threshold=5,
         young_age_fitness_boost=2.0,
-        old_age_threshold=20,
+        old_age_threshold=35,
         old_age_fitness_penalty=0.5,
         species_max_stagnation=30,
     )
@@ -142,12 +170,16 @@ def are_genomes_compatible_fn(genotype1: LSystemCPGHyperNEATGenotype,
                                      next_species_id)
 
     if do_recovery:
-        raise NotImplementedError('recovery not implemented')
         # loading a previous state of the experiment
         population.load_snapshot(gen_num)
         if gen_num >= 0:
             logger.info(f'Recovered snapshot {gen_num}, pop with {len(population.genus)} individuals')
 
+        # TODO partial recovery is not implemented, this is a substitute
+        has_offspring = False
+        next_robot_id = 1 + population.config.population_size + gen_num * population.config.offspring_size
+        population.next_robot_id = next_robot_id
+
         if has_offspring:
             raise NotImplementedError('partial recovery not implemented')
             recovered_individuals = population.load_partially_completed_generation(gen_num, population_size, offspring_size, next_robot_id)
diff --git a/pyrevolve/evolution/individual.py b/pyrevolve/evolution/individual.py
index 37ca523e54..d055cf34b1 100644
--- a/pyrevolve/evolution/individual.py
+++ b/pyrevolve/evolution/individual.py
@@ -6,6 +6,7 @@
     from typing import Optional, List
     from pyrevolve.revolve_bot import RevolveBot
     from pyrevolve.genotype import Genotype
+    from pyrevolve.evolution.speciation.species import Species
 
 
 class Individual:
@@ -38,15 +39,31 @@ def id(self) -> int:
             _id = self.genotype.id
         return _id
 
-    def export_genotype(self, folder) -> None:
+    def export_genotype(self, folder: str) -> None:
         self.genotype.export_genotype(os.path.join(folder, f'genotype_{self.phenotype.id}.txt'))
 
-    def export_phenotype(self, folder) -> None:
+    def export_phenotype(self, folder: str) -> None:
         if self.phenotype is None:
             self.develop()
         self.phenotype.save_file(os.path.join(folder, f'phenotype_{self.phenotype.id}.yaml'), conf_type='yaml')
 
-    def export_fitness(self, folder) -> None:
+    def export_phylogenetic_info(self, folder: str) -> None:
+        """
+        Export phylogenetic information
+        (parents and other possibly other information to build a phylogenetic tree)
+        :param folder: folder where to save the information
+        """
+        if self.parents is not None:
+            parents_ids: List[str] = [str(p.id) for p in self.parents]
+            parents_ids_str = ",".join(parents_ids)
+        else:
+            parents_ids_str = 'None'
+
+        filename = os.path.join(folder, f'parents_{self.id}.yaml')
+        with open(filename, 'w') as file:
+            file.write(f'parents:{parents_ids_str}')
+
+    def export_fitness_single_file(self, folder: str) -> None:
         """
         It's saving the fitness into a file. The fitness can be a floating point number or None
         :param folder: folder where to save the fitness
@@ -54,8 +71,9 @@ def export_fitness(self, folder) -> None:
         with open(os.path.join(folder, f'fitness_{self.id}.txt'), 'w') as f:
             f.write(str(self.fitness))
 
-    def export(self, folder) -> None:
+    def export(self, folder: str) -> None:
         self.export_genotype(folder)
+        self.export_phylogenetic_info(folder)
         self.export_phenotype(folder)
         self.export_fitness(folder)
 
diff --git a/pyrevolve/evolution/population/population.py b/pyrevolve/evolution/population/population.py
index 925fcccf6b..ad8d232225 100644
--- a/pyrevolve/evolution/population/population.py
+++ b/pyrevolve/evolution/population/population.py
@@ -3,13 +3,13 @@
 import os
 
 from pyrevolve.evolution.individual import Individual
-from pyrevolve.tol.manage import measures
 from pyrevolve.custom_logging.logger import logger
 from pyrevolve.evolution.population.population_config import PopulationConfig
 
 from typing import TYPE_CHECKING
 if TYPE_CHECKING:
-    from typing import List, Optional, AnyStr
+    from typing import List, Optional
+    from pyrevolve.evolution.speciation.species import Species
     from pyrevolve.tol.manage.measures import BehaviouralMeasurements
     from pyrevolve.util.supervisor.analyzer_queue import AnalyzerQueue, SimulatorQueue
 
@@ -44,10 +44,15 @@ def __init__(self,
         self.simulator_queue = simulator_queue
         self.next_robot_id = next_robot_id
 
-    def _new_individual(self, genotype):
+    def _new_individual(self,
+                        genotype,
+                        parents: Optional[List[Individual]] = None):
         individual = Individual(genotype)
         individual.develop()
         individual.phenotype.update_substrate()
+        if parents is not None:
+            individual.parents = parents
+
         self.config.experiment_management.export_genotype(individual)
         self.config.experiment_management.export_phenotype(individual)
         self.config.experiment_management.export_phenotype_images(individual)
diff --git a/pyrevolve/evolution/selection.py b/pyrevolve/evolution/selection.py
index 6cf888d6c0..6316fa1e37 100644
--- a/pyrevolve/evolution/selection.py
+++ b/pyrevolve/evolution/selection.py
@@ -32,9 +32,12 @@ def tournament_selection(population: List[Individual], k=2) -> Individual:
 
 def multiple_selection(population: List[Individual],
                        selection_size: int,
-                       selection_function: Callable[[List[Individual]], Individual]) -> List[Individual]:
+                       selection_function: Callable[[List[Individual]], Individual]
+                       ) -> List[Individual]:
     """
-    Perform selection on population of distinct group, can be used in the form parent selection or survival selection
+    Perform selection on population of distinct group, it can be used in the
+    form parent selection or survival selection.
+    It never selects the same individual more than once
     :param population: list of individuals where to select from
     :param selection_size: amount of individuals to select
     :param selection_function:
@@ -49,3 +52,22 @@ def multiple_selection(population: List[Individual],
                 selected_individuals.append(selected_individual)
                 new_individual = True
     return selected_individuals
+
+
+def multiple_selection_with_duplicates(population: List[Individual],
+                                       selection_size: int,
+                                       selection_function: Callable[[List[Individual]], Individual]
+                                       ) -> List[Individual]:
+    """
+    Perform selection on population of distinct group, it can be used in the
+    form parent selection or survival selection.
+    It can select the same individual more than once
+    :param population: list of individuals where to select from
+    :param selection_size: amount of individuals to select
+    :param selection_function:
+    """
+    selected_individuals = []
+    for _ in range(selection_size):
+        selected_individual = selection_function(population)
+        selected_individuals.append(selected_individual)
+    return selected_individuals
diff --git a/pyrevolve/evolution/speciation/genus.py b/pyrevolve/evolution/speciation/genus.py
index 73f49e4fef..fc62d6b0d6 100644
--- a/pyrevolve/evolution/speciation/genus.py
+++ b/pyrevolve/evolution/speciation/genus.py
@@ -127,18 +127,23 @@ async def next_generation(self,
         ##################################################################
         # MANAGE ORPHANS, POSSIBLY CREATE NEW SPECIES
         # recheck if other species can adopt the orphan individuals.
+        list_of_new_species = []
         for orphan in orphans:
             for species in new_species_collection:
                 if species.is_compatible(orphan, self.config):
                     species.append(orphan)
                     break
             else:
-                new_species_collection.add_species(Species([orphan], self._next_species_id))
+                new_species = Species([orphan], self._next_species_id)
+                new_species_collection.add_species(new_species)
+                list_of_new_species.append(new_species)
                 # add an entry for new species which does not have a previous iteration.
                 self._next_species_id += 1
 
         # Do a recount on the number of offspring per species.
-        offspring_amounts = self._count_offsprings(self.config.population_size - len(orphans))
+        new_species_size = sum(map(lambda species: len(species), list_of_new_species))
+        offspring_amounts = self._count_offsprings(self.config.population_size - new_species_size)
+        assert sum(offspring_amounts) == self.config.population_size - new_species_size
 
         ##################################################################
         # EVALUATE NEW INDIVIDUALS
@@ -248,10 +253,10 @@ def _correct_population_size(self, species_offspring_amount: List[int], missing_
         elif missing_offspring < 0:  # negative have excess individuals
             # remove missing number of individuals
             excess_offspring = -missing_offspring
-            excluded_list = set()
+            excluded_id_list = set()
 
             while excess_offspring > 0:
-                worst_species_index, species = self.species_collection.get_worst(1, excluded_list)
+                worst_species_index, species = self.species_collection.get_worst(1, excluded_id_list)
                 current_amount = species_offspring_amount[worst_species_index]
 
                 if current_amount > excess_offspring:
@@ -262,7 +267,7 @@ def _correct_population_size(self, species_offspring_amount: List[int], missing_
                     current_amount = 0
 
                 species_offspring_amount[worst_species_index] = current_amount
-                excluded_list.add(species)
+                excluded_id_list.add(species.id)
 
             assert excess_offspring == 0
 
diff --git a/pyrevolve/evolution/speciation/population_speciated.py b/pyrevolve/evolution/speciation/population_speciated.py
index 36d8204396..cfafd8ea85 100644
--- a/pyrevolve/evolution/speciation/population_speciated.py
+++ b/pyrevolve/evolution/speciation/population_speciated.py
@@ -34,6 +34,7 @@ async def initialize(self, recovered_individuals: Optional[List[Individual]] = N
         """
         Populates the population (individuals list) with Individual objects that contains their respective genotype.
         """
+        assert recovered_individuals is None
         individuals = []
 
         recovered_individuals = [] if recovered_individuals is None else recovered_individuals
@@ -91,6 +92,7 @@ def _generate_individual(self, individuals: List[Individual]) -> Individual:
             child = Individual(child_genotype)
         else:
             child = self.config.selection(individuals)
+            parents = [child]
 
         child.genotype.id = self.next_robot_id
         self.next_robot_id += 1
@@ -99,7 +101,7 @@ def _generate_individual(self, individuals: List[Individual]) -> Individual:
         child_genotype = self.config.mutation_operator(child.genotype, self.config.mutation_conf)
 
         # Create new individual
-        return self._new_individual(child_genotype)
+        return self._new_individual(child_genotype, parents)
 
     def load_snapshot(self, gen_num: int) -> None:
         """
@@ -121,3 +123,8 @@ def load_individual_fn(_id: int) -> Individual:
             if file.is_file() and file.name.endswith('.yaml'):
                 species = Species.Deserialize(file.path, loaded_individuals, load_individual_fn)
                 self.genus.species_collection.add_species(species)
+
+        n_loaded_individuals = count_individuals(self.genus.species_collection)
+        if n_loaded_individuals != self.config.population_size:
+            raise RuntimeError(f'The loaded population ({n_loaded_individuals}) '
+                               f'does not match the population size ({self.config.population_size})')
diff --git a/pyrevolve/evolution/speciation/population_speciated_config.py b/pyrevolve/evolution/speciation/population_speciated_config.py
index 257212afff..6c47105a4b 100644
--- a/pyrevolve/evolution/speciation/population_speciated_config.py
+++ b/pyrevolve/evolution/speciation/population_speciated_config.py
@@ -32,7 +32,7 @@ def __init__(self,
                  evaluation_time: float,
                  experiment_name: str,
                  experiment_management,
-                 are_genomes_compatible_fn: Callable[[Genotype, Genotype], bool],
+                 are_individuals_compatible_fn: Callable[[Individual, Individual], bool],
                  young_age_threshold: int = 5,
                  young_age_fitness_boost: float = 1.1,
                  old_age_threshold: int = 30,
@@ -60,7 +60,7 @@ def __init__(self,
         :param evaluation_time: duration of an experiment
         :param experiment_name: name for the folder of the current experiment
         :param experiment_management: object with methods for managing the current experiment
-        :param are_genomes_compatible_fn: function that determines if two genomes are compatible
+        :param are_individuals_compatible_fn: function that determines if two individuals are compatible
         :param young_age_threshold: define until what age (excluded) species are still young
             and need to be protected (with a fitness boost)
         :param young_age_fitness_boost: Fitness multiplier for young species.
@@ -88,7 +88,7 @@ def __init__(self,
                          experiment_name,
                          experiment_management,
                          offspring_size)
-        self.are_genomes_compatible = are_genomes_compatible_fn  # type: Callable[[Genotype, Genotype], bool]
+        self.are_individuals_compatible = are_individuals_compatible_fn  # type: Callable[[Individual, Individual], bool]
         self.young_age_threshold = young_age_threshold
         self.young_age_fitness_boost = young_age_fitness_boost
         self.old_age_threshold = old_age_threshold
diff --git a/pyrevolve/evolution/speciation/population_speciated_management.py b/pyrevolve/evolution/speciation/population_speciated_management.py
index 74da86402d..a440506f89 100644
--- a/pyrevolve/evolution/speciation/population_speciated_management.py
+++ b/pyrevolve/evolution/speciation/population_speciated_management.py
@@ -1,4 +1,4 @@
-from pyrevolve.evolution.selection import multiple_selection
+from pyrevolve.evolution.selection import multiple_selection, multiple_selection_with_duplicates
 
 
 def steady_state_speciated_population_management(old_individuals, new_individuals, number_of_individuals, selector):
@@ -7,7 +7,7 @@ def steady_state_speciated_population_management(old_individuals, new_individual
     # TODO old function: need parameter for ...
     selection_pool = old_individuals + new_individuals
 
-    return multiple_selection(selection_pool, number_of_individuals, selector)
+    return multiple_selection_with_duplicates(selection_pool, number_of_individuals, selector)
 
 
 def generational_population_speciated_management(old_individuals, new_individuals, number_of_individuals, selector):
diff --git a/pyrevolve/evolution/speciation/species.py b/pyrevolve/evolution/speciation/species.py
index de749df7d5..3199a52c76 100644
--- a/pyrevolve/evolution/speciation/species.py
+++ b/pyrevolve/evolution/speciation/species.py
@@ -54,7 +54,9 @@ def is_compatible(self, candidate: Individual, population_config: PopulationSpec
         :param population_config: config where to pick the `are genomes compatible` function
         :return: if the candidate individual is compatible or not
         """
-        return population_config.are_genomes_compatible(candidate.genotype, self._representative().genotype)
+        if self.empty():
+            return False
+        return population_config.are_individuals_compatible(candidate, self._representative())
 
     # TODO duplicate code with species collection best/worst function
     def get_best_fitness(self) -> float:
diff --git a/pyrevolve/evolution/speciation/species_collection.py b/pyrevolve/evolution/speciation/species_collection.py
index 5434727c2b..75c8dfc538 100644
--- a/pyrevolve/evolution/speciation/species_collection.py
+++ b/pyrevolve/evolution/speciation/species_collection.py
@@ -53,9 +53,9 @@ def _update_cache(self) -> None:
         assert len(self._collection) > 0
 
         # BEST
-        index = int(numpy.argmax(
-                [species.get_best_individual().fitness for species in self._collection]
-        ))
+        species_best_fitness = [species.get_best_individual().fitness for species in self._collection]
+        species_best_fitness = map(lambda f: -math.inf if f is None else f, species_best_fitness)
+        index = int(numpy.argmax(species_best_fitness))
         self._best = (index, self._collection[index])
 
         # cannot calculate WORST cache, because
@@ -76,18 +76,18 @@ def get_best(self) -> (int, Species):
 
     def get_worst(self,
                   minimal_size: int,
-                  exclude_list: Optional[Set[Species]] = None) -> (int, Species):
+                  exclude_id_list: Optional[Set[int]] = None) -> (int, Species):
         """
         Finds the worst species (based on the best fitness of that species)
         Crashes if there are no species with at least `minimal_size` individuals
 
         :param minimal_size: Species with less individuals than this will not be considered
-        :param exclude_list: Species in this list will be ignored
+        :param exclude_id_list: Species in this list will be ignored
         :return: the index and a reference to the worst species
         """
         assert len(self._collection) > 0
 
-        worst_species_index, worst_species = self._calculate_worst_fitness(minimal_size, exclude_list)
+        worst_species_index, worst_species = self._calculate_worst_fitness(minimal_size, exclude_id_list)
 
         assert worst_species_index != -1
         assert worst_species is not None
@@ -96,15 +96,15 @@ def get_worst(self,
 
     def _calculate_worst_fitness(self,
                                  minimal_size: int,
-                                 exclude_list: Optional[Set[Species]]) -> (int, Species):
+                                 exclude_id_list: Optional[Set[int]]) -> (int, Species):
         worst_species_index = -1
         worst_species_fitness = math.inf
         worst_species = None
 
         for i, species in enumerate(self._collection):
 
-            if exclude_list is not None \
-                    and species in exclude_list:
+            if exclude_id_list is not None \
+                    and species.id in exclude_id_list:
                 continue
 
             if len(species) < minimal_size:
@@ -113,6 +113,7 @@ def _calculate_worst_fitness(self,
                 # species_fitness = -math.inf
 
             species_fitness = species.get_best_fitness()
+            species_fitness = -math.inf if species_fitness is None else species_fitness
             if species_fitness < worst_species_fitness:
                 worst_species_fitness = species_fitness
                 worst_species = species
@@ -135,12 +136,12 @@ def cleanup(self) -> None:
         """
         Remove all empty species (cleanup routine for every case..)
         """
-        new_species = []
+        new_collection = []
         for species in self._collection:
             if not species.empty():
-                new_species.append(species)
+                new_collection.append(species)
 
-        self._collection = new_species
+        self._collection = new_collection
 
     def clear(self) -> None:
         self._collection.clear()
@@ -173,10 +174,10 @@ def update(self) -> None:
             old_best_species.age.reset_generations()
 
 
-def count_individuals(species_collection: Optional[SpeciesCollection] = None) -> int:
+def count_individuals(species_collection: SpeciesCollection) -> int:
     """
     Counts the number of individuals in the species_list.
-    :param species_list: if None, it will use self.species_list
+    :param species_collection: collection of species
     :return: the total number of individuals
     """
     # count the total number of individuals inside every species in the species_list
diff --git a/pyrevolve/experiment_management.py b/pyrevolve/experiment_management.py
index b673d3a6cd..8f9f84ce41 100644
--- a/pyrevolve/experiment_management.py
+++ b/pyrevolve/experiment_management.py
@@ -13,6 +13,7 @@
     from typing import List, AnyStr, Optional
     from pyrevolve.tol.manage.measures import BehaviouralMeasurements
     from pyrevolve.evolution.speciation.genus import Genus
+    from pyrevolve.evolution.speciation.species import Species
     from pyrevolve.evolution.population.population_config import PopulationConfig
 
 
@@ -28,6 +29,7 @@ def __init__(self, settings):
         self._experiment_folder: str = os.path.join(manager_folder, self.EXPERIMENT_FOLDER, self.settings.experiment_name, self.settings.run)
         self._data_folder: str = os.path.join(self._experiment_folder, self.DATA_FOLDER)
         self._genotype_folder: str = os.path.join(self.data_folder, 'genotypes')
+        self._phylogenetic_folder: str = os.path.join(self.data_folder, 'phylogeny')
         self._phenotype_folder: str = os.path.join(self.data_folder, 'phenotypes')
         self._phenotype_images_folder: str = os.path.join(self._phenotype_folder, 'images')
         self._fitness_file_path: str = os.path.join(self.data_folder, 'fitness.csv')
@@ -48,6 +50,7 @@ def create_exp_folders(self) -> None:
         os.makedirs(self.experiment_folder)
         os.mkdir(self.data_folder)
         os.mkdir(self._genotype_folder)
+        os.mkdir(self._phylogenetic_folder)
         os.mkdir(self._phenotype_folder)
         os.mkdir(self._descriptor_folder)
         os.mkdir(self._behavioural_desc_folder)
@@ -83,6 +86,7 @@ def export_genotype(self, individual: Individual) -> None:
         """
         if self.settings.recovery_enabled:
             individual.export_genotype(self._genotype_folder)
+            individual.export_phylogenetic_info(self._phylogenetic_folder)
 
     def export_phenotype(self, individual: Individual) -> None:
         """
@@ -152,12 +156,35 @@ def export_behavior_measures(self, _id: int, measures: BehaviouralMeasurements)
                 for key, val in measures.items():
                     f.write(f'{key} {val}\n')
 
-    def export_phenotype_images(self, individual: Individual, dirpath: Optional[str] = None) -> None:
+    def export_phenotype_images(self, individual: Individual, dirpath: Optional[str] = None, rename_if_present=False) -> None:
         dirpath = dirpath if dirpath is not None \
             else self._phenotype_images_folder
         try:
-            individual.phenotype.render_body(os.path.join(dirpath, f'body_{individual.phenotype.id}.png'))
-            individual.phenotype.render_brain(os.path.join(dirpath, f'brain_{individual.phenotype.id}.png'))
+            # -- Body image --
+            body_filename_part = os.path.join(dirpath, f'body_{individual.phenotype.id}')
+            if rename_if_present and os.path.exists(f'{body_filename_part}.png'):
+                counter = 1
+                while os.path.exists(f'{body_filename_part}_{counter}.png'):
+                    counter += 1
+                os.symlink(f'body_{individual.phenotype.id}.png',
+                           f'{body_filename_part}_{counter}.png',
+                           target_is_directory=False)
+            else:
+                # Write file
+                individual.phenotype.render_body(f'{body_filename_part}.png')
+
+            # -- Brain image --
+            brain_filename_part = os.path.join(dirpath, f'brain_{individual.phenotype.id}')
+            if rename_if_present and os.path.exists(f'{brain_filename_part}.png'):
+                counter = 1
+                while os.path.exists(f'{brain_filename_part}_{counter}.png'):
+                    counter += 1
+                os.symlink(f'brain_{individual.phenotype.id}.png',
+                           f'{brain_filename_part}_{counter}.png',
+                           target_is_directory=False)
+            else:
+                # Write file
+                individual.phenotype.render_brain(os.path.join(dirpath, f'{brain_filename_part}.png'))
         except Exception as e:
             logger.warning(f'Error rendering phenotype images: {e}')
 
@@ -193,7 +220,7 @@ def export_snapshots_species(self, genus: Genus, gen_num: int) -> None:
                 os.mkdir(species_folder)
                 species.serialize(species_on_disk)
                 for individual, _ in species.iter_individuals():
-                    self.export_phenotype_images(individual, species_folder)
+                    self.export_phenotype_images(individual, species_folder, rename_if_present=True)
 
     def experiment_is_new(self) -> bool:
         """
@@ -253,7 +280,7 @@ def read_recovery_state(self, population_size: int, offspring_size: int, species
                         if not species_on_disk.is_file():
                             continue
                         with open(species_on_disk.path) as file:
-                            species = yaml.load(file)
+                            species = yaml.load(file, Loader=yaml.SafeLoader)
                             n_exported_genomes += len(species['individuals_ids'])
                             species_id = species['id']
 
diff --git a/pyrevolve/genotype/lsystem_neat/lsystem_neat_genotype.py b/pyrevolve/genotype/lsystem_neat/lsystem_neat_genotype.py
index 778c95c1e1..85f4a6dde7 100644
--- a/pyrevolve/genotype/lsystem_neat/lsystem_neat_genotype.py
+++ b/pyrevolve/genotype/lsystem_neat/lsystem_neat_genotype.py
@@ -27,7 +27,7 @@ def __init__(self, conf: Optional[LSystemCPGHyperNEATGenotypeConfig] = None, rob
             self._brain_genome = None
         else:
             assert robot_id is not None
-            self._body_genome = random_initialization(conf.plasticoding, robot_id)
+            self._body_genome: Plasticoding = random_initialization(conf.plasticoding, robot_id)
             self._brain_genome = NeatBrainGenome(conf.neat, robot_id)
 
     @property
diff --git a/pyrevolve/genotype/plasticoding/alphabet.py b/pyrevolve/genotype/plasticoding/alphabet.py
index 0bc2ded3fc..17e4e43d64 100644
--- a/pyrevolve/genotype/plasticoding/alphabet.py
+++ b/pyrevolve/genotype/plasticoding/alphabet.py
@@ -1,7 +1,12 @@
+from __future__ import annotations
 from collections.abc import Iterable
 from enum import Enum
 from pyrevolve.custom_logging.logger import logger
 
+from typing import TYPE_CHECKING
+if TYPE_CHECKING:
+    from typing import List, Union
+
 INDEX_SYMBOL = 0
 INDEX_PARAMS = 1
 
@@ -43,7 +48,7 @@ class Alphabet(Enum):
     MOVE_REF_O = 'brainmoveTTS'
 
     @staticmethod
-    def modules(allow_vertical_brick: bool):
+    def modules(allow_vertical_brick: bool) -> List[Alphabet]:
         # this function MUST return the core always as the first element
         modules = [
             Alphabet.CORE_COMPONENT,
@@ -56,16 +61,16 @@ def modules(allow_vertical_brick: bool):
             modules.append(Alphabet.BLOCK_VERTICAL)
         return modules
 
-    def is_vertical_module(self):
+    def is_vertical_module(self) -> bool:
         return self is Alphabet.JOINT_VERTICAL \
                or self is Alphabet.BLOCK_VERTICAL
 
-    def is_joint(self):
+    def is_joint(self) -> bool:
         return self is Alphabet.JOINT_VERTICAL \
                or self is Alphabet.JOINT_HORIZONTAL
 
     @staticmethod
-    def morphology_mounting_commands():
+    def morphology_mounting_commands() -> List[Alphabet]:
         return [
             Alphabet.ADD_RIGHT,
             Alphabet.ADD_FRONT,
@@ -73,7 +78,9 @@ def morphology_mounting_commands():
         ]
 
     @staticmethod
-    def morphology_moving_commands(use_movement_commands: bool, use_rotation_commands: bool, use_movement_stack: bool):
+    def morphology_moving_commands(use_movement_commands: bool,
+                                   use_rotation_commands: bool,
+                                   use_movement_stack: bool) -> List[Alphabet]:
         commands = []
         if use_movement_commands:
             commands.append(Alphabet.MOVE_RIGHT)
@@ -93,7 +100,7 @@ def morphology_moving_commands(use_movement_commands: bool, use_rotation_command
         return commands
 
     @staticmethod
-    def controller_changing_commands():
+    def controller_changing_commands() -> List[Alphabet]:
         return [
             Alphabet.ADD_EDGE,
             Alphabet.MUTATE_EDGE,
@@ -104,15 +111,18 @@ def controller_changing_commands():
         ]
 
     @staticmethod
-    def controller_moving_commands():
+    def controller_moving_commands() -> List[Alphabet]:
         return [
             Alphabet.MOVE_REF_S,
             Alphabet.MOVE_REF_O,
         ]
 
     @staticmethod
-    def wordify(letters):
+    def wordify(letters: Union[Alphabet, List]
+                ) -> Union[(Alphabet, List), List[(Alphabet, List)]]:
         if isinstance(letters, Iterable):
             return [(a, []) for a in letters]
-        else:
+        elif isinstance(letters, Alphabet):
             return (letters, [])
+        else:
+            raise RuntimeError(f'Cannot wordify element of type {type(letters)}')
diff --git a/pyrevolve/genotype/plasticoding/initialization.py b/pyrevolve/genotype/plasticoding/initialization.py
index 7dec4e02f8..df6990d785 100644
--- a/pyrevolve/genotype/plasticoding/initialization.py
+++ b/pyrevolve/genotype/plasticoding/initialization.py
@@ -1,15 +1,19 @@
+from __future__ import annotations
+import random
 from pyrevolve.genotype.plasticoding.plasticoding import Alphabet
 from pyrevolve.genotype.plasticoding.plasticoding import Plasticoding
-import random
+
+from typing import TYPE_CHECKING
+if TYPE_CHECKING:
+    from typing import Dict, List
+    from pyrevolve.genotype.plasticoding.plasticoding import PlasticodingConfig
 
 
-def _generate_random_grammar(conf):
+def _generate_random_grammar(conf: PlasticodingConfig) -> Dict[Alphabet, List]:
     """
     Initializing a new genotype,
     :param conf: e_max_groups, maximum number of groups of symbols
-    :type conf: PlasticodingConfig
     :return: a random new Genome
-    :rtype: dictionary
     """
     s_segments = random.randint(1, conf.e_max_groups)
     grammar = {}
@@ -52,13 +56,13 @@ def _generate_random_grammar(conf):
     return grammar
 
 
-def random_initialization(conf, next_robot_id):
+def random_initialization(conf: PlasticodingConfig, _id: int) -> Plasticoding:
     """
     Initializing a random genotype.
-    :type conf: PlasticodingConfig
-    :return: a Genome
-    :rtype: Plasticoding
+    :param conf: Plasticoding genotype configuration
+    :param _id: id of the newly created genotype
+    :return: a Plasticoding Genome
     """
-    genotype = Plasticoding(conf, next_robot_id)
+    genotype = Plasticoding(conf, _id)
     genotype.grammar = _generate_random_grammar(conf)
-    return genotype
\ No newline at end of file
+    return genotype
diff --git a/pyrevolve/genotype/plasticoding/plasticoding.py b/pyrevolve/genotype/plasticoding/plasticoding.py
index 1fc0cb498c..61fe801dca 100644
--- a/pyrevolve/genotype/plasticoding/plasticoding.py
+++ b/pyrevolve/genotype/plasticoding/plasticoding.py
@@ -1,6 +1,7 @@
 from __future__ import annotations
 import random
 import math
+import copy
 
 from pyrevolve.custom_logging.logger import logger
 from pyrevolve.genotype import Genotype
@@ -10,8 +11,9 @@
 
 from typing import TYPE_CHECKING
 if TYPE_CHECKING:
-    from typing import Optional
+    from typing import Optional, Dict, List, Union
     from pyrevolve.genotype.plasticoding import PlasticodingConfig
+    from pyrevolve.revolve_bot import RevolveBot
 
 
 class Plasticoding(Genotype):
@@ -25,22 +27,34 @@ def __init__(self, conf: PlasticodingConfig, robot_id: Optional[int]):
         :param robot_id: unique id of the robot
         :type conf: PlasticodingConfig
         """
-        self.conf = conf
+        self.conf: PlasticodingConfig = conf
         assert robot_id is None or str(robot_id).isdigit()
-        self.id = int(robot_id) if robot_id is not None else -1
-        self.grammar = {}
+        self.id: int = int(robot_id) if robot_id is not None else -1
+        self.grammar: Dict[Alphabet, List] = {}
 
         # Auxiliary variables
-        self.valid = False
+        self.valid: bool = False
         self.intermediate_phenotype = None
-        self.phenotype = None
-
-    def load_genotype(self, genotype_file):
-        with open(genotype_file) as f:
+        self.phenotype: Optional[RevolveBot] = None
+
+    def clone(self):
+        # Cannot use deep clone for this genome, because it is bugged sometimes
+        _id = self.id if self.id >= 0 else None
+        other = Plasticoding(self.conf, _id)
+        other.grammar = {}
+        for key, value in self.grammar.items():
+            other.grammar[key] = copy.deepcopy(value)
+        other.valid = self.valid
+        # other.intermediate_phenotype = self.intermediate_phenotype
+        # other.phenotype = self.phenotype
+        return other
+
+    def load_genotype(self, genotype_filename: str) -> None:
+        with open(genotype_filename) as f:
             lines = f.readlines()
             self._load_genotype_from(lines)
 
-    def _load_genotype_from(self, lines):
+    def _load_genotype_from(self, lines: List[str]) -> None:
         for line in lines:
             line_array = line.split(' ')
             replaceable_symbol = Alphabet(line_array[0])
@@ -55,11 +69,11 @@ def _load_genotype_from(self, lines):
                     params = []
                 self.grammar[replaceable_symbol].append([symbol, params])
 
-    def export_genotype(self, filepath):
+    def export_genotype(self, filepath: str) -> None:
         with open(filepath, 'w+') as file:
             self._export_genotype_open_file(file)
 
-    def _export_genotype_open_file(self, file):
+    def _export_genotype_open_file(self, file) -> None:
         for key, rule in self.grammar.items():
             line = key.value + ' '
             for item_rule in range(0, len(rule)):
@@ -74,18 +88,20 @@ def _export_genotype_open_file(self, file):
                 line += symbol + ' '
             file.write(line + '\n')
 
-    def check_validity(self):
+    def check_validity(self) -> None:
         if self.phenotype._morphological_measurements.measurement_to_dict()['hinge_count'] > 0:
             self.valid = True
 
-    def develop(self):
+    def develop(self) -> RevolveBot:
         self.phenotype = GrammarExpander(self)\
             .expand_grammar(self.conf.i_iterations, self.conf.axiom_w)\
             .decode_sentence()
         return self.phenotype
 
     @staticmethod
-    def build_symbol(symbol, conf):
+    def build_symbol(symbol: Union[Alphabet, (Alphabet, List)],
+                     conf: PlasticodingConfig
+                     ) -> (Alphabet, List):
         """
         Adds params for alphabet symbols (when it applies).
         :return: