diff --git a/.github/workflows/ci-cd.yml b/.github/workflows/ci-cd.yml
index c999d27..43c103f 100644
--- a/.github/workflows/ci-cd.yml
+++ b/.github/workflows/ci-cd.yml
@@ -24,6 +24,24 @@ jobs:
             - name: Run Clippy
               run: cargo clippy --all-features -- -D warnings
     
+    examples:
+        runs-on: ubuntu-latest
+        name: Run Examples
+        permissions:
+            contents: read
+        steps:
+            - uses: actions/checkout@v4
+            - name: Run readme_ex example
+              run: cargo run --example readme_ex --features builtin,genrand
+            - name: Run crossover example
+              run: cargo run --example crossover --features crossover,genrand
+            - name: Run knockout example
+              run: cargo run --example knockout --features knockout,crossover,genrand
+            - name: Run derive example
+              run: cargo run --example derive --features derive
+            - name: Run speciation example
+              run: cargo run --example speciation --features speciation,genrand
+
     fmt:
         runs-on: ubuntu-latest
         name: Formatting Checks
diff --git a/Cargo.lock b/Cargo.lock
index 7904c9e..99c14bc 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -119,7 +119,7 @@ checksum = "d9c4f5dac5e15c24eb999c26181a6ca40b39fe946cbe4c263c7209467bc83af2"
 
 [[package]]
 name = "genetic-rs"
-version = "1.3.0"
+version = "1.4.0"
 dependencies = [
  "genetic-rs-common",
  "genetic-rs-macros",
@@ -128,7 +128,7 @@ dependencies = [
 
 [[package]]
 name = "genetic-rs-common"
-version = "1.3.0"
+version = "1.4.0"
 dependencies = [
  "itertools",
  "rand",
@@ -137,7 +137,7 @@ dependencies = [
 
 [[package]]
 name = "genetic-rs-macros"
-version = "1.3.0"
+version = "1.4.0"
 dependencies = [
  "darling",
  "genetic-rs-common",
diff --git a/Cargo.toml b/Cargo.toml
index d350add..972d04e 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -3,9 +3,13 @@ members = ["genetic-rs", "genetic-rs-common", "genetic-rs-macros"]
 resolver = "2"
 
 [workspace.package]
-version = "1.3.0"
+version = "1.4.0"
 authors = ["HyperCodec"]
 homepage = "https://github.com/hypercodec/genetic-rs"
 repository = "https://github.com/hypercodec/genetic-rs"
 license = "MIT"
 edition = "2021"
+
+[workspace.dependencies]
+genetic-rs-common = { path = "genetic-rs-common", version = "1.4.0", default-features = false }
+genetic-rs-macros = { path = "genetic-rs-macros", version = "1.4.0" }
diff --git a/README.md b/README.md
index 89644b8..cad7566 100644
--- a/README.md
+++ b/README.md
@@ -9,10 +9,14 @@ A small framework for managing genetic algorithms.
 ### Features
 First off, this crate comes with the `builtin`, `genrand`, `crossover`, `knockout`, and `speciation` features by default. If you want the simulation to be parallelized (which is most usecases), add the `rayon` feature. There are also some convenient macros with the `derive` feature.
 
-### How to Use
-> [!NOTE] 
-> If you are interested in implementing NEAT with this, or just want a more complex example, check out the [neat](https://crates.io/crates/neat) crate
+### Ecosystem
+This framework was created with a high degree of modularity in mind, allowing other crates to contribute to the ecosystem. Here's a list of some good crates:
+- [neat](https://crates.io/crates/neat) - Handles complex reproduction and prediction logic for the NEAT algorithm, allowing you to create AI simulations with ease. It also functions as a pretty good example for the more complex usecases of the crate's traits.
+- [genetic-rs-extras](https://crates.io/crates/genetic-rs-extras) - A companion crate with quality-of-life improvements and utility features.
+
+If you have a `genetic-rs`-based crate and you'd like it to be featured here, submit a PR or discussion post and I'll consider it.
 
+### How to Use
 Here's a simple genetic algorithm:
 
 ```rust
@@ -65,7 +69,7 @@ fn main() {
 }
 ```
 
-That is the minimal code for a working genetic algorithm on default features. You can [read the docs](https://docs.rs/genetic-rs) or [check the examples](/genetic-rs/examples/) for more complicated systems. I highly recommend looking into crossover reproduction, as it tends to produce better results than mitosis.
+That is the minimal code for a working genetic algorithm with just the `derive` feature (+ defaults). You can [read the docs](https://docs.rs/genetic-rs) or [check the examples](/genetic-rs/examples/) for more complicated systems. I highly recommend looking into crossover reproduction, as it tends to produce better results than mitosis.
 
 ### License
 This project falls under the `MIT` license.
diff --git a/genetic-rs-common/src/builtin/eliminator.rs b/genetic-rs-common/src/builtin/eliminator.rs
index b220913..9130359 100644
--- a/genetic-rs-common/src/builtin/eliminator.rs
+++ b/genetic-rs-common/src/builtin/eliminator.rs
@@ -1,5 +1,3 @@
-use std::ops::Not;
-
 use crate::Eliminator;
 use crate::FeatureBoundedGenome;
 
@@ -52,6 +50,16 @@ impl<G> FitnessObserver<G> for () {
     fn observe(&mut self, _fitnesses: &[(G, f32)]) {}
 }
 
+impl<F, G> FitnessObserver<G> for F
+where
+    F: FnMut(&[(G, f32)]),
+    G: FeatureBoundedGenome,
+{
+    fn observe(&mut self, fitnesses: &[(G, f32)]) {
+        (self)(fitnesses);
+    }
+}
+
 /// An observer that calls two observers in sequence.
 /// Created by [`FitnessObserver::layer`].
 pub struct LayeredObserver<G, A: FitnessObserver<G>, B: FitnessObserver<G>>(
@@ -71,6 +79,56 @@ where
     }
 }
 
+impl<G, A, B> Default for LayeredObserver<G, A, B>
+where
+    A: Default + FitnessObserver<G>,
+    B: Default + FitnessObserver<G>,
+{
+    fn default() -> Self {
+        Self(A::default(), B::default(), std::marker::PhantomData)
+    }
+}
+
+impl<G, A, B> Clone for LayeredObserver<G, A, B>
+where
+    A: Clone + FitnessObserver<G>,
+    B: Clone + FitnessObserver<G>,
+{
+    fn clone(&self) -> Self {
+        Self(self.0.clone(), self.1.clone(), std::marker::PhantomData)
+    }
+}
+
+impl<G, A, B> std::fmt::Debug for LayeredObserver<G, A, B>
+where
+    A: std::fmt::Debug + FitnessObserver<G>,
+    B: std::fmt::Debug + FitnessObserver<G>,
+{
+    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
+        f.debug_tuple("LayeredObserver")
+            .field(&self.0)
+            .field(&self.1)
+            .finish()
+    }
+}
+
+impl<G, A, B> PartialEq for LayeredObserver<G, A, B>
+where
+    A: PartialEq + FitnessObserver<G>,
+    B: PartialEq + FitnessObserver<G>,
+{
+    fn eq(&self, other: &Self) -> bool {
+        self.0 == other.0 && self.1 == other.1
+    }
+}
+
+impl<G, A, B> Eq for LayeredObserver<G, A, B>
+where
+    A: Eq + FitnessObserver<G>,
+    B: Eq + FitnessObserver<G>,
+{
+}
+
 #[cfg(not(feature = "rayon"))]
 #[doc(hidden)]
 pub trait FeatureBoundedFitnessObserver<G: FeatureBoundedGenome>: FitnessObserver<G> {}
@@ -330,7 +388,7 @@ mod knockout {
         }
     }
 
-    impl Not for KnockoutWinner {
+    impl std::ops::Not for KnockoutWinner {
         type Output = Self;
 
         fn not(self) -> Self::Output {
@@ -523,3 +581,150 @@ mod knockout {
 
 #[cfg(feature = "knockout")]
 pub use knockout::*;
+
+#[cfg(feature = "speciation")]
+mod speciation {
+    use crate::{prelude::*, speciation::SpeciatedPopulation};
+
+    /// An eliminator that attempts to preserve new experimental structures by dividing a genome's
+    /// fitness by the number of genomes in its species.
+    pub struct SpeciatedFitnessEliminator<
+        F: FeatureBoundedFitnessFn<G>,
+        G: Speciated + FeatureBoundedGenome,
+        O: FeatureBoundedFitnessObserver<G> = (),
+    > {
+        /// The divergence threshold used to determine whether a genome belongs in a species.
+        /// See [`SpeciatedPopulation::threshold`] for more info.
+        pub speciation_threshold: f32,
+
+        /// The inner fitness eliminator used to hold settings and such.
+        pub inner: FitnessEliminator<F, G, O>,
+
+        /// The context used to determine species membership.
+        /// This is necessary since the eliminator needs to calculate species membership to divide fitness by species size, and some genome types may require context to calculate divergence.
+        pub ctx: <G as Speciated>::Context,
+
+        _marker: std::marker::PhantomData<G>,
+    }
+
+    impl<F, G, O> SpeciatedFitnessEliminator<F, G, O>
+    where
+        F: FeatureBoundedFitnessFn<G>,
+        G: Speciated + FeatureBoundedGenome,
+        O: FeatureBoundedFitnessObserver<G>,
+    {
+        /// Creates a new [`SpeciatedFitnessEliminator`] with a given fitness function and thresholds.
+        pub fn new(
+            fitness_fn: F,
+            speciation_threshold: f32,
+            keep_threshold: f32,
+            observer: O,
+            ctx: <G as Speciated>::Context,
+        ) -> Self {
+            Self {
+                speciation_threshold,
+                inner: FitnessEliminator::new(fitness_fn, keep_threshold, observer),
+                ctx,
+                _marker: std::marker::PhantomData,
+            }
+        }
+
+        /// Creates a new [`SpeciatedFitnessEliminator`] from a regular [`FitnessEliminator`] and a speciation threshold.
+        /// Useful since you can use the builder for [`FitnessEliminator`] to construct the fitness function and observer, then convert it into a [`SpeciatedFitnessEliminator`] with this method.
+        pub fn from_fitness_eliminator(
+            fitness_eliminator: FitnessEliminator<F, G, O>,
+            speciation_threshold: f32,
+            ctx: <G as Speciated>::Context,
+        ) -> Self {
+            Self {
+                speciation_threshold,
+                inner: fitness_eliminator,
+                ctx,
+                _marker: std::marker::PhantomData,
+            }
+        }
+
+        /// Calculates the fitness of each genome, dividing by the number of genomes in its species, and sorts them by fitness.
+        /// Returns a vector of tuples containing the genome and its fitness score.
+        #[cfg(not(feature = "rayon"))]
+        pub fn calculate_and_sort(&self, genomes: Vec<G>) -> Vec<(G, f32)> {
+            let population =
+                SpeciatedPopulation::from_genomes(&genomes, self.speciation_threshold, &self.ctx);
+            let mut fitnesses = vec![0.0; genomes.len()];
+
+            for species in population.species() {
+                let len = species.len() as f32;
+                debug_assert!(len != 0.0);
+                for &index in species {
+                    let genome = &genomes[index];
+                    let fitness = self.inner.fitness_fn.fitness(genome);
+                    if fitness < 0.0 {
+                        fitnesses[index] = fitness * len;
+                    } else {
+                        fitnesses[index] = fitness / len;
+                    }
+                }
+            }
+
+            let mut fitnesses: Vec<(G, f32)> = genomes.into_iter().zip(fitnesses).collect();
+            fitnesses.sort_by(|(_a, afit), (_b, bfit)| bfit.partial_cmp(afit).unwrap());
+            fitnesses
+        }
+
+        /// Calculates the fitness of each genome, dividing by the number of genomes in its species, and sorts them by fitness.
+        /// Returns a vector of tuples containing the genome and its fitness score.
+        #[cfg(feature = "rayon")]
+        pub fn calculate_and_sort(&self, genomes: Vec<G>) -> Vec<(G, f32)> {
+            let population =
+                SpeciatedPopulation::from_genomes(&genomes, self.speciation_threshold, &self.ctx);
+
+            let mut fitnesses = vec![0.0; genomes.len()];
+
+            for species in population.species() {
+                let len = species.len() as f32;
+                debug_assert!(len != 0.0);
+                for &index in species {
+                    let genome = &genomes[index];
+                    let fitness = self.inner.fitness_fn.fitness(genome);
+                    if fitness < 0.0 {
+                        fitnesses[index] = fitness * len;
+                    } else {
+                        fitnesses[index] = fitness / len;
+                    }
+                }
+            }
+
+            let mut result: Vec<(G, f32)> = genomes.into_iter().zip(fitnesses).collect();
+            result.sort_by(|(_a, afit), (_b, bfit)| bfit.partial_cmp(afit).unwrap());
+            result
+        }
+    }
+
+    impl<F, G, O> Eliminator<G> for SpeciatedFitnessEliminator<F, G, O>
+    where
+        F: FeatureBoundedFitnessFn<G>,
+        G: Speciated + FeatureBoundedGenome,
+        O: FeatureBoundedFitnessObserver<G>,
+    {
+        #[cfg(not(feature = "rayon"))]
+        fn eliminate(&mut self, genomes: Vec<G>) -> Vec<G> {
+            let mut fitnesses = self.calculate_and_sort(genomes);
+            self.inner.observer.observe(&fitnesses);
+            let median_index = (fitnesses.len() as f32) * self.inner.threshold;
+            fitnesses.truncate(median_index as usize + 1);
+            fitnesses.into_iter().map(|(g, _)| g).collect()
+        }
+
+        #[cfg(feature = "rayon")]
+        fn eliminate(&mut self, genomes: Vec<G>) -> Vec<G> {
+            let mut fitnesses = self.calculate_and_sort(genomes);
+            self.inner.observer.observe(&fitnesses);
+            let median_index = (fitnesses.len() as f32) * self.inner.threshold;
+            fitnesses.truncate(median_index as usize + 1);
+            fitnesses.into_par_iter().map(|(g, _)| g).collect()
+        }
+    }
+}
+
+#[cfg(feature = "speciation")]
+pub use speciation::*;
diff --git a/genetic-rs-common/src/builtin/repopulator.rs b/genetic-rs-common/src/builtin/repopulator.rs
index 5f13994..0232698 100644
--- a/genetic-rs-common/src/builtin/repopulator.rs
+++ b/genetic-rs-common/src/builtin/repopulator.rs
@@ -84,6 +84,27 @@ where
     }
 }
 
+/// Helper trait for creating a new population from a parent genome. This is a shortcut for using [`MitosisRepopulator`]
+/// and is convenient for things like resuming a simulation from a deserialized elite.
+pub trait FromParent<G: Mitosis> {
+    /// Create a new population from a parent genome and a count of how many genomes to create.
+    /// The new population should be created by mutating the parent genome (i.e. [`Mitosis`]).
+    /// The function is essentially a shortcut for creating a new [`MitosisRepopulator`] and using
+    /// it to repopulate a new population from the parent genome.
+    fn from_parent(parent: G, count: usize, ctx: G::Context, rate: f32) -> Self;
+}
+
+impl<G: Mitosis> FromParent<G> for Vec<G> {
+    fn from_parent(parent: G, count: usize, ctx: G::Context, rate: f32) -> Self {
+        assert!(count > 0, "Count must be greater than 0");
+
+        let mut repopulator = MitosisRepopulator::new(rate, ctx);
+        let mut genomes = vec![parent];
+        repopulator.repopulate(&mut genomes, count);
+        genomes
+    }
+}
+
 #[cfg(feature = "crossover")]
 mod crossover {
     use rand::RngExt;
@@ -150,6 +171,16 @@ mod crossover {
             }
         }
     }
+
+    impl<G> Default for CrossoverRepopulator<G>
+    where
+        G: Crossover,
+        G::Context: Default,
+    {
+        fn default() -> Self {
+            Self::new(0.05, G::Context::default())
+        }
+    }
 }
 
 #[cfg(feature = "crossover")]
@@ -157,43 +188,82 @@ pub use crossover::*;
 
 #[cfg(feature = "speciation")]
 mod speciation {
-    use std::collections::HashMap;
-
     use rand::RngExt;
 
-    use super::*;
+    use crate::speciation::{Speciated, SpeciatedPopulation};
 
-    /// Used in speciated crossover nextgens. Allows for genomes to avoid crossover with ones that are too different.
-    pub trait Speciated {
-        /// The type used to distinguish
-        /// one genome's species from another.
-        type Species: Eq + std::hash::Hash; // I really don't like that we need `Eq` when `PartialEq` better fits the definiton.
+    use super::*;
 
-        /// Get/calculate this genome's species.
-        fn species(&self) -> Self::Species;
+    /// The action to take when a genome is found to be in a species by itself.
+    /// This can be used to prevent species from going extinct due to bad luck in crossover.
+    pub enum ActionIfIsolated {
+        /// Do nothing, allowing the species to go extinct if the genome is not fit enough.
+        /// If all species happen to be isolated (no species has >= 2 members), the
+        /// repopulator automatically falls back to [`CrossoverRepopulator`] to avoid
+        /// an infinite loop.
+        DoNothing,
+
+        /// Perform crossover between the genome and itself to create a new member of the species.
+        /// This can help prevent species from going extinct, but can also lead to less diversity in the population.
+        CrossoverSelf,
+
+        /// Perform crossover between the genome and a random member of the most similar species to create a new member of the species.
+        /// This can help prevent species from going extinct while maintaining more diversity than self-crossover, but can also lead to more computational overhead.
+        CrossoverSimilarSpecies,
+
+        /// Perform crossover between the genome and a random member of the entire population, ignoring species boundaries, to create a new member of the species.
+        /// This can help prevent species from going extinct, but can also contribute to more broken or dysfunctional species.
+        CrossoverRandom,
     }
 
     /// Repopulator that uses crossover reproduction to create new genomes, but only between genomes of the same species.
     pub struct SpeciatedCrossoverRepopulator<G: Crossover + Speciated> {
-        /// The inner crossover repopulator. This holds the settings for crossover operations,
-        /// but may also be called if [`allow_emergency_repr`][Self::allow_emergency_repr] is `true`.
-        pub crossover: CrossoverRepopulator<G>,
+        /// The inner crossover repopulator. This holds the settings for crossover operations.
+        pub inner: CrossoverRepopulator<G>,
+
+        /// The threshold used to determine if a genome belongs in a species.
+        /// See [`SpeciatedPopulation::threshold`] for more info.
+        pub speciation_threshold: f32,
 
-        /// Whether to allow genomes to reproduce across species boundaries
-        /// (effectively vanilla crossover)
-        /// in emergency situations where no genomes have compatible partners.
-        /// If disabled, the simulation will panic in such a situation.
-        pub allow_emergency_repr: bool,
+        /// The action to take when a genome is found to be in a species by itself.
+        pub action_if_isolated: ActionIfIsolated,
+
+        /// Additional context for speciation.
+        pub ctx: <G as Speciated>::Context,
 
         _marker: std::marker::PhantomData<G>,
     }
 
     impl<G: Crossover + Speciated> SpeciatedCrossoverRepopulator<G> {
         /// Creates a new [`SpeciatedCrossoverRepopulator`].
-        pub fn new(mutation_rate: f32, allow_emergency_repr: bool, ctx: G::Context) -> Self {
+        pub fn new(
+            mutation_rate: f32,
+            threshold: f32,
+            action_if_isolated: ActionIfIsolated,
+            crossover_ctx: <G as Crossover>::Context,
+            spec_ctx: <G as Speciated>::Context,
+        ) -> Self {
             Self {
-                crossover: CrossoverRepopulator::new(mutation_rate, ctx),
-                allow_emergency_repr,
+                inner: CrossoverRepopulator::new(mutation_rate, crossover_ctx),
+                ctx: spec_ctx,
+                speciation_threshold: threshold,
+                action_if_isolated,
+                _marker: std::marker::PhantomData,
+            }
+        }
+
+        /// Creates a new [`SpeciatedCrossoverRepopulator`] from an existing [`CrossoverRepopulator`], using the same mutation settings.
+        pub fn from_crossover(
+            inner: CrossoverRepopulator<G>,
+            threshold: f32,
+            action_if_isolated: ActionIfIsolated,
+            ctx: <G as Speciated>::Context,
+        ) -> Self {
+            Self {
+                inner,
+                ctx,
+                speciation_threshold: threshold,
+                action_if_isolated,
                 _marker: std::marker::PhantomData,
             }
         }
@@ -203,60 +273,124 @@ mod speciation {
     where
         G: Crossover + Speciated,
     {
-        // i'm still not really satisfied with this implementation,
-        // but it's better than the old one.
         fn repopulate(&mut self, genomes: &mut Vec<G>, target_size: usize) {
             let initial_size = genomes.len();
-            let mut rng = rand::rng();
-            let mut species: HashMap<<G as Speciated>::Species, Vec<&G>> = HashMap::new();
-
-            for genome in genomes.iter() {
-                let spec = genome.species();
-                species.entry(spec).or_insert_with(Vec::new).push(genome);
+            if initial_size >= target_size {
+                return;
             }
 
-            let mut species_iter = species.values();
-            let to_create = target_size - initial_size;
-            let mut new_genomes = Vec::with_capacity(to_create);
-
-            while new_genomes.len() < to_create {
-                if let Some(spec) = species_iter.next() {
-                    if spec.len() < 2 {
-                        continue;
-                    }
+            let mut rng = rand::rng();
+            let population =
+                SpeciatedPopulation::from_genomes(genomes, self.speciation_threshold, &self.ctx);
+
+            // if all species are isolated, we fall back to the inner crossover repopulator to avoid an infinite loop.
+            if (matches!(self.action_if_isolated, ActionIfIsolated::DoNothing)
+                && !population.species().iter().any(|s| s.len() >= 2))
+                || (matches!(self.action_if_isolated, ActionIfIsolated::CrossoverRandom)
+                    && initial_size == 1)
+            {
+                self.inner.repopulate(genomes, target_size);
+                return;
+            }
 
-                    for (i, &parent1) in spec.iter().enumerate() {
-                        let mut j = rng.random_range(1..spec.len());
-                        if j == i {
-                            j = 0;
+            let amount_to_make = target_size - initial_size;
+            let mut species_cycle = population.round_robin_enumerate();
+
+            let mut i = 0;
+            while i < amount_to_make {
+                let (species_i, genome_i) = species_cycle.next().unwrap();
+                let species = &population.species()[species_i];
+                let parent1 = &genomes[genome_i];
+                if species.len() < 2 {
+                    match self.action_if_isolated {
+                        ActionIfIsolated::DoNothing => continue,
+                        ActionIfIsolated::CrossoverSelf => {
+                            let child = parent1.crossover(
+                                parent1,
+                                &self.inner.ctx,
+                                self.inner.mutation_rate,
+                                &mut rng,
+                            );
+                            genomes.push(child);
+                            i += 1;
+                            continue;
                         }
-                        let parent2 = spec[j];
-
-                        new_genomes.push(parent1.crossover(
-                            parent2,
-                            &self.crossover.ctx,
-                            self.crossover.mutation_rate,
-                            &mut rng,
-                        ));
-                    }
-                } else {
-                    // reached the end, reset the iterator
-
-                    if new_genomes.is_empty() {
-                        // no genomes have compatible partners
-                        if self.allow_emergency_repr {
-                            self.crossover.repopulate(genomes, target_size);
-                            return;
-                        } else {
-                            panic!("no genomes with common species");
+                        ActionIfIsolated::CrossoverSimilarSpecies => {
+                            let mut best_species_i = 0;
+                            let mut best_divergence = f32::MAX;
+                            for (j, species) in population.species().iter().enumerate() {
+                                if j == species_i || species.is_empty() {
+                                    continue;
+                                }
+                                let representative = &genomes[species[0]];
+                                let divergence = parent1.divergence(representative, &self.ctx);
+                                if divergence < best_divergence {
+                                    best_divergence = divergence;
+                                    best_species_i = j;
+                                }
+                            }
+
+                            let best_species = &population.species()[best_species_i];
+                            let j = rng.random_range(0..best_species.len());
+                            let parent2 = &genomes[best_species[j]];
+                            let child = parent1.crossover(
+                                parent2,
+                                &self.inner.ctx,
+                                self.inner.mutation_rate,
+                                &mut rng,
+                            );
+                            genomes.push(child);
+                            i += 1;
+                            continue;
+                        }
+                        ActionIfIsolated::CrossoverRandom => {
+                            let mut j = rng.random_range(1..initial_size);
+                            if j == genome_i {
+                                j = 0;
+                            }
+
+                            let parent2 = &genomes[j];
+                            let child = parent1.crossover(
+                                parent2,
+                                &self.inner.ctx,
+                                self.inner.mutation_rate,
+                                &mut rng,
+                            );
+                            genomes.push(child);
+                            i += 1;
+                            continue;
                         }
                     }
+                }
 
-                    species_iter = species.values();
+                let mut j = rng.random_range(1..species.len());
+                if genome_i == species[j] {
+                    j = 0;
                 }
+                let parent2 = &genomes[species[j]];
+
+                let child =
+                    parent1.crossover(parent2, &self.inner.ctx, self.inner.mutation_rate, &mut rng);
+                genomes.push(child);
+
+                i += 1;
             }
+        }
+    }
 
-            genomes.extend(new_genomes);
+    impl<G> Default for SpeciatedCrossoverRepopulator<G>
+    where
+        G: Crossover + Speciated,
+        <G as Crossover>::Context: Default,
+        <G as Speciated>::Context: Default,
+    {
+        fn default() -> Self {
+            Self::from_crossover(
+                CrossoverRepopulator::default(),
+                0.1,
+                ActionIfIsolated::CrossoverSimilarSpecies,
+                <G as Speciated>::Context::default(),
+            )
         }
     }
 }
diff --git a/genetic-rs-common/src/lib.rs b/genetic-rs-common/src/lib.rs
index 05b0f2c..39e545c 100644
--- a/genetic-rs-common/src/lib.rs
+++ b/genetic-rs-common/src/lib.rs
@@ -8,6 +8,10 @@
 #[cfg(feature = "builtin")]
 pub mod builtin;
 
+/// Common speciation code used for speciated eliminators and repopulators.
+#[cfg(feature = "speciation")]
+pub mod speciation;
+
 /// Used to quickly import everything this crate has to offer.
 /// Simply add `use genetic_rs::prelude::*` to begin using this crate.
 pub mod prelude;
@@ -126,16 +130,7 @@ where
     fn gen_random(rng: &mut impl rand::Rng, amount: usize) -> Self;
 }
 
-/// Rayon version of the [`GenerateRandomCollection`] trait
-#[cfg(all(feature = "genrand", feature = "rayon"))]
-pub trait GenerateRandomCollectionParallel<T>
-where
-    T: GenerateRandom + Send,
-{
-    /// Generate a random collection of the inner objects with the given amount. Does not pass in rng like the sync counterpart.
-    fn par_gen_random(amount: usize) -> Self;
-}
-
+#[cfg(feature = "genrand")]
 impl<C, T> GenerateRandomCollection<T> for C
 where
     C: FromIterator<T>,
@@ -146,6 +141,16 @@ where
     }
 }
 
+/// Rayon version of the [`GenerateRandomCollection`] trait
+#[cfg(all(feature = "genrand", feature = "rayon"))]
+pub trait GenerateRandomCollectionParallel<T>
+where
+    T: GenerateRandom + Send,
+{
+    /// Generate a random collection of the inner objects with the given amount. Does not pass in rng like the sync counterpart.
+    fn par_gen_random(amount: usize) -> Self;
+}
+
 #[cfg(feature = "rayon")]
 impl<C, T> GenerateRandomCollectionParallel<T> for C
 where
diff --git a/genetic-rs-common/src/prelude.rs b/genetic-rs-common/src/prelude.rs
index bd3c19f..ade6ac9 100644
--- a/genetic-rs-common/src/prelude.rs
+++ b/genetic-rs-common/src/prelude.rs
@@ -5,4 +5,7 @@ pub use crate::*;
 #[cfg(feature = "builtin")]
 pub use crate::builtin::{eliminator::*, repopulator::*};
 
+#[cfg(feature = "speciation")]
+pub use crate::speciation::Speciated;
+
 pub use rand::prelude::*;
diff --git a/genetic-rs-common/src/speciation.rs b/genetic-rs-common/src/speciation.rs
new file mode 100644
index 0000000..402ab18
--- /dev/null
+++ b/genetic-rs-common/src/speciation.rs
@@ -0,0 +1,141 @@
+/// Trait that allows a genome to be separated into species.
+pub trait Speciated {
+    /// The context type for speciation. This can be used to provide additional information
+    /// that may be needed to calculate the divergence between two genomes.
+    type Context;
+
+    /// How structurally different this genome is to another genome. The higher the value, the more different they are.
+    /// When implementing this, you should consider only the major structural differences that could
+    /// take several generations to evolve an optimal solution. Do not consider minor differences that
+    /// are easily evolved in a single generation, such as a single weight change.
+    /// By common convention, this is typically a value between 0 and 1.
+    /// Additionally, it should also be symmetric, meaning that the divergence between genome A and genome B should be
+    /// the same as the divergence between genome B and genome A.
+    /// This value is used to separate genomes into a species.
+    fn divergence(&self, other: &Self, ctx: &Self::Context) -> f32;
+}
+
+/// A population of genomes that have been separated into species.
+pub struct SpeciatedPopulation {
+    /// The species in this population. Each species is a vector of indices into the original genome vector.
+    /// The first genome in a species is its representation (i.e. the one that gets compared to other genomes to determine
+    /// if they belong in the species).
+    /// Invariant: every inner `Vec` is non-empty.
+    species: Vec<Vec<usize>>,
+
+    /// The threshold used to determine if a genome belongs in a species. If the divergence between a genome and the representative genome
+    /// of a species is less than this threshold, then the genome belongs in that species.
+    pub threshold: f32,
+}
+
+impl SpeciatedPopulation {
+    /// Creates a new, empty [`SpeciatedPopulation`] with the given threshold.
+    pub fn new(threshold: f32) -> Self {
+        Self {
+            species: Vec::new(),
+            threshold,
+        }
+    }
+
+    /// Returns the species in this population.
+    /// Each inner slice is guaranteed to be non-empty.
+    pub fn species(&self) -> &[Vec<usize>] {
+        &self.species
+    }
+
+    /// Inserts a genome into the speciated population.
+    /// Returns whether a new species was created by this insertion.
+    pub fn insert_genome<G: Speciated>(
+        &mut self,
+        index: usize,
+        genomes: &[G],
+        ctx: &G::Context,
+    ) -> bool {
+        if let Some(species) = self.get_species_mut(index, genomes, ctx) {
+            species.push(index);
+            return false;
+        }
+        self.species.push(vec![index]);
+        true
+    }
+
+    /// Creates a new speciated population from a population of genomes.
+    /// Note that this can be O(n^2) worst case, but is typically much faster in practice,
+    /// especially if the genome structure doesn't mutate often.
+    pub fn from_genomes<G: Speciated>(population: &[G], threshold: f32, ctx: &G::Context) -> Self {
+        let mut speciated_population = SpeciatedPopulation::new(threshold);
+        for index in 0..population.len() {
+            speciated_population.insert_genome(index, population, ctx);
+        }
+        speciated_population
+    }
+
+    /// Gets the species that a genome belongs to, if any.
+    pub fn get_species<G: Speciated>(
+        &self,
+        index: usize,
+        genomes: &[G],
+        ctx: &G::Context,
+    ) -> Option<&Vec<usize>> {
+        let genome = &genomes[index];
+        self.species
+            .iter()
+            .find(|species| genome.divergence(&genomes[species[0]], ctx) < self.threshold)
+    }
+
+    /// Get a mutable reference to the species that a genome belongs to, if any.
+    pub fn get_species_mut<G: Speciated>(
+        &mut self,
+        index: usize,
+        genomes: &[G],
+        ctx: &G::Context,
+    ) -> Option<&mut Vec<usize>> {
+        let genome = &genomes[index];
+        self.species
+            .iter_mut()
+            .find(|species| genome.divergence(&genomes[species[0]], ctx) < self.threshold)
+    }
+
+    /// Round robin cyclic iterator over the genomes in this population.
+    /// Note that the index value is the index in the original genome vector, not the index in the species vector.
+    pub fn round_robin(&self) -> impl Iterator<Item = usize> + '_ {
+        let species = &self.species;
+        let mut idx_in_species = vec![0; self.species.len()];
+        let mut species_i = 0usize;
+
+        std::iter::from_fn(move || {
+            if species.is_empty() {
+                return None;
+            }
+
+            let cur_species = &species[species_i];
+            let genome_index = cur_species[idx_in_species[species_i]];
+            idx_in_species[species_i] = (idx_in_species[species_i] + 1) % cur_species.len();
+            species_i = (species_i + 1) % species.len();
+            Some(genome_index)
+        })
+    }
+
+    /// Round robin, but also returns the species index along with the genome index.
+    /// Note that the genome index value is the index in the original genome vector, not the index in the species vector.
+    /// However, the species index represents the index in [`Self::species`].
+    /// Format is (species_index, genome_index).
+    pub fn round_robin_enumerate(&self) -> impl Iterator<Item = (usize, usize)> + '_ {
+        let species = &self.species;
+        let mut idx_in_species = vec![0; self.species.len()];
+        let mut species_i = 0usize;
+
+        std::iter::from_fn(move || {
+            if species.is_empty() {
+                return None;
+            }
+
+            let cur_species = &species[species_i];
+            let cur_species_i = species_i;
+            let genome_index = cur_species[idx_in_species[species_i]];
+            idx_in_species[species_i] = (idx_in_species[species_i] + 1) % cur_species.len();
+            species_i = (species_i + 1) % species.len();
+            Some((cur_species_i, genome_index))
+        })
+    }
+}
diff --git a/genetic-rs-macros/Cargo.toml b/genetic-rs-macros/Cargo.toml
index 24d55b9..7f9aacc 100644
--- a/genetic-rs-macros/Cargo.toml
+++ b/genetic-rs-macros/Cargo.toml
@@ -28,7 +28,7 @@ unexpected_cfgs = { level = "allow", check-cfg = ["cfg(docsrs)"] }
 
 [dependencies]
 darling = "0.23.0"
-genetic-rs-common = { path = "../genetic-rs-common", version = "1.3.0" }
+genetic-rs-common.workspace = true
 proc-macro2 = "1.0.106"
 quote = "1.0.44"
 syn = "2.0.114"
diff --git a/genetic-rs/Cargo.toml b/genetic-rs/Cargo.toml
index 25d76c3..97fa327 100644
--- a/genetic-rs/Cargo.toml
+++ b/genetic-rs/Cargo.toml
@@ -22,8 +22,8 @@ rayon = ["genetic-rs-common/rayon"]
 derive = ["dep:genetic-rs-macros", "builtin"]
 
 [dependencies]
-genetic-rs-common = { path = "../genetic-rs-common", version = "1.3.0", default-features = false }
-genetic-rs-macros = { path = "../genetic-rs-macros", version = "1.3.0", optional = true }
+genetic-rs-common.workspace = true
+genetic-rs-macros = { workspace = true, optional = true }
 
 [package.metadata.docs.rs]
 all-features = true
@@ -35,18 +35,42 @@ unexpected_cfgs = { level = "allow", check-cfg = ["cfg(publish)", "cfg(docsrs)"]
 [dev-dependencies]
 rand = "0.10.0"
 
+[[example]]
+name = "readme_ex"
+required-features = ["builtin", "genrand"]
+
 [[example]]
 name = "crossover"
-required-features = ["crossover"]
+required-features = ["crossover", "genrand"]
+
+[[example]]
+name = "knockout"
+required-features = ["knockout", "crossover", "genrand"]
 
 [[example]]
 name = "speciation"
-required-features = ["speciation"]
+required-features = ["speciation", "genrand"]
 
 [[example]]
 name = "derive"
 required-features = ["derive"]
 
+[[test]]
+name = "simulation"
+required-features = ["builtin", "genrand"]
+
+[[test]]
+name = "knockout"
+required-features = ["knockout", "genrand"]
+
+[[test]]
+name = "repopulator"
+required-features = ["crossover", "genrand"]
+
+[[test]]
+name = "speciation"
+required-features = ["speciation", "genrand"]
+
 [[test]]
 name = "derive_macros"
 required-features = ["derive", "genrand", "crossover"]
diff --git a/genetic-rs/examples/speciation.rs b/genetic-rs/examples/speciation.rs
index 662ee45..a9f3635 100644
--- a/genetic-rs/examples/speciation.rs
+++ b/genetic-rs/examples/speciation.rs
@@ -1,26 +1,56 @@
 use genetic_rs::prelude::*;
 
+/// Max value of divergence() that will count two genomes as the same species.
+const SPECIATION_THRESHOLD: f32 = 0.1;
+
 #[derive(Clone, PartialEq, Debug)]
 struct MyGenome {
-    val1: f32,
-    val2: f32,
+    vals: Vec<f32>,
+}
+
+impl GenerateRandom for MyGenome {
+    fn gen_random(rng: &mut impl rand::Rng) -> Self {
+        Self {
+            vals: (0..10).map(|_| rng.random_range(-1.0..1.0)).collect(),
+        }
+    }
 }
 
 impl RandomlyMutable for MyGenome {
     type Context = ();
 
-    fn mutate(&mut self, _: &(), rate: f32, rng: &mut impl Rng) {
-        self.val1 += rng.random_range(-1.0..1.0) * rate;
-        self.val2 += rng.random_range(-1.0..1.0) * rate;
+    fn mutate(&mut self, _ctx: &Self::Context, rate: f32, rng: &mut impl rand::Rng) {
+        // internal mutation
+        for val in &mut self.vals {
+            if rng.random_bool(rate as f64) {
+                *val += rng.random_range(-0.1..0.1);
+            }
+        }
+
+        // structural mutations
+        if rng.random_bool(rate as f64) {
+            // add a new random gene
+            self.vals.push(rng.random_range(-1.0..1.0));
+        }
+        if self.vals.len() > 1 && rng.random_bool(rate as f64) {
+            // remove a random gene
+            let index = rng.random_range(0..self.vals.len());
+            self.vals.remove(index);
+        }
     }
 }
 
 impl Mitosis for MyGenome {
     type Context = ();
 
-    fn divide(&self, _: &(), rate: f32, rng: &mut impl Rng) -> Self {
+    fn divide(
+        &self,
+        ctx: &<Self as Mitosis>::Context,
+        rate: f32,
+        rng: &mut impl rand::Rng,
+    ) -> Self {
         let mut child = self.clone();
-        child.mutate(&(), rate, rng);
+        child.mutate(ctx, rate, rng);
         child
     }
 }
@@ -28,52 +58,96 @@ impl Mitosis for MyGenome {
 impl Crossover for MyGenome {
     type Context = ();
 
-    fn crossover(&self, other: &Self, _: &(), rate: f32, rng: &mut impl Rng) -> Self {
-        let mut child = Self {
-            val1: (self.val1 + other.val1) / 2.,
-            val2: (self.val2 + other.val2) / 2.,
-        };
-        child.mutate(&(), rate, rng);
+    fn crossover(
+        &self,
+        other: &Self,
+        ctx: &Self::Context,
+        rate: f32,
+        rng: &mut impl rand::Rng,
+    ) -> Self {
+        // even though they ideally shouldn't crossover if they aren't the same length
+        // because of speciation, we can still allow it in emergency scenarios to prevent extinction
+        // and encourage more diversity, at the cost of more broken or dysfunctional species
+        let mut child;
+        let smaller;
+        if self.vals.len() < other.vals.len() {
+            child = other.clone();
+            smaller = self;
+        } else {
+            child = self.clone();
+            smaller = other;
+        }
+
+        for i in 0..smaller.vals.len() {
+            if rng.random_bool(0.5) {
+                child.vals[i] = smaller.vals[i];
+            }
+        }
+
+        child.mutate(ctx, rate, rng);
+
         child
     }
 }
 
-#[derive(PartialEq, Eq, Hash)]
-struct MySpecies(u32);
-
 impl Speciated for MyGenome {
-    type Species = MySpecies;
+    type Context = ();
 
-    fn species(&self) -> Self::Species {
-        MySpecies((self.val1 + self.val2).round() as u32)
+    fn divergence(&self, other: &Self, _ctx: &Self::Context) -> f32 {
+        // distance in lengths divided by the larger length
+        let larger = self.vals.len().max(other.vals.len());
+        // can't be 0 because of how we set up mutate()
+        assert!(larger != 0);
+        let length_diff = (self.vals.len() as isize - other.vals.len() as isize).abs() as f32;
+        length_diff / larger as f32
     }
 }
 
-impl GenerateRandom for MyGenome {
-    fn gen_random(rng: &mut impl Rng) -> Self {
-        Self {
-            val1: rng.random_range(-1.0..1.0),
-            val2: rng.random_range(-1.0..1.0),
-        }
-    }
+fn fitness(genome: &MyGenome) -> f32 {
+    // the fitness function is just the sum of the genes, which encourages longer genomes.
+    // however, since the structure mutation can add negative values initially, speciation
+    // will help protect those genomes until they can mutate the internal state into positive values.
+    genome.vals.iter().sum()
 }
 
-fn fitness(g: &MyGenome) -> f32 {
-    // train to make val1 and val2 as different as possible
-    (g.val1 - g.val2).abs()
+fn print_fitnesses(fitnesses: &[(MyGenome, f32)]) {
+    // note that with SpeciatedFitnessEliminator,
+    // these values are divided by the number of genomes in the species if positive,
+    // multiplied if negative.
+    let hi = fitnesses[0].1;
+    let med = fitnesses[fitnesses.len() / 2].1;
+    let lo = fitnesses[fitnesses.len() - 1].1;
+
+    println!("hi: {hi} med: {med} lo: {lo}");
 }
 
 fn main() {
     let mut rng = rand::rng();
 
+    let fitness_eliminator = FitnessEliminator::builder()
+        .fitness_fn(fitness)
+        .observer(print_fitnesses)
+        .build_or_panic();
+
+    let crossover_rep = CrossoverRepopulator::default();
+
     let mut sim = GeneticSim::new(
         Vec::gen_random(&mut rng, 100),
-        FitnessEliminator::new_without_observer(fitness),
-        // 25% mutation rate, allow cross-species reproduction in emergency scenarios
-        SpeciatedCrossoverRepopulator::new(0.25, true, ()),
+        SpeciatedFitnessEliminator::from_fitness_eliminator(
+            fitness_eliminator,
+            SPECIATION_THRESHOLD,
+            (),
+        ),
+        SpeciatedCrossoverRepopulator::from_crossover(
+            crossover_rep,
+            SPECIATION_THRESHOLD,
+            ActionIfIsolated::CrossoverSimilarSpecies,
+            (),
+        ),
     );
 
     sim.perform_generations(100);
 
-    dbg!(sim.genomes);
+    // rerunning multiple times should show a lot of diversity in genome length
+    dbg!(&sim.genomes[0]);
 }
diff --git a/genetic-rs/tests/eliminator.rs b/genetic-rs/tests/eliminator.rs
index f4c2b26..8738c7e 100644
--- a/genetic-rs/tests/eliminator.rs
+++ b/genetic-rs/tests/eliminator.rs
@@ -12,6 +12,15 @@ impl<'a, G> FitnessObserver<G> for CountingObserver<'a> {
     }
 }
 
+#[derive(Debug, Default, Clone, PartialEq, Eq)]
+struct TrackingObserver(usize);
+
+impl<G> FitnessObserver<G> for TrackingObserver {
+    fn observe(&mut self, _fitnesses: &[(G, f32)]) {
+        self.0 += 1;
+    }
+}
+
 #[test]
 fn layered_observer_calls_both() {
     let count_a = Cell::new(0usize);
@@ -43,3 +52,50 @@ fn layered_observer_triple_layer() {
     assert_eq!(count_b.get(), 1);
     assert_eq!(count_c.get(), 1);
 }
+
+#[test]
+fn layered_observer_default() {
+    let layered: LayeredObserver<(), TrackingObserver, TrackingObserver> =
+        LayeredObserver::default();
+    assert_eq!(layered.0, TrackingObserver(0));
+    assert_eq!(layered.1, TrackingObserver(0));
+}
+
+#[test]
+fn layered_observer_clone() {
+    let original = TrackingObserver::default().layer(TrackingObserver::default());
+    let mut cloned: LayeredObserver<(), _, _> = original.clone();
+
+    let fitnesses: Vec<((), f32)> = vec![((), 1.0)];
+    cloned.observe(&fitnesses);
+
+    assert_eq!(cloned.0, TrackingObserver(1));
+    assert_eq!(cloned.1, TrackingObserver(1));
+    assert_eq!(original.0, TrackingObserver(0));
+    assert_eq!(original.1, TrackingObserver(0));
+}
+
+#[test]
+fn layered_observer_debug() {
+    let layered: LayeredObserver<(), TrackingObserver, TrackingObserver> =
+        TrackingObserver::default().layer(TrackingObserver::default());
+    let debug_str = format!("{:?}", layered);
+    assert!(debug_str.contains("LayeredObserver"));
+    assert!(debug_str.contains("TrackingObserver"));
+}
+
+#[test]
+fn layered_observer_partial_eq() {
+    let a: LayeredObserver<(), _, _> =
+        TrackingObserver::default().layer(TrackingObserver::default());
+    let b: LayeredObserver<(), _, _> =
+        TrackingObserver::default().layer(TrackingObserver::default());
+    assert_eq!(a, b);
+}
+
+#[test]
+fn layered_observer_partial_eq_not_equal() {
+    let a: LayeredObserver<(), _, _> = TrackingObserver(1).layer(TrackingObserver(0));
+    let b: LayeredObserver<(), _, _> = TrackingObserver(0).layer(TrackingObserver(0));
+    assert_ne!(a, b);
+}
diff --git a/genetic-rs/tests/knockout.rs b/genetic-rs/tests/knockout.rs
new file mode 100644
index 0000000..db950ef
--- /dev/null
+++ b/genetic-rs/tests/knockout.rs
@@ -0,0 +1,155 @@
+//! Integration tests for [`KnockoutEliminator`], [`KnockoutWinner`], and
+//! related types.
+
+use std::cmp::Ordering;
+
+use genetic_rs::prelude::*;
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Shared test genome
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[derive(Clone, Debug)]
+struct Genome(f32);
+
+impl GenerateRandom for Genome {
+    fn gen_random(rng: &mut impl rand::Rng) -> Self {
+        Self(rng.random())
+    }
+}
+
+/// Knockout function that keeps the genome with the smaller value.
+/// `a.total_cmp(&b)` returns `Less` when a < b, mapping to `KnockoutWinner::First` (a survives).
+fn smaller_wins(a: &Genome, b: &Genome) -> KnockoutWinner {
+    a.0.total_cmp(&b.0).into()
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// KnockoutWinner helpers
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// `!First == Second` and `!Second == First`.
+#[test]
+fn knockout_winner_not_impl() {
+    assert_eq!(!KnockoutWinner::First, KnockoutWinner::Second);
+    assert_eq!(!KnockoutWinner::Second, KnockoutWinner::First);
+}
+
+/// `From<Ordering>` must map correctly.
+#[test]
+fn knockout_winner_from_ordering() {
+    assert_eq!(KnockoutWinner::from(Ordering::Less), KnockoutWinner::First);
+    assert_eq!(KnockoutWinner::from(Ordering::Equal), KnockoutWinner::First);
+    assert_eq!(
+        KnockoutWinner::from(Ordering::Greater),
+        KnockoutWinner::Second
+    );
+}
+
+/// `Into<usize>`: `First` → 0, `Second` → 1.
+#[test]
+fn knockout_winner_into_usize() {
+    assert_eq!(usize::from(KnockoutWinner::First), 0usize);
+    assert_eq!(usize::from(KnockoutWinner::Second), 1usize);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// KnockoutEliminator output-size invariants
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// With an even-sized population the output must be exactly half.
+#[test]
+fn knockout_output_half_size_even_input() {
+    let genomes: Vec<Genome> = (0..10).map(|i| Genome(i as f32)).collect();
+    let mut elim = KnockoutEliminator::new(smaller_wins, ActionIfOdd::Panic);
+    let survivors = elim.eliminate(genomes);
+    assert_eq!(survivors.len(), 5);
+}
+
+/// `ActionIfOdd::KeepSingle` with 11 genomes: 1 bypasses contest + 5 winners = 6.
+#[test]
+fn knockout_action_if_odd_keep_single() {
+    let genomes: Vec<Genome> = (0..11).map(|i| Genome(i as f32)).collect();
+    let mut elim = KnockoutEliminator::new(smaller_wins, ActionIfOdd::KeepSingle);
+    let survivors = elim.eliminate(genomes);
+    assert_eq!(survivors.len(), 6);
+}
+
+/// `ActionIfOdd::DeleteSingle` with 11 genomes: 1 is discarded + 5 winners = 5.
+#[test]
+fn knockout_action_if_odd_delete_single() {
+    let genomes: Vec<Genome> = (0..11).map(|i| Genome(i as f32)).collect();
+    let mut elim = KnockoutEliminator::new(smaller_wins, ActionIfOdd::DeleteSingle);
+    let survivors = elim.eliminate(genomes);
+    assert_eq!(survivors.len(), 5);
+}
+
+/// `ActionIfOdd::Panic` must panic when given an odd number of genomes.
+#[test]
+#[should_panic]
+fn knockout_action_if_odd_panic_panics() {
+    let genomes: Vec<Genome> = (0..11).map(|i| Genome(i as f32)).collect();
+    let mut elim = KnockoutEliminator::new(smaller_wins, ActionIfOdd::Panic);
+    let _ = elim.eliminate(genomes);
+}
+
+/// With a population of 1, [`KnockoutEliminator`] must return it unchanged.
+#[test]
+fn knockout_single_genome_returns_unchanged() {
+    let genomes = vec![Genome(42.0)];
+    let mut elim = KnockoutEliminator::new(smaller_wins, ActionIfOdd::Panic);
+    let survivors = elim.eliminate(genomes);
+    assert_eq!(survivors.len(), 1);
+    assert!((survivors[0].0 - 42.0).abs() < 1e-6);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// KnockoutEliminator correctness
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// In each pair `(a, b)` the correct genome must survive based on the knockout fn.
+///
+/// With genomes [0, 1, 2, …, 9] and `smaller_wins`:
+/// - Pairs: (0,1), (2,3), (4,5), (6,7), (8,9)
+/// - `a < b` → `Less` → `First` → a survives → even-indexed values survive.
+#[test]
+fn knockout_correct_genome_survives() {
+    let genomes: Vec<Genome> = (0..10).map(|i| Genome(i as f32)).collect();
+    let mut elim = KnockoutEliminator::new(smaller_wins, ActionIfOdd::Panic);
+    let survivors = elim.eliminate(genomes);
+
+    for g in &survivors {
+        assert!(
+            (g.0 as i32) % 2 == 0,
+            "expected even (smaller) value to survive each pair, got {}",
+            g.0
+        );
+    }
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// FitnessKnockoutFn
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// [`FitnessKnockoutFn`] must delegate to the inner fitness function and
+/// produce a consistent result based on the `afit.total_cmp(&bfit)` comparison.
+///
+/// The `From<Ordering>` mapping is:
+/// - `Less`    → `First`  (a survives — a had lower fitness)
+/// - `Greater` → `Second` (b survives — b had lower fitness)
+///
+/// This tests the actual implementation behaviour rather than any particular
+/// notion of "higher fitness wins"; it guards against regressions in the
+/// delegation and ordering logic.
+#[test]
+fn fitness_knockout_fn_delegates_to_fitness_fn() {
+    let ko_fn = FitnessKnockoutFn::new(|g: &Genome| g.0);
+    let a = Genome(1.0);
+    let b = Genome(9.0);
+
+    // 1.0.total_cmp(&9.0) = Less → First (a survives)
+    assert_eq!(ko_fn.knockout(&a, &b), KnockoutWinner::First);
+
+    // 9.0.total_cmp(&1.0) = Greater → Second (b survives)
+    assert_eq!(ko_fn.knockout(&b, &a), KnockoutWinner::Second);
+}
diff --git a/genetic-rs/tests/repopulator.rs b/genetic-rs/tests/repopulator.rs
new file mode 100644
index 0000000..27235fd
--- /dev/null
+++ b/genetic-rs/tests/repopulator.rs
@@ -0,0 +1,156 @@
+//! Integration tests for [`MitosisRepopulator`], [`CrossoverRepopulator`],
+//! and the [`FromParent`] helper.
+
+use genetic_rs::prelude::*;
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Shared test genome
+// ─────────────────────────────────────────────────────────────────────────────
+
+#[derive(Clone, Debug, PartialEq)]
+struct Genome(f32);
+
+impl GenerateRandom for Genome {
+    fn gen_random(rng: &mut impl rand::Rng) -> Self {
+        Self(rng.random())
+    }
+}
+
+impl RandomlyMutable for Genome {
+    type Context = ();
+
+    fn mutate(&mut self, _: &(), rate: f32, rng: &mut impl rand::Rng) {
+        self.0 += rng.random::<f32>() * rate;
+    }
+}
+
+impl Mitosis for Genome {
+    type Context = ();
+
+    fn divide(&self, ctx: &(), rate: f32, rng: &mut impl rand::Rng) -> Self {
+        let mut child = self.clone();
+        child.mutate(ctx, rate, rng);
+        child
+    }
+}
+
+/// Deterministic crossover: the child is exactly the average of its parents.
+/// A zero mutation rate means no random component is added.
+impl Crossover for Genome {
+    type Context = ();
+
+    fn crossover(&self, other: &Self, _: &(), _rate: f32, _rng: &mut impl rand::Rng) -> Self {
+        Self((self.0 + other.0) / 2.0)
+    }
+}
+
+fn fitness(g: &Genome) -> f32 {
+    g.0
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// MitosisRepopulator
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// The repopulator must grow the population back to the target size.
+#[test]
+fn mitosis_repopulator_fills_to_target() {
+    let mut rng = rand::rng();
+    let mut genomes: Vec<Genome> = Vec::gen_random(&mut rng, 5);
+    MitosisRepopulator::new(0.1, ()).repopulate(&mut genomes, 20);
+    assert_eq!(genomes.len(), 20);
+}
+
+/// If the population is already at the target, it must not grow further.
+#[test]
+fn mitosis_repopulator_no_op_when_at_target() {
+    let mut rng = rand::rng();
+    let mut genomes: Vec<Genome> = Vec::gen_random(&mut rng, 10);
+    MitosisRepopulator::new(0.1, ()).repopulate(&mut genomes, 10);
+    assert_eq!(genomes.len(), 10);
+}
+
+/// At mutation rate 0.0 each child must be an exact clone of its parent.
+#[test]
+fn mitosis_zero_mutation_rate_produces_clone() {
+    let parent = Genome(42.0);
+    let mut rng = rand::rng();
+    let child = parent.divide(&(), 0.0, &mut rng);
+    assert_eq!(
+        child, parent,
+        "child at rate 0.0 must be identical to parent"
+    );
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// CrossoverRepopulator
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// The crossover repopulator must grow the population back to the target size.
+#[test]
+fn crossover_repopulator_fills_to_target() {
+    let mut rng = rand::rng();
+    let mut genomes: Vec<Genome> = Vec::gen_random(&mut rng, 5);
+    CrossoverRepopulator::new(0.1, ()).repopulate(&mut genomes, 20);
+    assert_eq!(genomes.len(), 20);
+}
+
+/// If the population is already at the target, it must not grow further.
+#[test]
+fn crossover_repopulator_no_op_when_at_target() {
+    let mut rng = rand::rng();
+    let mut genomes: Vec<Genome> = Vec::gen_random(&mut rng, 10);
+    CrossoverRepopulator::new(0.1, ()).repopulate(&mut genomes, 10);
+    assert_eq!(genomes.len(), 10);
+}
+
+/// Crossover child must equal the average of its parents when mutation rate is 0.
+#[test]
+fn crossover_child_is_average_of_parents() {
+    let p1 = Genome(2.0);
+    let p2 = Genome(4.0);
+    let mut rng = rand::rng();
+    let child = p1.crossover(&p2, &(), 0.0, &mut rng);
+    assert!(
+        (child.0 - 3.0).abs() < 1e-6,
+        "expected child value 3.0, got {}",
+        child.0
+    );
+}
+
+/// Population size stays constant when `GeneticSim` uses crossover over many generations.
+#[test]
+fn population_size_constant_with_crossover_sim() {
+    let mut rng = rand::rng();
+    let initial_size = 30;
+    let mut sim = GeneticSim::new(
+        Vec::<Genome>::gen_random(&mut rng, initial_size),
+        FitnessEliminator::new_without_observer(fitness),
+        CrossoverRepopulator::new(0.05, ()),
+    );
+    sim.perform_generations(50);
+    assert_eq!(sim.genomes.len(), initial_size);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// FromParent
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// [`Vec::from_parent`] must produce exactly the requested number of genomes.
+#[test]
+fn from_parent_creates_correct_count() {
+    let parent = Genome(1.0);
+    let population = Vec::<Genome>::from_parent(parent, 15, (), 0.0);
+    assert_eq!(population.len(), 15);
+}
+
+/// The original parent must be the first element of the returned population.
+#[test]
+fn from_parent_includes_original_parent() {
+    let parent = Genome(7.0);
+    let population = Vec::<Genome>::from_parent(parent.clone(), 5, (), 0.0);
+    assert_eq!(
+        population[0], parent,
+        "the first element must be the original parent"
+    );
+}
diff --git a/genetic-rs/tests/simulation.rs b/genetic-rs/tests/simulation.rs
new file mode 100644
index 0000000..12ed2c6
--- /dev/null
+++ b/genetic-rs/tests/simulation.rs
@@ -0,0 +1,233 @@
+//! Integration tests for [`GeneticSim`] core invariants, [`FitnessEliminator`] behaviour,
+//! and helper types such as [`FitnessObserver`].
+
+use genetic_rs::prelude::*;
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Shared test genome
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// A simple genome whose fitness is just its value.
+#[derive(Clone, Debug)]
+struct Genome(f32);
+
+impl GenerateRandom for Genome {
+    fn gen_random(rng: &mut impl rand::Rng) -> Self {
+        Self(rng.random())
+    }
+}
+
+impl RandomlyMutable for Genome {
+    type Context = ();
+
+    fn mutate(&mut self, _: &(), rate: f32, rng: &mut impl rand::Rng) {
+        self.0 += rng.random::<f32>() * rate;
+    }
+}
+
+impl Mitosis for Genome {
+    type Context = ();
+
+    fn divide(&self, ctx: &(), rate: f32, rng: &mut impl rand::Rng) -> Self {
+        let mut child = self.clone();
+        child.mutate(ctx, rate, rng);
+        child
+    }
+}
+
+fn fitness(g: &Genome) -> f32 {
+    g.0
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// GeneticSim population-size invariants
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// The population size must be identical before and after a single generation.
+#[test]
+fn population_size_preserved_after_next_generation() {
+    let mut rng = rand::rng();
+    let initial_size = 20;
+    let mut sim = GeneticSim::new(
+        Vec::<Genome>::gen_random(&mut rng, initial_size),
+        FitnessEliminator::new_without_observer(fitness),
+        MitosisRepopulator::new(0.1, ()),
+    );
+    sim.next_generation();
+    assert_eq!(sim.genomes.len(), initial_size);
+}
+
+/// The population size must remain constant over many generations.
+#[test]
+fn population_size_preserved_over_many_generations() {
+    let mut rng = rand::rng();
+    let initial_size = 50;
+    let mut sim = GeneticSim::new(
+        Vec::<Genome>::gen_random(&mut rng, initial_size),
+        FitnessEliminator::new_without_observer(fitness),
+        MitosisRepopulator::new(0.1, ()),
+    );
+    sim.perform_generations(100);
+    assert_eq!(sim.genomes.len(), initial_size);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// GenerateRandom / GenerateRandomCollection
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// [`Vec::gen_random`] must produce exactly the requested number of genomes.
+#[test]
+fn gen_random_collection_produces_correct_size() {
+    let mut rng = rand::rng();
+    for size in [1usize, 10, 100] {
+        let population: Vec<Genome> = Vec::gen_random(&mut rng, size);
+        assert_eq!(
+            population.len(),
+            size,
+            "expected {size} genomes, got {}",
+            population.len()
+        );
+    }
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// FitnessEliminator elimination behaviour
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// With the default threshold (0.5) and 10 genomes, exactly 6 must survive and
+/// all survivors must be in the top half by fitness.
+#[test]
+fn fitness_eliminator_keeps_top_by_fitness() {
+    // Genomes with distinct, known fitness values 0.0 … 9.0.
+    let genomes: Vec<Genome> = (0..10).map(|i| Genome(i as f32)).collect();
+    let mut eliminator = FitnessEliminator::new_without_observer(fitness);
+    let survivors = eliminator.eliminate(genomes);
+
+    // threshold=0.5, n=10 → floor(10*0.5) + 1 = 6 survivors.
+    assert_eq!(
+        survivors.len(),
+        6,
+        "expected 6 survivors with default threshold"
+    );
+
+    // Every survivor must be in the top half (fitness ≥ 4.0).
+    for g in &survivors {
+        assert!(
+            g.0 >= 4.0,
+            "genome with fitness {} survived but is below the top-half threshold",
+            g.0,
+        );
+    }
+}
+
+/// The genome with the highest fitness must always survive elimination.
+#[test]
+fn fitness_eliminator_highest_fitness_always_survives() {
+    let genomes: Vec<Genome> = (0..20).map(|i| Genome(i as f32)).collect();
+    let mut eliminator = FitnessEliminator::new_without_observer(fitness);
+    let survivors = eliminator.eliminate(genomes);
+    assert!(
+        survivors.iter().any(|g| g.0 == 19.0),
+        "the highest-fitness genome (19.0) must always survive"
+    );
+}
+
+/// The genome with the lowest fitness must always be eliminated.
+#[test]
+fn fitness_eliminator_lowest_fitness_never_survives() {
+    let genomes: Vec<Genome> = (0..20).map(|i| Genome(i as f32)).collect();
+    let mut eliminator = FitnessEliminator::new_without_observer(fitness);
+    let survivors = eliminator.eliminate(genomes);
+    assert!(
+        !survivors.iter().any(|g| g.0 == 0.0),
+        "the lowest-fitness genome (0.0) must always be eliminated"
+    );
+}
+
+/// [`FitnessEliminator::calculate_and_sort`] must return genomes sorted in
+/// descending order of fitness.
+#[test]
+fn fitness_eliminator_sorts_descending() {
+    let genomes = vec![Genome(3.0), Genome(1.0), Genome(4.0), Genome(2.0)];
+    let eliminator = FitnessEliminator::new_without_observer(fitness);
+    let sorted = eliminator.calculate_and_sort(genomes);
+
+    for window in sorted.windows(2) {
+        assert!(
+            window[0].1 >= window[1].1,
+            "fitness values are not in descending order: {} before {}",
+            window[0].1,
+            window[1].1,
+        );
+    }
+}
+
+/// A custom threshold controls the exact fraction of the population kept.
+#[test]
+fn fitness_eliminator_custom_threshold() {
+    // 10 genomes, threshold=0.3 → floor(10*0.3) + 1 = 4 survivors.
+    let genomes: Vec<Genome> = (0..10).map(|i| Genome(i as f32)).collect();
+    let mut eliminator = FitnessEliminator::new(fitness, 0.3, ());
+    let survivors = eliminator.eliminate(genomes);
+    assert_eq!(
+        survivors.len(),
+        4,
+        "expected 4 survivors with threshold=0.3"
+    );
+}
+
+/// [`FitnessEliminator::new`] must panic when the threshold is outside [0, 1].
+#[test]
+#[should_panic]
+fn fitness_eliminator_invalid_threshold_panics() {
+    FitnessEliminator::new(fitness, 1.5_f32, ());
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// FitnessEliminator builder
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// The builder pattern must produce an eliminator with the specified settings.
+#[test]
+fn fitness_eliminator_builder() {
+    let mut rng = rand::rng();
+    let mut eliminator: FitnessEliminator<_, Genome, ()> = FitnessEliminator::builder()
+        .fitness_fn(fitness)
+        .threshold(0.4)
+        .build();
+
+    let genomes: Vec<Genome> = Vec::gen_random(&mut rng, 10);
+    let _ = eliminator.eliminate(genomes);
+    assert!((eliminator.threshold - 0.4).abs() < 1e-6);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// FitnessObserver
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// A counting observer used below.
+struct Counter(usize);
+
+impl FitnessObserver<Genome> for Counter {
+    fn observe(&mut self, _: &[(Genome, f32)]) {
+        self.0 += 1;
+    }
+}
+
+/// The observer must be called exactly once per generation.
+#[test]
+fn observer_called_once_per_generation() {
+    let mut rng = rand::rng();
+    let eliminator = FitnessEliminator::new(fitness, 0.5, Counter(0));
+    let mut sim = GeneticSim::new(
+        Vec::<Genome>::gen_random(&mut rng, 10),
+        eliminator,
+        MitosisRepopulator::new(0.0, ()),
+    );
+
+    sim.perform_generations(7);
+    assert_eq!(
+        sim.eliminator.observer.0, 7,
+        "observer must be called once per generation"
+    );
+}
diff --git a/genetic-rs/tests/speciation.rs b/genetic-rs/tests/speciation.rs
new file mode 100644
index 0000000..751a7ca
--- /dev/null
+++ b/genetic-rs/tests/speciation.rs
@@ -0,0 +1,267 @@
+//! Integration tests for [`SpeciatedPopulation`], [`SpeciatedFitnessEliminator`],
+//! and [`SpeciatedCrossoverRepopulator`].
+
+use genetic_rs::prelude::*;
+use genetic_rs::speciation::SpeciatedPopulation;
+
+// ─────────────────────────────────────────────────────────────────────────────
+// Shared test genome
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// A genome whose species membership is fully determined by its integer `class`.
+/// Two genomes with the same class have divergence 0.0; different classes give 1.0.
+#[derive(Clone, Debug, PartialEq)]
+struct Genome {
+    class: i32,
+    val: f32,
+}
+
+impl Speciated for Genome {
+    type Context = ();
+
+    fn divergence(&self, other: &Self, _: &()) -> f32 {
+        if self.class == other.class {
+            0.0
+        } else {
+            1.0
+        }
+    }
+}
+
+impl GenerateRandom for Genome {
+    fn gen_random(rng: &mut impl rand::Rng) -> Self {
+        Self {
+            class: rng.random_range(0..3),
+            val: rng.random(),
+        }
+    }
+}
+
+impl RandomlyMutable for Genome {
+    type Context = ();
+
+    fn mutate(&mut self, _: &(), rate: f32, rng: &mut impl rand::Rng) {
+        self.val += rng.random::<f32>() * rate;
+    }
+}
+
+impl Crossover for Genome {
+    type Context = ();
+
+    fn crossover(&self, other: &Self, _: &(), _rate: f32, _rng: &mut impl rand::Rng) -> Self {
+        Self {
+            class: self.class,
+            val: (self.val + other.val) / 2.0,
+        }
+    }
+}
+
+fn fitness(g: &Genome) -> f32 {
+    g.val
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// SpeciatedPopulation — species grouping
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// All identical genomes must end up in a single species.
+#[test]
+fn identical_genomes_in_same_species() {
+    let genomes: Vec<Genome> = (0..5).map(|_| Genome { class: 0, val: 0.0 }).collect();
+    let pop = SpeciatedPopulation::from_genomes(&genomes, 0.5, &());
+    assert_eq!(pop.species().len(), 1, "expected a single species");
+    assert_eq!(
+        pop.species()[0].len(),
+        5,
+        "all genomes must belong to the single species"
+    );
+}
+
+/// Genomes of four distinct classes must form four separate species.
+#[test]
+fn different_class_genomes_in_different_species() {
+    let genomes: Vec<Genome> = (0..4).map(|i| Genome { class: i, val: 0.0 }).collect();
+    let pop = SpeciatedPopulation::from_genomes(&genomes, 0.5, &());
+    assert_eq!(pop.species().len(), 4);
+}
+
+/// With a threshold > 1.0 (larger than the max divergence) all genomes,
+/// regardless of class, must end up in the same species.
+#[test]
+fn high_threshold_groups_all_genomes() {
+    let genomes: Vec<Genome> = (0..6)
+        .map(|i| Genome {
+            class: i % 3,
+            val: 0.0,
+        })
+        .collect();
+    // Divergence is at most 1.0; with threshold 1.5 everything is "close enough".
+    let pop = SpeciatedPopulation::from_genomes(&genomes, 1.5, &());
+    assert_eq!(pop.species().len(), 1, "all genomes must be in one species");
+}
+
+/// Every genome index must appear in exactly one species.
+#[test]
+fn every_genome_index_appears_exactly_once() {
+    let n = 9usize;
+    let genomes: Vec<Genome> = (0..n)
+        .map(|i| Genome {
+            class: (i % 3) as i32,
+            val: i as f32,
+        })
+        .collect();
+    let pop = SpeciatedPopulation::from_genomes(&genomes, 0.5, &());
+
+    let mut seen = vec![false; n];
+    for species in pop.species() {
+        for &idx in species {
+            assert!(
+                !seen[idx],
+                "genome index {idx} appears in more than one species"
+            );
+            seen[idx] = true;
+        }
+    }
+    assert!(
+        seen.iter().all(|&v| v),
+        "not every genome index appeared in a species"
+    );
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// SpeciatedPopulation — insert_genome
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// Inserting a genome from a new class must create a new species.
+#[test]
+fn insert_genome_creates_new_species_for_novel_genome() {
+    let genomes = vec![
+        Genome { class: 0, val: 0.0 },
+        Genome { class: 1, val: 0.0 }, // divergence 1.0 > threshold 0.5 → new species
+    ];
+    let mut pop = SpeciatedPopulation::new(0.5);
+    pop.insert_genome(0, &genomes, &());
+    let created_new = pop.insert_genome(1, &genomes, &());
+    assert!(created_new, "expected a new species to be created");
+    assert_eq!(pop.species().len(), 2);
+}
+
+/// Inserting a genome from an existing class must join that species.
+#[test]
+fn insert_genome_joins_existing_species_for_similar_genome() {
+    let genomes = vec![
+        Genome { class: 0, val: 0.0 },
+        Genome { class: 0, val: 1.0 }, // divergence 0.0 < threshold 0.5 → joins species
+    ];
+    let mut pop = SpeciatedPopulation::new(0.5);
+    pop.insert_genome(0, &genomes, &());
+    let created_new = pop.insert_genome(1, &genomes, &());
+    assert!(
+        !created_new,
+        "must not create a new species for a similar genome"
+    );
+    assert_eq!(pop.species().len(), 1);
+    assert_eq!(pop.species()[0].len(), 2);
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// SpeciatedPopulation — round-robin iterator
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// The round-robin iterator must visit every genome index within one full cycle.
+/// With equal-sized species the cycle length equals the total number of genomes.
+#[test]
+fn round_robin_covers_all_genomes() {
+    // 3 classes × 2 genomes each = 6 genomes, 3 species of equal size.
+    let genomes: Vec<Genome> = (0..6)
+        .map(|i| Genome {
+            class: (i % 3) as i32,
+            val: i as f32,
+        })
+        .collect();
+    let pop = SpeciatedPopulation::from_genomes(&genomes, 0.5, &());
+
+    let mut seen = vec![false; 6];
+    for idx in pop.round_robin().take(6) {
+        seen[idx] = true;
+    }
+    assert!(
+        seen.iter().all(|&v| v),
+        "round_robin must cover all genome indices in one cycle"
+    );
+}
+
+/// The enumerate variant must yield matching (species_index, genome_index) pairs.
+#[test]
+fn round_robin_enumerate_species_index_is_valid() {
+    let genomes: Vec<Genome> = (0..6)
+        .map(|i| Genome {
+            class: (i % 3) as i32,
+            val: i as f32,
+        })
+        .collect();
+    let pop = SpeciatedPopulation::from_genomes(&genomes, 0.5, &());
+
+    for (species_i, genome_i) in pop.round_robin_enumerate().take(12) {
+        assert!(
+            pop.species()[species_i].contains(&genome_i),
+            "genome index {genome_i} is not in species {species_i}"
+        );
+    }
+}
+
+// ─────────────────────────────────────────────────────────────────────────────
+// SpeciatedFitnessEliminator — fitness scaling and population-size invariants
+// ─────────────────────────────────────────────────────────────────────────────
+
+/// The population size must remain constant when using [`SpeciatedFitnessEliminator`]
+/// together with [`SpeciatedCrossoverRepopulator`].
+#[test]
+fn speciated_sim_population_size_preserved() {
+    let mut rng = rand::rng();
+    let initial_size = 30;
+
+    let eliminator = SpeciatedFitnessEliminator::new(fitness, 0.5, 0.5, (), ());
+    let repopulator = SpeciatedCrossoverRepopulator::new(
+        0.1,
+        0.5,
+        ActionIfIsolated::CrossoverSimilarSpecies,
+        (),
+        (),
+    );
+
+    let mut sim = GeneticSim::new(
+        Vec::<Genome>::gen_random(&mut rng, initial_size),
+        eliminator,
+        repopulator,
+    );
+    sim.perform_generations(20);
+    assert_eq!(sim.genomes.len(), initial_size);
+}
+
+/// Speciation must give fitness bonuses to rare/isolated genomes:
+/// an isolated genome with *lower* raw fitness can out-compete a large species
+/// with higher raw fitness because the large species has its fitness divided by
+/// its size.
+///
+/// Setup:
+/// - 4 genomes of class 0 with val = 1.0  → adjusted fitness = 1.0 / 4 = 0.25
+/// - 1 genome  of class 1 with val = 0.5  → adjusted fitness = 0.5 / 1 = 0.5
+///
+/// With threshold = 0.5 and 5 genomes, 3 survive.
+/// Sorted adjusted fitness: [0.5, 0.25, 0.25, 0.25, 0.25]
+/// → the class-1 genome (raw 0.5) must survive despite the lower raw fitness.
+#[test]
+fn speciation_protects_rare_species() {
+    let mut class0_genomes: Vec<Genome> = (0..4).map(|_| Genome { class: 0, val: 1.0 }).collect();
+    let rare = Genome { class: 1, val: 0.5 };
+    class0_genomes.push(rare.clone());
+
+    let mut eliminator = SpeciatedFitnessEliminator::new(fitness, 0.5, 0.5, (), ());
+    let survivors = eliminator.eliminate(class0_genomes);
+
+    assert!(
+        survivors.iter().any(|g| g == &rare),
+        "the rare species genome must survive despite lower raw fitness"
+    );
+}