Repository: Helmholtz-AI-Energy/propulate (language: python)

Issue:
MPI communicator handling: Parallelism on the level of single individuals
Enable passing a specific MPI communicator to the evolution routine, defaulting to MPI_COMM_WORLD. Duplicate the communicator with MPI_Comm_dup beforehand.


Discussion / hints:
This is related to enabling parallelized evaluations, i.e., each evaluation of an individual is, e.g., a data-parallel training of a neural network using specific hyperparameters to be optimized. Requires one **"optimization" communicator** handling the optimization routine, with each rank corresponding to one worker evaluating an individual. Each of those optimization comm ranks (or workers) has another **"individual" communicator** attached which handles the parallelized evaluation of the individual itself. To communicate information between those two communicators, a worker's optimization comm rank must also be part of its individual comm as an "interface". This can be accomplished with splits or groups.

Migration and pollination between islands requires a communicator containing all workers from all islands, similar to what `COMM_WORLD` is with each individual corresponding to only one rank.

Fix patch:
diff --git a/README.md b/README.md
index 88ce0dc3..2b0ca68e 100644
--- a/README.md
+++ b/README.md
@@ -18,45 +18,204 @@
 
 ## What `Propulate` can do for you
 
-`Propulate` is an HPC-tailored software for solving optimization problems in parallel. It is openly accessible and easy to use. Compared to a widely used competitor, `Propulate` is consistently faster - at least an order of magnitude for a set of typical benchmarks - and in some cases even more accurate.
-
-Inspired by biology, `Propulate` borrows mechanisms from biological evolution, such as selection, recombination, and mutation. Evolution begins with a population of solution candidates, each with randomly initialized genes. It is an iterative "survival of the fittest" process where the population at each iteration can be viewed as a generation. For each generation, the fitness of each candidate in the population is evaluated. The genes of the fittest candidates are incorporated in the next generation.
-
-Like in nature, `Propulate` does not wait for all compute units to finish the evaluation of the current generation. Instead, the compute units communicate the currently available information and use that to breed the next candidate immediately. This avoids waiting idly for other units and thus a load imbalance.
-Each unit is responsible for evaluating a single candidate. The result is a fitness level corresponding with that candidate’s genes, allowing us to compare and rank all candidates. This information is sent to other compute units as soon as it becomes available.
-When a unit is finished evaluating a candidate and communicating the resulting fitness, it breeds the candidate for the next generation using the fitness values of all candidates it evaluated and received from other units so far. 
+`Propulate` is an HPC-tailored software for solving optimization problems in parallel. It is openly accessible and easy 
+to use. Compared to a widely used competitor, `Propulate` is consistently faster - at least an order of magnitude for a 
+set of typical benchmarks - and in some cases even more accurate.
+
+Inspired by biology, `Propulate` borrows mechanisms from biological evolution, such as selection, recombination, and 
+mutation. Evolution begins with a population of solution candidates, each with randomly initialized genes. It is an 
+iterative "survival of the fittest" process where the population at each iteration can be viewed as a generation. For 
+each generation, the fitness of each candidate in the population is evaluated. The genes of the fittest candidates are 
+incorporated in the next generation.
+
+Like in nature, `Propulate` does not wait for all compute units to finish the evaluation of the current generation. 
+Instead, the compute units communicate the currently available information and use that to breed the next candidate 
+immediately. This avoids waiting idly for other units and thus a load imbalance.
+Each unit is responsible for evaluating a single candidate. The result is a fitness level corresponding with that 
+candidate’s genes, allowing us to compare and rank all candidates. This information is sent to other compute units as 
+soon as it becomes available.
+When a unit is finished evaluating a candidate and communicating the resulting fitness, it breeds the candidate for the 
+next generation using the fitness values of all candidates it evaluated and received from other units so far. 
 
 `Propulate` can be used for hyperparameter optimization and neural architecture search. 
-It was already successfully applied in several accepted scientific publications. Applications include grid load forecasting, remote sensing, and structural molecular biology.
+It was already successfully applied in several accepted scientific publications. Applications include grid load 
+forecasting, remote sensing, and structural molecular biology.
 
 ## In more technical terms
 
-``Propulate`` is a massively parallel evolutionary hyperparameter optimizer based on the island model with asynchronous propagation of populations and asynchronous migration.
-In contrast to classical GAs, ``Propulate`` maintains a continuous population of already evaluated individuals with a softened notion of the typically strictly separated, discrete generations.
+``Propulate`` is a massively parallel evolutionary hyperparameter optimizer based on the island model with asynchronous 
+propagation of populations and asynchronous migration.
+In contrast to classical GAs, ``Propulate`` maintains a continuous population of already evaluated individuals with a 
+softened notion of the typically strictly separated, discrete generations.
 Our contributions include:
 - A novel parallel genetic algorithm based on a fully asynchronized island model with independently processing workers.
-- Massive parallelism by asynchronous propagation of continuous populations and migration via efficient communication using the message passing interface.
+- Massive parallelism by asynchronous propagation of continuous populations and migration via efficient communication 
+- using the message passing interface.
 - Optimized use efficiency of parallel hardware by minimizing idle times in distributed computing environments.
 
 To be more efficient, the generations are less well separated than they usually are in evolutionary algorithms.
-New individuals are generated from a pool of currently active, already evaluated individuals that may be from any generation.
+New individuals are generated from a pool of currently active, already evaluated individuals that may be from any 
+generation.
 Individuals may be removed from the breeding population based on different criteria.
 
 You can find the corresponding publication [here](https://doi.org/10.1007/978-3-031-32041-5_6):  
->Taubert, O. *et al.* (2023). Massively Parallel Genetic Optimization Through Asynchronous Propagation of Populations. In: Bhatele, A., Hammond, J., Baboulin, M., Kruse, C. (eds) High Performance Computing. ISC High Performance 2023. Lecture Notes in Computer Science, vol 13948. Springer, Cham. [doi.org/10.1007/978-3-031-32041-5_6](https://doi.org/10.1007/978-3-031-32041-5_6)
+> Taubert, O. *et al.* (2023). Massively Parallel Genetic Optimization Through Asynchronous Propagation of Populations. 
+> In: Bhatele, A., Hammond, J., Baboulin, M., Kruse, C. (eds) High Performance Computing. ISC High Performance 2023. 
+> Lecture Notes in Computer Science, vol 13948. Springer, Cham. 
+> [doi.org/10.1007/978-3-031-32041-5_6](https://doi.org/10.1007/978-3-031-32041-5_6)
 
 ## Documentation
 
-Check out the full documentation at [https://propulate.readthedocs.io/](https://propulate.readthedocs.io/) :rocket:! Here you can 
-find installation instructions, tutorials, theoretical background, and API references.
+Check out the full documentation at [https://propulate.readthedocs.io/](https://propulate.readthedocs.io/) :rocket:! Here you can find installation 
+instructions, tutorials, theoretical background, and API references.
+
+**:point_right: If you have any questions or run into any challenges while using `Propulate`, don't hesitate to post an 
+[issue](https://github.com/Helmholtz-AI-Energy/propulate/issues) :bookmark:, reach out via [GitHub 
+discussions](https://github.com/Helmholtz-AI-Energy/propulate/discussions) :octocat:, or contact us directly via e-mail 
+:email: to [propulate@lists.kit.edu](mailto:propulate@lists.kit.edu).** 
 
 ## Installation
 
 You can install the latest stable release from PyPI: ``pip install propulate``  
 If you need the latest updates, you can also install ``Propulate`` directly from the master branch at you own risk.
 Pull and run ``pip install -e .`` or ``python setup.py develop``.  
-``Propulate`` depends on ``mpi4py`` and requires an MPI implementation under the hood. Currently, it is only tested with 
-OpenMPI.
+``Propulate`` depends on [``mpi4py``](https://mpi4py.readthedocs.io/en/stable/) and requires an MPI implementation under 
+the hood. Currently, it is only tested with [OpenMPI](https://www.open-mpi.org/).
+
+## Quickstart
+*Below, you can find a quick recipe for how to use `Propulate` in general. Check out the official 
+[ReadTheDocs](https://propulate.readthedocs.io/en/latest/tut_propulator.html) documentation for more detailed tutorials 
+and explanations.*
+
+Let's minimize the sphere function $f_\text{sphere}\left(x,y\right)=x^2 +y^2$ with `Propulate` as a quick example. The 
+minimum is at $\left(x, y\right)=\left(0,0\right)$ at the orange star.
+![](./docs/images/sphere.png)
+First, we need to define the key ingredients that define our optimization problem:
+- The **search space** of the parameters to be optimized as a `Python` dictionary. `Propulate` can handle three different
+  parameter types:
+    - A tuple of `float` for a continuous parameter, e.g., `{"learning_rate": (0.0001, 0.01)}`
+    - A tuple of `int` for an ordinal parameter, e.g., `{"conv_layers": (2, 10)}`
+    - A tuple of `str` for a categorical parameter, e.g., `{"activation": ("relu", "sigmoid", "tanh")}`
+  
+  Thus, an exemplary search space might look like this:
+  ```python
+  search_space = {
+      "learning_rate": (0.0001, 0.01),  # Search a continuous space between 0.0001 and 0.01.
+      "num_layers": (2, 10),  # Search the integer space between 2 and 10 (inclusive).
+      "activation": ("relu", "sigmoid", "tanh"),  # Search the categorical space with the specified possibilities.
+  }
+  ```
+
+  The sphere function has two continuous parameters, $x$ and $y$, and we consider $x,y\in\left[-5.12,5.12\right]$. The 
+  search space in our example thus looks like this:
+  ```python
+  limits = {
+      "x": (-5.12, 5.12),
+      "y": (-5.12, 5.12)
+  } 
+  ```
+- The **loss function**. This is the function we want to minimize in order to find the best parameters. It can be any 
+  `Python` function that
+  - takes a set of parameters as a `Python` dictionary as an input.
+  - returns a scalar loss value that determines how good the tested parameter set is.
+  
+  In this example, the loss function whose minimum we want to find is the sphere function:
+  ```python
+  def sphere(params: Dict[str, float]) -> float:
+    """
+    Sphere function: continuous, convex, separable, differentiable, unimodal
+
+    Input domain: -5.12 <= x, y <= 5.12
+    Global minimum 0 at (x, y) = (0, 0)
+
+    Parameters
+    ----------
+    params: Dict[str, float]
+        The function parameters.
+  
+    Returns
+    -------
+    float
+        The function value.
+    """
+    return numpy.sum(numpy.array(list(params.values())) ** 2).item()
+  ```
+Next, we need to define the **evolutionary operator** or propagator that we want to use to breed new individuals during the 
+optimization process. `Propulate` provides a reasonable default propagator via a utility function:
+```python
+# Set up logger for Propulate optimization.
+propulate.set_logger_config()
+# Set up separate random number generator for Propulate optimization. DO NOT USE SOMEWHERE ELSE!
+rng = random.Random(
+    <your-random-seed> + mpi4py.MPI.COMM_WORLD.rank
+)
+# Set up evolutionary operator.
+propagator = propulate.get_default_propagator(
+    pop_size=config.pop_size,  # The breeding population size
+    limits=limits,  # The search-space limits
+    rng=rng,  # Random number generator
+)
+```
+We also need to set up the asynchronous parallel evolutionary **optimizer**, that is a so-called ``Propulator`` instance:
+```python
+# Set up Propulator performing actual optimization.
+propulator = propulate.Propulator(
+    loss_fn=sphere,
+    propagator=propagator,
+    rng=rng,
+    generations=config.generations,
+    checkpoint_path=config.checkpoint,
+)
+```
+Now we can run the actual optimization. Overall, ``generations * mpi4py.MPI.COMM_WORLD.size`` evaluations will be 
+performed:
+```python
+# Run optimization and print summary of results.
+propulator.propulate()
+propulator.summarize()
+```
+The output should look something like this:
+```text
+#################################################
+# PROPULATE: Parallel Propagator of Populations #
+#################################################
+
+[2024-03-12 14:37:01,374][propulate.propulator][INFO] - No valid checkpoint file given. Initializing population randomly...
+[2024-03-12 14:37:01,374][propulate.propulator][INFO] - Island 0 has 4 workers.
+[2024-03-12 14:37:01,374][propulate.propulator][INFO] - Island 0 Worker 0: In generation 0...
+[2024-03-12 14:37:01,374][propulate.propulator][INFO] - Island 0 Worker 3: In generation 0...
+[2024-03-12 14:37:01,374][propulate.propulator][INFO] - Island 0 Worker 2: In generation 0...
+[2024-03-12 14:37:01,374][propulate.propulator][INFO] - Island 0 Worker 1: In generation 0...
+[2024-03-12 14:37:01,377][propulate.propulator][INFO] - Island 0 Worker 3: In generation 10...
+[2024-03-12 14:37:01,377][propulate.propulator][INFO] - Island 0 Worker 1: In generation 10...
+[2024-03-12 14:37:01,378][propulate.propulator][INFO] - Island 0 Worker 0: In generation 10...
+[2024-03-12 14:37:01,378][propulate.propulator][INFO] - Island 0 Worker 2: In generation 10...
+
+...
+[2024-03-12 14:37:02,197][propulate.propulator][INFO] - Island 0 Worker 1: In generation 960...
+[2024-03-12 14:37:02,206][propulate.propulator][INFO] - Island 0 Worker 2: In generation 990...
+[2024-03-12 14:37:02,206][propulate.propulator][INFO] - Island 0 Worker 1: In generation 970...
+[2024-03-12 14:37:02,215][propulate.propulator][INFO] - Island 0 Worker 1: In generation 980...
+[2024-03-12 14:37:02,224][propulate.propulator][INFO] - Island 0 Worker 1: In generation 990...
+[2024-03-12 14:37:02,232][propulate.propulator][INFO] - OPTIMIZATION DONE.
+NEXT: Final checks for incoming messages...
+[2024-03-12 14:37:02,244][propulate.propulator][INFO] - ###########
+# SUMMARY #
+###########
+Number of currently active individuals is 4000.
+Expected overall number of evaluations is 4000.
+[2024-03-12 14:37:03,703][propulate.propulator][INFO] - Top 1 result(s) on island 0:
+(1): [{'a': '2.91E-3', 'b': '-3.05E-3'}, loss 1.78E-5, island 0, worker 0, generation 956]
+```
+### Let's get your hands dirty
+Do the following to run the [example script](https://github.com/Helmholtz-AI-Energy/propulate/blob/master/tutorials/propulator_example.py):
+
+- Make sure you have a working MPI installation on your machine.
+- If you have not already done this, create a fresh virtual environment with ``Python``: ``$ python3 -m venv best-venv-ever``
+- Activate it: ``$ source best-venv-ever/bin/activate``
+- Upgrade ``pip``: ``$ pip install --upgrade pip``
+- Install ``Propulate``: ``$ pip install propulate``
+- Run the example script ``propulator_example.py``: ``$ mpirun --use-hwthread-cpus python propulator_example.py``
 
 ## Acknowledgments
 *This work is supported by the Helmholtz AI platform grant.*
diff --git a/coverage.svg b/coverage.svg
index 9d027c7d..ae1c3500 100644
--- a/coverage.svg
+++ b/coverage.svg
@@ -15,7 +15,7 @@
     <g fill="#fff" text-anchor="middle" font-family="DejaVu Sans,Verdana,Geneva,sans-serif" font-size="11">
         <text x="31.5" y="15" fill="#010101" fill-opacity=".3">coverage</text>
         <text x="31.5" y="14">coverage</text>
-        <text x="80" y="15" fill="#010101" fill-opacity=".3">49%</text>
-        <text x="80" y="14">49%</text>
+        <text x="80" y="15" fill="#010101" fill-opacity=".3">54%</text>
+        <text x="80" y="14">54%</text>
     </g>
 </svg>
diff --git a/docs/cmaes.rst b/docs/cmaes.rst
index 1d651201..b6644c52 100644
--- a/docs/cmaes.rst
+++ b/docs/cmaes.rst
@@ -4,5 +4,4 @@ Covariance Matrix Adaption Evolution Strategy
 =============================================
 
 .. note::
-   Covariance Matrix Adaption Evolution Strategy (CMA-ES) is under active development and not yet provided in the latest
-   version of ``Propulate``. It will be available soon.
\ No newline at end of file
+   The documentation for Covariance Matrix Adaption Evolution Strategy (CMA-ES) will be available soon.
\ No newline at end of file
diff --git a/docs/conf.py b/docs/conf.py
index d6f84618..08485e5d 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -36,7 +36,7 @@
 ]
 
 autoapi_dirs = ["../propulate"]
-autoapi_python_class_content = "both"
+autoapi_python_class_content = "init"
 
 # Add any paths that contain templates here, relative to this directory.
 templates_path = ["_templates"]
diff --git a/docs/index.rst b/docs/index.rst
index 9291e5ce..21186afd 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -103,4 +103,4 @@ Indices and Tables
 
 * :ref:`genindex`
 * :ref:`modindex`
-* :ref:`search`
\ No newline at end of file
+* :ref:`search`
diff --git a/docs/pso.rst b/docs/pso.rst
index 0fd69c2f..c393557c 100644
--- a/docs/pso.rst
+++ b/docs/pso.rst
@@ -4,5 +4,4 @@ Particle Swarm Optimization
 ===========================
 
 .. note::
-   Particle swarm optimization (PSO) is under active development and not yet provided in the latest version of
-   ``Propulate``. It will be available soon.
\ No newline at end of file
+   The documentation for particle swarm optimization (PSO) will be available soon.
\ No newline at end of file
diff --git a/docs/quickstart.rst b/docs/quickstart.rst
index 943970e5..2795f526 100644
--- a/docs/quickstart.rst
+++ b/docs/quickstart.rst
@@ -14,7 +14,51 @@ separate ``Python3`` virtual environment:
     $ pip install propulate
 
 This will install the ``Propulate`` optimizer functionality along with some basic example scripts to get you started.
-You can find those scripts in the Github repository's `subfolder`_ ``scripts/``.
+You can find those scripts in the Github repository's `subfolder`_ ``tutorials/``.
+
+To quickly test whether you installation was successful and everything works as expected, check out the minimum
+working example provided below:
+
+.. code-block:: python
+
+    """Minimum working example showing how to use Propulate."""
+    import propulate
+    from mpi4py import MPI
+    import random
+
+    # Set the communicator and the optimization parameters
+    comm = MPI.COMM_WORLD
+    rng = random.Random(MPI.COMM_WORLD.rank)
+    population_size = comm.size * 2
+    generations = 100
+    checkpoint = "./propulate_checkpoints"
+    propulate.utils.set_logger_config()
+
+    # Define the function to minimize and the search space, e.g., a 2D sphere function on (-5.12, 5.12)^2.
+    def loss_fn(params):
+        """Loss function to minimize."""
+        return params["x"] ** 2 + params["y"] ** 2
+
+    limits = {"x": (-5.12, 5.12), "y": (-5.12, 5.12)}
+
+    # Initialize the propagator and propulator with default parameters.
+    propagator = propulate.utils.get_default_propagator(
+        pop_size=population_size,
+        limits=limits,
+        rng=rng
+    )
+    propulator = propulate.Propulator(
+        loss_fn=loss_fn,
+        propagator=propagator,
+        rng=rng,
+        island_comm=comm,
+        generations=generations,
+        checkpoint_path=checkpoint,
+    )
+
+    # Run optimization and get summary of results.
+    propulator.propulate()
+    propulator.summarize()
 
 ``Propulate`` has more functionality and configurable options available, to accomodate various use cases and workflows.
 To learn more about how to use ``Propulate`` and adapt it to your needs, check out the :ref:`usage` section of the
@@ -22,4 +66,4 @@ documentation |:raised_hands:|.
 
 
 .. Links
-.. _subfolder: https://github.com/Helmholtz-AI-Energy/propulate/tree/master/scripts
\ No newline at end of file
+.. _subfolder: https://github.com/Helmholtz-AI-Energy/propulate/tree/master/tutorials
\ No newline at end of file
diff --git a/docs/tut_hpo.rst b/docs/tut_hpo.rst
index 47213b06..4c7f0c30 100644
--- a/docs/tut_hpo.rst
+++ b/docs/tut_hpo.rst
@@ -31,9 +31,11 @@ Thus, our search space dictionary looks as follows:
 
 .. code-block:: python
 
-    limits = {"conv_layers": (2, 10),  # number of convolutional layers, int for ordinal
-              "activation": ("relu", "sigmoid", "tanh"),  # activation function, str for categorical
-              "learning_rate": (0.01, 0.0001)}  # learning rate, float for continuous
+    limits = {
+        "conv_layers": (2, 10),  # Number of convolutional layers, int for ordinal
+        "activation": ("relu", "sigmoid", "tanh"),  # Activation function, str for categorical
+        "learning_rate": (0.01, 0.0001),  # Learning rate, float for continuous
+    }
 
 When tuning an ML model's hyperparameters, evaluating an individual during the optimization corresponds to training
 a neural network instance using a specific combination of hyperparameters to be optimized. In addition, we need some
@@ -47,33 +49,65 @@ Below, we show you how to do this using the example of the most important code s
 the hyperparameters we want to optimize are highlighted in pink. We start with defining the neural network which looks like this:
 
 .. code-block:: python
-    :emphasize-lines: 5,6,7,26,28,29,33-35,106,107,120
+    :emphasize-lines: 38-40,59,61,65-77
 
     class Net(LightningModule):
-        """Neural network class."""
+        """
+        Neural network class.
+
+        Attributes
+        ----------
+        best_accuracy : float
+            The model's best validation accuracy.
+        conv_layers : torch.nn.modules.container.Sequential
+            The model's convolutional layers.
+        fc : torch.nn.modules.linear.Linear
+            The model's fully connected layers.
+        loss_fn : torch.nn.modules.loss
+            The loss function used for training the model.
+        lr : float
+            The learning rate.
+        train_acc : torchmetrics.classification.accuracy.Accuracy
+            The accuracy metric used for evaluating model performance on the training dataset.
+        val_acc : torchmetrics.classification.accuracy.Accuracy
+            The accuracy metric used for evaluating model performance on the validation dataset.
+
+        Methods
+        -------
+        forward()
+            The forward pass.
+        training_step()
+            Calculate loss for training step in Lightning train loop.
+        validation_step()
+            Calculate loss for validation step in Lightning validation loop during training.
+        configure_optimizers()
+            Configure the optimizer.
+        on_validation_epoch_end()
+            Calculate and store the model's validation accuracy after each epoch.
+        """
 
         def __init__(
-                self,
-                conv_layers: int,
-                activation: torch.nn.modules.activation,
-                lr: float,
-                loss_fn: torch.nn.modules.loss,
+            self,
+            conv_layers: int,
+            activation: torch.nn.modules.activation,
+            lr: float,
+            loss_fn: torch.nn.modules.loss,
         ) -> None:
             """
             Set up neural network.
 
             Parameters
             ----------
-            conv_layers: int
-                         number of convolutional layers
-            activation: torch.nn.modules.activation
-                        activation function to use
-            lr: float
-                learning rate
-            loss_fn: torch.nn.modules.loss
-                     loss function
+            conv_layers : int
+                The number of convolutional layers.
+            activation : torch.nn.modules.activation
+                The activation function to use.
+            lr : float
+                The learning rate.
+            loss_fn : torch.nn.modules.loss
+                The loss function.
             """
-            super(Net, self).__init__()
+            super().__init__()
 
             self.lr = lr  # Set learning rate.
             self.loss_fn = loss_fn  # Set the loss function used for training the model.
@@ -106,13 +140,13 @@ the hyperparameters we want to optimize are highlighted in pink. We start with d
 
             Parameters
             ----------
-            x: torch.Tensor
-               data sample
+            x : torch.Tensor
+               The data sample.
 
             Returns
             -------
             torch.Tensor
-                The model's predictions for input data sample
+                The model's predictions for input data sample.
             """
             b, c, w, h = x.size()
             x = self.conv_layers(x)
@@ -128,15 +162,15 @@ the hyperparameters we want to optimize are highlighted in pink. We start with d
 
             Parameters
             ----------
-            batch: Tuple[torch.Tensor, torch.Tensor]
-                   input batch
-            batch_idx: int
-                       batch index
+            batch : Tuple[torch.Tensor, torch.Tensor]
+                The input batch.
+            batch_idx : int
+                Its batch index.
 
             Returns
             -------
             torch.Tensor
-                training loss for input batch
+                The training loss for this input batch.
             """
             x, y = batch
             pred = self(x)
@@ -154,15 +188,15 @@ the hyperparameters we want to optimize are highlighted in pink. We start with d
 
             Parameters
             ----------
-            batch: Tuple[torch.Tensor, torch.Tensor]
-                   current batch
-            batch_idx: int
-                       batch index
+            batch : Tuple[torch.Tensor, torch.Tensor]
+                The current batch
+            batch_idx : int
+                The batch index.
 
             Returns
             -------
             torch.Tensor
-                validation loss for input batch
+                The validation loss for the input batch.
             """
             x, y = batch
             pred = self(x)
@@ -174,19 +208,17 @@ the hyperparameters we want to optimize are highlighted in pink. We start with d
 
         def configure_optimizers(self) -> torch.optim.SGD:
             """
-            Configure optimizer.
+            Configure the optimizer.
 
             Returns
             -------
             torch.optim.sgd.SGD
-                stochastic gradient descent optimizer
+                A stochastic gradient descent optimizer.
             """
             return torch.optim.SGD(self.parameters(), lr=self.lr)
 
         def on_validation_epoch_end(self):
-            """
-            Calculate and store the model's validation accuracy after each epoch.
-            """
+            """Calculate and store the model's validation accuracy after each epoch."""
             val_acc_val = self.val_acc.compute()
             self.val_acc.reset()
             if val_acc_val > self.best_accuracy:
@@ -202,42 +234,51 @@ We also need some helper function to load the MNIST data:
 
         Parameters
         ----------
-        batch_size: int
-                    batch size
+        batch_size : int
+            The batch size.
 
         Returns
         -------
-        DataLoader
-            training dataloader
-        DataLoader
-            validation dataloader
+        torch.utils.data.DataLoader
+            The training dataloader.
+        torch.utils.data.DataLoader
+            The validation dataloader.
         """
-        ...
-
         data_transform = Compose([ToTensor(), Normalize((0.1307,), (0.3081,))])
+        num_workers = NUM_WORKERS
+        log.info(f"Use {num_workers} workers in dataloader.")
 
-        if MPI.COMM_WORLD.Get_rank() == 0:  # Only root downloads data.
+        if MPI.COMM_WORLD.rank == 0:  # Only root downloads data.
             train_loader = DataLoader(
                 dataset=MNIST(
                     download=True, root=".", transform=data_transform, train=True
                 ),  # Use MNIST training dataset.
                 batch_size=batch_size,  # Batch size
+                num_workers=num_workers,
+                pin_memory=True,
+                persistent_workers=True,
                 shuffle=True,  # Shuffle data.
             )
 
-        MPI.COMM_WORLD.Barrier()
-        if MPI.COMM_WORLD.Get_rank() != 0:
+        MPI.COMM_WORLD.barrier()
+        if MPI.COMM_WORLD.rank != 0:
             train_loader = DataLoader(
                 dataset=MNIST(
                     download=False, root=".", transform=data_transform, train=True
                 ),  # Use MNIST training dataset.
                 batch_size=batch_size,  # Batch size
+                num_workers=num_workers,
+                pin_memory=True,
+                persistent_workers=True,
                 shuffle=True,  # Shuffle data.
             )
         val_loader = DataLoader(
             dataset=MNIST(
                 download=False, root=".", transform=data_transform, train=False
             ),  # Use MNIST testing dataset.
+            num_workers=num_workers,
+            pin_memory=True,
+            persistent_workers=True,
             batch_size=1,  # Batch size
             shuffle=False,  # Do not shuffle data.
         )
@@ -247,53 +288,65 @@ Now we are ready to set up the ``Propulate`` loss function that is minimized dur
 order to find the best hyperparameters for our model:
 
 .. code-block:: python
-    :emphasize-lines: 18-20, 24-26, 29, 30, 48
+    :emphasize-lines: 16-18, 22-26, 32-35, 60
 
-    def ind_loss(
-            params: Dict[str, Union[int, float, str]]
-    ) -> float:
+    def ind_loss(params: Dict[str, Union[int, float, str]]) -> float:
         """
-        Loss function for evolutionary optimization with Propulate.
-        We minimize the model's negative validation accuracy.
+        Loss function for evolutionary optimization with Propulate. Minimize the model's negative validation accuracy.
 
         Parameters
         ----------
-        params: dict[str, int | float | str]]
+        params : Dict[str, int | float | str]
+            The hyperparameters to be optimized evolutionarily.
 
         Returns
         -------
         float
-            The trained model's negative validation accuracy
+            The trained model's negative validation accuracy.
         """
         # Extract hyperparameter combination to test from input dictionary.
         conv_layers = params["conv_layers"]  # Number of convolutional layers
         activation = params["activation"]  # Activation function
         lr = params["lr"]  # Learning rate
 
-        epochs = 2  # Number of epochs to train
+        epochs = 100
 
-        # Define the activation function mapping.
-        activations = {"relu": nn.ReLU, "sigmoid": nn.Sigmoid, "tanh": nn.Tanh}
+        activations = {
+            "relu": nn.ReLU,
+            "sigmoid": nn.Sigmoid,
+            "tanh": nn.Tanh,
+        }  # Define activation function mapping.
         activation = activations[activation]  # Get activation function.
-        loss_fn = torch.nn.CrossEntropyLoss()  # Use cross-entropy loss for multi-class classification.
+        loss_fn = (
+            torch.nn.CrossEntropyLoss()
+        )  # Use cross-entropy loss for multi-class classification.
 
-        model = Net(conv_layers, activation, lr, loss_fn)  # Set up neural network with specified hyperparameters.
+        model = Net(
+            conv_layers, activation, lr, loss_fn
+        )  # Set up neural network with specified hyperparameters.
         model.best_accuracy = 0.0  # Initialize the model's best validation accuracy.
 
-        train_loader, val_loader = get_data_loaders(batch_size=8)  # Get training and validation data loaders.
+        train_loader, val_loader = get_data_loaders(
+            batch_size=8
+        )  # Get training and validation data loaders.
+
+        tb_logger = loggers.TensorBoardLogger(
+            save_dir=log_path + "/lightning_logs"
+        )  # Get tensor board logger.
 
         # Under the hood, the Lightning Trainer handles the training loop details.
-        trainer = Trainer(max_epochs=epochs,  # Stop training once this number of epochs is reached.
-                          accelerator="gpu",  # Pass accelerator type.
-                          devices=[  # Devices to train on
-                              MPI.COMM_WORLD.Get_rank() % GPUS_PER_NODE
-                                  ],
-                          enable_progress_bar=False,  # Disable progress bar.
-                          )
+        trainer = Trainer(
+            max_epochs=epochs,  # Stop training once this number of epochs is reached.
+            accelerator="gpu",  # Pass accelerator type.
+            devices=[MPI.COMM_WORLD.rank % GPUS_PER_NODE],  # Devices to train on
+            callbacks=[EarlyStopping(monitor="val_loss", mode="min")],
+            enable_progress_bar=False,  # Disable progress bar.
+            logger=tb_logger,  # Logger
+        )
         trainer.fit(  # Run full model training optimization routine.
             model=model,  # Model to train
             train_dataloaders=train_loader,  # Dataloader for training samples
-            val_dataloaders=val_loader  # Dataloader for validation samples
+            val_dataloaders=val_loader,  # Dataloader for validation samples
         )
         # Return negative best validation accuracy as an individual's loss.
         return -model.best_accuracy.item()
@@ -303,37 +356,53 @@ asynchronous island model (``Islands``) which takes care of the actual genetic o
 
 .. code-block:: python
 
-    num_generations = 3  # Number of generations
-    pop_size = 2 * MPI.COMM_WORLD.size  # Breeding population size
-    limits = {
-        "conv_layers": (2, 10),
-        "activation": ("relu", "sigmoid", "tanh"),
-        "lr": (0.01, 0.0001),
-    }  # Search space
-    rng = random.Random(
-        MPI.COMM_WORLD.rank
-    )  # Set up separate random number generator for evolutionary optimizer.
-    propagator = get_default_propagator(  # Get default evolutionary operator.
-        pop_size=pop_size,  # Breeding population size
-        limits=limits,  # Search space
-        mate_prob=0.7,  # Crossover probability
-        mut_prob=0.4,  # Mutation probability
-        random_prob=0.1,  # Random-initialization probability
-        rng=rng,  # Random number generator for evolutionary optimizer
-    )
-    islands = Islands(  # Set up island model.
-        loss_fn=ind_loss,  # Loss function to optimize
-        propagator=propagator,  # Evolutionary operator
-        rng=rng,  # Random number generator
-        generations=num_generations,  # Number of generations per worker
-        num_islands=1,  # Number of islands
-        checkpoint_path=log_path,
-    )
-    islands.evolve(  # Run evolutionary optimization.
-        top_n=1,  # Print top-n best individuals on each island in summary.
-        logging_interval=1,  # Logging interval
-        debug=2,  # Verbosity level
-    )
+    if __name__ == "__main__":
+        comm = MPI.COMM_WORLD
+        num_generations = 10  # Number of generations
+        pop_size = 2 * comm.size  # Breeding population size
+        limits = {
+            "conv_layers": (2, 10),
+            "activation": ("relu", "sigmoid", "tanh"),
+            "lr": (0.01, 0.0001),
+        }  # Define search space.
+        rng = random.Random(
+            comm.rank
+        )  # Set up separate random number generator for evolutionary optimizer.
+        propagator = get_default_propagator(  # Get default evolutionary operator.
+            pop_size=pop_size,  # Breeding population size
+            limits=limits,  # Search space
+            crossover_prob=0.7,  # Crossover probability
+            mutation_prob=0.4,  # Mutation probability
+            random_init_prob=0.1,  # Random-initialization probability
+            rng=rng,  # Separate random number generator for Propulate optimization
+        )
+
+        # Set up separate logger for Propulate optimization.
+        set_logger_config(
+            level=logging.INFO,  # Logging level
+            log_file=f"{log_path}/{pathlib.Path(__file__).stem}.log",  # Logging path
+            log_to_stdout=True,  # Print log on stdout.
+            log_rank=False,  # Do not prepend MPI rank to logging messages.
+            colors=True,  # Use colors.
+        )
+
+        # Set up propulator performing actual optimization.
+        propulator = Propulator(
+            loss_fn=ind_loss,  # Loss function to optimize
+            propagator=propagator,  # Evolutionary operator
+            rng=rng,  # Random number generator
+            island_comm=comm,  # Communicator
+            generations=num_generations,  # Number of generations per worker
+            checkpoint_path=log_path,  # Path to save checkpoints to
+        )
+
+        # Run optimization and print summary of results.
+        propulator.propulate(
+            logging_interval=1, debug=2  # Logging interval and verbosity level
+        )
+        propulator.summarize(
+            top_n=1, debug=2  # Print top-n best individuals on each island in summary.
+        )
 
 After additionally installing ``torch``, ``lightning``, and ``tensorboard``, we can use the the virtual environment from
 the previous example to run a search with four GPUs on a single node:
diff --git a/docs/tut_islands.rst b/docs/tut_islands.rst
index c968bcca..02a7c6ed 100644
--- a/docs/tut_islands.rst
+++ b/docs/tut_islands.rst
@@ -29,10 +29,12 @@ island. With ``num_migrants = 1``, this is the default behaviour in ``Propulate`
 
 .. code-block:: python
 
-    migration_topology = config.num_migrants * np.ones(  # Set up fully connected migration topology.
+    # Set up fully connected migration topology.
+    migration_topology = config.num_migrants * np.ones(
         (config.num_islands, config.num_islands),
         dtype=int)
-    np.fill_diagonal(migration_topology, 0)  # An island does not send migrants to itself.
+    # An island does not send migrants to itself.
+    np.fill_diagonal(migration_topology, 0)
 
 Next, we set up the island model itself using the ``Islands`` class. In addition to the ``Propulator`` arguments defining
 the islands' internal asynchronous optimization process, ``Islands`` takes all migration-relevant arguments
@@ -40,7 +42,8 @@ to configure the islands and the migration between them:
 
 .. code-block:: python
 
-    islands = Islands(  # Set up island model.
+    # Set up the island model.
+    islands = Islands(
         loss_fn=sphere,  # Loss function
         propagator=propagator,  # Evolutionary operator
         rng=rng,  # Separate random number generator for the optimization process
@@ -67,7 +70,7 @@ Now we are ready to run the optimization:
     islands.evolve(  # Run actual optimization.
         top_n=config.top_n, # Number of best individuals to print in the summary.
         logging_interval=config.logging_int, # Logging interval
-        debug=config.verbosity)  # Verbosity level
+        debug=config.verbosity)  # Debug level
 
 
 .. note::
@@ -78,7 +81,7 @@ You can run the example script ``islands_example.py``:
 
 .. code-block:: console
 
-    $ mpirun --use-hwthread-cpus python propulator_example.py
+    $ mpirun --use-hwthread-cpus python islands_example.py
 
 With ten MPI ranks and two islands with five workers each, the output looks like this:
 
@@ -103,49 +106,43 @@ With ten MPI ranks and two islands with five workers each, the output looks like
     ⠀        ⢶⣤⣀⣀⣠⣴⠶⠛⠋⠙⠻⣦⣄⣀⣀⣠⣤⣴⠶⠶⣦⣄⣀⣀⣠⣤⣤⡶⠀
             ⠀⠀⠈⠉⠉⠉⠀⠀⠀⠀⠀⠀⠀⠉⠉⠉⠉⠀⠀⠀⠀⠀⠉⠉⠉⠉⠀⠀⠀⠀
 
-    Worker distribution [0 0 0 0 0 1 1 1 1 1] with island counts [5 5] and island displacements [0 5].
-    Migration topology [[0 1]
+    [2024-03-13 12:30:35,048][propulate.islands][INFO] - Worker distribution [0 0 0 0 0 1 1 1 1 1] with island counts [5 5] and island displacements [0 5].
+    [2024-03-13 12:30:35,049][propulate.islands][INFO] - Migration topology [[0 1]
      [1 0]] has shape (2, 2).
-    NOTE: Island migration probability 0.9 results in per-rank migration probability 0.18.
+    [2024-03-13 12:30:35,049][propulate.islands][INFO] - NOTE: Island migration probability 0.9 results in per-rank migration probability 0.18.
     Starting parallel optimization process.
-    NOTE: No valid checkpoint file given. Initializing population randomly...
-    Island 1 has 5 workers.
-    No pollination.
-    NOTE: No valid checkpoint file given. Initializing population randomly...
-    Island 0 has 5 workers.
-    Island 1 Worker 1: In generation 0...
-    Island 1 Worker 2: In generation 0...
-    Island 0 Worker 2: In generation 0...
-    Island 0 Worker 3: In generation 0...
-    Island 0 Worker 0: In generation 0...
-    Island 1 Worker 0: In generation 0...
-    Island 1 Worker 3: In generation 0...
-    Island 0 Worker 4: In generation 0...
-    Island 1 Worker 4: In generation 0...
-    Island 0 Worker 1: In generation 0...
-
+    [2024-03-13 12:30:35,049][propulate.islands][INFO] - Use island model with real migration.
+    [2024-03-13 12:30:35,049][propulate.propulator][INFO] - No valid checkpoint file given. Initializing population randomly...
+    [2024-03-13 12:30:35,049][propulate.migrator][INFO] - Island 0 has 5 workers.
+    [2024-03-13 12:30:35,049][propulate.migrator][INFO] - Island 0 Worker 0: In generation 0...
+    [2024-03-13 12:30:35,049][propulate.propulator][INFO] - No valid checkpoint file given. Initializing population randomly...
+    [2024-03-13 12:30:35,049][propulate.migrator][INFO] - Island 1 has 5 workers.
+    [2024-03-13 12:30:35,049][propulate.migrator][INFO] - Island 1 Worker 0: In generation 0...
+    [2024-03-13 12:30:35,049][propulate.migrator][INFO] - Island 0 Worker 2: In generation 0...
+    [2024-03-13 12:30:35,049][propulate.migrator][INFO] - Island 1 Worker 3: In generation 0...
+    [2024-03-13 12:30:35,049][propulate.migrator][INFO] - Island 0 Worker 3: In generation 0...
+    [2024-03-13 12:30:35,049][propulate.migrator][INFO] - Island 1 Worker 1: In generation 0...
+    [2024-03-13 12:30:35,049][propulate.migrator][INFO] - Island 1 Worker 2: In generation 0...
+    [2024-03-13 12:30:35,049][propulate.migrator][INFO] - Island 0 Worker 4: In generation 0...
+    [2024-03-13 12:30:35,049][propulate.migrator][INFO] - Island 1 Worker 4: In generation 0...
+    [2024-03-13 12:30:35,049][propulate.migrator][INFO] - Island 0 Worker 1: In generation 0...
     ...
-
-    Island 0 Worker 0: In generation 980...
-    Island 1 Worker 4: In generation 990...
-    Island 1 Worker 2: In generation 940...
-    Island 0 Worker 0: In generation 990...
-    Island 1 Worker 2: In generation 950...
-    Island 1 Worker 2: In generation 960...
-    Island 1 Worker 2: In generation 970...
-    Island 1 Worker 2: In generation 980...
-    Island 1 Worker 2: In generation 990...
-    OPTIMIZATION DONE.
-    NEXT: Final checks for incoming messages...
-
-    ###########
+    [2024-03-13 12:30:42,928][propulate.migrator][INFO] - Island 0 Worker 1: In generation 990...
+    [2024-03-13 12:30:42,959][propulate.migrator][INFO] - Island 0 Worker 0: In generation 950...
+    [2024-03-13 12:30:42,980][propulate.migrator][INFO] - Island 0 Worker 0: In generation 960...
+    [2024-03-13 12:30:43,010][propulate.migrator][INFO] - Island 0 Worker 0: In generation 970...
+    [2024-03-13 12:30:43,038][propulate.migrator][INFO] - Island 0 Worker 0: In generation 980...
+    [2024-03-13 12:30:43,092][propulate.migrator][INFO] - Island 0 Worker 0: In generation 990...
+    [2024-03-13 12:30:43,121][propulate.migrator][INFO] - OPTIMIZATION DONE.
+    [2024-03-13 12:30:43,121][propulate.migrator][INFO] - NEXT: Final checks for incoming messages...
+    [2024-03-13 12:30:43,292][propulate.propulator][INFO] - ###########
     # SUMMARY #
     ###########
-
     Number of currently active individuals is 10000.
     Expected overall number of evaluations is 10000.
-    Top 1 result(s) on island 1:
-    (1): [{'a': '-1.84E-4', 'b': '9.84E-4'}, loss 1.00E-6, island 1, worker 0, generation 248]
+    [2024-03-13 12:30:46,536][propulate.propulator][INFO] - Top 1 result(s) on island 1:
+    (1): [{'a': '-2.83E-4', 'b': '1.04E-3'}, loss 1.16E-6, island 0, worker 3, generation 901]
+
+    [2024-03-13 12:30:46,611][propulate.propulator][INFO] - Top 1 result(s) on island 0:
+    (1): [{'a': '-2.83E-4', 'b': '1.04E-3'}, loss 1.16E-6, island 0, worker 3, generation 901]
 
-    Top 1 result(s) on island 0:
-    (1): [{'a': '-1.84E-4', 'b': '9.84E-4'}, loss 1.00E-6, island 1, worker 0, generation 248]
diff --git a/docs/tut_multi_rank_worker.rst b/docs/tut_multi_rank_worker.rst
new file mode 100644
index 00000000..2db2fabc
--- /dev/null
+++ b/docs/tut_multi_rank_worker.rst
@@ -0,0 +1,14 @@
+.. _tut_multi_rank_worker:
+
+Using Propulate with Multi-Rank Workers
+=======================================
+.. note::
+
+   You can find the corresponding ``Python`` script here:
+   https://github.com/Helmholtz-AI-Energy/propulate/blob/master/tutorials/multi_rank_workers_example.py
+
+In addition to the already explained functionality, ``Propulate`` enables using multi-rank workers for an internally
+parallelized evaluation of the loss function. This is useful for, e.g., data-parallel training of neural networks during
+the hyperparameter optimization, where each individual network is trained on multiple GPUs.
+
+A more detailed explanation of the tutorial will be available soon |:rocket:|.
\ No newline at end of file
diff --git a/docs/tut_propulator.rst b/docs/tut_propulator.rst
index eb7f513c..d7e1d0fb 100644
--- a/docs/tut_propulator.rst
+++ b/docs/tut_propulator.rst
@@ -7,10 +7,11 @@ Evolutionary Optimization of a Mathematical Function
    You can find the corresponding ``Python`` script here:
    https://github.com/Helmholtz-AI-Energy/propulate/blob/master/tutorials/propulator_example.py
 
-The basic optimization mechanism in ``Propulate`` is that of Darwinian evolution, i.e.,
-beneficial traits are selected, recombined, and mutated to breed more fit individuals.
-To show you how ``Propulate`` works, we use its *basic asynchronous evolutionary optimizer* to minimize
-two-dimensional mathematical functions. Let us consider the sphere function:
+The basic optimization mechanism in ``Propulate`` is that of Darwinian evolution, i.e., beneficial traits are selected,
+recombined, and mutated to breed more fit individuals.
+To show you how ``Propulate`` works, we use its *basic asynchronous evolutionary optimizer* to minimize two-dimensional
+mathematical functions.
+Let us consider the sphere function:
 
 .. math::
     f_\mathrm{sphere}\left(x,y\right)=x^2+y^2
@@ -43,17 +44,21 @@ As the very first step, we need to define the key ingredients that define the op
 
   .. code-block:: python
 
-    limits = {"learning_rate": (0.001, 0.01),
-              "conv_layers": (2, 10),
-              "activation": ("relu", "sigmoid", "tanh")}
+    limits = {
+        "learning_rate": (0.001, 0.01),
+        "conv_layers": (2, 10),
+        "activation": ("relu", "sigmoid", "tanh")
+    }
 
   The sphere function has two continuous parameters, :math:`x` and :math:`y`, and we consider
   :math:`x,y \in\left[-5.12, 5.12\right]`. The search space in our example thus looks like this:
 
   .. code-block:: python
 
-    limits = {"x": (-5.12, 5.12),
-              "y": (-5.12, 5.12)}
+    limits = {
+        "x": (-5.12, 5.12),
+        "y": (-5.12, 5.12)
+    }
 
 * The fitness or *loss function* (also known as the objective function). This is the function we want to optimize in order
   to find the best parameters. It can be any ``Python`` function with the following characteristics:
@@ -74,40 +79,54 @@ As the very first step, we need to define the key ingredients that define the op
 
     def sphere(params: Dict[str, float]) -> float:
         """
-        Sphere function: continuous, convex, separable, differentiable, unimodal
+        Sphere function: continuous, convex, separable, differentiable, unimodal.
 
         Input domain: -5.12 <= x, y <= 5.12
         Global minimum 0 at (x, y) = (0, 0)
 
         Parameters
         ----------
-        params: dict[str, float]
-                function parameters
+        params: Dict[str, float]
+            The function parameters.
+            
         Returns
         -------
         float
-            function value
+            The function value.
         """
-        return numpy.sum(numpy.array(list(params.values())) ** 2)
-
+        return numpy.sum(numpy.array(list(params.values())) ** 2).item()
 
 Next, we need to define the evolutionary operator or propagator that we want to use to breed new individuals during the
-optimization process. ``Propulate`` provides a reasonable default propagator via a utility function, ``get_default_propagator``,
-that serves as a good start for the most optimization problems. You can adapt its hyperparameters, such as crossover and mutation
-probability, as you wish. In the example script, you can pass those hyperparameters as command-line options (this is the
-``config`` in the code snippet below) or just use the default values. You also need to pass a separate random number
-generator that is used exclusively in the evolutionary optimization process (and not in the objective function):
+optimization process. ``Propulate`` provides a reasonable default propagator via a utility function,
+``get_default_propagator``, that serves as a good start for the most optimization problems. You can adapt its
+hyperparameters, such as crossover and mutation probability, as you wish. In the example script, you can pass those
+hyperparameters as command-line options (this is the ``config`` in the code snippet below) or just use the default
+values. You also need to pass a separate random number generator that is used exclusively in the evolutionary
+optimization process (and not in the objective function).
+In addition, you can adapt the separate logger used to track the ``Propulate`` optimization with the utility function
+``set_logger_config`` as shown below:
 
 .. code-block:: python
 
-    rng = random.Random(config.seed+MPI.COMM_WORLD.rank)  # Separate random number generator for optimization.
+    # Set up separate logger for Propulate optimization.
+    propulate.set_logger_config(
+        level=config.logging_level,  # Logging level
+        log_file=f"{config.checkpoint}/{pathlib.Path(__file__).stem}.log",  # Logging path
+        log_to_stdout=True,  # Print log on stdout.
+        log_rank=False,  # Do not prepend MPI rank to logging messages.
+        colors=True,  # Use colors.
+    )
+    rng = random.Random(
+        config.seed + MPI.COMM_WORLD.rank
+    )  # Separate random number generator for optimization.
     propagator = propulate.utils.get_default_propagator(  # Get default evolutionary operator.
         pop_size=config.pop_size,  # Breeding pool size
         limits=limits,  # Search-space limits
-        mate_prob=config.crossover_probability,  # Crossover probability
-        mut_prob=config.mutation_probability,  # Mutation probability
-        random_prob=config.random_init_probability,  # Random-initialization probability
-        rng=rng)  # Random number generator for the optimization process
+        crossover_prob=config.crossover_probability,  # Crossover probability
+        mutation_prob=config.mutation_probability,  # Mutation probability
+        random_init_prob=config.random_init_probability,  # Random-initialization probability
+        rng=rng  # Random number generator for the optimization process
+    )
 
 We also need to set up the actual evolutionary optimizer, that is a so-called ``Propulator`` instance. This will handle the
 parallel asynchronous optimization process for us:
@@ -117,10 +136,10 @@ parallel asynchronous optimization process for us:
     propulator = Propulator(  # Set up propulator performing actual optimization.
         loss_fn=sphere,  # Loss function to minimize
         propagator=propagator,  # Evolutionary operator
-        comm=MPI.COMM_WORLD,  # Communicator
+        rng=rng,  # Random number generator for optimization process
         generations=config.generations,  # Number of generations
-        checkpoint_path=config.checkpoint,  # Checkpoint path
-        rng=rng)  # Random number generator for optimization process
+        checkpoint_path=config.checkpoint  # Checkpoint path
+    )
 
 Now it's time to run the actual optimization. Overall, ``generations * MPI.COMM_WORLD.size`` evaluations will be performed:
 
@@ -138,36 +157,32 @@ The output looks like this:
     # PROPULATE: Parallel Propagator of Populations #
     #################################################
 
-    NOTE: No valid checkpoint file given. Initializing population randomly...
-    Island 0 has 4 workers.
-    Island 0 Worker 0: In generation 0...
-    Island 0 Worker 2: In generation 0...
-    Island 0 Worker 1: In generation 0...
-    Island 0 Worker 3: In generation 0...
-    Island 0 Worker 0: In generation 10...
-    Island 0 Worker 1: In generation 10...
-    Island 0 Worker 3: In generation 10...
-    Island 0 Worker 2: In generation 10...
-    Island 0 Worker 0: In generation 20...
-    Island 0 Worker 1: In generation 20...
-    Island 0 Worker 2: In generation 20...
-    Island 0 Worker 3: In generation 20...
+    [2024-03-12 14:37:01,374][propulate.propulator][INFO] - No valid checkpoint file given. Initializing population randomly...
+    [2024-03-12 14:37:01,374][propulate.propulator][INFO] - Island 0 has 4 workers.
+    [2024-03-12 14:37:01,374][propulate.propulator][INFO] - Island 0 Worker 0: In generation 0...
+    [2024-03-12 14:37:01,374][propulate.propulator][INFO] - Island 0 Worker 3: In generation 0...
+    [2024-03-12 14:37:01,374][propulate.propulator][INFO] - Island 0 Worker 2: In generation 0...
+    [2024-03-12 14:37:01,374][propulate.propulator][INFO] - Island 0 Worker 1: In generation 0...
+    [2024-03-12 14:37:01,377][propulate.propulator][INFO] - Island 0 Worker 3: In generation 10...
+    [2024-03-12 14:37:01,377][propulate.propulator][INFO] - Island 0 Worker 1: In generation 10...
+    [2024-03-12 14:37:01,378][propulate.propulator][INFO] - Island 0 Worker 0: In generation 10...
+    [2024-03-12 14:37:01,378][propulate.propulator][INFO] - Island 0 Worker 2: In generation 10...
 
     ...
-
-    Island 0 Worker 0: In generation 990...
-    Island 0 Worker 1: In generation 990...
-    OPTIMIZATION DONE.
+    [2024-03-12 14:37:02,197][propulate.propulator][INFO] - Island 0 Worker 1: In generation 960...
+    [2024-03-12 14:37:02,206][propulate.propulator][INFO] - Island 0 Worker 2: In generation 990...
+    [2024-03-12 14:37:02,206][propulate.propulator][INFO] - Island 0 Worker 1: In generation 970...
+    [2024-03-12 14:37:02,215][propulate.propulator][INFO] - Island 0 Worker 1: In generation 980...
+    [2024-03-12 14:37:02,224][propulate.propulator][INFO] - Island 0 Worker 1: In generation 990...
+    [2024-03-12 14:37:02,232][propulate.propulator][INFO] - OPTIMIZATION DONE.
     NEXT: Final checks for incoming messages...
-
-    ###########
+    [2024-03-12 14:37:02,244][propulate.propulator][INFO] - ###########
     # SUMMARY #
     ###########
-
     Number of currently active individuals is 4000.
     Expected overall number of evaluations is 4000.
-    Top 1 result(s) on island 0:
-    (1): [{'a': '3.90E-3', 'b': '1.16E-3'}, loss 1.66E-5, island 0, worker 2, generation 739]
+    [2024-03-12 14:37:03,703][propulate.propulator][INFO] - Top 1 result(s) on island 0:
+    (1): [{'a': '2.91E-3', 'b': '-3.05E-3'}, loss 1.78E-5, island 0, worker 0, generation 956]
 
 Let's Get Your Hands Dirty (At Least a Bit)
 -------------------------------------------
diff --git a/docs/usage.rst b/docs/usage.rst
index 9af87f9d..876b4b6c 100644
--- a/docs/usage.rst
+++ b/docs/usage.rst
@@ -12,6 +12,7 @@ If you want to use ``Propulate`` for your own applications, you can use these sc
    tut_propulator
    tut_islands
    tut_hpo
+   tut_multi_rank_worker
 
 .. Links
 .. _Github: https://github.com/Helmholtz-AI-Energy/propulate
\ No newline at end of file
diff --git a/propulate/islands.py b/propulate/islands.py
index 67b8e4aa..a4cc042c 100644
--- a/propulate/islands.py
+++ b/propulate/islands.py
@@ -2,7 +2,7 @@
 import platform
 import random
 from pathlib import Path
-from typing import Callable, Union, List, Type
+from typing import Callable, List, Optional, Type, Union
 
 import numpy as np
 from mpi4py import MPI
@@ -24,6 +24,16 @@ class Islands:
     migrants to be included in the target islands' populations. With worse performing islands typically receiving
     candidates from better performing ones, islands communicate genetic information competitively, thus increasing
     diversity among the subpopulations compared to panmictic models.
+
+    Attributes
+    ----------
+    propulator : propulate.Propulator
+        The Propulator instance performing the actual optimization.
+
+    Methods
+    -------
+    evolve()
+        Run Propulate optimization and summarize results.
     """
 
     def __init__(
@@ -33,67 +43,64 @@ def __init__(
         rng: random.Random,
         generations: int = 0,
         num_islands: int = 1,
-        island_sizes: np.ndarray = None,
-        migration_topology: np.ndarray = None,
+        island_sizes: Optional[np.ndarray] = None,
+        migration_topology: Optional[np.ndarray] = None,
         migration_probability: float = 0.0,
         emigration_propagator: Type[Propagator] = SelectMin,
         immigration_propagator: Type[Propagator] = SelectMax,
         pollination: bool = True,
         checkpoint_path: Union[str, Path] = Path("./"),
+        ranks_per_worker: int = 1,
     ) -> None:
         """
-        Initialize island model with given parameters.
+        Initialize an island model with the given parameters.
 
         Parameters
         ----------
-        loss_fn: Callable
-                 loss function to be minimized
-        propagator: propulate.propagators.Propagator
-                    propagator to apply for breeding
-        rng: random.Random
-             random number generator
-        generations: int
-                     number of generations
-        num_islands: int
-                     number of separate, equally sized evolutionary islands (ignored if ``island_sizes`` is not None)
-                     (differences +-1 possible due to load balancing)
-        island_sizes: numpy.ndarray
-                      array with numbers of workers for each island (heterogeneous case)
-        migration_topology: numpy.ndarray
-                            2D matrix where entry (i,j) specifies how many individuals are sent
-                            by island i to island j
-                            (int: absolute number, float: relative fraction of population)
-        migration_probability: float
-                               probability of migration after each generation
-        emigration_propagator: type[propulate.propagators.Propagator]
-                               emigration propagator, i.e., how to choose individuals for emigration
-                               that are sent to destination island.
-                               Should be some kind of selection operator.
-        immigration_propagator: type[propulate.propagators.Propagator]
-                                immigration propagator, i.e., how to choose individuals on target island
-                                to be replaced by immigrants.
-                                Should be some kind of selection operator.
-        pollination: bool
-                     If True, copies of emigrants are sent, otherwise, emigrants are removed from
-                     original island.
-        checkpoint_path: Union[Path, str]
-                         Path where checkpoints are loaded from and stored.
+        loss_fn : Callable
+            The loss function to be minimized.
+        propagator : propulate.propagators.Propagator
+            The propagator, i.e., evolutionary operator, to apply for breeding.
+        rng : random.Random
+             The separate random number generator for the Propulate optimization.
+        generations : int, optional
+            The number of generations. Default is 0.
+        num_islands : int, optional
+            The number of separate equally sized evolutionary islands. Ignored if ``island_sizes`` is not ``None``.
+            Differences of +-1 are possible due to load balancing. Default is 1.
+        island_sizes : numpy.ndarray[int], optional
+            An array with numbers of workers for each island (heterogeneous case).
+        migration_topology : numpy.ndarray[int], optional
+            A two-dimensional matrix where entry (i,j) specifies how many individuals are sent by island i to island j
+        migration_probability : float, optional
+            The probability of migration after each generation.
+        emigration_propagator : Type[propulate.propagators.Propagator], optional
+            The emigration propagator, i.e., how to choose individuals for emigration that are sent to the destination
+            island. Should be some kind of selection operator. Default is ``SelectMin``.
+        immigration_propagator : Type[propulate.propagators.Propagator], optional
+            The immigration propagator, i.e., how to choose individuals on the target island to be replaced by the
+            immigrants. Should be some kind of selection operator. Default is ``SelectMax``.
+        pollination : bool, optional
+            If True, copies of emigrants are sent; otherwise, emigrants are removed from the original island.
+            Default is True.
+        checkpoint_path : pathlib.Path | str, optional
+            The path where checkpoints are loaded from and stored to. Default is the current working directory.
+        ranks_per_worker : int, optional
+            The number of ranks per worker. Default is 1.
 
         Raises
         ------
         ValueError
-            If specified number of islands is smaller than 1.
-        ValueError
-            If number of workers in custom worker distribution does not equal overall number of processors.
-        ValueError
-            If custom migration topology has the wrong shape.
-        ValueError
-            If migration probability is not within [0, 1].
+            If the overall number of ranks is not evenly divisible by the requested number of ranks per worker.
+            If the specified number of islands is smaller than 1.
+            If the number of workers in the custom worker distribution does not equal overall number of ranks.
+            If a custom migration topology has the wrong shape.
+            If the migration probability is not within [0, 1].
         """
-        # Set up communicator.
-        rank, size = MPI.COMM_WORLD.rank, MPI.COMM_WORLD.size
+        # Set up full world communicator.
+        full_world_rank, full_world_size = MPI.COMM_WORLD.rank, MPI.COMM_WORLD.size
 
-        if rank == 0:
+        if full_world_rank == 0:
             print(
                 "#################################################\n"
                 "# PROPULATE: Parallel Propagator of Populations #\n"
@@ -117,94 +124,146 @@ def __init__(
                     "        ⠀⠀⠈⠉⠉⠉⠀⠀⠀⠀⠀⠀⠀⠉⠉⠉⠉⠀⠀⠀⠀⠀⠉⠉⠉⠉⠀⠀⠀⠀\n"
                 )
 
-        # Homogeneous case with equal island sizes (differences of +-1 possible due to load balancing).
-        if island_sizes is None:
-            if num_islands < 1:
-                raise ValueError(
-                    f"Invalid number of evolutionary islands, needs to be >= 1 but was {num_islands}."
-                )
-            base_size = size // num_islands  # Base number of workers of each island
-            remainder = (
-                size % num_islands
-            )  # Number of remaining workers to be distributed
-            island_sizes = base_size * np.ones(num_islands, dtype=int)
-            island_sizes[
-                :remainder
-            ] += 1  # Distribute remaining workers equally for balanced load.
-
-        # Heterogeneous case with user-defined island sizes.
-        if np.sum(island_sizes) != size:
+        # Split the full world communicator into (multi rank) workers.
+        if full_world_size % ranks_per_worker != 0:
             raise ValueError(
-                f"There should be COMM_WORLD.size = {size} workers but only {np.sum(island_sizes)} were specified."
+                f"The number of overall ranks, i.e., {full_world_size}, should be evenly divisible by "
+                f"the number of ranks per worker, i.e., {ranks_per_worker}."
+            )
+        worker_idx = (
+            full_world_rank // ranks_per_worker
+        )  # This is the same for full world ranks belonging to the same worker.
+        if ranks_per_worker > 1:
+            # Create new communicators by splitting MPI.COMM_WORLD into group of sub-communicators based on
+            # input values `color` and `key`. `color` determines to which new communicator each processes will belong.
+            # `key` determines the ordering (rank) within each new communicator.
+            worker_sub_comm = MPI.COMM_WORLD.Split(
+                color=worker_idx, key=full_world_rank
             )
-        num_islands = island_sizes.size  # Determine number of islands.
 
-        #  Set up intra-island communicator for communication within each island.
-        intra_color = np.concatenate(
-            [idx * np.ones(el, dtype=int) for idx, el in enumerate(island_sizes)]
-        ).ravel()
-        island_idx = intra_color[
-            rank
-        ]  # Determine island index (which is also each rank's intra color).
-        intra_key = rank
+        else:
+            worker_sub_comm = MPI.COMM_SELF
 
-        # Determine displacements as positions of unique elements, where # unique elements equals number of islands.
-        _, island_displs = np.unique(intra_color, return_index=True)
+        # Create the Propulate world communicator, consisting of rank 0 of each worker's sub communicator.
+        worker_root_ranks = [
+            rank
+            for rank in list(range(full_world_size))
+            if rank % ranks_per_worker == 0
+        ]
+        propulate_world_group = MPI.COMM_WORLD.group.Incl(worker_root_ranks)
+        propulate_world_comm = MPI.COMM_WORLD.Create_group(propulate_world_group)
 
-        if rank == 0:
-            log.info(
-                f"Worker distribution {intra_color} with island counts "
-                f"{island_sizes} and island displacements {island_displs}."
+        # Make sure that the Propulate world communicator is only defined on rank 0 of each worker's sub communicator.
+        # Only those ranks are involved in the actual Propulate optimization logic and need to know about the related
+        # variables and settings.
+        if full_world_rank in worker_root_ranks:
+            propulate_world_rank, propulate_world_size = (
+                propulate_world_comm.rank,
+                propulate_world_comm.size,
             )
+        else:
+            propulate_world_comm = None
+
+        # Set up Propulate optimization logic only on ranks that are part of the Propulate world communicator.
+        if propulate_world_comm is not None:
+            # Homogeneous case with equal island sizes (differences of +-1 possible due to load balancing).
+            if island_sizes is None:
+                if num_islands < 1:
+                    raise ValueError(
+                        f"Invalid number of evolutionary islands, needs to be >= 1 but was {num_islands}."
+                    )
+                base_size = (
+                    propulate_world_size // num_islands
+                )  # Base number of workers of each island
+                remainder = (
+                    propulate_world_size % num_islands
+                )  # Number of remaining workers to be distributed
+                island_sizes = base_size * np.ones(num_islands, dtype=int)
+                island_sizes[
+                    :remainder
+                ] += 1  # Distribute remaining workers equally for balanced load.
+
+            # Heterogeneous case with user-defined island sizes.
+            if np.sum(island_sizes) != propulate_world_size:
+                raise ValueError(
+                    f"There should be COMM_WORLD.size = {propulate_world_size} workers "
+                    f"but only {np.sum(island_sizes)} were specified."
+                )
+            num_islands = island_sizes.size  # Determine number of islands.
+
+            #  Set up intra-island communicator for communication within each island.
+            island_colors = np.concatenate(
+                [idx * np.ones(el, dtype=int) for idx, el in enumerate(island_sizes)]
+            ).ravel()
+            island_idx = island_colors[
+                propulate_world_rank
+            ]  # Determine island index (which is also each rank's intra color).
+            island_key = propulate_world_rank
 
-        # Create new communicators by splitting MPI.COMM_WORLD into group of sub-communicators based on
-        # input values `color` and `key`. `color` determines to which new communicator each processes will belong.
-        # `key` determines the ordering (rank) within each new communicator.
-        comm_intra = MPI.COMM_WORLD.Split(color=island_idx, key=intra_key)
+            # Determine displacements as positions of unique elements, where # unique elements equals number of islands.
+            _, island_displs = np.unique(island_colors, return_index=True)
 
-        # Set up migration topology.
-        if migration_topology is None:
-            migration_topology = np.ones((num_islands, num_islands), dtype=int)
-            np.fill_diagonal(migration_topology, 0)  # No island self-talk.
-            if rank == 0:
+            if full_world_rank == 0:
                 log.info(
-                    "NOTE: No migration topology given, using fully connected top-1 topology."
+                    f"Worker distribution {island_colors} with island counts "
+                    f"{island_sizes} and island displacements {island_displs}."
                 )
 
-        if rank == 0:
-            log.info(
-                f"Migration topology {migration_topology} has shape {migration_topology.shape}."
-            )
+            island_comm = propulate_world_comm.Split(color=island_idx, key=island_key)
 
-        if migration_topology.shape != (num_islands, num_islands):
-            raise ValueError(
-                f"Migration topology must be a quadratic matrix of size "
-                f"{island_displs.size} x {island_displs.size} but has shape {migration_topology.shape}."
-            )
+            # Set up migration topology.
+            if migration_topology is None:
+                migration_topology = np.ones((num_islands, num_islands), dtype=int)
+                np.fill_diagonal(migration_topology, 0)  # No island self-talk.
+                if full_world_rank == 0:
+                    log.info(
+                        "NOTE: No migration topology given, using fully connected top-1 topology."
+                    )
 
-        if migration_probability > 1.0:
-            raise ValueError(
-                f"Migration probability must be in [0, 1] but was set to {migration_probability}."
-            )
-        migration_prob_rank = migration_probability / comm_intra.size
+            if full_world_rank == 0:
+                log.info(
+                    f"Migration topology {migration_topology} has shape {migration_topology.shape}."
+                )
 
-        if rank == 0:
-            log.info(
-                f"NOTE: Island migration probability {migration_probability} "
-                f"results in per-rank migration probability {migration_prob_rank}.\n"
-                "Starting parallel optimization process."
-            )
+            if migration_topology.shape != (num_islands, num_islands):
+                raise ValueError(
+                    f"Migration topology must be a quadratic matrix of size "
+                    f"{island_displs.size} x {island_displs.size} but has shape {migration_topology.shape}."
+                )
+
+            if migration_probability > 1.0:
+                raise ValueError(
+                    f"Migration probability must be in [0, 1] but was set to {migration_probability}."
+                )
+            migration_prob_rank = migration_probability / island_comm.size
+
+            if full_world_rank == 0:
+                log.info(
+                    f"NOTE: Island migration probability {migration_probability} "
+                    f"results in per-rank migration probability {migration_prob_rank}.\n"
+                    "Starting parallel optimization process."
+                )
+        else:
+            migration_prob_rank = None
+            island_displs = None
+            island_idx = None
+            island_comm = None
+            emigration_propagator = None
+            immigration_propagator = None
 
         MPI.COMM_WORLD.barrier()
         # Set up one Propulator for each island.
         if pollination is False:
-            if rank == 0:
+            if full_world_rank == 0:
                 log.info("Use island model with real migration.")
             self.propulator = Migrator(
                 loss_fn=loss_fn,
                 propagator=propagator,
+                rng=rng,
                 island_idx=island_idx,
-                comm=comm_intra,
+                island_comm=island_comm,
+                propulate_comm=propulate_world_comm,
+                worker_sub_comm=worker_sub_comm,
                 generations=generations,
                 checkpoint_path=checkpoint_path,
                 migration_topology=migration_topology,
@@ -212,16 +271,18 @@ def __init__(
                 emigration_propagator=emigration_propagator,
                 island_displs=island_displs,
                 island_counts=island_sizes,
-                rng=rng,
             )
         else:
-            if rank == 0:
+            if full_world_rank == 0:
                 log.info("Use island model with pollination.")
             self.propulator = Pollinator(
                 loss_fn=loss_fn,
                 propagator=propagator,
+                rng=rng,
                 island_idx=island_idx,
-                comm=comm_intra,
+                island_comm=island_comm,
+                propulate_comm=propulate_world_comm,
+                worker_sub_comm=worker_sub_comm,
                 generations=generations,
                 checkpoint_path=checkpoint_path,
                 migration_topology=migration_topology,
@@ -230,7 +291,6 @@ def __init__(
                 immigration_propagator=immigration_propagator,
                 island_displs=island_displs,
                 island_counts=island_sizes,
-                rng=rng,
             )
 
     def _run(
@@ -241,17 +301,17 @@ def _run(
 
         Parameters
         ----------
-        top_n: int
-               number of best results to report
-        logging_interval: int
-                          logging interval
-        debug: int
-               verbosity/debug level
+        top_n : int
+            The number of best results to report.
+        logging_interval : int
+            The logging interval.
+        debug : int
+            The debug level.
 
         Returns
         -------
-        list[list[Individual] | Individual]
-            top-n best individuals on each island
+        List[List[propulate.Individual] | propulate.Individual]
+            The top-n best individuals on each island.
         """
         self.propulator.propulate(logging_interval, debug)
         return self.propulator.summarize(top_n, debug)
@@ -260,20 +320,20 @@ def evolve(
         self, top_n: int = 3, logging_interval: int = 10, debug: int = 1
     ) -> List[Union[List[Individual], Individual]]:
         """
-        Run Propulate optimization.
+        Run Propulate optimization and summarize results.
 
         Parameters
         ----------
-        top_n: int
-               number of best results to report
-        logging_interval: int
-                          logging interval
-        debug: int
-               verbosity/debug level
+        top_n : int, optional
+            The number of best results to report. Default is 3.
+        logging_interval : int, optional
+            The logging interval. Default is 10.
+        debug : int
+            The debug level; 0 - silent; 1 - moderate, 2 - noisy (debug mode). Default is 1.
 
         Returns
         -------
-        list[list[Individual] | Individual]
-            top-n best individuals on each island
+        List[List[propulate.Individual] | propulate.Individual]
+            The top-n best individuals on each island.
         """
         return self._run(top_n, logging_interval, debug)
diff --git a/propulate/migrator.py b/propulate/migrator.py
index e16cefc9..121aeb9d 100644
--- a/propulate/migrator.py
+++ b/propulate/migrator.py
@@ -2,7 +2,7 @@
 import random
 import logging
 from pathlib import Path
-from typing import Callable, Union, Type
+from typing import Callable, Optional, Type, Union
 
 import numpy as np
 from mpi4py import MPI
@@ -21,63 +21,89 @@ class Migrator(Propulator):
 
     Individuals can only exist on one evolutionary island at a time, i.e., they are removed
     (i.e. deactivated for breeding) from the sending island upon emigration.
+
+    Attributes
+    ----------
+    emigrated : List[propulate.population.Individual]
+        A list of emigrated individuals to be deactivated on the sending island.
+
+    Methods
+    -------
+    propulate()
+        Run asynchronous evolutionary optimization routine with actual migration.
+
+    Notes
+    -----
+    The ``Migrator`` class inherits all methods and attributes from the ``Propulator`` class.
+
+    See Also
+    --------
+    :class:`Propulator` : The parent class.
     """
 
     def __init__(
         self,
         loss_fn: Callable,
         propagator: Propagator,
+        rng: random.Random,
         island_idx: int = 0,
-        comm: MPI.Comm = MPI.COMM_WORLD,
+        island_comm: MPI.Comm = MPI.COMM_WORLD,
+        propulate_comm: MPI.Comm = MPI.COMM_WORLD,
+        worker_sub_comm: MPI.Comm = MPI.COMM_SELF,
         generations: int = -1,
         checkpoint_path: Union[str, Path] = Path("./"),
-        migration_topology: np.ndarray = None,
+        migration_topology: Optional[np.ndarray] = None,
         migration_prob: float = 0.0,
         emigration_propagator: Type[Propagator] = SelectMin,
-        island_displs: np.ndarray = None,
-        island_counts: np.ndarray = None,
-        rng: random.Random = None,
+        island_displs: Optional[np.ndarray] = None,
+        island_counts: Optional[np.ndarray] = None,
     ) -> None:
         """
         Initialize ``Migrator`` with given parameters.
 
         Parameters
         ----------
-        loss_fn: Callable
-                 loss function to be minimized
+        loss_fn : Callable
+            The loss function to be minimized.
         propagator: propulate.propagators.Propagator
-                    propagator to apply for breeding
-        island_idx: int
-                    index of island
-        comm: MPI.Comm
-              intra-island communicator
-        generations: int
-                     number of generations to run
-        checkpoint_path: Union[Path, str]
-                         Path where checkpoints are loaded from and stored.
-        migration_topology: numpy.ndarray
-                            2D matrix where entry (i,j) specifies how many
-                            individuals are sent by island i to island j
-        migration_prob: float
-                        per-worker migration probability
-        emigration_propagator: type[propulate.propagators.Propagator]
-                               emigration propagator, i.e., how to choose individuals
-                               for emigration that are sent to destination island.
-                               Should be some kind of selection operator.
-        island_displs: numpy.ndarray
-                    array with MPI.COMM_WORLD rank of each island's worker 0
-                    Element i specifies MPI.COMM_WORLD rank of worker 0 on island with index i.
-        island_counts: numpy.ndarray
-                       array with number of workers per island
-                       Element i specifies number of workers on island with index i.
-        rng: random.Random
-             random number generator
+            The propagator to apply for breeding.
+        rng : random.Random
+            The separate random number generator for the Propulate optimization.
+        island_idx: int, optional
+            The island's index. Default is 0.
+        island_comm: MPI.Comm, optional
+            The intra-island communicator for communication within that island. Default is ``MPI.COMM_WORLD``.
+        propulate_comm : MPI.Comm, optional
+            The Propulate world communicator, consisting of rank 0 of each worker's sub communicator.
+            Default is ``MPI.COMM_WORLD``.
+        worker_sub_comm : MPI.Comm, optional
+            The sub communicator for each (multi rank) worker. Default is ``MPI.COMM_SELF``.
+        generations : int, optional
+            The number of generations to run. Default is -1, i.e., run into wall-clock time limit.
+        checkpoint_path : pathlib.Path | str, optional
+            The path where the checkpoints are loaded from and stored. Default is current working directory.
+        migration_topology : numpy.ndarray, optional
+            The migration topology, i.e., a 2D matrix where entry (i,j) specifies how many individuals are sent by
+            island i to island j.
+        migration_prob : float, optional
+            The per-worker migration probability. Default is 0.0.
+        emigration_propagator : Type[propulate.propagators.Propagator]
+            The emigration propagator, i.e., how to choose individuals for emigration that are sent to the destination
+            island. Should be some kind of selection operator. Default is ``SelectMin``.
+        island_displs : numpy.ndarray, optional
+            An array with ``propulate_comm`` rank of each island's worker 0. Element i specifies the rank of worker 0 on
+            island with index i in the Propulate communicator.
+        island_counts : numpy.ndarray, optional
+            An array with the number of workers per island. Element i specifies the number of workers on island i.
         """
         super().__init__(
             loss_fn,
             propagator,
+            rng,
             island_idx,
-            comm,
+            island_comm,
+            propulate_comm,
+            worker_sub_comm,
             generations,
             checkpoint_path,
             migration_topology,
@@ -85,16 +111,16 @@ def __init__(
             emigration_propagator,
             island_displs,
             island_counts,
-            rng,
         )
         # Set class attributes.
-        self.emigrated = []  # emigrated individuals to be deactivated on sending island
+        self.emigrated = []  # Emigrated individuals to be deactivated on sending island
 
     def _send_emigrants(self) -> None:
-        """
-        Perform migration, i.e. island sends individuals out to other islands.
-        """
-        log_string = f"Island {self.island_idx} Worker {self.comm.rank} Generation {self.generation}: EMIGRATION\n"
+        """Perform migration, i.e. island sends individuals out to other islands."""
+        log_string = (
+            f"Island {self.island_idx} Worker {self.island_comm.rank} "
+            f"Generation {self.generation}: EMIGRATION\n"
+        )
         # Determine relevant line of migration topology.
         to_migrate = self.migration_topology[self.island_idx, :]
         num_emigrants = np.sum(
@@ -103,7 +129,7 @@ def _send_emigrants(self) -> None:
         eligible_emigrants = [
             ind
             for ind in self.population
-            if ind.active and ind.current == self.comm.rank
+            if ind.active and ind.current == self.island_comm.rank
         ]
 
         # Only perform migration if overall number of emigrants to be sent
@@ -122,7 +148,7 @@ def _send_emigrants(self) -> None:
             for target_island, offspring in enumerate(to_migrate):
                 if offspring == 0:
                     continue
-                # Determine MPI.COMM_WORLD ranks of workers on target island.
+                # Determine self.propulate_comm ranks of workers on target island.
                 displ = self.island_displs[target_island]
                 count = self.island_counts[target_island]
                 dest_island = np.arange(displ, displ + count)
@@ -136,11 +162,11 @@ def _send_emigrants(self) -> None:
 
                 # Deactivate emigrants on sending island (true migration).
                 for r in range(
-                    self.comm.size
+                    self.island_comm.size
                 ):  # Send emigrants to other intra-island workers for deactivation.
-                    if r == self.comm.rank:
+                    if r == self.island_comm.rank:
                         continue  # No self-talk.
-                    self.comm.send(
+                    self.island_comm.send(
                         copy.deepcopy(emigrants), dest=r, tag=SYNCHRONIZATION_TAG
                     )
                     log_string += (
@@ -153,8 +179,10 @@ def _send_emigrants(self) -> None:
                 # Determine new responsible worker on target island.
                 for ind in departing:
                     ind.current = self.rng.randrange(0, count)
-                for r in dest_island:  # Loop over MPI.COMM_WORLD destination ranks.
-                    MPI.COMM_WORLD.send(
+                for (
+                    r
+                ) in dest_island:  # Loop over self.propulate_comm destination ranks.
+                    self.propulate_comm.send(
                         copy.deepcopy(departing), dest=r, tag=MIGRATION_TAG
                     )
                     log_string += (
@@ -190,7 +218,7 @@ def _send_emigrants(self) -> None:
 
         else:
             log.debug(
-                f"Island {self.island_idx} worker {self.comm.rank} generation {self.generation}: \n"
+                f"Island {self.island_idx} worker {self.island_comm.rank} generation {self.generation}: \n"
                 f"Population size {len(eligible_emigrants)} too small "
                 f"to select {num_emigrants} migrants."
             )
@@ -204,16 +232,19 @@ def _receive_immigrants(self) -> None:
         RuntimeError
             If identical immigrant is already active on target island for real migration.
         """
-        log_string = f"Island {self.island_idx} Worker {self.comm.rank} Generation {self.generation}: IMMIGRATION\n"
+        log_string = (
+            f"Island {self.island_idx} Worker {self.island_comm.rank} "
+            f"Generation {self.generation}: IMMIGRATION\n"
+        )
         probe_migrants = True
         while probe_migrants:
             stat = MPI.Status()
-            probe_migrants = MPI.COMM_WORLD.iprobe(
+            probe_migrants = self.propulate_comm.iprobe(
                 source=MPI.ANY_SOURCE, tag=MIGRATION_TAG, status=stat
             )
             log_string += f"Immigrant(s) to receive?...{probe_migrants}\n"
             if probe_migrants:
-                immigrants = MPI.COMM_WORLD.recv(
+                immigrants = self.propulate_comm.recv(
                     source=stat.Get_source(), tag=MIGRATION_TAG
                 )
                 log_string += (
@@ -282,9 +313,9 @@ def _check_emigrants_to_deactivate(self) -> bool:
                         break
 
                 log.info(
-                    f"Island {self.island_idx} Worker {self.comm.rank} Generation {self.generation}:\n"
+                    f"Island {self.island_idx} Worker {self.island_comm.rank} Generation {self.generation}:\n"
                     f"Currently in emigrated: {emigrant}\n"
-                    f"Island {self.island_idx} Worker {self.comm.rank} Generation {self.generation}: "
+                    f"Island {self.island_idx} Worker {self.island_comm.rank} Generation {self.generation}: "
                     f"Currently in population: {existing_ind}\nEquivalence check: {existing_ind[0] == emigrant} "
                     f"{compare_traits} {existing_ind[0].loss == self.emigrated[idx].loss} "
                     f"{existing_ind[0].active == emigrant.active} {existing_ind[0].current == emigrant.current} "
@@ -296,20 +327,21 @@ def _check_emigrants_to_deactivate(self) -> bool:
         return check
 
     def _deactivate_emigrants(self) -> None:
-        """
-        Check for and possibly receive emigrants from other intra-island workers to be deactivated.
-        """
-        log_string = f"Island {self.island_idx} Worker {self.comm.rank} Generation {self.generation}: DEACTIVATION\n"
+        """Check for and possibly receive emigrants from other intra-island workers to be deactivated."""
+        log_string = (
+            f"Island {self.island_idx} Worker {self.island_comm.rank} "
+            f"Generation {self.generation}: DEACTIVATION\n"
+        )
         probe_sync = True
         while probe_sync:
             stat = MPI.Status()
-            probe_sync = self.comm.iprobe(
+            probe_sync = self.island_comm.iprobe(
                 source=MPI.ANY_SOURCE, tag=SYNCHRONIZATION_TAG, status=stat
             )
             log_string += f"Emigrants from others to be deactivated to be received?...{probe_sync}\n"
             if probe_sync:
                 # Receive new emigrants.
-                new_emigrants = self.comm.recv(
+                new_emigrants = self.island_comm.recv(
                     source=stat.Get_source(), tag=SYNCHRONIZATION_TAG
                 )
                 # Add new emigrants to list of emigrants to be deactivated.
@@ -356,35 +388,45 @@ def _deactivate_emigrants(self) -> None:
         )
         log.debug(log_string)
 
-    def _work(self, logging_interval: int, debug: int):
+    def _work(self, logging_interval: int = 10, debug: int = 1):
         """
         Execute evolutionary algorithm using island model with real migration in parallel.
 
         Parameters
         ----------
-        logging_interval: int
-                          logging interval
-        debug: int
-               verbosity/debug level; 0 - silent; 1 - moderate, 2 - noisy (debug mode)
+        logging_interval : int, optional
+            The logging interval. Default is 10.
+        debug : int
+            The debug level; 0 - silent; 1 - moderate, 2 - noisy (debug mode). Default is 1.
 
         Raises
         ------
         ValueError
             If any individuals are left that should have been deactivated before (only for debug > 0).
         """
+        if self.worker_sub_comm != MPI.COMM_SELF:
+            self.generation = self.worker_sub_comm.bcast(self.generation, root=0)
+        if self.propulate_comm is None:
+            while self.generations <= -1 or self.generation < self.generations:
+                # Breed and evaluate individual.
+                # print(self.generation)
+
+                self._evaluate_individual()
+                self.generation += 1
+            return
 
-        if self.comm.rank == 0:
-            log.info(f"Island {self.island_idx} has {self.comm.size} workers.")
+        if self.island_comm.rank == 0:
+            log.info(f"Island {self.island_idx} has {self.island_comm.size} workers.")
 
-        dump = True if self.comm.rank == 0 else False
+        dump = True if self.island_comm.rank == 0 else False
         migration = True if self.migration_prob > 0 else False
-        MPI.COMM_WORLD.barrier()
+        self.propulate_comm.barrier()
 
         # Loop over generations.
         while self.generations <= -1 or self.generation < self.generations:
             if self.generation % int(logging_interval) == 0:
                 log.info(
-                    f"Island {self.island_idx} Worker {self.comm.rank}: In generation {self.generation}..."
+                    f"Island {self.island_idx} Worker {self.island_comm.rank}: In generation {self.generation}..."
                 )
 
             # Breed and evaluate individual.
@@ -422,20 +464,20 @@ def _work(self, logging_interval: int, debug: int):
         # so that each of them holds the complete final population and the found optimum
         # irrespective of the order they finished.
 
-        MPI.COMM_WORLD.barrier()
-        if MPI.COMM_WORLD.rank == 0:
+        self.propulate_comm.barrier()
+        if self.propulate_comm.rank == 0:
             log.info("OPTIMIZATION DONE.")
             log.info("NEXT: Final checks for incoming messages...")
-        MPI.COMM_WORLD.barrier()
+        self.propulate_comm.barrier()
 
         # Final check for incoming individuals evaluated by other intra-island workers.
         self._receive_intra_island_individuals()
-        MPI.COMM_WORLD.barrier()
+        self.propulate_comm.barrier()
 
         if migration:
             # Final check for incoming individuals from other islands.
             self._receive_immigrants()
-            MPI.COMM_WORLD.barrier()
+            self.propulate_comm.barrier()
 
             # Emigration: Final check for emigrants from other intra-island workers to be deactivated.
             self._deactivate_emigrants()
@@ -443,10 +485,10 @@ def _work(self, logging_interval: int, debug: int):
             if debug > 0:
                 check = self._check_emigrants_to_deactivate()
                 assert check is False
-                MPI.COMM_WORLD.barrier()
+                self.propulate_comm.barrier()
                 if len(self.emigrated) > 0:
                     log.info(
-                        f"Island {self.island_idx} Worker {self.comm.rank} Generation {self.generation}: "
+                        f"Island {self.island_idx} Worker {self.island_comm.rank} Generation {self.generation}: "
                         f"Finally {len(self.emigrated)} individual(s) in emigrated: {self.emigrated}:\n"
                         f"{self.population}"
                     )
@@ -456,11 +498,11 @@ def _work(self, logging_interval: int, debug: int):
                             "There should not be any individuals left that need to be deactivated."
                         )
 
-            MPI.COMM_WORLD.barrier()
+            self.propulate_comm.barrier()
 
         # Final checkpointing on rank 0.
-        if self.comm.rank == 0:
+        if self.island_comm.rank == 0:
             self._dump_final_checkpoint()  # Dump checkpoint.
-        MPI.COMM_WORLD.barrier()
+        self.propulate_comm.barrier()
         _ = self._determine_worker_dumping_next()
-        MPI.COMM_WORLD.barrier()
+        self.propulate_comm.barrier()
diff --git a/propulate/pollinator.py b/propulate/pollinator.py
index 059e03d4..fa3a6ac5 100644
--- a/propulate/pollinator.py
+++ b/propulate/pollinator.py
@@ -3,7 +3,7 @@
 import random
 import time
 from pathlib import Path
-from typing import Callable, Union, Tuple, List, Type
+from typing import Callable, List, Optional, Tuple, Type, Union
 
 import numpy as np
 from mpi4py import MPI
@@ -20,72 +20,97 @@ class Pollinator(Propulator):
     """
     Parallel propagator of populations with pollination.
 
-    Individuals can actively exist on multiple evolutionary islands at a time, i.e.,
-    copies of emigrants are sent out and emigrating individuals are not deactivated on
-    sending island for breeding. Instead, immigrants replace individuals on the target
-    island according to an immigration policy set by the immigration propagator.
+    Individuals can actively exist on multiple evolutionary islands at a time, i.e., copies of emigrants are sent out
+    and emigrating individuals are not deactivated on sending island for breeding. Instead, immigrants replace
+    individuals on the target island according to an immigration policy set by the immigration propagator.
+
+    Attributes
+    ----------
+    immigration_propagator : Type[propulate.propagators.Propagator]
+        The immigration propagator, i.e., how to choose individuals to be replaced by immigrants on target island.
+    replaced : List[propulate.population.Individual]
+        The individuals to be replaced by the immigrants.
+
+    Methods
+    -------
+    propulate()
+        Run asynchronous evolutionary optimization routine with pollination.
+
+    Notes
+    -----
+    The ``Pollinator`` class inherits all methods and attributes from the ``Propulator`` class.
+
+    See Also
+    --------
+    :class:`Propulator` : The parent class.
     """
 
     def __init__(
         self,
         loss_fn: Callable,
         propagator: Propagator,
+        rng: random.Random,
         island_idx: int = 0,
-        comm: MPI.Comm = MPI.COMM_WORLD,
+        island_comm: MPI.Comm = MPI.COMM_WORLD,
+        propulate_comm: MPI.Comm = MPI.COMM_WORLD,
+        worker_sub_comm: MPI.Comm = MPI.COMM_SELF,
         generations: int = 0,
         checkpoint_path: Union[Path, str] = Path("./"),
-        migration_topology: np.ndarray = None,
+        migration_topology: Optional[np.ndarray] = None,
         migration_prob: float = 0.0,
         emigration_propagator: Type[Propagator] = SelectMin,
         immigration_propagator: Type[Propagator] = SelectMax,
-        island_displs: np.ndarray = None,
-        island_counts: np.ndarray = None,
-        rng: random.Random = None,
+        island_displs: Optional[np.ndarray] = None,
+        island_counts: Optional[np.ndarray] = None,
     ) -> None:
         """
         Initialize ``Pollinator`` with given parameters.
 
         Parameters
         ----------
-        loss_fn: Callable
-                 loss function to be minimized
-        propagator: propulate.propagators.Propagator
-                    propagator to apply for breeding
-        island_idx: int
-                    index of island
-        comm: MPI.Comm
-              intra-island communicator
-        generations: int
-                     number of generations to run
-        checkpoint_path: Path
-                         path where checkpoints are loaded from and stored.
-        migration_topology: numpy.ndarray
-                            2D matrix where entry (i,j) specifies how many
-                            individuals are sent by island i to island j
-        migration_prob: float
-                        per-worker migration probability
-        emigration_propagator: type[propulate.propagators.Propagator]
-                               emigration propagator, i.e., how to choose individuals
-                               for emigration that are sent to destination island.
-                               Should be some kind of selection operator.
-        immigration_propagator: type[propulate.propagators.Propagator]
-                                immigration propagator, i.e., how to choose individuals
-                                to be replaced by immigrants on target island.
-                                Should be some kind of selection operator.
-        island_displs: numpy.ndarray
-                    array with MPI.COMM_WORLD rank of each island's worker 0
-                    Element i specifies MPI.COMM_WORLD rank of worker 0 on island with index i.
-        island_counts: numpy.ndarray
-                       array with number of workers per island
-                       Element i specifies number of workers on island with index i.
-        rng: random.Random
-             random number generator
+        loss_fn : Callable
+            The loss function to be minimized.
+        propagator : propulate.propagators.Propagator
+            The propagator to apply for breeding.
+        rng : random.Random
+            The separate random number generator for the Propulate optimization.
+        island_idx : int, optional
+            The island index. Default is 0.
+        island_comm : MPI.Comm, optional
+            The intra-island communicator for communication within this island. Default is ``MPI.COMM_WORLD``.
+        propulate_comm : MPI.Comm, optional
+            The Propulate world communicator, consisting of rank 0 of each worker's sub communicator.
+            Default is ``MPI.COMM_WORLD``.
+        worker_sub_comm : MPI.Comm, optional
+            The sub communicator for each (multi rank) worker. Default is ``MPI.COMM_SELF``.
+        generations : int, optional
+            The number of generations to run. Default is -1, i.e., run into wall-clock time limit.
+        checkpoint_path : pathlib.Path, optional
+            The path where checkpoints are loaded from and stored. Default is current working directory.
+        migration_topology : numpy.ndarray, optional
+            The 2D matrix where entry (i,j) specifies how many individuals are sent by island i to island j.
+        migration_prob : float, optional
+            The per-worker migration probability. Default is 0.0.
+        emigration_propagator : Type[propulate.propagators.Propagator], optional
+            The emigration propagator, i.e., how to choose individuals for emigration that are sent to the destination
+            island. Should be some kind of selection operator. Default is ``SelectMin``.
+        immigration_propagator : Type[propulate.propagators.Propagator], optional
+            The immigration propagator, i.e., how to choose individuals to be replaced by immigrants on a target island.
+            Should be some kind of selection operator. Default is ``SelectMax``.
+        island_displs : numpy.ndarray, optional
+            An array with ``propulate_comm`` rank of each island's worker 0. Element i specifies the rank of worker 0 on
+            island with index i in the Propulate communicator.
+        island_counts : numpy.ndarray, optional
+            An array with the number of workers per island. Element i specifies the number of workers on island i.
         """
         super().__init__(
             loss_fn,
             propagator,
+            rng,
             island_idx,
-            comm,
+            island_comm,
+            propulate_comm,
+            worker_sub_comm,
             generations,
             checkpoint_path,
             migration_topology,
@@ -93,35 +118,35 @@ def __init__(
             emigration_propagator,
             island_displs,
             island_counts,
-            rng,
         )
         # Set class attributes.
-        self.immigration_propagator = immigration_propagator  # immigration propagator
-        self.replaced = []  # individuals to be replaced by immigrants
+        self.immigration_propagator = immigration_propagator  # Immigration propagator
+        self.replaced = []  # Individuals to be replaced by immigrants
 
     def _send_emigrants(self) -> None:
-        """
-        Perform migration, i.e. island sends individuals out to other islands.
-        """
-        log_string = f"Island {self.island_idx} Worker {self.comm.rank} Generation {self.generation}: EMIGRATION\n"
+        """Perform migration, i.e. island sends individuals out to other islands."""
+        log_string = (
+            f"Island {self.island_idx} Worker {self.island_comm.rank} "
+            f"Generation {self.generation}: EMIGRATION\n"
+        )
         # Determine relevant line of migration topology.
         to_migrate = self.migration_topology[self.island_idx, :]
         num_emigrants = np.amax(
             to_migrate
         )  # Determine maximum number of emigrants to be sent out at once.
-        # NOTE For pollination, emigration-responsible worker not necessary as emigrating individuals are
-        # not deactivated and copies are allowed.
+        # For pollination, an emigration-responsible worker is not necessary as emigrating individuals are not
+        # deactivated and copies are allowed.
         # All active individuals are eligible emigrants.
         eligible_emigrants, _ = self._get_active_individuals()
 
-        # Only perform migration if maximum number of emigrants to be sent
-        # out at once is smaller than current number of eligible emigrants.
+        # Only perform migration if maximum number of emigrants to be sent out at once is smaller than current number
+        # of eligible emigrants.
         if num_emigrants <= len(eligible_emigrants):
             # Loop through relevant part of migration topology.
             for target_island, offspring in enumerate(to_migrate):
                 if offspring == 0:
                     continue
-                # Determine MPI.COMM_WORLD ranks of workers on target island.
+                # Determine ranks of workers on target island in the Propulate world communicator.
                 displ = self.island_displs[target_island]
                 count = self.island_counts[target_island]
                 dest_island = np.arange(displ, displ + count)
@@ -149,8 +174,8 @@ def _send_emigrants(self) -> None:
                     log_string += (
                         f"{ind} with migration history {ind.migration_history}\n"
                     )
-                for r in dest_island:  # Loop through MPI.COMM_WORLD destination ranks.
-                    MPI.COMM_WORLD.send(
+                for r in dest_island:  # Loop through Propulate world destination ranks.
+                    self.propulate_comm.send(
                         copy.deepcopy(departing), dest=r, tag=MIGRATION_TAG
                     )
                     log_string += (
@@ -166,26 +191,27 @@ def _send_emigrants(self) -> None:
 
         else:
             log.debug(
-                f"Island {self.island_idx} Worker {self.comm.rank} Generation {self.generation}: "
+                f"Island {self.island_idx} Worker {self.island_comm.rank} Generation {self.generation}: "
                 f"Population size {len(eligible_emigrants)} too small "
                 f"to select {num_emigrants} migrants."
             )
 
     def _receive_immigrants(self) -> None:
-        """
-        Check for and possibly receive immigrants send by other islands.
-        """
+        """Check for and possibly receive immigrants send by other islands."""
         replace_num = 0
-        log_string = f"Island {self.island_idx} Worker {self.comm.rank} Generation {self.generation}: IMMIGRATION\n"
+        log_string = (
+            f"Island {self.island_idx} Worker {self.island_comm.rank} "
+            f"Generation {self.generation}: IMMIGRATION\n"
+        )
         probe_migrants = True
         while probe_migrants:
             stat = MPI.Status()
-            probe_migrants = MPI.COMM_WORLD.iprobe(
+            probe_migrants = self.propulate_comm.iprobe(
                 source=MPI.ANY_SOURCE, tag=MIGRATION_TAG, status=stat
             )
             log_string += f"Immigrant(s) to receive?...{probe_migrants}\n"
             if probe_migrants:
-                immigrants = MPI.COMM_WORLD.recv(
+                immigrants = self.propulate_comm.recv(
                     source=stat.Get_source(), tag=MIGRATION_TAG
                 )
                 log_string += (
@@ -202,7 +228,7 @@ def _receive_immigrants(self) -> None:
                     )  # Append immigrant to population.
 
                     replace_num = 0
-                    if self.comm.rank == immigrant.current:
+                    if self.island_comm.rank == immigrant.current:
                         replace_num += 1
                     log_string += (
                         f"Responsible for choosing {replace_num} individual(s) "
@@ -213,11 +239,10 @@ def _receive_immigrants(self) -> None:
                 # cannot choose the same individual independently for replacement and thus deactivation.
                 if replace_num > 0:
                     # From current population, choose `replace_num` individuals to be replaced.
-                    # eligible_for_replacement = [ind for ind in self.population[:-len(immigrants)] if ind.active \
                     eligible_for_replacement = [
                         ind
                         for ind in self.population
-                        if ind.active and ind.current == self.comm.rank
+                        if ind.active and ind.current == self.island_comm.rank
                     ]
 
                     immigrator = self.immigration_propagator(
@@ -228,10 +253,10 @@ def _receive_immigrants(self) -> None:
                     )  # Choose individual to be replaced by immigrant.
 
                     # Send individuals to be replaced to other intra-island workers for deactivation.
-                    for r in range(self.comm.size):
-                        if r == self.comm.rank:
+                    for r in range(self.island_comm.size):
+                        if r == self.island_comm.rank:
                             continue  # No self-talk.
-                        self.comm.send(
+                        self.island_comm.send(
                             copy.deepcopy(to_replace), dest=r, tag=SYNCHRONIZATION_TAG
                         )
                         log_string += (
@@ -245,27 +270,25 @@ def _receive_immigrants(self) -> None:
                         individual.active = False
 
         _, num_active = self._get_active_individuals()
-        log_string += (
-            f"After immigration: {num_active}/{len(self.population)} active.\n"
-        )
+        log_string += f"After immigration: {num_active}/{len(self.population)} active."
         log.debug(log_string)
 
     def _deactivate_replaced_individuals(self) -> None:
-        """
-        Check for and possibly receive individuals from other intra-island workers to be deactivated
-        because of immigration.
-        """
-        log_string = f"Island {self.island_idx} Worker {self.comm.rank} Generation {self.generation}: REPLACEMENT\n"
+        """Check for and receive individuals from other intra-island workers to be deactivated due to immigration."""
+        log_string = (
+            f"Island {self.island_idx} Worker {self.island_comm.rank} "
+            f"Generation {self.generation}: REPLACEMENT\n"
+        )
         probe_sync = True
         while probe_sync:
             stat = MPI.Status()
-            probe_sync = self.comm.iprobe(
+            probe_sync = self.island_comm.iprobe(
                 source=MPI.ANY_SOURCE, tag=SYNCHRONIZATION_TAG, status=stat
             )
             log_string += f"Individual(s) to replace...{probe_sync}\n"
             if probe_sync:
                 # Receive new individuals.
-                to_replace = self.comm.recv(
+                to_replace = self.island_comm.recv(
                     source=stat.Get_source(), tag=SYNCHRONIZATION_TAG
                 )
                 # Add new emigrants to list of emigrants to be deactivated.
@@ -309,27 +332,26 @@ def _deactivate_replaced_individuals(self) -> None:
         log.debug(log_string)
 
     def _check_for_duplicates(
-        self, active: bool, debug: int
+        self, active: bool, debug: int = 1
     ) -> Tuple[List[List[Union[Individual, int]]], List[Individual]]:
         """
         Check for duplicates in current population.
 
-        For pollination, duplicates are allowed as emigrants are sent as copies
-        and not deactivated on sending island.
+        For pollination, duplicates are allowed as emigrants are sent as copies and not deactivated on sending island.
 
         Parameters
         ----------
-        active: bool
-                Whether to consider active individuals (True) or all individuals (False)
-        debug: int
-               verbosity/debug level; 0 - silent; 1 - moderate, 2 - noisy (debug mode)
+        active : bool
+            Whether to consider active individuals (True) or all individuals (False).
+        debug : int, optional
+            The debug level; 0 - silent; 1 - moderate, 2 - noisy (debug mode). Default is 1.
 
         Returns
         -------
-        list[list[propulate.individual.Individual | int]]
-            individuals and their occurrences
-        list[propulate.individual.Individual]
-            unique individuals in population
+        List[List[propulate.individual.Individual | int]]
+            The individuals and their occurrences.
+        List[propulate.individual.Individual]
+            All unique individuals in the population.
         """
         if active:
             population, _ = self._get_active_individuals()
@@ -351,36 +373,44 @@ def _check_for_duplicates(
             if not considered:
                 num_copies = population.count(individual)
                 log.debug(
-                    f"Island {self.island_idx} Worker {self.comm.rank} Generation {self.generation}: "
+                    f"Island {self.island_idx} Worker {self.island_comm.rank} Generation {self.generation}: "
                     f"{individual} occurs {num_copies} time(s)."
                 )
                 unique_inds.append(individual)
                 occurrences.append([individual, num_copies])
         return occurrences, unique_inds
 
-    def _work(self, logging_interval: int, debug: int):
+    def _work(self, logging_interval: int = 10, debug: int = 1):
         """
         Execute evolutionary algorithm using island model with pollination in parallel.
 
         Parameters
         ----------
-        logging_interval: int
-                          logging interval
-        debug: int
-               verbosity/debug level; 0 - silent; 1 - moderate, 2 - noisy (debug mode)
+        logging_interval : int, optional
+            The logging interval. Default is 10.
+        debug : int, optional
+            The debug level; 0 - silent; 1 - moderate, 2 - noisy (debug mode). Default is 1.
         """
-        if self.comm.rank == 0:
-            log.info(f"Island {self.island_idx} has {self.comm.size} workers.")
-
-        dump = True if self.comm.rank == 0 else False
+        if self.worker_sub_comm != MPI.COMM_SELF:
+            self.generation = self.worker_sub_comm.bcast(self.generation, root=0)
+        if self.propulate_comm is None:
+            while self.generations <= -1 or self.generation < self.generations:
+                # Breed and evaluate individual.
+                self._evaluate_individual()
+                self.generation += 1
+            return
+        if self.island_comm.rank == 0:
+            log.info(f"Island {self.island_idx} has {self.island_comm.size} workers.")
+
+        dump = True if self.island_comm.rank == 0 else False
         migration = True if self.migration_prob > 0 else False
-        MPI.COMM_WORLD.barrier()
+        self.propulate_comm.barrier()
 
         # Loop over generations.
         while self.generations <= -1 or self.generation < self.generations:
             if debug == 1 and self.generation % int(logging_interval) == 0:
                 log.info(
-                    f"Island {self.island_idx} Worker {self.comm.rank}: In generation {self.generation}..."
+                    f"Island {self.island_idx} Worker {self.island_comm.rank}: In generation {self.generation}..."
                 )
 
             # Breed and evaluate individual.
@@ -414,37 +444,35 @@ def _work(self, logging_interval: int, debug: int):
         # so that each of them holds the complete final population and the found optimum
         # irrespective of the order they finished.
 
-        MPI.COMM_WORLD.barrier()
-        if MPI.COMM_WORLD.rank == 0:
-            log.info("OPTIMIZATION DONE.")
-            log.info("NEXT: Final checks for incoming messages...")
-        MPI.COMM_WORLD.barrier()
+        self.propulate_comm.barrier()
+        if self.propulate_comm.rank == 0:
+            log.info("OPTIMIZATION DONE.\nNEXT: Final checks for incoming messages...")
+        self.propulate_comm.barrier()
 
         # Final check for incoming individuals evaluated by other intra-island workers.
         self._receive_intra_island_individuals()
-        MPI.COMM_WORLD.barrier()
+        self.propulate_comm.barrier()
 
         if migration:
             # Final check for incoming individuals from other islands.
             self._receive_immigrants()
-            MPI.COMM_WORLD.barrier()
+            self.propulate_comm.barrier()
 
             # Immigration: Final check for individuals replaced by other intra-island workers to be deactivated.
             self._deactivate_replaced_individuals()
-            MPI.COMM_WORLD.barrier()
+            self.propulate_comm.barrier()
 
             if len(self.replaced) > 0:
                 log.info(
-                    f"Island {self.island_idx} Worker {self.comm.rank} Generation {self.generation}: "
-                    f"Finally {len(self.replaced)} individual(s) in replaced: {self.replaced}:\n"
-                    f"{self.population}"
+                    f"Island {self.island_idx} Worker {self.island_comm.rank} Generation {self.generation}: "
+                    f"Finally {len(self.replaced)} individual(s) in replaced: {self.replaced}:\n{self.population}"
                 )
                 self._deactivate_replaced_individuals()
-            MPI.COMM_WORLD.barrier()
+            self.propulate_comm.barrier()
 
         # Final checkpointing on rank 0.
-        if self.comm.rank == 0:
+        if self.island_comm.rank == 0:
             self._dump_final_checkpoint()  # Dump checkpoint.
-        MPI.COMM_WORLD.barrier()
+        self.propulate_comm.barrier()
         _ = self._determine_worker_dumping_next()
-        MPI.COMM_WORLD.barrier()
+        self.propulate_comm.barrier()
diff --git a/propulate/propagators/__init__.py b/propulate/propagators/__init__.py
index d852105d..1293e952 100644
--- a/propulate/propagators/__init__.py
+++ b/propulate/propagators/__init__.py
@@ -1,6 +1,4 @@
-"""
-This package bundles all classes that are used as propagators in Propulate's optimization routine.
-"""
+"""This package bundles all classes that are used as propagators in Propulate's optimization routine."""
 
 from .base import (
     Propagator,
@@ -17,9 +15,9 @@
     PointMutation,
     RandomPointMutation,
     IntervalMutationNormal,
-    MateUniform,
-    MateMultiple,
-    MateSigmoid,
+    CrossoverUniform,
+    CrossoverMultiple,
+    CrossoverSigmoid,
 )
 
 from .pso import (
@@ -50,9 +48,9 @@
     "PointMutation",
     "RandomPointMutation",
     "IntervalMutationNormal",
-    "MateUniform",
-    "MateMultiple",
-    "MateSigmoid",
+    "CrossoverUniform",
+    "CrossoverMultiple",
+    "CrossoverSigmoid",
     "BasicPSO",
     "VelocityClampingPSO",
     "ConstrictionPSO",
diff --git a/propulate/propagators/base.py b/propulate/propagators/base.py
index 7b98397d..a5402731 100644
--- a/propulate/propagators/base.py
+++ b/propulate/propagators/base.py
@@ -1,5 +1,5 @@
 import random
-from typing import List, Dict, Union, Tuple
+from typing import List, Dict, Optional, Tuple, Union
 
 from ..population import Individual
 
@@ -11,14 +11,14 @@ def _check_compatible(out1: int, in2: int) -> bool:
     Parameters
     ----------
     out1: int
-        number of output individuals returned by first propagator
+        The number of output individuals returned by the first propagator.
     in2: int
-        number of input individuals taken by second propagator
+        The number of input individuals taken by the second propagator.
 
     Returns
     -------
     bool
-        True if propagators can be stacked, False if not.
+        True if the input propagators can be stacked: False if not.
     """
     return out1 == in2 or in2 == -1
 
@@ -28,22 +28,36 @@ class Propagator:
     Abstract base class for all propagators, i.e., evolutionary operators.
 
     A propagator takes a collection of individuals and uses them to breed a new collection of individuals.
+
+    Attributes
+    ----------
+    offspring : int
+        The number of output individuals to breed.
+    parents : int
+        The number of input individuals to use as parents.
+    rng : random.Random
+        The separate random number generator for the Propulate optimization.
+
+    Methods
+    -------
+    __call__()
+        Apply the propagator.
     """
 
     def __init__(
-        self, parents: int = 0, offspring: int = 0, rng: random.Random = None
+        self, parents: int = 0, offspring: int = 0, rng: Optional[random.Random] = None
     ) -> None:
         """
         Initialize a propagator with given parameters.
 
         Parameters
         ----------
-        parents: int
-            number of input individuals (-1 for any)
-        offspring: int
-            number of output individuals
-        rng: random.Random
-            random number generator
+        parents : int, optional
+            The number of input individuals (-1 for any). Default is 0 for abstract base class.
+        offspring : int, optional
+            The number of output individuals to breed. Default is 0 for abstract base class.
+        rng : random.Random, optional
+            The separate random number generator for the Propulate optimization.
 
         Raises
         ------
@@ -54,6 +68,8 @@ def __init__(
             raise ValueError("Propagator has to sire more than 0 offspring.")
         self.offspring = offspring  # Number of offspring individuals to breed
         self.parents = parents  # Number of parent individuals
+        if rng is None:
+            rng = random.Random()
         self.rng = rng  # Random number generator
 
     def __call__(self, inds: List[Individual]) -> Union[List[Individual], Individual]:
@@ -62,13 +78,13 @@ def __call__(self, inds: List[Individual]) -> Union[List[Individual], Individual
 
         Parameters
         ----------
-        inds: list[propulate.individual.Individual]
-            input individuals the propagator is applied to
+        inds : List[propulate.Individual]
+            The input individuals the propagator is applied to.
 
         Returns
         -------
-        list[Individual] | Individual
-            individual(s) bred by applying the propagator
+        List[propulate.Individual] | propulate.Individual
+            The individual(s) bred by applying the propagator.
             While this abstract base class method actually returns ``None``, each concrete child class of ``Propagator``
             should return an ``Individual`` instance or a list of them.
 
@@ -85,6 +101,19 @@ class Stochastic(Propagator):
     Apply a propagator with a given probability.
 
     If the propagator is not applied, the output still has to adhere to the defined number of offspring.
+
+    Attributes
+    ----------
+    probability : float
+        The probability of applying the propagator.
+
+    Notes
+    -----
+    The ``Stochastic`` class inherits all methods and attributes from the ``Propagator`` class.
+
+    See Also
+    --------
+    :class:`Propagator` : The parent class.
     """
 
     def __init__(
@@ -99,21 +128,21 @@ def __init__(
 
         Parameters
         ----------
-        parents: int
-            number of input individuals (-1 for any)
-        offspring: int
-            number of output individuals
-        probability: float
-            probability of application
-        rng: random.Random
-            random number generator
+        parents : int, optional
+            The number of input individuals (-1 for any). Default is 0.
+        offspring : int, optional
+            The number of output individuals. Default is 0.
+        probability : float, optional
+            The probability of application. Default is 0.0.
+        rng : random.Random, optional
+            The separate random number generator for the Propulate optimization.
 
         Raises
         ------
         ValueError
             If the number of offspring to breed is zero.
         """
-        super(Stochastic, self).__init__(parents, offspring, rng)
+        super().__init__(parents, offspring, rng)
         self.probability = probability
         if offspring == 0:
             raise ValueError("Propagator has to sire more than 0 offspring.")
@@ -125,6 +154,23 @@ class Conditional(Propagator):
 
     If the population consists of the specified number of individuals required for breeding (or more),
     a different propagator is applied than if not.
+
+    Attributes
+    ----------
+    pop_size : int
+        The breeding population size.
+    true_prop : propulate.propagators.Propagator
+        The propagator applied if the current population's size is greater than ``pop_size``.
+    false_prop : propulate.propagators.Propagator
+        The propagator applied if the current population's size equals at least ``pop_size``.
+
+    Notes
+    -----
+    The ``Conditional`` class inherits all methods and attributes from the ``Propagator`` class.
+
+    See Also
+    --------
+    :class:`Propagator` : The parent class.
     """
 
     def __init__(
@@ -136,22 +182,22 @@ def __init__(
         offspring: int = -1,
     ) -> None:
         """
-        Initialize the conditional propagator.
+        Initialize a conditional propagator.
 
         Parameters
         ----------
-        pop_size: int
-            breeding population size
-        true_prop: propulate.propagators.Propagator
-            propagator applied if size of current population >= pop_size.
-        false_prop: propulate.propagators.Propagator
-            propagator applied if size of current population < pop_size.
-        parents: int
-            number of input individuals (-1 for any)
-        offspring: int
-            number of output individuals
+        pop_size : int
+            The breeding population size.
+        true_prop : propulate.propagators.Propagator
+            The propagator applied if the current population's size equals at least ``pop_size``.
+        false_prop : propulate.propagators.Propagator
+            The propagator applied if the current population's size is less than ``pop_size``.
+        parents : int, optional
+            The number of input individuals (-1 for any). Default is -1.
+        offspring : int
+            The number of output individuals to breed. Default is -1.
         """
-        super(Conditional, self).__init__(parents, offspring)
+        super().__init__(parents, offspring)
         self.pop_size = pop_size
         self.true_prop = true_prop
         self.false_prop = false_prop
@@ -162,46 +208,59 @@ def __call__(self, inds: List[Individual]) -> List[Individual]:
 
         Parameters
         ----------
-        inds: list[propulate.individual.Individual]
-            input individuals the propagator is applied to
+        inds : List[propulate.Individual]
+            The input individuals the propagator is applied to.
 
         Returns
         -------
-        list[propulate.individual.Individual]
-            output individuals returned by the conditional propagator
+        List[propulate.Individual]
+            The output individuals returned by the conditional propagator.
         """
         if (
             len(inds) >= self.pop_size
-        ):  # If number of evaluated individuals >= pop_size apply true_prop.
+        ):  # If number of evaluated individuals >= `pop_size`, apply `true_prop`.
             return self.true_prop(inds)
-        else:  # Else apply false_prop.
+        else:  # Else apply `false_prop`.
             return self.false_prop(inds)
 
 
 class Compose(Propagator):
     """
     Stack propagators together sequentially for successive application.
+
+    Attributes
+    ----------
+    propagators : List[propulate.propagators.Propagator]
+        The propagators to be stacked together.
+
+    Notes
+    -----
+    The ``Compose`` class inherits all methods and attributes from the ``Propagator`` class.
+
+    See Also
+    --------
+    :class:`Propagator` : The parent class.
     """
 
     def __init__(self, propagators: List[Propagator]) -> None:
         """
-        Initialize composed propagator.
+        Initialize a composed propagator.
 
         Parameters
         ----------
-        propagators: list[propulate.propagators.Propagator]
-            propagators to be stacked together sequentially
+        propagators : List[propulate.propagators.Propagator]
+            The propagators to be stacked together sequentially.
 
         Raises
         ------
         ValueError
-            If propagators to stack are incompatible in terms of number of input and output individuals.
+            If the propagators to stack are incompatible in terms of number of input and output individuals.
         """
         if len(propagators) < 1:
             raise ValueError(
                 f"Not enough propagators given ({len(propagators)}). At least 1 is required."
             )
-        super(Compose, self).__init__(propagators[0].parents, propagators[-1].offspring)
+        super().__init__(propagators[0].parents, propagators[-1].offspring)
         for i in range(len(propagators) - 1):
             # Check compatibility of consecutive propagators in terms of number of parents + offsprings.
             if not _check_compatible(
@@ -220,17 +279,17 @@ def __init__(self, propagators: List[Propagator]) -> None:
 
     def __call__(self, inds: List[Individual]) -> List[Individual]:
         """
-        Apply composed propagator.
+        Apply the composed propagator.
 
         Parameters
         ----------
-        inds: list[propulate.individual.Individual]
-            input individuals the propagator is applied to
+        inds : List[propulate.Individual]
+            The input individuals the propagator is applied to.
 
         Returns
         -------
-        list[propulate.individual.Individual]
-            output individuals after application of propagator
+        List[propulate.Individual]
+            The output individuals after application of the propagator.
         """
         for p in self.propagators:
             inds = p(inds)
@@ -239,8 +298,15 @@ def __call__(self, inds: List[Individual]) -> List[Individual]:
 
 class SelectMin(Propagator):
     """
-    Select specified number of best performing individuals as evaluated by their losses.
-    i.e., those individuals with minimum losses.
+    Select a specified number of best performing individuals in terms of their losses, i.e., with the smallest losses.
+
+    Notes
+    -----
+    The ``SelectMin`` class inherits all methods and attributes from the ``Propagator`` class.
+
+    See Also
+    --------
+    :class:`Propagator` : The parent class.
     """
 
     def __init__(
@@ -248,28 +314,28 @@ def __init__(
         offspring: int,
     ) -> None:
         """
-        Initialize elitist selection propagator.
+        Initialize an elitist selection propagator.
 
         Parameters
         ----------
-        offspring: int
-            number of offsprings (individuals to be selected)
+        offspring : int
+            The number of offspring (individuals to be selected).
         """
-        super(SelectMin, self).__init__(-1, offspring)
+        super().__init__(-1, offspring)
 
     def __call__(self, inds: List[Individual]) -> List[Individual]:
         """
-        Apply elitist-selection propagator.
+        Apply the elitist-selection propagator.
 
         Parameters
         ----------
-        inds: list[propulate.individual.Individual]
-            input individuals the propagator is applied to
+        inds : List[propulate.Individual]
+            The input individuals the propagator is applied to.
 
         Returns
         -------
-        list[propulate.individual.Individual]
-            selected output individuals after application of the propagator
+        List[propulate.Individual]
+            The selected output individuals after application of the propagator.
 
         Raises
         ------
@@ -288,8 +354,15 @@ def __call__(self, inds: List[Individual]) -> List[Individual]:
 
 class SelectMax(Propagator):
     """
-    Select specified number of worst performing individuals as evaluated by their losses,
-    i.e., those individuals with maximum losses.
+    Select a specified number of worst performing individuals in terms of their losses, i.e., with the greatest losses.
+
+    Notes
+    -----
+    The ``SelectMax`` class inherits all methods and attributes from the ``Propagator`` class.
+
+    See Also
+    --------
+    :class:`Propagator` : The parent class.
     """
 
     def __init__(
@@ -297,28 +370,28 @@ def __init__(
         offspring: int,
     ) -> None:
         """
-        Initialize anti-elitist propagator.
+        Initialize an anti-elitist propagator.
 
         Parameters
         ----------
-        offspring: int
-            number of offspring (individuals to be selected)
+        offspring : int
+            The number of offspring (individuals to be selected).
         """
-        super(SelectMax, self).__init__(-1, offspring)
+        super().__init__(-1, offspring)
 
     def __call__(self, inds: List[Individual]) -> List[Individual]:
         """
-        Apply anti-elitist-selection propagator.
+        Apply the anti-elitist-selection propagator.
 
         Parameters
         ----------
-        inds: list[propulate.individual.Individual]
-            individuals the propagator is applied to
+        inds : List[propulate.Individual]
+            The individuals the propagator is applied to.
 
         Returns
         -------
-        list[propulate.individual.Individual]
-            selected individuals after application of the propagator
+        List[propulate.Individual]
+            The selected individuals after application of the propagator.
 
         Raises
         ------
@@ -337,35 +410,43 @@ def __call__(self, inds: List[Individual]) -> List[Individual]:
 
 class SelectUniform(Propagator):
     """
-    Select specified number of individuals randomly.
+    Select a specified number of individuals randomly.
+
+    Notes
+    -----
+    The ``SelectUniform`` class inherits all methods and attributes from the ``Propagator`` class.
+
+    See Also
+    --------
+    :class:`Propagator` : The parent class.
     """
 
-    def __init__(self, offspring: int, rng: random.Random = None) -> None:
+    def __init__(self, offspring: int, rng: Optional[random.Random] = None) -> None:
         """
-        Initialize random-selection propagator.
+        Initialize a random-selection propagator.
 
         Parameters
         ----------
-        offspring: int
-            number of offspring (individuals to be selected)
-        rng: random.Random
-            random number generator
+        offspring : int
+            The number of offspring (individuals to be selected).
+        rng : random.Random, optional
+            The separate random number generator for the Propulate optimization.
         """
-        super(SelectUniform, self).__init__(-1, offspring, rng)
+        super().__init__(-1, offspring, rng)
 
     def __call__(self, inds: List[Individual]) -> List[Individual]:
         """
-        Apply uniform-selection propagator.
+        Apply the uniform-selection propagator.
 
         Parameters
         ----------
-        inds: list[propulate.individual.Individual]
-            individuals the propagator is applied to
+        inds : List[propulate.Individual]
+            The individuals the propagator is applied to.
 
         Returns
         -------
-        list[propulate.individual.Individual]
-            selected individuals after application of propagator
+        List[propulate.Individual]
+            The selected individuals after application of the propagator.
 
         Raises
         ------
@@ -384,7 +465,20 @@ def __call__(self, inds: List[Individual]) -> List[Individual]:
 # TODO parents should be fixed to one NOTE see utils reason why it is not right now
 class InitUniform(Stochastic):
     """
-    Initialize individual by uniformly sampling specified limits for each trait.
+    Initialize an individual by uniformly sampling the specified limits for each trait.
+
+    Attributes
+    ----------
+    limits : Dict[str, Tuple[float, float]] | Dict[str, Tuple[int, int]] | Dict[str, Tuple[str, ...]]
+        The search space, i.e., the limits of (hyper-)parameters to be optimized.
+
+    Notes
+    -----
+    The ``InitUniform`` class inherits all methods and attributes from the ``Stochastic`` class.
+
+    See Also
+    --------
+    :class:`Stochastic` : The parent class.
     """
 
     def __init__(
@@ -396,38 +490,38 @@ def __init__(
         ],
         parents: int = 0,
         probability: float = 1.0,
-        rng: random.Random = None,
+        rng: Optional[random.Random] = None,
     ) -> None:
         """
-        Initialize random-initialization propagator.
+        Initialize a random-initialization propagator.
 
         Parameters
         ----------
-        limits: dict[str, tuple[float, float]] | dict[str, tuple[int, int]] | dict[str, tuple[str, ...]]
-            search space, i.e., limits of (hyper-)parameters to be optimized
-        parents: int
-            number of parents
-        probability: float
-            probability of application
-        rng: random.Random
-            random number generator
+        limits : Dict[str, Tuple[float, float]] | Dict[str, Tuple[int, int]] | Dict[str, Tuple[str, ...]]
+            The search space, i.e., the limits of (hyper-)parameters to be optimized.
+        parents : int, optional
+            The number of parents. Default is 0.
+        probability : float, optional
+            The probability of application. Default is 1.0.
+        rng : random.Random, optional
+            The separate random number generator for the Propulate optimization.
         """
-        super(InitUniform, self).__init__(parents, 1, probability, rng)
+        super().__init__(parents, 1, probability, rng)
         self.limits = limits
 
     def __call__(self, *inds: Individual) -> Individual:
         """
-        Apply uniform-initialization propagator.
+        Apply the uniform-initialization propagator.
 
         Parameters
         ----------
-        inds: propulate.individual.Individual
-            individuals the propagator is applied to
+        inds : propulate.Individual
+            The individuals the propagator is applied to.
 
         Returns
         -------
-        propulate.individual.Individual
-            output individual after application of propagator
+        propulate.Individual
+            The output individual after application of the propagator.
 
         Raises
         ------
diff --git a/propulate/propagators/cmaes.py b/propulate/propagators/cmaes.py
index efdb4ed2..c3107202 100644
--- a/propulate/propagators/cmaes.py
+++ b/propulate/propagators/cmaes.py
@@ -1,5 +1,5 @@
 import random
-from typing import List, Dict, Tuple
+from typing import List, Dict, Optional, Tuple
 
 import numpy as np
 
@@ -10,11 +10,77 @@
 class CMAParameter:
     """
     Handle and store all basic/active CMA-related constants/variables and strategy parameters.
+
+    Attributes
+    ----------
+    b_matrix : numpy.ndarray
+        The B matrix in the covariance matrix decomposition.
+    c_1 : float
+        The learning rate for the rank-one update of the covariance matrix update.
+    c_c : float
+        The decay rate for evolution path for the rank-one update of the covariance matrix.
+    c_mu : float
+        The learning rate for the rank-mu update of the covariance matrix update.
+    c_sigma : float
+        A step-size control parameter.
+    chi_n : float
+        The expectation value of ||N(0,I)||.
+    covariance_inv_sqrt : numpy.ndarray
+        Square root of the inverse of the covariance matrix: C^-1/2 = B*D^(-1)*B^T
+    covariance_matrix : numpy.ndarray
+        The covariance matrix.
+    condition_limit : float
+        The maximum allowed condition of the covariance matrix to ensure numerical stability.
+    constant_trace : bool
+        Whether to keep the trace (sum of diagonal elements) of ``self.covariance_matrix`` constant.
+    count_eval : int
+        The number of individuals evaluated.
+    d_matrix : numpy.ndarray
+        The D matrix in the covariance matrix decomposition.
+    d_sigma : float
+        A step-size control parameter.
+    eigen_eval : int
+        The Number of individuals evaluated when the covariance matrix was last decomposed into B and D.
+    exploration : bool
+        If True decompose covariance matrix for each generation (worse runtime, less exploitation, more
+        ``decompose_in_each_generation``); else decompose covariance matrix only after a certain number of
+        individuals evaluated (better runtime, more exploitation, less ``decompose_in_each_generation``).
+    lambd : int
+        The number of individuals considered for each generation.
+    limits : Dict[str, float]
+        The limits of the search space.
+    mean : numpy.ndarray
+        The distribution's mean.
+    mu : int
+        The number of positive recombination weights.
+    mu_eff : float
+        The variance effective selection mass.
+    old_mean : numpy.ndarray
+        The mean of the last generation.
+    p_c : numpy.ndarray
+        A dynamic strategy variable.
+    p_sigma : numpy.ndarray
+        A dynamic strategy variable.
+    problem_dimension : int
+        The number of dimensions in the search space.
+    sigma : float
+        The distribution's standard deviation.
+    weights : numpy.ndarray
+        The recombination weights.
+
+    Methods
+    -------
+    update_mean()
+        Update mean and old mean.
+    update_covariance_matrix()
+        Update the covariance matrix.
+    mahalanobis_norm()
+        Compute the Mahalanobis distance.
     """
 
     def __init__(
         self,
-        lamb: int,
+        lambd: int,
         mu: int,
         problem_dimension: int,
         weights: np.ndarray,
@@ -23,6 +89,7 @@ def __init__(
         c_1: float,
         c_mu: float,
         limits: Dict,
+        initial_mean: np.ndarray,
         exploration: bool,
     ) -> None:
         """
@@ -30,24 +97,26 @@ def __init__(
 
         Parameters
         ----------
-        lamb : int
-            number of individuals considered for each generation
+        lambd : int
+            The number of individuals considered for each generation.
         mu : int
-            number of positive recombination weights
+            The number of positive recombination weights.
         problem_dimension : int
-            number of dimensions in the search space
+            The number of dimensions in the search space.
         weights : numpy.ndarray
-            recombination weights
+            The recombination weights.
         mu_eff : float
-            variance effective selection mass
+            The variance effective selection mass.
         c_c : float
-            decay rate for evolution path for the rank-one update of the covariance matrix
+            The decay rate for the evolution path for the rank-one update of the covariance matrix.
         c_1 : float
-            learning rate for the rank-one update of the covariance matrix update
+            The learning rate for the rank-one update of the covariance matrix update.
         c_mu : float
-            learning rate for the rank-mu update of the covariance matrix update
+            The learning rate for the rank-mu update of the covariance matrix update.
         limits : dict
-            limits of search space
+            The limits of the search space.
+        initial_mean : np.ndarray
+            The initial mean of the distribution.
         exploration : bool
             If True decompose covariance matrix for each generation (worse runtime, less exploitation, more
             ``decompose_in_each_generation``); else decompose covariance matrix only after a certain number of
@@ -55,7 +124,7 @@ def __init__(
         """
         self.problem_dimension = problem_dimension
         self.limits = limits
-        self.lamb = lamb
+        self.lambd = lambd
         self.mu = mu
         self.weights = weights
         self.mu_eff = mu_eff
@@ -76,32 +145,30 @@ def __init__(
         self.p_c = np.zeros((problem_dimension, 1))
 
         # Prevent equal eigenvalues, hack from https://github.com/CMA-ES/pycma/blob/development/cma/sampler.py
-        self.co_matrix = np.diag(
+        self.covariance_matrix = np.diag(
             np.ones(problem_dimension)
             * np.exp(
                 (1e-4 / self.problem_dimension) * np.arange(self.problem_dimension)
             )
         )
         self.b_matrix = np.eye(self.problem_dimension)
-        # Assume ``self.co_matrix`` to be initialized as a diagonal matrix.
-        self.d_matrix = np.diag(self.co_matrix) ** 0.5
+        # Assume ``self.covariance_matrix`` to be initialized as a diagonal matrix.
+        self.d_matrix = np.diag(self.covariance_matrix) ** 0.5
         # Sort eigenvalues in ascending order
         indices_eig = self.d_matrix.argsort()
         self.d_matrix = self.d_matrix[indices_eig]
         self.b_matrix = self.b_matrix[:, indices_eig]
         # Square root of the inverse of the covariance matrix: C^-1/2 = B*D^(-1)*B^T
-        self.co_inv_sqrt = (
+        self.covariance_inv_sqrt = (
             self.b_matrix @ np.diag(self.d_matrix ** (-1)) @ self.b_matrix.T
         )
         # Maximum allowed condition of the covariance matrix to ensure numerical stability
         self.condition_limit = 1e5 - 1
-        # Whether to keep the trace (sum of diagonal elements) of ``self.co_matrix`` constant.
+        # Whether to keep the trace (sum of diagonal elements) of ``self.covariance_matrix`` constant.
         self.constant_trace = False
 
         # Use this initial mean when using multiple islands.
-        self.mean = np.array(
-            [[np.random.uniform(*limits[limit]) for limit in limits]]
-        ).reshape((problem_dimension, 1))
+        self.mean = initial_mean
         # 0.3 instead of 0.2 is also often used for greater initial step size
         self.sigma = 0.2 * (
             (max(max(limits[i]) for i in limits)) - min(min(limits[i]) for i in limits)
@@ -111,110 +178,81 @@ def __init__(
         self.old_mean = None
         self.exploration = exploration
 
-        # Number of individuals evaluated when the covariance matrix was last decomposed into B and D
+        # Number of individuals evaluated when the covariance matrix was last decomposed into B and D.
         self.eigen_eval = 0
         # Number of individuals evaluated
         self.count_eval = 0
 
         # Expectation value of ||N(0,I)||
-        self.chiN = problem_dimension**0.5 * (
+        self.chi_n = problem_dimension**0.5 * (
             1 - 1.0 / (4 * problem_dimension) + 1.0 / (21 * problem_dimension**2)
         )
 
-    def set_mean(self, new_mean: np.ndarray) -> None:
+    def update_mean(self, new_mean: np.ndarray) -> None:
         """
-        Setter for mean property. Updates the old mean as well.
+        Update mean and old mean property.
 
         Parameters
         ----------
         new_mean : numpy.ndarray
-            new mean
+            The new mean.
         """
         self.old_mean = self.mean
         self.mean = new_mean
 
-    def set_p_sigma(self, new_p_sigma: np.ndarray) -> None:
-        """
-        Setter for evolution path of step-size adaptation.
-
-        Parameters
-        ----------
-        new_p_sigma : numpy.ndarray
-            new evolution path
+    def update_covariance_matrix(self, new_covariance_matrix: np.ndarray) -> None:
         """
-        self.p_sigma = new_p_sigma
+        Update the covariance matrix.
 
-    def set_p_c(self, new_p_c: np.ndarray) -> None:
-        """
-        Setter for evolution path of covariance matrix adaptation.
+        Computes new values for ``b_matrix``, ``d_matrix``, and ``covariance_inv_sqrt``. Decomposition of
+        ``covariance_matrix`` is O(n^3), hence the possibility of lazy updating ``b_matrix`` and ``d_matrix``.
 
         Parameters
         ----------
-        new_p_c : numpy.ndarray
-            evolution path
-        """
-        self.p_c = new_p_c
-
-    def set_sigma(self, new_sigma: float) -> None:
+        new_covariance_matrix : numpy.ndarray
+            The new covariance matrix.
         """
-        Setter for step-size.
-
-        Parameters
-        ----------
-        new_sigma : float
-            step-size
-        """
-        self.sigma = new_sigma
-
-    def set_co_matrix(self, new_co_matrix: np.ndarray) -> None:
-        """
-        Setter for the covariance matrix.
-
-        Computes new values for ``b_matrix``, ``d_matrix``, and ``co_inv_sqrt``. Decomposition of ``co_matrix`` is
-        O(n^3), hence the possibility of lazy updating ``b_matrix`` and ``d_matrix``.
-
-        Parameters
-        ----------
-        new_co_matrix : numpy.ndarray
-            new covariance matrix
-        """
-        # Update b and d matrix and co_inv_sqrt only after certain number of evaluations to ensure 0(n^2).
-        # Also trade-off decompose_in_each_generation or not.
+        # Update b and d matrix and covariance_inv_sqrt only after certain number of evaluations to ensure 0(n^2).
+        # Also, trade-off decompose_in_each_generation or not.
         if self.exploration or (
             self.count_eval - self.eigen_eval
-            > self.lamb / (self.c_1 + self.c_mu) / self.problem_dimension / 10
+            > self.lambd / (self.c_1 + self.c_mu) / self.problem_dimension / 10
         ):
             self.eigen_eval = self.count_eval
-            self._decompose_co_matrix(new_co_matrix)
-            self.co_inv_sqrt = (
+            self._decompose_co_matrix(new_covariance_matrix)
+            self.covariance_inv_sqrt = (
                 self.b_matrix @ np.diag(self.d_matrix ** (-1)) @ self.b_matrix.T
             )
             # Ensure symmetry.
-            self.co_inv_sqrt = (self.co_inv_sqrt + self.co_inv_sqrt.T) / 2
+            self.covariance_inv_sqrt = (
+                self.covariance_inv_sqrt + self.covariance_inv_sqrt.T
+            ) / 2
 
     def _decompose_co_matrix(self, new_co_matrix: np.ndarray) -> None:
         """
-        Eigendecomposition of the covariance matrix into eigenvalues (d_matrix) and eigenvectors (columns of b_matrix)
+        Eigen-decomposition of the covariance matrix into eigenvalues (d_matrix) and eigenvectors (columns of b_matrix).
+
         Parameters
         ----------
-        new_co_matrix: the new covariance matrix that should be decomposed
+        new_co_matrix : numpy.ndarray
+            The new covariance matrix that should be decomposed.
         """
         # Enforce symmetry.
-        self.co_matrix = np.triu(new_co_matrix) + np.triu(new_co_matrix, 1).T
+        self.covariance_matrix = np.triu(new_co_matrix) + np.triu(new_co_matrix, 1).T
         d_matrix_old = self.d_matrix
         try:
-            self.d_matrix, self.b_matrix = np.linalg.eigh(self.co_matrix)
-            if any(self.d_matrix <= 0):
-                # Covariance matrix eigen decomposition failed, consider reformulating objective function.
+            self.d_matrix, self.b_matrix = np.linalg.eigh(self.covariance_matrix)
+            if np.any(self.d_matrix <= 0):
+                # Covariance matrix eigen-decomposition failed, consider reformulating objective function.
                 raise ValueError("Covariance matrix not positive definite.")
         except Exception as _:
-            # Add min(eigenvalues(self.co_matrix_old)) to diag(self.co_matrix) and try again
+            # Add min(eigenvalues(self.co_matrix_old)) to diag(self.covariance_matrix) and try again
             min_eig_old = min(d_matrix_old) ** 2
             for i in range(self.problem_dimension):
-                self.co_matrix[i, i] += min_eig_old
+                self.covariance_matrix[i, i] += min_eig_old
             # Replace eigenvalues with standard deviations
             self.d_matrix = (d_matrix_old**2 + min_eig_old) ** 0.5
-            self._decompose_co_matrix(self.co_matrix)
+            self._decompose_co_matrix(self.covariance_matrix)
         else:
             assert all(np.isfinite(self.d_matrix))
             self._sort_b_d_matrix()
@@ -223,8 +261,8 @@ def _decompose_co_matrix(self, new_co_matrix: np.ndarray) -> None:
             if self.constant_trace:
                 s = 1 / np.mean(
                     self.d_matrix
-                )  # normalize co_matrix to control overall magnitude
-                self.co_matrix *= s
+                )  # normalize covariance_matrix to control overall magnitude
+                self.covariance_matrix *= s
                 self.d_matrix *= s
             self.d_matrix **= 0.5
 
@@ -238,23 +276,21 @@ def _limit_condition(self, limit: float) -> None:
         Parameters
         ----------
         limit : float
-            threshold for the condition of the matrix
+            The threshold for the condition of the matrix.
         """
         # Check if condition number of matrix is too big.
         if (self.d_matrix[-1] / self.d_matrix[0]) ** 2 > limit:
             eps = (self.d_matrix[-1] ** 2 - limit * self.d_matrix[0] ** 2) / (limit - 1)
             for i in range(self.problem_dimension):
                 # Decrease ratio of largest to smallest eigenvalue, absolute difference remains.
-                self.co_matrix[i, i] += eps
+                self.covariance_matrix[i, i] += eps
             # Eigenvalues are definitely positive now.
             self.d_matrix **= 2
             self.d_matrix += eps
             self.d_matrix **= 0.5
 
     def _sort_b_d_matrix(self) -> None:
-        """
-        Sort columns of ``b_matrix`` and ``d_matrix`` according to the eigenvalues in ``d_matrix``.
-        """
+        """Sort columns of ``b_matrix`` and ``d_matrix`` according to the eigenvalues in ``d_matrix``."""
         indices_eig = np.argsort(self.d_matrix)
         self.d_matrix = self.d_matrix[indices_eig]
         self.b_matrix = self.b_matrix[:, indices_eig]
@@ -265,37 +301,52 @@ def _sort_b_d_matrix(self) -> None:
 
     def mahalanobis_norm(self, dx: np.ndarray) -> np.ndarray:
         """
-        Compute the Mahalanobis distance by using C^(-1/2) and the difference vector of a point to the mean of a
-        distribution.
+        Compute the Mahalanobis distance using C^(-1/2) and the difference vector of a point to the distribution's mean.
 
         Parameters
         ----------
         dx : numpy.ndarray
-            difference vector
+            The difference vector.
 
         Returns
         -------
         numpy.ndarray
-            resulting Mahalanobis distance
+            The resulting Mahalanobis distance.
         """
-        return np.linalg.norm(np.dot(self.co_inv_sqrt, dx))
+        return np.linalg.norm(np.dot(self.covariance_inv_sqrt, dx))
 
 
 class CMAAdapter:
     """
     Abstract base class for the adaption of strategy parameters of CMA-ES.
+
     Strategy class from the viewpoint of the strategy design pattern.
+
+
+    Methods
+    -------
+    update_mean()
+        Abstract method for updating of mean in CMA-ES variants.
+    update_step_size()
+        Update step-size in CMA-ES variants.
+    update_covariance_matrix()
+        Abstract method for the adaptation of the covariance matrix of CMA-ES variants.
+    compute_weights()
+        Abstract method for computing the recombination weights of a CMA-ES variant.
+    compute_learning_rates()
+        Compute the learning rates for the CMA-variants.
     """
 
     def update_mean(self, par: CMAParameter, arx: np.ndarray) -> None:
         """
         Abstract method for updating of mean in CMA-ES variants.
+
         Parameters
         ----------
         par : CMAParameter
-            parameter object of the CMA-ES propagation
+            The parameter object of the CMA-ES propagation.
         arx : numpy.ndarray
-            individuals of the distribution
+            The individuals of the distribution.
 
         Raises
         ------
@@ -312,33 +363,27 @@ def update_step_size(par: CMAParameter) -> None:
         Parameters
         ----------
         par : CMAParameter
-            parameter object of the CMA-ES propagation
-        """
-        par.set_p_sigma(
-            (1 - par.c_sigma) * par.p_sigma
-            + np.sqrt(par.c_sigma * (2 - par.c_sigma) * par.mu_eff)
-            * par.co_inv_sqrt
-            @ (par.mean - par.old_mean)
-            / par.sigma
-        )
-        par.set_sigma(
-            par.sigma
-            * np.exp(
-                (par.c_sigma / par.d_sigma)
-                * (np.linalg.norm(par.p_sigma, ord=2) / par.chiN - 1)
-            )
+            The parameter object of the CMA-ES propagation.
+        """
+        par.p_sigma = (1 - par.c_sigma) * par.p_sigma + np.sqrt(
+            par.c_sigma * (2 - par.c_sigma) * par.mu_eff
+        ) * par.covariance_inv_sqrt @ (par.mean - par.old_mean) / par.sigma
+
+        par.sigma = par.sigma * np.exp(
+            (par.c_sigma / par.d_sigma)
+            * (np.linalg.norm(par.p_sigma, ord=2) / par.chi_n - 1)
         )
 
     def update_covariance_matrix(self, par: CMAParameter, arx: np.ndarray) -> None:
         """
-        Abstract method for the adaption of the covariance matrix of CMA-ES variants.
+        Abstract method for the adaptation of the covariance matrix of CMA-ES variants.
 
         Parameters
         ----------
         par : CMAParameter
-            parameter object of the CMA-ES propagation
+            The parameter object of the CMA-ES propagation.
         arx : numpy.ndarray
-            individuals of the distribution
+            The individuals of the distribution.
 
         Raises
         ------
@@ -356,11 +401,11 @@ def compute_weights(
         Parameters
         ----------
         mu : int
-            number of positive recombination weights
+            The number of positive recombination weights.
         lamb : int
-            number of individuals considered for each generation
+            The number of individuals considered for each generation.
         problem_dimension : int
-            number of dimensions in the search space
+            The number of dimensions in the search space.
 
         Returns
         -------
@@ -384,14 +429,18 @@ def compute_learning_rates(
         Parameters
         ----------
         mu_eff : float
-            variance effective selection mass
+            The variance effective selection mass.
         problem_dimension : int
-            number of dimensions in the search space
+            The number of dimensions in the search space.
 
         Returns
         -------
-        tuple[float, float, float]
-            tuple of c_c, c_1, c_mu
+        float
+            The decay rate for evolution path for the rank-one update of the covariance matrix, ``c_c``.
+        float
+            The learning rate for the rank-one update of the covariance matrix update, ``c_1``.
+        float
+            The learning rate for the rank-mu update of the covariance matrix update, ``c_mu``.
         """
         c_c = (4 + mu_eff / problem_dimension) / (
             problem_dimension + 4 + 2 * mu_eff / problem_dimension
@@ -406,8 +455,17 @@ def compute_learning_rates(
 
 class BasicCMA(CMAAdapter):
     """
-    Adaption of strategy parameters of CMA-ES according to the original CMA-ES algorithm. Concrete strategy class from
-    the viewpoint of the strategy design pattern.
+    Adaption of strategy parameters of CMA-ES according to the original CMA-ES algorithm.
+
+    Concrete strategy class from the viewpoint of the strategy design pattern.
+
+    Notes
+    -----
+    The ``BasicCMA`` class inherits all methods and attributes from the ``CMAAdapter`` class.
+
+    See Also
+    --------
+    :class:`CMAAdapter` : The parent class.
     """
 
     def compute_weights(
@@ -419,16 +477,24 @@ def compute_weights(
         Parameters
         ----------
         mu : int
-            number of positive recombination weights
+            The number of positive recombination weights
         lamb : int
-            number of individuals considered for each generation
+            The number of individuals considered for each generation
         problem_dimension : int
-            number of dimensions in the search space
+            The number of dimensions in the search space
 
         Returns
         -------
-        tuple[np.ndarray, float, float, float, float]
-            tuple of the weights, mu_eff, c_1, c_c and c_mu.
+        numpy.ndarray
+            The weights.
+        float
+            The variance effective selection mass, ``mu_eff``.
+        float
+            The learning rate for the rank-one update of the covariance matrix update, ``c_1``.
+        float
+            The decay rate for evolution path for the rank-one update of the covariance matrix, ``c_c``.
+        float
+            The learning rate for the rank-mu update of the covariance matrix update, ``c_mu``.
         """
         weights = np.log(mu + 0.5) - np.log(np.arange(1, mu + 1))
         weights /= np.sum(weights)
@@ -443,12 +509,12 @@ def update_mean(self, par: CMAParameter, arx: np.ndarray) -> None:
         Parameters
         ----------
         par : CMAParameter
-            parameter object of the CMA-ES propagation
+            The parameter object of the CMA-ES propagation.
         arx : numpy.ndarray
-            individuals of the distribution
+            The individuals of the distribution.
         """
         # Matrix vector multiplication (reshape weights to column vector)
-        par.set_mean(arx @ par.weights.reshape(-1, 1))
+        par.update_mean(arx @ par.weights.reshape(-1, 1))
 
     def update_covariance_matrix(self, par: CMAParameter, arx: np.ndarray) -> None:
         """
@@ -457,36 +523,33 @@ def update_covariance_matrix(self, par: CMAParameter, arx: np.ndarray) -> None:
         Parameters
         ----------
         par : CMAParameter
-            parameter object of the CMA-ES propagation
+            The parameter object of the CMA-ES propagation.
         arx : numpy.ndarray
-            individuals of the distribution
+            The individuals of the distribution.
         """
         # Turn off rank-one accumulation when sigma increases quickly.
         h_sig = np.sum(par.p_sigma**2) / (
-            1 - (1 - par.c_sigma) ** (2 * (par.count_eval / par.lamb))
+            1 - (1 - par.c_sigma) ** (2 * (par.count_eval / par.lambd))
         ) / par.problem_dimension < 2 + 4.0 / (par.problem_dimension + 1)
         # Update evolution path.
-        par.set_p_c(
-            (1 - par.c_c) * par.p_c
-            + h_sig
-            * np.sqrt(par.c_c * (2 - par.c_c) * par.mu_eff)
-            * (par.mean - par.old_mean)
-            / par.sigma
-        )
+        par.p_c = (1 - par.c_c) * par.p_c + h_sig * np.sqrt(
+            par.c_c * (2 - par.c_c) * par.mu_eff
+        ) * (par.mean - par.old_mean) / par.sigma
+
         # Use ``h_sig`` to the power of two (unlike in paper) for the variance loss from ``h_sig``.
         ar_tmp = (1 / par.sigma) * (
             arx[:, : par.mu] - np.tile(par.old_mean, (1, par.mu))
         )
         new_co_matrix = (
-            (1 - par.c_1 - par.c_mu) * par.co_matrix
+            (1 - par.c_1 - par.c_mu) * par.covariance_matrix
             + par.c_1
             * (
                 par.p_c @ par.p_c.T
-                + (1 - h_sig) * par.c_c * (2 - par.c_c) * par.co_matrix
+                + (1 - h_sig) * par.c_c * (2 - par.c_c) * par.covariance_matrix
             )
             + par.c_mu * ar_tmp @ (par.weights * ar_tmp).T
         )
-        par.set_co_matrix(new_co_matrix)
+        par.update_covariance_matrix(new_co_matrix)
 
 
 class ActiveCMA(CMAAdapter):
@@ -496,6 +559,14 @@ class ActiveCMA(CMAAdapter):
     Differently from the original CMA-ES algorithm, active CMA-ES uses negative recombination weights (only for the
     covariance matrix adaptation) for individuals with relatively low fitness.
     Concrete strategy class from the viewpoint of the strategy design pattern.
+
+    Notes
+    -----
+    The ``ActiveCMA`` class inherits all methods and attributes from the ``CMAAdapter`` class.
+
+    See Also
+    --------
+    :class:`CMAAdapter` : The parent class.
     """
 
     def compute_weights(
@@ -507,23 +578,31 @@ def compute_weights(
         Parameters
         ----------
         mu : int
-            number of positive recombination weights
+            The number of positive recombination weights.
         lamb : int
-            number of individuals considered for each generation
+            The number of individuals considered for each generation.
         problem_dimension : int
-            number of dimensions in the search space
+            The number of dimensions in the search space.
 
         Returns
         -------
-        tuple[np.ndarray, float, float, float, float]
-            tuple of the weights, mu_eff, c_1, c_c and c_mu
+        numpy.ndarray
+            The weights.
+        float
+            The variance effective selection mass, ``mu_eff``.
+        float
+            The learning rate for the rank-one update of the covariance matrix update, ``c_1``.
+        float
+            The decay rate for evolution path for the rank-one update of the covariance matrix, ``c_c``.
+        float
+            The learning rate for the rank-mu update of the covariance matrix update, ``c_mu``.
         """
         weights_preliminary = np.log(lamb / 2 + 0.5) - np.log(np.arange(1, lamb + 1))
         mu_eff = np.sum(weights_preliminary[:mu]) ** 2 / np.sum(
             weights_preliminary[:mu] ** 2
         )
         c_c, c_1, c_mu = ActiveCMA.compute_learning_rates(mu_eff, problem_dimension)
-        # Now compute final weights.
+        # Compute final weights.
         mu_eff_minus = np.sum(weights_preliminary[mu:]) ** 2 / np.sum(
             weights_preliminary[mu:] ** 2
         )
@@ -546,13 +625,13 @@ def update_mean(self, par: CMAParameter, arx: np.ndarray) -> None:
         Parameters
         ----------
         par : CMAParameter
-            parameter object of the CMA-ES propagation
+            The parameter object of the CMA-ES propagation.
         arx : numpy.ndarray
-            individuals of the distribution
+            The individuals of the distribution.
         """
         # Matrix vector multiplication (reshape weights to column vector)
         # Only consider positive weights.
-        par.set_mean(arx @ par.weights[: par.mu].reshape(-1, 1))
+        par.update_mean(arx @ par.weights[: par.mu].reshape(-1, 1))
 
     def update_covariance_matrix(self, par: CMAParameter, arx: np.ndarray) -> None:
         """
@@ -561,23 +640,20 @@ def update_covariance_matrix(self, par: CMAParameter, arx: np.ndarray) -> None:
         Parameters
         ----------
         par : CMAParameter
-            parameter object of the CMA-ES propagation
+            The parameter object of the CMA-ES propagation.
         arx : numpy.ndarray
-            individuals of the distribution
+            The individuals of the distribution.
         """
         # Turn off rank-one accumulation when sigma increases quickly.
         h_sig = np.sum(par.p_sigma**2) / (
-            1 - (1 - par.c_sigma) ** (2 * (par.count_eval / par.lamb))
+            1 - (1 - par.c_sigma) ** (2 * (par.count_eval / par.lambd))
         ) / par.problem_dimension < 2 + 4.0 / (par.problem_dimension + 1)
         # Update evolution path.
-        par.set_p_c(
-            (1 - par.c_c) * par.p_c
-            + h_sig
-            * np.sqrt(par.c_c * (2 - par.c_c) * par.mu_eff)
-            * (par.mean - par.old_mean)
-            / par.sigma
-        )
-        weights_circle = np.zeros((par.lamb,))
+        par.p_c = (1 - par.c_c) * par.p_c + h_sig * np.sqrt(
+            par.c_c * (2 - par.c_c) * par.mu_eff
+        ) * (par.mean - par.old_mean) / par.sigma
+
+        weights_circle = np.zeros((par.lambd,))
         for i, w_i in enumerate(par.weights):
             # Guarantee positive definiteness.
             weights_circle[i] = w_i
@@ -591,17 +667,17 @@ def update_covariance_matrix(self, par: CMAParameter, arx: np.ndarray) -> None:
                     ** 2
                 )
         # Use ``h_sig`` to the power of two (unlike in paper) for the variance loss from ``h_sig``.
-        ar_tmp = (1 / par.sigma) * (arx - np.tile(par.old_mean, (1, par.lamb)))
+        ar_tmp = (1 / par.sigma) * (arx - np.tile(par.old_mean, (1, par.lambd)))
         new_co_matrix = (
-            (1 - par.c_1 - par.c_mu) * par.co_matrix
+            (1 - par.c_1 - par.c_mu) * par.covariance_matrix
             + par.c_1
             * (
                 par.p_c @ par.p_c.T
-                + (1 - h_sig) * par.c_c * (2 - par.c_c) * par.co_matrix
+                + (1 - h_sig) * par.c_c * (2 - par.c_c) * par.covariance_matrix
             )
             + par.c_mu * ar_tmp @ (weights_circle * ar_tmp).T
         )
-        par.set_co_matrix(new_co_matrix)
+        par.update_covariance_matrix(new_co_matrix)
 
 
 class CMAPropagator(Propagator):
@@ -610,17 +686,45 @@ class CMAPropagator(Propagator):
 
     Uses ``CMAAdapter`` to adapt strategy parameters like mean, step-size, and covariance matrix and stores them in a
     ``CMAParameter`` object. The context class from the viewpoint of the strategy design pattern.
+
+    Attributes
+    ----------
+    adapter : CMAAdapter
+        The adaptation strategy of CMA-ES.    par
+    pool_size : int
+        The size of the pool of individuals pre-selected before selecting the best from this pool.
+    select_single_best : SelectMin
+        Selection operator to select the best individual.
+    select_from_pool : SelectUniform
+        Select randomly from breeding pool.
+    select_pool : SelectMin
+        Select the best individuals from available ones as breeding pool.
+    select_worst : SelectMax
+        Select a specified number of the worst individuals.
+    select_worst_all_time : bool
+        If True, use the worst individuals for negative recombination weights in active CMA-ES, else use the worst
+        (lambda - mu) individuals of the best lambda individuals.
+    rng : random.Random
+        The separate random number generator for the Propulate optimization.
+
+    Notes
+    -----
+    The ``CMAPropagator`` class inherits all methods and attributes from the ``Propagator`` class.
+
+    See Also
+    --------
+    :class:`Propagator` : The parent class.
     """
 
     def __init__(
         self,
         adapter: CMAAdapter,
         limits: Dict,
-        rng: random.Random,
         decompose_in_each_generation: bool = False,
         select_worst_all_time: bool = False,
-        pop_size: int = None,
+        pop_size: Optional[int] = None,
         pool_size: int = 3,
+        rng: Optional[random.Random] = None,
     ) -> None:
         """
         Instantiate a CMA-ES propagator.
@@ -628,44 +732,51 @@ def __init__(
         Parameters
         ----------
         adapter : CMAAdapter
-            adaption strategy of CMA-ES
-        limits : dict
-            limits of the search space
-        rng: random.Random
-            random number generator
-        decompose_in_each_generation : bool
+            The adaptation strategy of CMA-ES.
+        limits : Dict[str, float]
+            The limits of the search space.
+        decompose_in_each_generation : bool, optional
             If True, decompose covariance matrix for each generation (worse runtime, less exploitation, more
             exploration); else decompose covariance matrix only after a certain number of individuals evaluated
-            (better runtime, more exploitation, less exploration)
-        select_worst_all_time : bool
+            (better runtime, more exploitation, less exploration). Default is False.
+        select_worst_all_time : bool, optional
             If True, use the worst individuals for negative recombination weights in active CMA-ES, else use the worst
-            (lambda - mu) individuals of the best lambda individuals. If BasicCMA is used, the given value is irrelevant
-            regarding functionality.
-        pop_size : int
-            number of individuals to be considered in each generation
-        pool_size : iny
-            size of the pool of individuals preselected before selecting the best of this pool
+            (lambda - mu) individuals of the best lambda individuals. If ``BasicCMA`` is used, the given value is
+            irrelevant regarding functionality. Default is False.
+        pop_size : int, optional
+            The number of individuals to be considered in each generation.
+        pool_size : int, optional
+            The size of the pool of individuals pre-selected before selecting the best from this pool. Default is 3.
+        rng: random.Random, optional
+            The separate random number generator for the Propulate optimization.
         """
         self.adapter = adapter
         problem_dimension = len(limits)
         # Number of individuals considered for each generation
-        lamb = (
+        lambd = (
             pop_size if pop_size else 4 + int(np.floor(3 * np.log(problem_dimension)))
         )
-        super(CMAPropagator, self).__init__(lamb, 1)
-
+        super().__init__(lambd, 1, rng=rng)
+        self.numpy_rng = np.random.default_rng(
+            seed=self.rng.randint(a=0, b=np.iinfo(np.int32).max)
+        )
         # Number of positive recombination weights
-        mu = lamb // 2
-        self.select_worst = SelectMax(lamb - mu)
+        mu = lambd // 2
+        self.select_worst = SelectMax(lambd - mu)
         self.select_worst_all_time = select_worst_all_time
 
         # CMA-ES variant specific weights and learning rates
         weights, mu_eff, c_c, c_1, c_mu = adapter.compute_weights(
-            mu, lamb, problem_dimension
+            mu, lambd, problem_dimension
         )
 
+        initial_mean = np.array(
+            [[self.numpy_rng.uniform(*limits[limit]) for limit in limits]]
+        ).reshape((problem_dimension, 1))
+        # 0.3 instead of 0.2 is also often used for greater initial step size
+
         self.par = CMAParameter(
-            lamb,
+            lambd,
             mu,
             problem_dimension,
             weights,
@@ -674,12 +785,13 @@ def __init__(
             c_1,
             c_mu,
             limits,
+            initial_mean,
             decompose_in_each_generation,
         )
         self.pool_size = int(pool_size) if int(pool_size) >= 1 else 3
-        self.select_pool = SelectMin(self.pool_size * lamb)
+        self.select_pool = SelectMin(self.pool_size * lambd)
         self.select_from_pool = SelectUniform(mu - 1, rng=rng)
-        self.select_best_1 = SelectMin(1)
+        self.select_single_best = SelectMin(1)
 
     def __call__(self, inds: List[Individual]) -> Individual:
         """
@@ -690,13 +802,13 @@ def __call__(self, inds: List[Individual]) -> Individual:
 
         Parameters
         ----------
-        inds: list[Individual]
-            individuals available
+        inds: List[propulate.Individual]
+            Available individuals.
 
         Returns
         -------
-        new_ind : Individual
-            newly sampled individual
+        new_ind : propulate.Individual
+            The newly sampled individual.
         """
         num_inds = len(inds)
         # Add individuals from different workers to ``eval_count``.
@@ -704,9 +816,9 @@ def __call__(self, inds: List[Individual]) -> Individual:
         # Sample new individual.
         new_ind = self._sample_cma()
         # Check if ``len(inds)`` >= or < ``pool_size * lambda`` and make sample or sample + update.
-        if num_inds >= self.pool_size * self.par.lamb:
+        if num_inds >= self.pool_size * self.par.lambd:
             inds_pooled = self.select_pool(inds)
-            best = self.select_best_1(inds_pooled)
+            best = self.select_single_best(inds_pooled)
             if not self.select_worst_all_time:
                 worst = self.select_worst(inds_pooled)
             else:
@@ -729,17 +841,17 @@ def __call__(self, inds: List[Individual]) -> Individual:
 
     def _transform_individuals_to_matrix(self, inds: List[Individual]) -> np.ndarray:
         """
-        Take a list of individuals and transform it to numpy matrix for easier subsequent computation.
+        Take a list of individuals and transform it to a numpy array for easier subsequent computation.
 
         Parameters
         ----------
-        inds : list[Individual]
-            list of individuals
+        inds : list[propulate.Individual]
+            The list of individuals.
 
         Returns
         -------
         arx : numpy.ndarray
-            Array of shape (problem_dimension, len(inds))
+            Array of shape [problem_dimension, len(inds)].
         """
         arx = np.zeros((self.par.problem_dimension, len(inds)))
         for k, ind in enumerate(inds):
@@ -753,12 +865,11 @@ def _sample_cma(self) -> Individual:
 
         Returns
         -------
-        new_ind : Individual
-            the newly sampled individual
+        new_ind : propulate.Individual
+            The newly sampled individual.
         """
-        new_x = None
         # Generate new offspring
-        random_vector = np.random.randn(self.par.problem_dimension, 1)
+        random_vector = self.numpy_rng.standard_normal((self.par.problem_dimension, 1))
         try:
             new_x = self.par.mean + self.par.sigma * self.par.b_matrix @ (
                 self.par.d_matrix * random_vector
@@ -774,58 +885,3 @@ def _sample_cma(self) -> Individual:
         for i, (dim, _) in enumerate(self.par.limits.items()):
             new_ind[dim] = new_x[i, 0]
         return new_ind
-
-    def get_mean(self) -> np.ndarray:
-        """
-        Getter for mean attribute.
-
-        Returns
-        -------
-        mean : numpy.ndarray
-            current CMA-ES mean of the best mu individuals
-        """
-        return self.par.mean
-
-    def get_sigma(self) -> float:
-        """
-        Getter for step size.
-
-        Returns
-        -------
-        sigma : float
-            current step-size
-        """
-        return self.par.sigma
-
-    def get_co_matrix(self) -> np.ndarray:
-        """
-        Getter for covariance matrix.
-
-        Returns
-        -------
-        co_matrix : numpy.ndarray
-            current covariance matrix
-        """
-        return self.par.co_matrix
-
-    def get_evolution_path_sigma(self) -> np.ndarray:
-        """
-        Getter for evolution path of step-size adaption.
-
-        Returns
-        -------
-        p_sigma : numpy.ndarray
-            evolution path for step-size adaption
-        """
-        return self.par.p_sigma
-
-    def get_evolution_path_co_matrix(self) -> np.ndarray:
-        """
-        Getter for evolution path of covariance matrix adaption.
-
-        Returns
-        -------
-        p_c : numpy.ndarray
-            evolution path for covariance matrix adaption
-        """
-        return self.par.p_c
diff --git a/propulate/propagators/ga.py b/propulate/propagators/ga.py
index 3d7bbc70..4f136142 100644
--- a/propulate/propagators/ga.py
+++ b/propulate/propagators/ga.py
@@ -1,6 +1,6 @@
 import copy
 import random
-from typing import List, Dict, Union, Tuple
+from typing import Dict, List, Optional, Tuple, Union
 
 import numpy as np
 
@@ -11,6 +11,21 @@
 class PointMutation(Stochastic):
     """
     Point-mutate given number of traits with given probability.
+
+    Attributes
+    ----------
+    points : int
+        The number of points to mutate.
+    limits : Dict[str, Tuple[float, float]] | Dict[str, Tuple[int, int]] | Dict[str, Tuple[str, ...]]
+        The search space, i.e., the limits of (hyper-)parameters to be optimized.
+
+    Notes
+    -----
+    The ``PointMutation`` class inherits all methods and attributes from the ``Stochastic`` class.
+
+    See Also
+    --------
+    :class:`Stochastic` : The parent class.
     """
 
     def __init__(
@@ -22,28 +37,28 @@ def __init__(
         ],
         points: int = 1,
         probability: float = 1.0,
-        rng: random.Random = None,
+        rng: Optional[random.Random] = None,
     ) -> None:
         """
         Initialize point-mutation propagator.
 
         Parameters
         ----------
-        limits: dict[str, tuple[float, float]] | dict[str, tuple[int, int]] | dict[str, tuple[str, ...]]
-            limits of (hyper-)parameters to be optimized
-        points: int
-            number of points to mutate
-        probability: float
-            probability of application
-        rng: random.Random
-            random number generator
+        limits : Dict[str, Tuple[float, float]] | Dict[str, Tuple[int, int]] | Dict[str, Tuple[str, ...]]
+            The search space, i.e., the limits of the (hyper-)parameters to be optimized.
+        points : int, optional
+            The number of points to mutate. Default is 1.
+        probability : float, optional
+            The probability of application. Default is 1.0.
+        rng : random.Random, optional
+            The separate random number generator of the Propulate optimization.
 
         Raises
         ------
         ValueError
             If the requested number of points to mutate is greater than the number of traits.
         """
-        super(PointMutation, self).__init__(1, 1, probability, rng)
+        super().__init__(1, 1, probability, rng)
         self.points = points
         self.limits = limits
         if len(limits) < points:
@@ -53,17 +68,17 @@ def __init__(
 
     def __call__(self, ind: Individual) -> Individual:
         """
-        Apply point-mutation propagator.
+        Apply the point-mutation propagator.
 
         Parameters
         ----------
-        ind: propulate.individual.Individual
-            individual the propagator is applied to
+        ind : propulate.Individual
+            The individual the propagator is applied to.
 
         Returns
         -------
-        propulate.individual.Individual
-            possibly point-mutated individual after application of propagator
+        propulate.Individual
+            The possibly point-mutated individual after application of the propagator.
         """
         if (
             self.rng.random() < self.probability
@@ -85,13 +100,29 @@ def __call__(self, ind: Individual) -> Individual:
                 elif isinstance(ind[i], str):
                     # Return random element from non-empty sequence.
                     ind[i] = self.rng.choice(self.limits[i])
-
         return ind  # Return point-mutated individual.
 
 
 class RandomPointMutation(Stochastic):
     """
     Point-mutate random number of traits with given probability.
+
+    Attributes
+    ----------
+    min_points : int
+        The minimum number of points to mutate.
+    max_points : int
+        The maximum number of points to mutate.
+    limits : Dict[str, Tuple[float, float]] | Dict[str, Tuple[int, int]] | Dict[str, Tuple[str, ...]]
+        The search space, i.e., the limits of (hyper-)parameters to be optimized.
+
+    Notes
+    -----
+    The ``RandomPointMutation`` class inherits all methods and attributes from the ``Stochastic`` class.
+
+    See Also
+    --------
+    :class:`Stochastic` : The parent class.
     """
 
     def __init__(
@@ -104,34 +135,32 @@ def __init__(
         min_points: int = 1,
         max_points: int = 1,
         probability: float = 1.0,
-        rng: random.Random = None,
+        rng: Optional[random.Random] = None,
     ) -> None:
         """
         Initialize random point-mutation propagator.
 
         Parameters
         ----------
-        limits: dict[str, tuple[float, float]] | dict[str, tuple[int, int]] | dict[str, tuple[str, ...]]
-            limits of parameters to optimize, i.e., search space
-        min_points: int
-            minimum number of points to mutate
-        max_points: int
-            maximum number of points to mutate
-        probability: float
-            probability of application
-        rng: random.Random
+        limits : Dict[str, Tuple[float, float]] | Dict[str, Tuple[int, int]] | Dict[str, Tuple[str, ...]]
+            The limits of the parameters to optimize, i.e., the search space.
+        min_points : int, optional
+            minimum number of points to mutate. Default is 1.
+        max_points : int, optional
+            maximum number of points to mutate. Default is 1.
+        probability : float, optional
+            probability of application. Default is 1.0.
+        rng : random.Random, optional
             random number generator
 
         Raises
         ------
         ValueError
             If no or a negative number of points shall be mutated.
-        ValueError
             If there are fewer traits than requested number of points to mutate.
-        ValueError
             If the requested minimum number of points to mutate is greater than the requested maximum number.
         """
-        super(RandomPointMutation, self).__init__(1, 1, probability, rng)
+        super().__init__(1, 1, probability, rng)
         if min_points <= 0:
             raise ValueError(
                 f"Minimum number of points to mutate must be > 0 but was {min_points}."
@@ -151,17 +180,17 @@ def __init__(
 
     def __call__(self, ind: Individual) -> Individual:
         """
-        Apply random-point-mutation propagator.
+        Apply the random-point-mutation propagator.
 
         Parameters
         ----------
-        ind: propulate.population.Individual
-            individual the propagator is applied to
+        ind : propulate.Individual
+            The individual the propagator is applied to.
 
         Returns
         -------
-        propulate.population.Individual
-            possibly point-mutated individual after application of propagator
+        propulate.Individual
+            The possibly point-mutated individual after application of the propagator.
         """
         if (
             self.rng.random() < self.probability
@@ -190,7 +219,24 @@ def __call__(self, ind: Individual) -> Individual:
 
 class IntervalMutationNormal(Stochastic):
     """
-    Mutate given number of traits according to Gaussian distribution around current value with given probability.
+    Mutate a number of traits according to a Gaussian distribution around the current value with the given probability.
+
+    Attributes
+    ----------
+    sigma_factor : float
+        The scaling factor for the interval width to obtain standard deviation.
+    points : int
+        The number of points to mutate.
+    limits : Dict[str, Tuple[float, float]] | Dict[str, Tuple[int, int]] | Dict[str, Tuple[str, ...]]
+        The search space, i.e., the limits of (hyper-)parameters to be optimized.
+
+    Notes
+    -----
+    The ``IntervalMutationNormal`` class inherits all methods and attributes from the ``Stochastic`` class.
+
+    See Also
+    --------
+    :class:`Stochastic` : The parent class.
     """
 
     def __init__(
@@ -203,31 +249,31 @@ def __init__(
         sigma_factor: float = 0.1,
         points: int = 1,
         probability: float = 1.0,
-        rng: random.Random = None,
+        rng: Optional[random.Random] = None,
     ) -> None:
         """
         Initialize interval-mutation propagator.
 
         Parameters
         ----------
-        limits: dict[str, tuple[float, float]] | dict[str, tuple[int, int]] | dict[str, tuple[str, ...]]
-            limits of (hyper-)parameters to be optimized, i.e., search space
-        sigma_factor: float
-            scaling factor for interval width to obtain standard deviation
-        points: int
-            number of points to mutate
-        probability: float
-            probability of application
-        rng: random.Random
-            random number generator
+        limits : Dict[str, Tuple[float, float]] | Dict[str, Tuple[int, int]] | Dict[str, Tuple[str, ...]]
+            The limits of the (hyper-)parameters to be optimized, i.e., the search space.
+        sigma_factor : float, optional
+            The scaling factor for the interval width to obtain the standard deviation. Default is 0.1.
+        points : int, optional
+            The number of points to mutate. Default is 1.
+        probability : float, optional
+            The probability of application, Default is 1.0
+        rng : random.Random, optional
+            The separate random number generator for the Propulate optimization.
 
         Raises
         ------
         ValueError
             If the individuals has fewer continuous traits than the requested number of points to mutate.
         """
-        super(IntervalMutationNormal, self).__init__(1, 1, probability, rng)
-        self.points = points  # number of traits to point-mutate
+        super().__init__(1, 1, probability, rng)
+        self.points = points  # Number of traits to point-mutate
         self.limits = limits
         self.sigma_factor = sigma_factor
         n_interval_traits = len([x for x in limits if isinstance(limits[x][0], float)])
@@ -238,17 +284,17 @@ def __init__(
 
     def __call__(self, ind: Individual) -> Individual:
         """
-        Apply interval-mutation propagator.
+        Apply the interval-mutation propagator.
 
         Parameters
         ----------
-        ind: propulate.individual.Individual
-            input individual the propagator is applied to
+        ind : propulate.Individual
+            The input individual the propagator is applied to.
 
         Returns
         -------
-        propulate.individual.Individual
-            possibly interval-mutated output individual after application of propagator
+        propulate.Individual
+            The possibly interval-mutated output individual after application of the propagator.
         """
         if (
             self.rng.random() < self.probability
@@ -277,56 +323,67 @@ def __call__(self, ind: Individual) -> Individual:
         return ind  # Return point-mutated individual.
 
 
-class MateUniform(Stochastic):  # uniform crossover
+class CrossoverUniform(Stochastic):  # uniform crossover
     """
     Generate new individual by uniform crossover of two parents with specified relative parent contribution.
+
+    Attributes
+    ----------
+    rel_parent_contrib : float
+        The relative parent contribution with respect to the first parent.
+
+    Notes
+    -----
+    The ``CrossoverUniform`` class inherits all methods and attributes from the ``Stochastic`` class.
+
+    See Also
+    --------
+    :class:`Stochastic` : The parent class.
     """
 
     def __init__(
         self,
-        rel_parent_contrib: float = 0.5,
+        relative_parent_contribution: float = 0.5,
         probability: float = 1.0,
-        rng: random.Random = None,
+        rng: Optional[random.Random] = None,
     ) -> None:
         """
         Initialize uniform crossover propagator.
 
         Parameters
         ----------
-        rel_parent_contrib: float
-            relative parent contribution w.r.t. first parent
-        probability: float
-            probability of application
-        rng: random.Random
-            random number generator
+        relative_parent_contribution : float, optional
+            The relative parent contribution with respect to the first parent. Default is 0.5.
+        probability: float, optional
+            The probability of application. Default is 1.0.
+        rng: random.Random, optional
+            The separate random number generator for the Propulate optimization.
 
         Raises
         ------
         ValueError
             If the relative parent contribution is not within [0, 1].
         """
-        super(MateUniform, self).__init__(
-            2, 1, probability, rng
-        )  # Breed 1 offspring from 2 parents.
-        if rel_parent_contrib <= 0 or rel_parent_contrib >= 1:
+        super().__init__(2, 1, probability, rng)  # Breed 1 offspring from 2 parents.
+        if relative_parent_contribution <= 0 or relative_parent_contribution >= 1:
             raise ValueError(
-                f"Relative parent contribution must be within (0, 1) but was {rel_parent_contrib}."
+                f"Relative parent contribution must be within (0, 1) but was {relative_parent_contribution}."
             )
-        self.rel_parent_contrib = rel_parent_contrib
+        self.rel_parent_contrib = relative_parent_contribution
 
     def __call__(self, inds: List[Individual]) -> Individual:
         """
-        Apply uniform-crossover propagator.
+        Apply the uniform-crossover propagator.
 
         Parameters
         ----------
-        inds: List[propulate.individual.Individual]
-            individuals the propagator is applied to
+        inds : List[propulate.Individual]
+            The individuals the propagator is applied to.
 
         Returns
         -------
-        propulate.individual.Individual
-            possibly cross-bred individual after application of propagator
+        propulate.Individual
+            The possibly cross-bred individual after application of the propagator.
         """
         ind = copy.deepcopy(inds[0])  # Consider 1st parent.
         ind.loss = None  # Initialize individual's loss attribute.
@@ -340,43 +397,54 @@ def __call__(self, inds: List[Individual]) -> Individual:
         return ind  # Return offspring.
 
 
-class MateMultiple(Stochastic):  # uniform crossover
+class CrossoverMultiple(Stochastic):  # uniform crossover
     """
     Breed new individual by uniform crossover of multiple parents.
+
+    Notes
+    -----
+    The ``CrossoverMultiple`` class inherits all methods and attributes from the ``Stochastic`` class.
+
+    See Also
+    --------
+    :class:`Stochastic` : The parent class.
     """
 
     def __init__(
-        self, parents: int = -1, probability: float = 1.0, rng: random.Random = None
+        self,
+        parents: int = -1,
+        probability: float = 1.0,
+        rng: Optional[random.Random] = None,
     ) -> None:
         """
-        Initialize multiple-crossover propagator.
+        Initialize a multiple-crossover propagator.
 
         Parameters
         ----------
-        probability: float
-            probability of application
-        parents: int
-            number of parents
-        rng: random.Random
-            random number generator
+        probability : float, optional
+            The probability of application. Default is 1.0.
+        parents: int, optional
+            The number of parents (not used) here. Default is -1.
+        rng : random.Random, optional
+            The separate random number generator for the Propulate optimization.
         """
-        super(MateMultiple, self).__init__(
+        super().__init__(
             parents, 1, probability, rng
         )  # Breed 1 offspring from specified number of parents.
 
     def __call__(self, inds: List[Individual]) -> Individual:
         """
-        Apply multiple-crossover propagator.
+        Apply the multi-crossover propagator.
 
         Parameters
         ----------
-        inds: list[propulate.individual.Individual]
-            individuals the propagator is applied to
+        inds : list[propulate.Individual]
+            The individuals the propagator is applied to.
 
         Returns
         -------
-        propulate.individual.Individual
-            possibly cross-bred individual after application of propagator
+        propulate.Individual
+            The possibly cross-bred individual after application of propagator
         """
         ind = copy.deepcopy(inds[0])  # Consider 1st parent.
         ind.loss = None  # Initialize individual's loss attribute.
@@ -389,50 +457,63 @@ def __call__(self, inds: List[Individual]) -> Individual:
         return ind  # Return offspring.
 
 
-class MateSigmoid(Stochastic):
+class CrossoverSigmoid(Stochastic):
     """
     Generate new individual by crossover of two parents according to Boltzmann sigmoid probability.
 
     Consider two parent individuals with fitness values f1 and f2. Let f1 <= f2. For each trait,
     the better parent's value is accepted with the probability sigmoid(- (f1-f2) / temperature).
+
+    Attributes
+    ----------
+    temperature : float
+        The temperature hyperparameter of the Boltzmann distribution.
+
+    Notes
+    -----
+    The ``CrossoverSigmoid`` class inherits all methods and attributes from the ``Stochastic`` class.
+
+    See Also
+    --------
+    :class:`Stochastic` : The parent class.
     """
 
     def __init__(
         self,
         temperature: float = 1.0,
         probability: float = 1.0,
-        rng: random.Random = None,
+        rng: Optional[random.Random] = None,
     ) -> None:
         """
-        Initialize sigmoid-crossover propagator.
+        Initialize a sigmoid-crossover propagator.
 
         Parameters
         ----------
-        temperature: float
-            temperature for Boltzmann factor in sigmoid probability
-        probability: float
-            probability of application
-        rng: random.Random
-            random number generator
+        temperature : float, optional
+            The temperature in the Boltzmann factor of the sigmoid probability. Default is 1.0.
+        probability : float, optional
+            The probability of application. Default is 1.0.
+        rng : random.Random, optional
+            The separate random number generator for the Propulate optimization.
         """
-        super(MateSigmoid, self).__init__(
+        super(CrossoverSigmoid, self).__init__(
             2, 1, probability, rng
         )  # Breed 1 offspring from 2 parents.
         self.temperature = temperature
 
     def __call__(self, inds: List[Individual]) -> Individual:
         """
-        Apply sigmoid-crossover propagator.
+        Apply the sigmoid-crossover propagator.
 
         Parameters
         ----------
-        inds: list[propulate.individual.Individual]
-            individuals the propagator is applied to
+        inds : List[propulate.Individual]
+            The individuals the propagator is applied to.
 
         Returns
         -------
-        propulate.individual.Individual
-            possibly cross-bred individual after application of propagator
+        propulate.Individual
+            The possibly cross-bred individual after application of the propagator.
         """
         ind = copy.deepcopy(inds[0])  # Consider 1st parent.
         ind.loss = None  # Initialize individual's loss attribute.
diff --git a/propulate/propagators/pso.py b/propulate/propagators/pso.py
index 6dce4b86..8fe1999f 100644
--- a/propulate/propagators/pso.py
+++ b/propulate/propagators/pso.py
@@ -24,6 +24,31 @@ class BasicPSO(Propagator):
 
     This variant was first proposed in 1998 by Y. Shi and R. Eberhart, "A modified particle swarm optimizer"
     https://doi.org/10.1109/ICEC.1998.699146.
+
+    Attributes
+    ----------
+    c_cognitive : float
+        The constant cognitive factor for scaling the distance to the particle's personal best value.
+    c_social : float
+        The constant social factor for scaling the distance to the swarm's global best value.
+    inertia : float
+        The inertia weight.
+    limits : Dict[str, Tuple[float, float]]
+        The borders of the continuous search domain.
+    limits_as_array : numpy.ndarray
+        The limits converted to a numpy array.
+    rank : int
+        The global rank of the worker the propagator is living on.
+    rng : random.Random
+        The separate random number generator for introducing non-linearity.
+
+    Notes
+    -----
+    The ``BasicPSO`` class inherits all methods and attributes from the ``Propagator`` class.
+
+    See Also
+    --------
+    :class:`Propagator` : The parent class.
     """
 
     def __init__(
@@ -45,17 +70,17 @@ def __init__(
         Parameters
         ----------
         inertia : float
-            inertia weight.
+            The inertia weight.
         c_cognitive : float
-            constant cognitive factor for scaling the distance to the particle's personal best value
+            The constant cognitive factor for scaling the distance to the particle's personal best value.
         c_social : float
-            constant social factor for scaling the distance to the swarm's global best value
+            The constant social factor for scaling the distance to the swarm's global best value.
         rank : int
-            global rank of the worker the propagator is living on
+            The global rank of the worker the propagator is living on.
         limits : Dict[str, Tuple[float, float]]
-            borders of the continuous search domain
+            The borders of the continuous search domain.
         rng : random.Random
-            random number generator for introducing non-linearity
+            The separate random number generator for introducing non-linearity.
         """
         super().__init__(parents=-1, offspring=1)
         self.c_social = c_social
@@ -75,16 +100,16 @@ def __call__(self, individuals: List[Individual]) -> Particle:
 
         Parameters
         ----------
-        individuals: List[Individual]
-            list of individuals that must at least contain one individual that belongs to the propagator
+        individuals : List[propulate.Individual]
+            A list of individuals that must at least contain one individual that belongs to the propagator.
             This list is used to calculate personal and global best of the particle and the swarm,
             respectively, and then to update the particle based on the retrieved results. Individuals that
             cannot be used as ``Particle`` objects are converted to particles first.
 
         Returns
         -------
-        propulate.population.Particle
-            updated particle
+        propulate.Particle
+            The updated particle.
         """
         old_p, p_best, g_best = self._prepare_data(individuals)
 
@@ -106,12 +131,12 @@ def _prepare_data(
 
         Parameters
         ----------
-        individuals : List[Individual]
-            ``Individual`` objects that shall be used as data basis for a PSO update step
+        individuals : List[propulate.Individual]
+            ``Individual`` objects that shall be used as data basis for a PSO update step.
 
         Returns
         -------
-        Tuple[propulate.population.Particle, propulate.population.Particle, propulate.population.Particle]
+        Tuple[propulate.Particle, propulate.Particle, propulate.Particle]
             The following particles in this very order:
             1.  old_p: the current particle to be updated now
             2.  p_best: the personal best value of this particle
@@ -162,16 +187,16 @@ def _make_new_particle(
         Parameters
         ----------
         position : np.ndarray
-            position of the particle to be created
+            The position of the particle to be created.
         velocity : np.ndarray
-            velocity of the particle to be created
+            The velocity of the particle to be created.
         generation : int
-            generation of the new particle
+            The generation of the new particle.
 
         Returns
         -------
-        propulate.population.Particle
-            new ``Particle`` object resulting from the PSO update step
+        propulate.Particle
+            The new ``Particle`` object resulting from the PSO update step.
         """
         new_p = Particle(position, velocity, generation, self.rank)
         for i, k in enumerate(self.limits):
@@ -187,6 +212,19 @@ class VelocityClampingPSO(BasicPSO):
     to determine each parameter's maximum velocity value relative to its search-space limits.
 
     Based on these values, the velocities of the particles are cut down to a reasonable value.
+
+    Attributes
+    ----------
+    v_cap : numpy.ndarray
+        The clamping velocities for each dimension of the search space.
+
+    Notes
+    -----
+    The ``VelocityClampingPSO`` class inherits all methods and attributes from the ``BasicPSO`` class.
+
+    See Also
+    --------
+    :class:`BasicPSO` : The parent class.
     """
 
     def __init__(
@@ -205,19 +243,19 @@ def __init__(
         Parameters
         ----------
         inertia : float
-            inertia factor
+            The inertia factor.
         c_cognitive : float
-            constant cognitive factor for scaling the distance to the particle's personal best value
+            The constant cognitive factor for scaling the distance to the particle's personal best value.
         c_social : float
-            constant social factor for scaling the distance to the swarm's global best value
+            The constant social factor for scaling the distance to the swarm's global best value.
         rank : int
-            global rank of the worker the propagator is living on
+            The global rank of the worker the propagator is living on.
         limits : Dict[str, Tuple[float, float]]
-            borders of the continuous search domain
+            The borders of the continuous search domain.
         rng : random.Random
-            random number generator for introducing non-linearity
+            The separate random number generator for introducing non-linearity.
         v_limits : Union[float, np.ndarray]
-            clamping factor to be multiplied with the clamping limit in order to reduce it further
+            The clamping factor to be multiplied with the clamping limit in order to reduce it further.
             Should be in (0, 1). If this parameter has float type, it is applied to all dimensions of the search
             domain; else, each of its elements is applied to the corresponding dimension of the search domain.
         """
@@ -236,16 +274,16 @@ def __call__(self, individuals: List[Individual]) -> Particle:
 
         Parameters
         ----------
-        individuals: List[Individual]
-            list of individuals that must at least contain one individual that belongs to the propagator
+        individuals : List[propulate.Individual]
+            The list of individuals that must at least contain one individual that belongs to the propagator.
             This list is used to calculate personal and global best of the particle and the swarm,
             respectively, and then to update the particle based on the retrieved results. Individuals that
             cannot be used as ``Particle`` objects are converted to particles first.
 
         Returns
         -------
-        propulate.population.Particle
-            updated particle
+        propulate.Particle
+            The updated particle.
         """
         old_p, p_best, g_best = self._prepare_data(individuals)
 
@@ -272,6 +310,14 @@ class ConstrictionPSO(BasicPSO):
     The constriction factor is calculated from cognitive and social factors and thus no separate hyperparameter.
 
     This propagator runs on ``Particle`` objects.
+
+    Notes
+    -----
+    The ``ConstrictionPSO`` class inherits all methods and attributes from the ``BasicPSO`` class.
+
+    See Also
+    --------
+    :class:`BasicPSO` : The parent class.
     """
 
     def __init__(
@@ -290,17 +336,17 @@ def __init__(
         Parameters
         ----------
         c_cognitive : float
-            constant cognitive factor for scaling the distance to the particle's personal best value
+            The constant cognitive factor for scaling the distance to the particle's personal best value.
             *Has to sum up with ``c_social`` to a number greater than 4!*
         c_social : float
-            constant social factor for scaling the distance to the swarm's global best value
+            The constant social factor for scaling the distance to the swarm's global best value.
             *Has to sum up with ``c_cognitive`` to a number greater than 4!*
         rank : int
-            global rank of the worker the propagator is living on
+            The global rank of the worker the propagator is living on.
         limits : Dict[str, Tuple[float, float]]
-            borders of the continuous search domain
+            The borders of the continuous search domain.
         rng : random.Random
-            random number generator for introducing non-linearity
+            The random number generator for introducing non-linearity.
 
         Raises
         ------
@@ -324,16 +370,16 @@ def __call__(self, individuals: List[Individual]) -> Particle:
 
         Parameters
         ----------
-        individuals: List[Individual]
-            list of individuals that must at least contain one individual that belongs to the propagator
+        individuals: List[propulate.Individual]
+            A list of individuals that must at least contain one individual that belongs to the propagator.
             This list is used to calculate personal and global best of the particle and the swarm,
             respectively, and then to update the particle based on the retrieved results. Individuals that
             cannot be used as ``Particle`` objects are converted to particles first.
 
         Returns
         -------
-        propulate.population.Particle
-            updated particle
+        propulate.Particle
+            The updated particle.
         """
         old_p, p_best, g_best = self._prepare_data(individuals)
 
@@ -363,9 +409,13 @@ class CanonicalPSO(ConstrictionPSO):
     Proceedings of the 2000 Congress on Evolutionary Computation. CEC00 (Cat. No.00TH8512), La Jolla, CA, USA, 2000,
     pp. 84-88 vol.1, https://10.1109/CEC.2000.870279.
 
+    Notes
+    -----
+    The ``CanonicalPSO`` class inherits all methods and attributes from the ``ConstrictionPSO`` class.
+
     See Also
     --------
-    :class: ``ConstrictionPSO`` : Parent class
+    :class:`ConstrictionPSO` : The parent class.
     """
 
     def __init__(
@@ -386,15 +436,15 @@ def __init__(
         Parameters
         ----------
         c_cognitive : float
-            constant cognitive factor for scaling the distance to the particle's personal best value
+            The constant cognitive factor for scaling the distance to the particle's personal best value.
         c_social : float
-            constant social factor to scaling the distance to the swarm's global best value
+            The constant social factor to scaling the distance to the swarm's global best value.
         rank : int
-            global rank of the worker the propagator is living on
+            The global rank of the worker the propagator is living on.
         limits : Dict[str, Tuple[float, float]]
-            borders of the continuous search domain
+            The borders of the continuous search domain.
         rng : random.Random
-            random number generator for introducing non-linearity
+            The random number generator for introducing non-linearity.
         """
         super().__init__(c_cognitive, c_social, rank, limits, rng)
         x_min, x_max = self.limits_as_array
@@ -410,16 +460,16 @@ def __call__(self, individuals: List[Individual]) -> Particle:
 
         Parameters
         ----------
-        individuals: List[Individual]
-            list of individuals that must at least contain one individual that belongs to the propagator
+        individuals : List[propulate.Individual]
+            The list of individuals that must at least contain one individual that belongs to the propagator.
             This list is used to calculate personal and global best of the particle and the swarm,
             respectively, and then to update the particle based on the retrieved results. Individuals that
             cannot be used as ``Particle`` objects are converted to particles first.
 
         Returns
         -------
-        propulate.population.Particle
-            update particle
+        propulate.Particle
+            The updated particle.
         """
         # Abuse Constriction's update rule, so I don't have to rewrite it.
         victim = super().__call__(individuals)
@@ -435,6 +485,25 @@ def __call__(self, individuals: List[Individual]) -> Particle:
 class InitUniformPSO(Stochastic):
     """
     Initialize ``Particle`` by uniformly sampling specified limits for each trait.
+
+    Attributes
+    ----------
+    limits : Dict[str, Tuple[float, float]]
+        The borders of the continuous search domain.
+    limits_as_array : numpy.ndarray
+        The limits converted to a numpy array.
+    v_limits : float | np.ndarray
+        The multiplicative constant to reduce initial random velocity values.
+    rank : int
+        The global rank of the worker the propagator is living on.
+
+    Notes
+    -----
+    The ``InitUniformPSO`` class inherits all methods and attributes from the ``Stochastic`` class.
+
+    See Also
+    --------
+    :class:`Stochastic` : The parent class.
     """
 
     def __init__(
@@ -453,51 +522,54 @@ def __init__(
 
         Parameters
         ----------
-        limits : dict[str, tuple[float, float]]
-            limits of the search space, i.e., limits of (hyper-)parameters to be optimized
+        limits : Dict[str, Tuple[float, float]]
+            The limits of the search space, i.e., the limits of (hyper-)parameters to be optimized.
         rank : int
-            rank of the worker in MPI.COMM_WORLD
-        parents : int
-            number of input individuals (-1 for any)
-        probability : float
-            probability of creating a completely new individual
-        rng : random.Random
-            random number generator
-        v_init_limit : float | np.ndarray
-            multiplicative constant to reduce initial random velocity values
+            The rank of the worker in the Propulate communicator
+        parents : int, optional
+            The number of input individuals (-1 for any). Default is 0.
+        probability : float, optional
+            The probability of creating a completely new individual. Default is 1.0.
+        rng : random.Random, optional
+            The separate random number generator for the Propulate optimization.
+        v_init_limit : float | np.ndarray, optional
+            The multiplicative constant to reduce initial random velocity values. Default is 0.1.
         """
         super().__init__(parents, 1, probability, rng)
         self.limits = limits
-        self.laa = np.array(list(limits.values())).T
+        self.limits_as_array = np.array(list(limits.values())).T
         if isinstance(v_init_limit, np.ndarray):
-            assert v_init_limit.shape[-1] == self.laa.shape[-1]
+            assert v_init_limit.shape[-1] == self.limits_as_array.shape[-1]
         self.v_limits = v_init_limit
         self.rank = rank
 
     def __call__(self, individuals: List[Individual]) -> Particle:
         """
-        Apply uniform-initialization propagator.
+        Apply the uniform-initialization propagator.
 
         Parameters
         ----------
-        individuals : list[Individual]
-            individuals the propagator is applied to
+        individuals : List[propulate.Individual]
+            The individuals the propagator is applied to.
 
         Returns
         -------
-        propulate.population.Particle
-            one particle object
+        propulate.Particle
+            A single particle object.
         """
         if (
             len(individuals) == 0 or self.rng.random() < self.probability
         ):  # Apply only with specified `probability`.
             position = np.array(
-                [self.rng.uniform(*self.laa[..., i]) for i in range(self.laa.shape[-1])]
+                [
+                    self.rng.uniform(*self.limits_as_array[..., i])
+                    for i in range(self.limits_as_array.shape[-1])
+                ]
             )
             velocity = np.array(
                 [
-                    self.rng.uniform(*(self.v_limits * self.laa)[..., i])
-                    for i in range(self.laa.shape[-1])
+                    self.rng.uniform(*(self.v_limits * self.limits_as_array)[..., i])
+                    for i in range(self.limits_as_array.shape[-1])
                 ]
             )
 
@@ -532,6 +604,27 @@ class StatelessPSO(Propagator):
     As this propagator works without velocity, there is also no inertia weight used.
 
     It uses only classes provided by vanilla Propulate.
+
+    Attributes
+    ----------
+    c_cognitive : float
+        The constant cognitive factor for scaling the distance to the particle's personal best value.
+    c_social : float
+        The constant social factor for scaling the distance to the swarm's global best value.
+    limits : Dict[str, Tuple[float, float]]
+        The borders of the continuous search domain.
+    rank : int
+        The global rank of the worker the propagator is living on.
+    rng : random.Random
+        The separate random number generator for introducing non-linearity.
+
+    Notes
+    -----
+    The ``StatelessPSO`` class inherits all methods and attributes from the ``Propagator`` class.
+
+    See Also
+    --------
+    :class:`Propagator` : The parent class.
     """
 
     def __init__(
@@ -548,15 +641,15 @@ def __init__(
         Parameters
         ----------
         c_cognitive : float
-            constant cognitive factor for scaling individual's personal best value
+            The constant cognitive factor for scaling the individual's personal best value.
         c_social : float
-            constant social factor for scaling swarm's global best value
+            The constant social factor for scaling the swarm's global best value.
         rank : int
-            global rank of the worker the propagator is living on
-        limits : dict[str, tuple[float, float]
-            borders of the continuous search domain
+            The global rank of the worker the propagator is living on.
+        limits : Dict[str, Tuple[float, float]
+            The borders of the continuous search domain.
         rng : random.Random
-            random number generator required for non-linearity of update
+            The random number generator required for non-linearity of the update.
         """
         super().__init__(parents=-1, offspring=1)
         self.c_social = c_social
@@ -567,17 +660,17 @@ def __init__(
 
     def __call__(self, individuals: List[Individual]) -> Individual:
         """
-        Apply standard PSO update without inertia and old velocity.
+        Apply the standard PSO update without inertia and old velocity.
 
         Parameters
         ----------
-        individuals : list[Individual]
-            individuals are used as data basis for the PSO update
+        individuals : List[propulate.Individual]
+            The individuals used as data basis for the PSO update.
 
         Returns
         -------
-        propulate.population.Individual
-            updated individual
+        propulate.Individual
+            The updated individual.
 
         Raises
         ------
diff --git a/propulate/propulator.py b/propulate/propulator.py
index 63f0714c..4cfad97f 100644
--- a/propulate/propulator.py
+++ b/propulate/propulator.py
@@ -6,7 +6,7 @@
 import time
 from operator import attrgetter
 from pathlib import Path
-from typing import Callable, Union, List, Tuple, Type
+from typing import Callable, List, Optional, Tuple, Type, Union
 
 import deepdiff
 import numpy as np
@@ -26,80 +26,136 @@ class Propulator:
     This class provides Propulate's basic asynchronous population-based optimization routine (without an island model).
     At the same time, it serves as a base class of ``Migrator`` and ``Pollinator``, which implement an asynchronous
     island model on top of the asynchronous base optimizer with real migration and pollination, respectively.
+
+    Attributes
+    ----------
+    checkpoint_path : str | pathlib.Path
+        The path where checkpoints are loaded from and stored to.
+    emigration_propagator : Type[propulate.Propagator]
+        The emigration propagator, i.e., how to choose individuals for emigration that are sent to the destination
+        island. Should be some kind of selection operator.
+    generation : int
+        The current generation.
+    generations : int
+        The overall number of generations.
+    island_comm : MPI.Comm
+        The intra-island communicator for communication within the island.
+    island_counts : numpy.ndarray
+        The numbers of workers for each island.
+    island_displs : numpy.ndarray
+        The island displacements.
+    island_idx : int
+        The island's index.
+    loss_fn : Callable
+        The loss function to be minimized.
+    migration_prob : float
+        The migration probability.
+    migration_topology : np.ndarray
+        The migration topology.
+    population : List[propulate.Individual]
+        The population list of individuals on that island.
+    propagator : propulate.Propagator
+        The evolutionary operator.
+    propulate_comm : MPI.Comm
+        The Propulate world communicator, consisting of rank 0 of each worker's sub communicator.
+    rng : random.Random
+        The separate random number generator for the Propulate optimization.
+    worker_sub_comm : MPI.Comm
+        The worker's internal communicator for parallelized evaluation of single individuals.
+
+    Methods
+    -------
+    propulate()
+        Run asynchronous evolutionary optimization routine without migration or pollination.
+    summarize()
+        Get top-n results from Propulate optimization.
     """
 
     def __init__(
         self,
         loss_fn: Callable,
         propagator: Propagator,
+        rng: random.Random,
         island_idx: int = 0,
-        comm: MPI.Comm = MPI.COMM_WORLD,
+        island_comm: MPI.Comm = MPI.COMM_WORLD,
+        propulate_comm: MPI.Comm = MPI.COMM_WORLD,
+        worker_sub_comm: MPI.Comm = MPI.COMM_SELF,
         generations: int = -1,
         checkpoint_path: Union[str, Path] = Path("./"),
-        migration_topology: np.ndarray = None,
+        migration_topology: Optional[np.ndarray] = None,
         migration_prob: float = 0.0,
         emigration_propagator: Type[Propagator] = SelectMin,
-        island_displs: np.ndarray = None,
-        island_counts: np.ndarray = None,
-        rng: random.Random = None,
+        island_displs: Optional[np.ndarray] = None,
+        island_counts: Optional[np.ndarray] = None,
     ) -> None:
         """
         Initialize Propulator with given parameters.
 
         Parameters
         ----------
-        loss_fn: Callable
-                 loss function to be minimized
-        propagator: propulate.propagators.Propagator
-                    propagator to apply for breeding
-        island_idx: int
-                    index of island
-        comm: MPI.Comm
-              intra-island communicator
-        generations: int
-                     number of generations to run
-        checkpoint_path: Union[Path, str]
-                         Path where checkpoints are loaded from and stored.
-        migration_topology: numpy.ndarray
-                            2D matrix where entry (i,j) specifies how many
-                            individuals are sent by island i to island j
-        migration_prob: float
-                        per-worker migration probability
-        emigration_propagator: type[propulate.propagators.Propagator]
-                               emigration propagator, i.e., how to choose individuals
-                               for emigration that are sent to destination island.
-                               Should be some kind of selection operator.
-        island_displs: numpy.ndarray
-                    array with MPI.COMM_WORLD rank of each island's worker 0
-                    Element i specifies MPI.COMM_WORLD rank of worker 0 on island with index i.
-        island_counts: numpy.ndarray
-                       array with number of workers per island
-                       Element i specifies number of workers on island with index i.
-        rng: random.Random
-             random number generator
+        loss_fn : Callable
+            The loss function to be minimized.
+        propagator : propulate.propagators.Propagator
+            The propagator to apply for breeding.
+        rng : random.Random
+            The separate random number generator for the Propulate optimization.
+        island_idx : int, optional
+            The island's index. Default is 0.
+        island_comm : MPI.Comm, optional
+            The intra-island communicator. Default is ``MPI.COMM_WORLD``.
+        propulate_comm : MPI.Comm, optional
+            The Propulate world communicator, consisting of rank 0 of each worker's sub communicator.
+            Default is ``MPI.COMM_WORLD``.
+        worker_sub_comm : MPI.Comm, optional
+            The sub communicator for each (multi rank) worker. Default is ``MPI.COMM_SELF``.
+        generations : int, optional
+            The number of generations to run. Default is -1, i.e., run into wall-clock time limit.
+        checkpoint_path : pathlib.Path | str, optional
+            The path where checkpoints are loaded from and stored. Default is current working directory.
+        migration_topology : numpy.ndarray, optional
+            The migration topology, i.e., a 2D matrix where entry (i,j) specifies how many individuals are sent by
+            island i to island j.
+        migration_prob : float, optional
+            The per-worker migration probability. Default is 0.0.
+        emigration_propagator : Type[propulate.propagators.Propagator], optional
+            The emigration propagator, i.e., how to choose individuals for emigration that are sent to destination
+            island. Should be some kind of selection operator. Default is ``SelectMin``.
+        island_displs : numpy.ndarray, optional
+            An array with ``propulate_comm`` rank of each island's worker 0. Element i specifies the rank of worker 0 on
+            island with index i in the Propulate communicator.
+        island_counts : numpy.ndarray, optional
+            An array with the number of workers per island. Element i specifies the number of workers on island with
+            index i.
         """
         # Set class attributes.
-        self.loss_fn = loss_fn  # callable loss function
-        self.propagator = propagator  # evolutionary propagator
+        self.loss_fn = loss_fn  # Callable loss function
+        self.propagator = propagator  # Evolutionary propagator
         if generations == 0:  # If number of iterations requested == 0.
             if MPI.COMM_WORLD.rank == 0:
                 log.info("Requested number of generations is zero...[RETURN]")
             return
         self.generations = (
-            generations  # number of generations (evaluations per individual)
+            generations  # Number of generations (evaluations per individual)
+        )
+        self.generation = 0  # Current generation not yet evaluated
+        self.island_idx = island_idx  # Island index
+        self.island_comm = island_comm  # Intra-island communicator
+        self.propulate_comm = propulate_comm  # Propulate world communicator
+        self.worker_sub_comm = (
+            worker_sub_comm  # Sub communicator for each (multi rank) worker
         )
-        self.generation = 0  # current generation not yet evaluated
-        self.island_idx = island_idx  # island index
-        self.comm = comm  # intra-island communicator
-        self.checkpoint_path = Path(checkpoint_path)  # checkpoint path
+        if self.propulate_comm is None:  # Exit early for sub-worker only ranks.
+            # These ranks are not used for anything aside from the calculation of the user-defined loss function.
+            return
+        self.checkpoint_path = Path(checkpoint_path)  # Checkpoint path
         self.checkpoint_path.mkdir(parents=True, exist_ok=True)
-        self.migration_prob = migration_prob  # per-rank migration probability
-        self.migration_topology = migration_topology  # migration topology
+        self.migration_prob = migration_prob  # Per-rank migration probability
+        self.migration_topology = migration_topology  # Migration topology
         self.island_displs = (
-            island_displs  # MPI.COMM_WORLD rank of each island's worker 0
+            island_displs  # Propulate world rank of each island's worker
         )
-        self.island_counts = island_counts  # number of workers on each island
-        self.emigration_propagator = emigration_propagator  # emigration propagator
+        self.island_counts = island_counts  # Number of workers on each island
+        self.emigration_propagator = emigration_propagator  # Emigration propagator
         self.rng = rng
 
         # Load initial population of evaluated individuals from checkpoint if exists.
@@ -116,25 +172,25 @@ def __init__(
                             [
                                 x.generation
                                 for x in self.population
-                                if x.rank == self.comm.rank
+                                if x.rank == self.island_comm.rank
                             ]
                         )
                         + 1
-                    )
-                    if self.comm.rank == 0:
+                    )  # Determine generation to be evaluated next from population checkpoint.
+                    if self.island_comm.rank == 0:
                         log.info(
                             "Valid checkpoint file found. "
                             f"Resuming from generation {self.generation} of loaded population..."
                         )
                 except OSError:
                     self.population = []
-                    if self.comm.rank == 0:
+                    if self.island_comm.rank == 0:
                         log.info(
                             "No valid checkpoint file. Initializing population randomly..."
                         )
         else:
             self.population = []
-            if self.comm.rank == 0:
+            if self.island_comm.rank == 0:
                 log.info(
                     "No valid checkpoint file given. Initializing population randomly..."
                 )
@@ -145,10 +201,10 @@ def propulate(self, logging_interval: int = 10, debug: int = 1) -> None:
 
         Parameters
         ----------
-        logging_interval: int
-                          Print each worker's progress every ``logging_interval`` th generation.
-        debug: int
-               verbosity/debug level; 0 - silent; 1 - moderate, 2 - noisy (debug mode)
+        logging_interval : int, optional
+            Print each worker's progress every ``logging_interval``-th generation. Default is 10.
+        debug : int, optional
+            The debug level; 0 - silent; 1 - moderate, 2 - noisy (debug mode). Default is 1.
         """
         self._work(logging_interval, debug)
 
@@ -158,13 +214,12 @@ def _get_active_individuals(self) -> Tuple[List[Individual], int]:
 
         Returns
         -------
-        list[propulate.individual.Individual]
-            currently active individuals in population
+        List[propulate.Individual]
+            All active individuals in the current population.
         int
-            number of currently active individuals
+            The number of currently active individuals.
         """
         active_pop = [ind for ind in self.population if ind.active]
-
         return active_pop, len(active_pop)
 
     def _breed(self) -> Individual:
@@ -173,53 +228,66 @@ def _breed(self) -> Individual:
 
         Returns
         -------
-        propulate.individual.Individual
-            newly bred individual
-        """
-        active_pop, _ = self._get_active_individuals()
-        ind = self.propagator(
-            active_pop
-        )  # Breed new individual from active population.
-        ind.generation = self.generation  # Set generation.
-        ind.rank = self.comm.rank  # Set worker rank.
-        ind.active = True  # If True, individual is active for breeding.
-        ind.island = self.island_idx  # Set birth island.
-        ind.current = self.comm.rank  # Set worker responsible for migration.
-        ind.migration_steps = 0  # Set number of migration steps performed.
-        ind.migration_history = str(self.island_idx)
+        propulate.Individual
+            The newly bred individual.
+        """
+        if (
+            self.propulate_comm is not None
+        ):  # Only processes in the Propulate world communicator, consisting of rank 0 of each worker's sub
+            # communicator, are involved in the actual optimization routine.
+            active_pop, _ = self._get_active_individuals()
+            ind = self.propagator(
+                active_pop
+            )  # Breed new individual from active population.
+            ind.generation = self.generation  # Set generation.
+            ind.rank = self.island_comm.rank  # Set worker rank.
+            ind.active = True  # If True, individual is active for breeding.
+            ind.island = self.island_idx  # Set birth island.
+            ind.current = self.island_comm.rank  # Set worker responsible for migration.
+            ind.migration_steps = 0  # Set number of migration steps performed.
+            ind.migration_history = str(self.island_idx)
+        else:  # The other processes do not breed themselves.
+            ind = None
+
+        if (
+            self.worker_sub_comm != MPI.COMM_SELF
+        ):  # Broadcast newly bred individual to all internal ranks of a worker from rank 0,
+            # which is also part of the Propulate comm.
+            ind = self.worker_sub_comm.bcast(obj=ind, root=0)
 
         return ind  # Return new individual.
 
     def _evaluate_individual(self) -> None:
-        """
-        Breed and evaluate individual.
-        """
+        """Breed and evaluate individual."""
         ind = self._breed()  # Breed new individual.
         start_time = time.time()  # Start evaluation timer.
-        ind.loss = self.loss_fn(ind)  # Evaluate its loss.
+        if self.worker_sub_comm != MPI.COMM_SELF:
+            ind.loss = self.loss_fn(ind, self.worker_sub_comm)  # Evaluate its loss.
+        else:
+            ind.loss = self.loss_fn(ind)
+        if self.propulate_comm is None:
+            return
         ind.evaltime = time.time()  # Stop evaluation timer.
         ind.evalperiod = ind.evaltime - start_time  # Calculate evaluation duration.
         self.population.append(
             ind
         )  # Add evaluated individual to worker-local population.
         log.debug(
-            f"Island {self.island_idx} Worker {self.comm.rank} Generation {self.generation}: BREEDING\n"
-            f"Bred and evaluated individual {ind}.\n"
+            f"Island {self.island_idx} Worker {self.island_comm.rank} Generation {self.generation}: BREEDING\n"
+            f"Bred and evaluated individual {ind}."
         )
-
         # Tell other workers in own island about results to synchronize their populations.
-        for r in range(self.comm.size):  # Loop over ranks in intra-island communicator.
-            if r == self.comm.rank:
+        for r in range(
+            self.island_comm.size
+        ):  # Loop over ranks in intra-island communicator.
+            if r == self.island_comm.rank:
                 continue  # No self-talk.
-            self.comm.send(copy.deepcopy(ind), dest=r, tag=INDIVIDUAL_TAG)
+            self.island_comm.send(copy.deepcopy(ind), dest=r, tag=INDIVIDUAL_TAG)
 
     def _receive_intra_island_individuals(self) -> None:
-        """
-        Check for and possibly receive incoming individuals
-        evaluated by other workers within own island.
-        """
+        """Check for and possibly receive incoming individuals evaluated by other workers within own island."""
         log_string = (
-            f"Island {self.island_idx} Worker {self.comm.rank} Generation {self.generation}: "
+            f"Island {self.island_idx} Worker {self.island_comm.rank} Generation {self.generation}: "
             f"INTRA-ISLAND SYNCHRONIZATION\n"
         )
         probe_ind = True
@@ -227,7 +295,7 @@ def _receive_intra_island_individuals(self) -> None:
             stat = (
                 MPI.Status()
             )  # Retrieve status of reception operation, including source and tag.
-            probe_ind = self.comm.iprobe(
+            probe_ind = self.island_comm.iprobe(
                 source=MPI.ANY_SOURCE, tag=INDIVIDUAL_TAG, status=stat
             )
             # If True, continue checking for incoming messages. Tells whether message corresponding
@@ -236,20 +304,22 @@ def _receive_intra_island_individuals(self) -> None:
             log_string += f"Incoming individual to receive?...{probe_ind}\n"
             if probe_ind:
                 # Receive individual and add it to own population.
-                ind_temp = self.comm.recv(source=stat.Get_source(), tag=INDIVIDUAL_TAG)
+                ind_temp = self.island_comm.recv(
+                    source=stat.Get_source(), tag=INDIVIDUAL_TAG
+                )
                 self.population.append(
                     ind_temp
                 )  # Add received individual to own worker-local population.
                 log_string += f"Added individual {ind_temp} from W{stat.Get_source()} to own population.\n"
         _, n_active = self._get_active_individuals()
         log_string += (
-            f"After probing within island: {n_active}/{len(self.population)} active.\n"
+            f"After probing within island: {n_active}/{len(self.population)} active."
         )
         log.debug(log_string)
 
     def _send_emigrants(self) -> None:
         """
-        Perform migration, i.e. island sends individuals out to other islands.
+        Perform migration, i.e., island sends individuals out to other islands.
 
         Raises
         ------
@@ -277,8 +347,8 @@ def _get_unique_individuals(self) -> List[Individual]:
 
         Returns
         -------
-        list[propulate.individual.Individual]
-            unique individuals
+        List[propulate.Individual]
+            All unique individuals in the current population.
         """
         unique_inds = []
         for individual in self.population:
@@ -301,8 +371,8 @@ def _check_intra_island_synchronization(
 
         Parameters
         ----------
-        populations: list[list[propulate.individual.Individual]]
-                     list of islands' sorted population lists
+        populations : List[List[propulate.Individual]]
+            A list of all islands' sorted population lists.
 
         Returns
         -------
@@ -317,22 +387,22 @@ def _check_intra_island_synchronization(
             if len(difference) == 0:
                 continue
             log.info(
-                f"Island {self.island_idx} Worker {self.comm.rank}: Population not synchronized:\n"
-                f"{difference}\n"
+                f"Island {self.island_idx} Worker {self.island_comm.rank}: Population not synchronized:\n"
+                f"{difference}"
             )
             synchronized = False
         return synchronized
 
-    def _work(self, logging_interval: int, debug: int):
+    def _work(self, logging_interval: int = 10, debug: int = -1):
         """
         Execute evolutionary algorithm in parallel.
 
         Parameters
         ----------
-        logging_interval: int
-                          logging interval
-        debug: int
-               verbosity/debug level; 0 - silent; 1 - moderate, 2 - noisy (debug mode)
+        logging_interval : int, optional
+            The logging interval. Default is 10.
+        debug : int, optional
+            The debug level; 0 - silent; 1 - moderate, 2 - noisy (debug mode). Default is 1.
 
         Raises
         ------
@@ -340,18 +410,26 @@ def _work(self, logging_interval: int, debug: int):
             If any individuals are left that should have been deactivated before (only for debug > 0).
 
         """
+        if self.worker_sub_comm != MPI.COMM_SELF:
+            self.generation = self.worker_sub_comm.bcast(self.generation, root=0)
+        if self.propulate_comm is None:
+            while self.generations <= -1 or self.generation < self.generations:
+                # Breed and evaluate individual.
+                self._evaluate_individual()
+                self.generation += 1
+            return
 
-        if self.comm.rank == 0:
-            log.info(f"Island {self.island_idx} has {self.comm.size} workers.")
+        if self.island_comm.rank == 0:
+            log.info(f"Island {self.island_idx} has {self.island_comm.size} workers.")
 
-        dump = True if self.comm.rank == 0 else False
-        MPI.COMM_WORLD.barrier()
+        dump = True if self.island_comm.rank == 0 else False
+        self.propulate_comm.barrier()
 
         # Loop over generations.
         while self.generations <= -1 or self.generation < self.generations:
             if self.generation % int(logging_interval) == 0:
                 log.info(
-                    f"Island {self.island_idx} Worker {self.comm.rank}: In generation {self.generation}..."
+                    f"Island {self.island_idx} Worker {self.island_comm.rank}: In generation {self.generation}..."
                 )
 
             # Breed and evaluate individual.
@@ -375,29 +453,26 @@ def _work(self, logging_interval: int, debug: int):
         # so that each of them holds the complete final population and the found optimum
         # irrespective of the order they finished.
 
-        MPI.COMM_WORLD.barrier()
-        if MPI.COMM_WORLD.rank == 0:
-            log.info("OPTIMIZATION DONE.")
-            log.info("NEXT: Final checks for incoming messages...")
-        MPI.COMM_WORLD.barrier()
+        self.propulate_comm.barrier()
+        if self.propulate_comm.rank == 0:
+            log.info("OPTIMIZATION DONE.\nNEXT: Final checks for incoming messages...")
+        self.propulate_comm.barrier()
 
         # Final check for incoming individuals evaluated by other intra-island workers.
         self._receive_intra_island_individuals()
-        MPI.COMM_WORLD.barrier()
+        self.propulate_comm.barrier()
 
         # Final checkpointing on rank 0.
-        if self.comm.rank == 0:
+        if self.island_comm.rank == 0:
             self._dump_final_checkpoint()  # Dump checkpoint.
-        MPI.COMM_WORLD.barrier()
+        self.propulate_comm.barrier()
         _ = self._determine_worker_dumping_next()
-        MPI.COMM_WORLD.barrier()
+        self.propulate_comm.barrier()
 
     def _dump_checkpoint(self):
-        """
-        Dump checkpoint to file.
-        """
+        """Dump checkpoint to file."""
         log.debug(
-            f"Island {self.island_idx} Worker {self.comm.rank} Generation {self.generation}: "
+            f"Island {self.island_idx} Worker {self.island_comm.rank} Generation {self.generation}: "
             f"Dumping checkpoint..."
         )
         save_ckpt_file = self.checkpoint_path / f"island_{self.island_idx}_ckpt.pickle"
@@ -409,28 +484,30 @@ def _dump_checkpoint(self):
         with open(save_ckpt_file, "wb") as f:
             pickle.dump(self.population, f)
 
-        dest = self.comm.rank + 1 if self.comm.rank + 1 < self.comm.size else 0
-        self.comm.send(True, dest=dest, tag=DUMP_TAG)
+        dest = (
+            self.island_comm.rank + 1
+            if self.island_comm.rank + 1 < self.island_comm.size
+            else 0
+        )
+        self.island_comm.send(True, dest=dest, tag=DUMP_TAG)
 
     def _determine_worker_dumping_next(self):
-        """
-        Determine the worker who dumps the checkpoint in the next generation.
-        """
+        """Determine the worker who dumps the checkpoint in the next generation."""
         dump = False
         stat = MPI.Status()
-        probe_dump = self.comm.iprobe(source=MPI.ANY_SOURCE, tag=DUMP_TAG, status=stat)
+        probe_dump = self.island_comm.iprobe(
+            source=MPI.ANY_SOURCE, tag=DUMP_TAG, status=stat
+        )
         if probe_dump:
-            dump = self.comm.recv(source=stat.Get_source(), tag=DUMP_TAG)
+            dump = self.island_comm.recv(source=stat.Get_source(), tag=DUMP_TAG)
             log.debug(
-                f"Island {self.island_idx} Worker {self.comm.rank} Generation {self.generation}: "
+                f"Island {self.island_idx} Worker {self.island_comm.rank} Generation {self.generation}: "
                 f"Going to dump next: {dump}. Before: Worker {stat.Get_source()}"
             )
         return dump
 
     def _dump_final_checkpoint(self):
-        """
-        Dump final checkpoint.
-        """
+        """Dump final checkpoint."""
         save_ckpt_file = self.checkpoint_path / f"island_{self.island_idx}_ckpt.pickle"
         if os.path.isfile(save_ckpt_file):
             try:
@@ -441,7 +518,7 @@ def _dump_final_checkpoint(self):
                 pickle.dump(self.population, f)
 
     def _check_for_duplicates(
-        self, active: bool, debug: int
+        self, active: bool, debug: int = 1
     ) -> Tuple[List[List[Union[Individual, int]]], List[Individual]]:
         """
         Check for duplicates in current population.
@@ -451,17 +528,17 @@ def _check_for_duplicates(
 
         Parameters
         ----------
-        active: bool
-                Whether to consider active individuals (True) or all individuals (False)
-        debug: int
-               verbosity/debug level; 0 - silent; 1 - moderate, 2 - noisy (debug mode)
+        active : bool
+            Whether to consider active individuals (True) or all individuals (False).
+        debug : int, optional
+            The debug level; 0 - silent; 1 - moderate, 2 - noisy (debug mode). Default is 1.
 
         Returns
         -------
-        list[list[propulate.individual.Individual | int]]
-            individuals and their occurrences
-        list[propulate.individual.Individual]
-            unique individuals in population
+        List[List[propulate.Individual | int]]
+            The individuals and their occurrences.
+        List[propulate.Individual]
+            The unique individuals in the population.
         """
         if active:
             population, _ = self._get_active_individuals()
@@ -478,7 +555,7 @@ def _check_for_duplicates(
             if not considered:
                 num_copies = population.count(individual)
                 log.debug(
-                    f"Island {self.island_idx} Worker {self.comm.rank} Generation {self.generation}: "
+                    f"Island {self.island_idx} Worker {self.island_comm.rank} Generation {self.generation}: "
                     f"{individual} occurs {num_copies} time(s)."
                 )
                 unique_inds.append(individual)
@@ -487,47 +564,47 @@ def _check_for_duplicates(
 
     def summarize(
         self, top_n: int = 1, debug: int = 1
-    ) -> List[Union[List[Individual], Individual]]:
+    ) -> Union[List[Union[List[Individual], Individual]], None]:
         """
-        Get top-n results from propulate optimization.
+        Get top-n results from Propulate optimization.
 
         Parameters
         ----------
-        top_n: int
-               number of best results to report
-        debug: int
-               verbosity/debug level; 0 - silent; 1 - moderate, 2 - noisy (debug mode)
+        top_n : int, optional
+            The number of best results to report. Default is 1.
+        debug : int, optional
+            The debug level; 0 - silent; 1 - moderate, 2 - noisy (debug mode). Default is 1.
 
         Returns
         -------
-        list[list[Individual] | Individual]
-            top-n best individuals on each island
+        List[List[propulate.Individual] | propulate.Individual]
+            The top-n best individuals on each island.
         """
+        if self.propulate_comm is None:
+            return
         active_pop, num_active = self._get_active_individuals()
         assert np.all(
-            np.array(self.comm.allgather(num_active), dtype=int) == num_active
+            np.array(self.island_comm.allgather(num_active), dtype=int) == num_active
         )
         if self.island_counts is not None:
             num_active = int(
-                MPI.COMM_WORLD.allreduce(
+                self.propulate_comm.allreduce(
                     num_active / self.island_counts[self.island_idx]
                 )
             )
 
-        MPI.COMM_WORLD.barrier()
-        if MPI.COMM_WORLD.rank == 0:
-            log.info("###########")
-            log.info("# SUMMARY #")
-            log.info("###########")
-            log.info(f"Number of currently active individuals is {num_active}. ")
+        self.propulate_comm.barrier()
+        if self.propulate_comm.rank == 0:
             log.info(
-                f"Expected overall number of evaluations is {self.generations*MPI.COMM_WORLD.size}."
+                "###########\n# SUMMARY #\n###########\n"
+                f"Number of currently active individuals is {num_active}.\n"
+                f"Expected overall number of evaluations is {self.generations*self.propulate_comm.size}."
             )
         # Only double-check number of occurrences of each individual for DEBUG level 2.
         if debug == 2:
-            populations = self.comm.gather(self.population, root=0)
+            populations = self.island_comm.gather(self.population, root=0)
             occurrences, _ = self._check_for_duplicates(True, debug)
-            if self.comm.rank == 0:
+            if self.island_comm.rank == 0:
                 if self._check_intra_island_synchronization(populations):
                     log.info(
                         f"Island {self.island_idx}: Populations among workers synchronized."
@@ -540,18 +617,18 @@ def summarize(
                     f"Island {self.island_idx}: {len(active_pop)}/{len(self.population)} "
                     f"individuals active ({len(occurrences)} unique)"
                 )
-        MPI.COMM_WORLD.barrier()
+        self.propulate_comm.barrier()
         if debug == 0:
             best = min(self.population, key=attrgetter("loss"))
-            if self.comm.rank == 0:
-                log.info(f"Top result on island {self.island_idx}: {best}\n")
+            if self.island_comm.rank == 0:
+                log.info(f"Top result on island {self.island_idx}: {best}")
         else:
             unique_pop = self._get_unique_individuals()
             unique_pop.sort(key=lambda x: x.loss)
             best = unique_pop[:top_n]
-            if self.comm.rank == 0:
+            if self.island_comm.rank == 0:
                 res_str = f"Top {top_n} result(s) on island {self.island_idx}:\n"
                 for i in range(top_n):
                     res_str += f"({i+1}): {unique_pop[i]}\n"
                 log.info(res_str)
-        return MPI.COMM_WORLD.allgather(best)
+        return self.propulate_comm.allgather(best)
diff --git a/propulate/utils.py b/propulate/utils.py
index 0815f8a8..6efaa117 100644
--- a/propulate/utils.py
+++ b/propulate/utils.py
@@ -3,7 +3,7 @@
 import random
 import sys
 from pathlib import Path
-from typing import Dict, Union, Tuple
+from typing import Dict, Optional, Union, Tuple
 
 import colorlog
 import numpy as np
@@ -15,7 +15,7 @@
     Conditional,
     InitUniform,
     IntervalMutationNormal,
-    MateUniform,
+    CrossoverUniform,
     PointMutation,
     Propagator,
     SelectMin,
@@ -30,31 +30,32 @@ def get_default_propagator(
         Dict[str, Tuple[int, int]],
         Dict[str, Tuple[str, ...]],
     ],
-    mate_prob: float,
-    mut_prob: float,
-    random_prob: float,
+    crossover_prob: float = 0.7,
+    mutation_prob: float = 0.4,
+    random_init_prob: float = 0.1,
     sigma_factor: float = 0.05,
-    rng: random.Random = None,
+    rng: Optional[random.Random] = None,
 ) -> Propagator:
     """
     Get Propulate's default evolutionary optimization propagator.
 
     Parameters
     ----------
-    pop_size: int
-              number of individuals in breeding population
-    limits: dict
-            (hyper-)parameters to be optimized, i.e., search space
-    mate_prob: float
-               uniform-crossover probability
-    mut_prob: float
-              point-mutation probability
-    random_prob: float
-                 random-initialization probability
-    sigma_factor: float
-                  scaling factor for obtaining std from search-space boundaries for interval mutation
-    rng: random.Random
-         random number generator
+    pop_size : int
+        The number of individuals in the breeding population.
+    limits : Dict[str, Tuple[float, float]] | Dict[str, Tuple[int, int]] | Dict[str, Tuple[str, ...]]
+        The (hyper-)parameters to be optimized, i.e., the search space.
+    crossover_prob : float, optional
+        The uniform-crossover probability. Default is 0.7.
+    mutation_prob : float, optional
+        The point-mutation probability. Default is 0.4.
+    random_init_prob : float, optional
+        The random-initialization probability. Default is 0.1.
+    sigma_factor : float
+        The scaling factor for obtaining the standard deviation from the search-space boundaries for interval mutation.
+        Default is 0.05.
+    rng : random.Random, optional
+        The separate random number generator for the Propulate optimization.
 
     Returns
     -------
@@ -68,12 +69,12 @@ def get_default_propagator(
             [  # Compose propagator out of basic evolutionary operators with Compose(...).
                 SelectMin(pop_size),
                 SelectUniform(offspring=2, rng=rng),
-                MateUniform(mate_prob, rng=rng),
-                PointMutation(limits, probability=mut_prob, rng=rng),
+                CrossoverUniform(crossover_prob, rng=rng),
+                PointMutation(limits, probability=mutation_prob, rng=rng),
                 IntervalMutationNormal(
                     limits, sigma_factor=sigma_factor, probability=1.0, rng=rng
                 ),
-                InitUniform(limits, parents=1, probability=random_prob, rng=rng),
+                InitUniform(limits, parents=1, probability=random_init_prob, rng=rng),
             ]
         )
     else:
@@ -81,9 +82,9 @@ def get_default_propagator(
             [  # Compose propagator out of basic evolutionary operators with Compose(...).
                 SelectMin(pop_size),
                 SelectUniform(offspring=2, rng=rng),
-                MateUniform(mate_prob, rng=rng),
-                PointMutation(limits, probability=mut_prob, rng=rng),
-                InitUniform(limits, parents=1, probability=random_prob, rng=rng),
+                CrossoverUniform(crossover_prob, rng=rng),
+                PointMutation(limits, probability=mutation_prob, rng=rng),
+                InitUniform(limits, parents=1, probability=random_init_prob, rng=rng),
             ]
         )
 
@@ -91,37 +92,31 @@ def get_default_propagator(
     propagator = Conditional(
         pop_size, propagator, init
     )  # Initialize random if population size < specified `pop_size`.
-
     return propagator
 
 
 def set_logger_config(
     level: int = logging.INFO,
-    log_file: Union[str, Path] = None,
+    log_file: Optional[Union[str, Path]] = None,
     log_to_stdout: bool = True,
     log_rank: bool = False,
     colors: bool = True,
 ) -> None:
     """
-    Set up the logger. Should only need to be done once.
-    Generally, logging should only be done on the master rank.
+    Set up the logger. Should only need to be done once. Generally, logging should only be done on the master rank.
 
     Parameters
     ----------
-    level: logging.INFO, ...
-           default level for logging
-           default: INFO
-    log_file: str, Path
-              file to save the log to
-              default: None
-    log_to_stdout: bool
-                   flag indicating if the log should be printed on stdout
-                   default: True
-    log_rank: bool
-              flag for prepending the MPI rank to the logging message
-    colors: bool
-            flag for using colored logs
-            default: True
+    level : int
+        The default level for logging. Default is ``logging.INFO``.
+    log_file : str | Path, optional
+        The file to save the log to.
+    log_to_stdout : bool
+        A flag indicating if the log should be printed on stdout. Default is True.
+    log_rank : bool
+        A flag for prepending the MPI rank to the logging message. Default is False.
+    colors : bool
+        A flag for using colored logs. Default is True.
     """
     rank = f"{MPI.COMM_WORLD.Get_rank()}:" if log_rank else ""
     # Get base logger for Propulate.
@@ -157,7 +152,6 @@ def set_logger_config(
         file_handler.setFormatter(simple_formatter)
         base_logger.addHandler(file_handler)
     base_logger.setLevel(level)
-    return
 
 
 def make_particle(individual: Individual) -> Particle:
@@ -166,13 +160,13 @@ def make_particle(individual: Individual) -> Particle:
 
     Parameters
     ----------
-    individual: Individual
-                Individual to be converted to a particle
+    individual : propulate.Individual
+        An individual to be converted to a particle.
 
     Returns
     --------
-    Particle
-        Converted individual
+    propulate.Particle
+        The converted individual.
     """
     p = Particle(generation=individual.generation)
     p.position = np.zeros(len(individual))
diff --git a/tutorials/cmaes_example.py b/tutorials/cmaes_example.py
index 64ec2d66..aca7765f 100644
--- a/tutorials/cmaes_example.py
+++ b/tutorials/cmaes_example.py
@@ -1,14 +1,13 @@
-#!/usr/bin/env python3
+"""Simple example script using CMA-ES"""
 import random
-import argparse
-import logging
+import pathlib
 
 from mpi4py import MPI
 
 from propulate import Propulator
 from propulate.propagators import BasicCMA, ActiveCMA, CMAPropagator
 from propulate.utils import set_logger_config
-from function_benchmark import get_function_search_space
+from function_benchmark import get_function_search_space, parse_arguments
 
 
 if __name__ == "__main__":
@@ -21,51 +20,12 @@
             "#################################################\n"
         )
 
-    parser = argparse.ArgumentParser(
-        prog="Simple Propulator example",
-        description="Set up and run a basic Propulator optimization of mathematical functions.",
-    )
-    parser.add_argument(  # Function to optimize
-        "-f",
-        "--function",
-        type=str,
-        choices=[
-            "bukin",
-            "eggcrate",
-            "himmelblau",
-            "keane",
-            "leon",
-            "rastrigin",
-            "schwefel",
-            "sphere",
-            "step",
-            "rosenbrock",
-            "quartic",
-            "bisphere",
-            "birastrigin",
-            "griewank",
-        ],
-        default="sphere",
-    )
-    parser.add_argument(
-        "-g", "--generations", type=int, default=1000
-    )  # Number of generations
-    parser.add_argument(
-        "-s", "--seed", type=int, default=0
-    )  # Seed for Propulate random number generator
-    parser.add_argument("-a", "--adapter", type=str, default="basic")
-    parser.add_argument("-v", "--verbosity", type=int, default=1)  # Verbosity level
-    parser.add_argument(
-        "-ckpt", "--checkpoint", type=str, default="./"
-    )  # Path for loading and writing checkpoints.
-    parser.add_argument("-t", "--top_n", type=int, default=1)
-    parser.add_argument("-l", "--logging_int", type=int, default=10)
-    config = parser.parse_args()
+    config, _ = parse_arguments(comm)
 
     # Set up separate logger for Propulate optimization.
     set_logger_config(
-        level=logging.INFO,  # logging level
-        log_file=f"{config.checkpoint}/propulator.log",  # logging path
+        level=config.logging_level,  # Logging level
+        log_file=f"{config.checkpoint}/{pathlib.Path(__file__).stem}.log",  # Logging path
         log_to_stdout=True,  # Print log on stdout.
         log_rank=False,  # Do not prepend MPI rank to logging messages.
         colors=True,  # Use colors.
@@ -92,12 +52,14 @@
     propulator = Propulator(
         loss_fn=function,
         propagator=propagator,
-        comm=comm,
+        rng=rng,
+        island_comm=comm,
         generations=config.generations,
         checkpoint_path=config.checkpoint,
-        rng=rng,
     )
 
     # Run optimization and print summary of results.
-    propulator.propulate(logging_interval=config.logging_int, debug=config.verbosity)
+    propulator.propulate(
+        logging_interval=config.logging_interval, debug=config.verbosity
+    )
     propulator.summarize(top_n=config.top_n, debug=config.verbosity)
diff --git a/tutorials/function_benchmark.py b/tutorials/function_benchmark.py
index 6a381910..8a704d40 100644
--- a/tutorials/function_benchmark.py
+++ b/tutorials/function_benchmark.py
@@ -1,4 +1,9 @@
+"""Benchmark function module."""
+import argparse
+import logging
 from typing import Callable, Dict, Tuple
+
+from mpi4py import MPI
 import numpy as np
 
 
@@ -11,13 +16,13 @@ def rosenbrock(params: Dict[str, float]) -> float:
 
     Parameters
     ----------
-    params: dict[str, float]
-            function parameters
+    params : Dict[str, float]
+        The function parameters.
 
     Returns
     -------
     float
-        function value
+        The function value.
     """
     params = np.array(list(params.values()))
     return 100 * (params[0] ** 2 - params[1]) ** 2 + (1 - params[0]) ** 2
@@ -37,16 +42,16 @@ def step(params: Dict[str, float]) -> float:
 
     Parameters
     ----------
-    params: dict[str, float]
-            function parameters
+    params : Dict[str, float]
+        The function parameters.
 
     Returns
     -------
     float
-        function value
+        The function value.
     """
     params = np.array(list(params.values()))
-    return np.sum(params.astype(int), dtype=float)
+    return np.sum(params.astype(int), dtype=float).item()
 
 
 def quartic(params: Dict[str, float]) -> float:
@@ -62,13 +67,13 @@ def quartic(params: Dict[str, float]) -> float:
 
     Parameters
     ----------
-    params: dict[str, float]
-            function parameters
+    params : Dict[str, float]
+        The function parameters.
 
     Returns
     -------
     float
-        function value
+        The function value.
     """
     params = np.array(list(params.values()))
     idx = np.arange(1, len(params) + 1)
@@ -90,13 +95,13 @@ def rastrigin(params: Dict[str, float]) -> float:
 
     Parameters
     ----------
-    params: dict[str, float]
-            function parameters
+    params : Dict[str, float]
+        The function parameters.
 
     Returns
     -------
     float
-        function value
+        The function value.
     """
     a = 10.0
     params = np.array(list(params.values()))
@@ -107,7 +112,7 @@ def griewank(params: Dict[str, float]) -> float:
     """
     Griewank function.
 
-    Griewank's product creates sub-populations strongly codependent to parallel GAs, while the summation produces a
+    Griewank's product creates subpopulations strongly codependent to parallel GAs, while the summation produces a
     parabola. Its local optima lie above parabola level but decrease with increasing dimensions, i.e., the larger the
     search range, the flatter the function.
 
@@ -117,13 +122,13 @@ def griewank(params: Dict[str, float]) -> float:
 
     Parameters
     ----------
-    params: dict[str, float]
-            function parameters
+    params : Dict[str, float]
+        The function parameters.
 
     Returns
     -------
     float
-        function value
+        The function value.
     """
     params = np.array(list(params.values()))
     idx = np.arange(1, len(params) + 1)
@@ -142,13 +147,13 @@ def schwefel(params: Dict[str, float]) -> float:
 
     Parameters
     ----------
-    params: dict[str, float]
-            function parameters
+    params : Dict[str, float]
+        The function parameters.
 
     Returns
     -------
     float
-        function value
+        The function value.
     """
     v = 418.982887
     params = np.array(list(params.values()))
@@ -175,13 +180,13 @@ def bisphere(params: Dict[str, float]) -> float:
 
     Parameters
     ----------
-    params: dict[str, float]
-            function parameters
+    params : Dict[str, float]
+        The function parameters.
 
     Returns
     -------
     float
-        function value
+        The function value.
     """
     params = np.array(list(params.values()))
     n = len(params)
@@ -189,7 +194,9 @@ def bisphere(params: Dict[str, float]) -> float:
     s = 1 - np.sqrt(1 / (2 * np.sqrt(n + 20) - 8.2))
     mu1 = 2.5
     mu2 = -np.sqrt((mu1**2 - d) / s)
-    return min(np.sum((params - mu1) ** 2), d * n + s * np.sum((params - mu2) ** 2))
+    return min(
+        np.sum((params - mu1) ** 2), d * n + s * np.sum((params - mu2) ** 2)
+    ).item()
 
 
 def birastrigin(params: Dict[str, float]) -> float:
@@ -210,13 +217,13 @@ def birastrigin(params: Dict[str, float]) -> float:
 
     Parameters
     ----------
-    params: dict[str, float]
-            function parameters
+    params : Dict[str, float]
+        The function parameters.
 
     Returns
     -------
     float
-        function value
+        The function value.
     """
     params = np.array(list(params.values()))
     n = len(params)
@@ -238,13 +245,13 @@ def bukin_n6(params: Dict[str, float]) -> float:
 
     Parameters
     ----------
-    params: dict[str, float]
-            function parameters
+    params : Dict[str, float]
+        The function parameters.
 
     Returns
     -------
     float
-        function value
+        The function value.
     """
     params = np.array(list(params.values()))
     return 100 * np.sqrt(np.abs(params[1] - 0.01 * params[0] ** 2)) + 0.01 * np.abs(
@@ -261,13 +268,13 @@ def egg_crate(params: Dict[str, float]) -> float:
 
     Parameters
     ----------
-    params: dict[str, float]
-            function parameters
+    params : Dict[str, float]
+        The function parameters.
 
     Returns
     -------
     float
-        function value
+        The function value.
     """
     params = np.array(list(params.values()))
     return (
@@ -286,13 +293,13 @@ def himmelblau(params: Dict[str, float]) -> float:
 
     Parameters
     ----------
-    params: dict[str, float]
-            function parameters
+    params : Dict[str, float]
+        The function parameters.
 
     Returns
     -------
     float
-        function value
+        The function value.
     """
     params = np.array(list(params.values()))
     return (params[0] ** 2 + params[1] - 11) ** 2 + (
@@ -309,13 +316,13 @@ def keane(params: Dict[str, float]) -> float:
 
     Parameters
     ------
-    params: dict[str, float]
-            function parameters
+    params : Dict[str, float]
+        The function parameters.
 
     Returns
     -------
     float
-        function value
+        The function value.
     """
     params = np.array(list(params.values()))
     return (
@@ -334,13 +341,13 @@ def leon(params: Dict[str, float]) -> float:
 
     Parameters
     ----------
-    params: dict[str, float]
-            function parameters
+    params : Dict[str, float]
+        The function parameters.
 
     Returns
     -------
     float
-        function value
+        The function value.
     """
     params = np.array(list(params.values()))
     return 100 * (params[1] - params[0] ** 3) ** 2 + (1 - params[0]) ** 2
@@ -355,14 +362,15 @@ def sphere(params: Dict[str, float]) -> float:
 
     Parameters
     ----------
-    params: dict[str, float]
-            function parameters
+    params : Dict[str, float]
+        The function parameters.
+
     Returns
     -------
     float
-        function value
+        The function value.
     """
-    return np.sum(np.array(list(params.values())) ** 2)
+    return np.sum(np.array(list(params.values())) ** 2).item()
 
 
 def get_function_search_space(
@@ -374,14 +382,14 @@ def get_function_search_space(
     Parameters
     ----------
     fname: str
-           function name
+        The function name.
 
     Returns
     -------
     Callable
-        function
-    dict[str, tuple[float, float]]
-        search space
+        The callable function.
+    Dict[str, tuple[float, float]]
+        The search space.
     """
     if fname == "bukin":
         function = bukin_n6
@@ -593,3 +601,127 @@ def get_function_search_space(
         ValueError(f"Function {fname} undefined...exiting")
 
     return function, limits
+
+
+def parse_arguments(
+    propulate_comm: MPI.Comm = MPI.COMM_WORLD,
+) -> Tuple[argparse.Namespace, Dict[str, bool]]:
+    """
+    Set up argument parser for Propulate optimization of simple mathematical functions.
+
+    Parameters
+    ----------
+    propulate_comm : MPI.Comm, optional
+        The communicator used to run the Propulate optimization. Default is ``MPI.COMM_WORLD``.
+
+    Returns
+    -------
+    Namespace
+        The namespace of all parsed arguments.
+    Dict[str, bool]
+        A dictionary logging if one of the PSO hyperparameters was actually set. Only relevant for PSO.
+    """
+    parser = argparse.ArgumentParser(
+        prog="Simple Propulator example",
+        description="Set up and run a basic Propulator optimization of mathematical functions.",
+    )
+    parser.add_argument(  # Function to optimize
+        "--function",
+        type=str,
+        choices=[
+            "bukin",
+            "eggcrate",
+            "himmelblau",
+            "keane",
+            "leon",
+            "rastrigin",
+            "schwefel",
+            "sphere",
+            "step",
+            "rosenbrock",
+            "quartic",
+            "bisphere",
+            "birastrigin",
+            "griewank",
+        ],
+        default="sphere",
+    )
+    parser.add_argument(
+        "--generations", type=int, default=1000
+    )  # Number of generations
+    parser.add_argument(
+        "--seed", type=int, default=0
+    )  # Seed for Propulate random number generator
+    parser.add_argument("--verbosity", type=int, default=1)  # Verbosity level
+    parser.add_argument(
+        "--checkpoint", type=str, default="./"
+    )  # Path for loading and writing checkpoints.
+    parser.add_argument(
+        "--pop_size", type=int, default=2 * propulate_comm.size
+    )  # Breeding pool size
+    parser.add_argument(
+        "--crossover_probability", type=float, default=0.7
+    )  # Crossover probability
+    parser.add_argument(
+        "--mutation_probability", type=float, default=0.4
+    )  # Mutation probability
+    parser.add_argument("--random_init_probability", type=float, default=0.1)
+    parser.add_argument("--top_n", type=int, default=1)
+    parser.add_argument("--logging_interval", type=int, default=10)
+    parser.add_argument("--logging_level", type=int, default=logging.INFO)
+
+    # -------- Island-model specific arguments (ignored if not needed) --------
+    parser.add_argument(
+        "--num_islands", type=int, default=2
+    )  # Number of separate evolutionary islands
+    parser.add_argument(
+        "--migration_probability", type=float, default=0.9
+    )  # Migration probability
+    parser.add_argument("--num_migrants", type=int, default=1)
+    parser.add_argument("--pollination", action="store_true")
+
+    # -------- PSO-specific arguments (ignored if not needed) --------
+    parser.add_argument(
+        "--variant",
+        type=str,
+        choices=["Basic", "VelocityClamping", "Constriction", "Canonical"],
+        default="Basic",
+    )  # PSO variant to run
+    hp_set: Dict[str, bool] = {
+        "inertia": False,
+        "cognitive": False,
+        "social": False,
+    }
+
+    class ParamSettingCatcher(argparse.Action):
+        """
+        This class extends ``argparse``'s ``Action`` class in order to allow for an action that logs if one of the PSO
+        hyperparameters was actually set.
+        """
+
+        def __call__(self, parser, namespace, values, option_string=None):
+            hp_set[self.dest] = True
+            super().__call__(parser, namespace, values, option_string)
+
+    parser.add_argument(
+        "--inertia", type=float, default=0.729, action=ParamSettingCatcher
+    )  # Inertia weight
+    parser.add_argument(
+        "--cognitive", type=float, default=1.49445, action=ParamSettingCatcher
+    )  # Cognitive factor
+    parser.add_argument(
+        "--social", type=float, default=1.49445, action=ParamSettingCatcher
+    )  # Social factor
+    parser.add_argument(
+        "--clamping_factor", type=float, default=0.6
+    )  # Velocity clamping factor
+
+    # -------- CMA-ES specific arguments (ignored if not needed)
+    parser.add_argument("--adapter", type=str, default="basic")
+
+    # -------- Multi-rank worker specific arguments (ignored if not needed)
+    parser.add_argument(
+        "--ranks_per_worker", type=int, default=2
+    )  # Number of sub ranks that each worker will use
+
+    return parser.parse_args(), hp_set
diff --git a/tutorials/islands_example.py b/tutorials/islands_example.py
index e6b595c7..71720391 100755
--- a/tutorials/islands_example.py
+++ b/tutorials/islands_example.py
@@ -1,6 +1,5 @@
-#!/usr/bin/env python3
-import argparse
-import logging
+"""Simple island model example script."""
+import pathlib
 import random
 
 import numpy as np
@@ -9,75 +8,18 @@
 from propulate import Islands
 from propulate.propagators import SelectMin, SelectMax
 from propulate.utils import get_default_propagator, set_logger_config
-from function_benchmark import get_function_search_space
+from function_benchmark import get_function_search_space, parse_arguments
 
 if __name__ == "__main__":
     comm = MPI.COMM_WORLD
 
-    parser = argparse.ArgumentParser(
-        prog="Simple Propulator example",
-        description="Set up and run a basic Propulator optimization of mathematical functions.",
-    )
-    parser.add_argument(  # Function to optimize
-        "-f",
-        "--function",
-        type=str,
-        choices=[
-            "bukin",
-            "eggcrate",
-            "himmelblau",
-            "keane",
-            "leon",
-            "rastrigin",
-            "schwefel",
-            "sphere",
-            "step",
-            "rosenbrock",
-            "quartic",
-            "bisphere",
-            "birastrigin",
-            "griewank",
-        ],
-        default="sphere",
-    )
-    parser.add_argument(
-        "-g", "--generations", type=int, default=1000
-    )  # Number of generations
-    parser.add_argument(
-        "-s", "--seed", type=int, default=0
-    )  # Seed for Propulate random number generator
-    parser.add_argument("-v", "--verbosity", type=int, default=1)  # Verbosity level
-    parser.add_argument(
-        "-ckpt", "--checkpoint", type=str, default="./"
-    )  # Path for loading and writing checkpoints.
-    parser.add_argument(
-        "-p", "--pop_size", type=int, default=2 * comm.size
-    )  # Breeding pool size
-    parser.add_argument(
-        "-cp", "--crossover_probability", type=float, default=0.7
-    )  # Crossover probability
-    parser.add_argument(
-        "-mp", "--mutation_probability", type=float, default=0.4
-    )  # Mutation probability
-    parser.add_argument("-rp", "--random_init_probability", type=float, default=0.1)
-    parser.add_argument(
-        "-i", "--num_islands", type=int, default=2
-    )  # Number of separate evolutionary islands
-    parser.add_argument(
-        "-migp", "--migration_probability", type=float, default=0.9
-    )  # Migration probability
-    parser.add_argument("-m", "--num_migrants", type=int, default=1)
-    parser.add_argument("-pln", "--pollination", action="store_true")
-    parser.add_argument(
-        "-t", "--top_n", type=int, default=1
-    )  # Print top-n best individuals on each island in summary.
-    parser.add_argument("-l", "--logging_int", type=int, default=10)  # Logging interval
-    config = parser.parse_args()
+    # Parse command-line arguments.
+    config, _ = parse_arguments(comm)
 
     # Set up separate logger for Propulate optimization.
     set_logger_config(
-        level=logging.INFO,  # logging level
-        log_file=f"{config.checkpoint}/islands.log",  # logging path
+        level=config.logging_level,  # Logging level
+        log_file=f"{config.checkpoint}/{pathlib.Path(__file__).stem}.log",  # Logging path
         log_to_stdout=True,  # Print log on stdout.
         log_rank=False,  # Do not prepend MPI rank to logging messages.
         colors=True,  # Use colors.
@@ -94,10 +36,10 @@
     propagator = get_default_propagator(  # Get default evolutionary operator.
         pop_size=config.pop_size,  # Breeding pool size
         limits=limits,  # Search-space limits
-        mate_prob=config.crossover_probability,  # Crossover probability
-        mut_prob=config.mutation_probability,  # Mutation probability
-        random_prob=config.random_init_probability,  # Random-initialization probability
-        rng=rng,  # Random number generator
+        crossover_prob=config.crossover_probability,  # Crossover probability
+        mutation_prob=config.mutation_probability,  # Mutation probability
+        random_init_prob=config.random_init_probability,  # Random-initialization probability
+        rng=rng,  # Separate random number generator for Propulate optimization
     )
 
     # Set up migration topology.
@@ -128,5 +70,7 @@
 
     # Run actual optimization.
     islands.evolve(
-        top_n=config.top_n, logging_interval=config.logging_int, debug=config.verbosity
+        top_n=config.top_n,
+        logging_interval=config.logging_interval,
+        debug=config.verbosity,
     )
diff --git a/tutorials/minimum_working_example.py b/tutorials/minimum_working_example.py
new file mode 100644
index 00000000..1b5e081f
--- /dev/null
+++ b/tutorials/minimum_working_example.py
@@ -0,0 +1,38 @@
+"""Minimum working example showing how to use Propulate."""
+import propulate
+from mpi4py import MPI
+import random
+
+# Set the communicator and the optimization parameters
+comm = MPI.COMM_WORLD
+rng = random.Random(MPI.COMM_WORLD.rank)
+population_size = comm.size * 2
+generations = 100
+checkpoint = "./propulate_checkpoints"
+propulate.utils.set_logger_config()
+
+
+# Define the function to minimize and the search space, e.g., a 2D sphere function on (-5.12, 5.12)^2.
+def loss_fn(params):
+    """Loss function to minimize."""
+    return params["x"] ** 2 + params["y"] ** 2
+
+
+limits = {"x": (-5.12, 5.12), "y": (-5.12, 5.12)}
+
+# Initialize the propagator and propulator with default parameters.
+propagator = propulate.utils.get_default_propagator(
+    pop_size=population_size, limits=limits, rng=rng
+)
+propulator = propulate.Propulator(
+    loss_fn=loss_fn,
+    propagator=propagator,
+    rng=rng,
+    island_comm=comm,
+    generations=generations,
+    checkpoint_path=checkpoint,
+)
+
+# Run optimization and get summary of results.
+propulator.propulate()
+propulator.summarize()
diff --git a/tutorials/multi_rank_workers_example.py b/tutorials/multi_rank_workers_example.py
new file mode 100644
index 00000000..446d46f0
--- /dev/null
+++ b/tutorials/multi_rank_workers_example.py
@@ -0,0 +1,117 @@
+"""
+Multi-rank workers in Propulate using the example of a parallel sphere function.
+Tested with 8 processes overall, 2 islands, and 2 ranks per worker, where each worker calculates one of the squared
+terms in the (in this case) two-dimensional sphere function. In general, the parallel sphere function's dimension
+should equal the number of ranks per worker.
+"""
+import random
+import pathlib
+from typing import Dict
+
+import numpy as np
+from mpi4py import MPI
+
+from propulate import Islands
+from propulate.propagators import SelectMin, SelectMax
+from propulate.utils import get_default_propagator, set_logger_config
+from function_benchmark import parse_arguments
+
+
+def parallel_sphere(params: Dict[str, float], comm: MPI.Comm = MPI.COMM_SELF) -> float:
+    """
+    Parallel sphere function to showcase using multi-rank workers in Propulate.
+
+    Sphere function: continuous, convex, separable, differentiable, unimodal
+
+    Input domain: -5.12 <= x, y <= 5.12
+    Global minimum 0 at (x, y) = (0, 0)
+
+    Parameters
+    ----------
+    params : Dict[str, float]
+        The function parameters.
+    comm : MPI.Comm
+        The communicator of the worker.
+
+    Returns
+    -------
+    float
+        The function value.
+    """
+    if comm != MPI.COMM_SELF:
+        term = (
+            list(params.values())[comm.rank] ** 2
+        )  # Each rank squares one of the inputs.
+        return comm.allreduce(term)  # Return the sum over all squared inputs.
+    else:
+        return np.sum(np.array(list(params.values())) ** 2).item()
+
+
+if __name__ == "__main__":
+    full_world_comm = MPI.COMM_WORLD  # Get full world communicator.
+
+    config, _ = parse_arguments()
+
+    # Set up separate logger for Propulate optimization.
+    set_logger_config(
+        level=config.logging_level,  # Logging level
+        log_file=f"{config.checkpoint}/{pathlib.Path(__file__).stem}.log",  # Logging path
+        log_to_stdout=True,  # Print log on stdout.
+        log_rank=False,  # Do not prepend MPI rank to logging messages.
+        colors=True,  # Use colors.
+    )
+
+    rng = random.Random(
+        config.seed + full_world_comm.rank
+    )  # Separate random number generator for optimization.
+    # Set callable function + search-space limits.
+    limits = {
+        "a": (-5.12, 5.12),
+        "b": (-5.12, 5.12),
+    }
+    function = parallel_sphere
+
+    # Set up evolutionary operator.
+    propagator = get_default_propagator(  # Get default evolutionary operator.
+        pop_size=config.pop_size,  # Breeding pool size
+        limits=limits,  # Search-space limits
+        crossover_prob=config.crossover_probability,  # Crossover probability
+        mutation_prob=config.mutation_probability,  # Mutation probability
+        random_init_prob=config.random_init_probability,  # Random-initialization probability
+        rng=rng,  # Separate random number generator for Propulate optimization
+    )
+
+    # Set up migration topology.
+    migration_topology = (
+        config.num_migrants
+        * np.ones(  # Set up fully connected migration topology.
+            (config.num_islands, config.num_islands), dtype=int
+        )
+    )
+    np.fill_diagonal(
+        migration_topology, 0
+    )  # An island does not send migrants to itself.
+
+    # Set up island model.
+    islands = Islands(
+        loss_fn=function,  # Loss function to be minimized
+        propagator=propagator,  # Propagator, i.e., evolutionary operator to be used
+        rng=rng,  # Separate random number generator for Propulate optimization
+        generations=config.generations,  # Overall number of generations
+        num_islands=config.num_islands,  # Number of islands
+        migration_topology=migration_topology,  # Migration topology
+        migration_probability=config.migration_probability,  # Migration probability
+        emigration_propagator=SelectMin,  # Selection policy for emigrants
+        immigration_propagator=SelectMax,  # Selection policy for immigrants
+        pollination=config.pollination,  # Whether to use pollination or migration
+        checkpoint_path=config.checkpoint,  # Checkpoint path
+        # ----- SPECIFIC FOR MULTI-RANK UCS ----
+        ranks_per_worker=config.ranks_per_worker,  # Number of ranks per (multi rank) worker
+    )
+
+    # Run actual optimization.
+    islands.evolve(
+        top_n=config.top_n,  # Print top-n best individuals on each island in summary.
+        logging_interval=config.logging_interval,  # Logging interval
+        debug=config.verbosity,  # Debug level
+    )
diff --git a/tutorials/propulator_example.py b/tutorials/propulator_example.py
index 58b05453..f0d13cf3 100755
--- a/tutorials/propulator_example.py
+++ b/tutorials/propulator_example.py
@@ -1,13 +1,11 @@
-#!/usr/bin/env python3
+"""Simple Propulator example script using the default genetic propagator."""
 import random
-import argparse
-import logging
+import pathlib
 
 from mpi4py import MPI
-
 from propulate import Propulator
 from propulate.utils import get_default_propagator, set_logger_config
-from function_benchmark import get_function_search_space
+from function_benchmark import get_function_search_space, parse_arguments
 
 if __name__ == "__main__":
     comm = MPI.COMM_WORLD
@@ -18,61 +16,13 @@
             "# PROPULATE: Parallel Propagator of Populations #\n"
             "#################################################\n"
         )
-
-    parser = argparse.ArgumentParser(
-        prog="Simple Propulator example",
-        description="Set up and run a basic Propulator optimization of mathematical functions.",
-    )
-    parser.add_argument(  # Function to optimize
-        "-f",
-        "--function",
-        type=str,
-        choices=[
-            "bukin",
-            "eggcrate",
-            "himmelblau",
-            "keane",
-            "leon",
-            "rastrigin",
-            "schwefel",
-            "sphere",
-            "step",
-            "rosenbrock",
-            "quartic",
-            "bisphere",
-            "birastrigin",
-            "griewank",
-        ],
-        default="sphere",
-    )
-    parser.add_argument(
-        "-g", "--generations", type=int, default=1000
-    )  # Number of generations
-    parser.add_argument(
-        "-s", "--seed", type=int, default=0
-    )  # Seed for Propulate random number generator
-    parser.add_argument("-v", "--verbosity", type=int, default=1)  # Verbosity level
-    parser.add_argument(
-        "-ckpt", "--checkpoint", type=str, default="./"
-    )  # Path for loading and writing checkpoints.
-    parser.add_argument(
-        "-p", "--pop_size", type=int, default=2 * comm.size
-    )  # Breeding pool size
-    parser.add_argument(
-        "-cp", "--crossover_probability", type=float, default=0.7
-    )  # Crossover probability
-    parser.add_argument(
-        "-mp", "--mutation_probability", type=float, default=0.4
-    )  # Mutation probability
-    parser.add_argument("-rp", "--random_init_probability", type=float, default=0.1)
-    parser.add_argument("-t", "--top_n", type=int, default=1)
-    parser.add_argument("-l", "--logging_int", type=int, default=10)
-    config = parser.parse_args()
+    # Parse command-line arguments.
+    config, _ = parse_arguments(comm)
 
     # Set up separate logger for Propulate optimization.
     set_logger_config(
-        level=logging.INFO,  # logging level
-        log_file=f"{config.checkpoint}/propulator.log",  # logging path
+        level=config.logging_level,  # Logging level
+        log_file=f"{config.checkpoint}/{pathlib.Path(__file__).stem}.log",  # Logging path
         log_to_stdout=True,  # Print log on stdout.
         log_rank=False,  # Do not prepend MPI rank to logging messages.
         colors=True,  # Use colors.
@@ -84,27 +34,32 @@
     function, limits = get_function_search_space(
         config.function
     )  # Get callable function + search-space limits.
-
     # Set up evolutionary operator.
     propagator = get_default_propagator(  # Get default evolutionary operator.
         pop_size=config.pop_size,  # Breeding pool size
         limits=limits,  # Search-space limits
-        mate_prob=config.crossover_probability,  # Crossover probability
-        mut_prob=config.mutation_probability,  # Mutation probability
-        random_prob=config.random_init_probability,  # Random-initialization probability
-        rng=rng,  # Random number generator
+        crossover_prob=config.crossover_probability,  # Crossover probability
+        mutation_prob=config.mutation_probability,  # Mutation probability
+        random_init_prob=config.random_init_probability,  # Random-initialization probability
+        rng=rng,  # Separate random number generator for Propulate optimization
     )
 
     # Set up propulator performing actual optimization.
     propulator = Propulator(
         loss_fn=function,
         propagator=propagator,
-        comm=comm,
+        rng=rng,
+        island_comm=comm,
         generations=config.generations,
         checkpoint_path=config.checkpoint,
-        rng=rng,
     )
 
     # Run optimization and print summary of results.
-    propulator.propulate(logging_interval=config.logging_int, debug=config.verbosity)
-    propulator.summarize(top_n=config.top_n, debug=config.verbosity)
+    propulator.propulate(
+        logging_interval=config.logging_interval,  # Logging interval
+        debug=config.verbosity,  # Verbosity level
+    )
+    propulator.summarize(
+        top_n=config.top_n,  # Print top-n best individuals on each island in summary.
+        debug=config.verbosity,  # Debug level
+    )
diff --git a/tutorials/pso_example.py b/tutorials/pso_example.py
index 5f0249a9..5f70f576 100644
--- a/tutorials/pso_example.py
+++ b/tutorials/pso_example.py
@@ -1,12 +1,9 @@
 """
 This files contains an example use case for the PSO propagators. Here, you can choose between benchmark functions and
-optimize them.
-
-The example shows, how to set up Propulate in order to use it with PSO.
+optimize them. The example shows how to set up Propulate in order to use it with PSO.
 """
-import argparse
+import pathlib
 import random
-from typing import Dict
 
 from mpi4py import MPI
 
@@ -19,7 +16,7 @@
     CanonicalPSO,
     InitUniformPSO,
 )
-from tutorials.function_benchmark import get_function_search_space
+from function_benchmark import get_function_search_space, parse_arguments
 
 if __name__ == "__main__":
     comm = MPI.COMM_WORLD
@@ -31,91 +28,15 @@
             "#################################################\n"
         )
 
-    parser = argparse.ArgumentParser(
-        prog="Simple PSO example",
-        description="Set up and run a basic particle swarm optimization of mathematical functions.",
-    )
-    parser.add_argument(  # Function to optimize
-        "-f",
-        "--function",
-        type=str,
-        choices=[
-            "bukin",
-            "eggcrate",
-            "himmelblau",
-            "keane",
-            "leon",
-            "rastrigin",
-            "schwefel",
-            "sphere",
-            "step",
-            "rosenbrock",
-            "quartic",
-            "bisphere",
-            "birastrigin",
-            "griewank",
-        ],
-        default="sphere",
-    )
-    parser.add_argument(
-        "-g", "--generations", type=int, default=1000
-    )  # Number of generations
-    parser.add_argument(
-        "-s", "--seed", type=int, default=0
-    )  # Seed for Propulate random number generator
-    parser.add_argument(
-        "-v", "--verbosity", type=int, default=1, choices=range(6)
-    )  # Verbosity level
-    parser.add_argument(
-        "-ckpt", "--checkpoint", type=str, default="./"
-    )  # Path for loading and writing checkpoints.
-    parser.add_argument(
-        "-p", "--pop_size", type=int, default=2 * comm.size
-    )  # Breeding pool size
-    parser.add_argument(
-        "-var",
-        "--variant",
-        type=str,
-        choices=["Basic", "VelocityClamping", "Constriction", "Canonical"],
-        default="Basic",
-    )  # PSO variant to run
-
-    hp_set: Dict[str, bool] = {
-        "inertia": False,
-        "cognitive": False,
-        "social": False,
-    }
-
-    class ParamSettingCatcher(argparse.Action):
-        """
-        This class extends argparse's Action class in order to allow for an action, that logs, if one of the PSO HP
-        was actually set.
-        """
-
-        def __call__(self, parser, namespace, values, option_string=None):
-            hp_set[self.dest] = True
-            super().__call__(parser, namespace, values, option_string)
-
-    parser.add_argument(
-        "--inertia", type=float, default=0.729, action=ParamSettingCatcher
-    )  # Inertia weight
-    parser.add_argument(
-        "--cognitive", type=float, default=1.49445, action=ParamSettingCatcher
-    )  # Cognitive factor
-    parser.add_argument(
-        "--social", type=float, default=1.49445, action=ParamSettingCatcher
-    )  # Social factor
-    parser.add_argument(
-        "--clamping_factor", type=float, default=0.6
-    )  # Clamping factor for velocity clamping
-    parser.add_argument("-t", "--top_n", type=int, default=1)
-    parser.add_argument("-l", "--logging_int", type=int, default=20)
-    config = parser.parse_args()
+    config, hp_set = parse_arguments(comm)
 
     # Set up separate logger for Propulate optimization.
     set_logger_config(
-        level=config.logging_int,  # logging level
-        log_file=f"{config.checkpoint}/propulator.log",  # logging path
+        level=config.logging_level,  # Logging level
+        log_file=f"{config.checkpoint}/{pathlib.Path(__file__).stem}.log",  # Logging path
+        log_to_stdout=True,  # Print log on stdout.
+        log_rank=False,  # Do not prepend MPI rank to logging messages.
+        colors=True,  # Use colors.
     )
 
     rng = random.Random(
@@ -167,10 +88,12 @@ def __call__(self, parser, namespace, values, option_string=None):
     propulator = Propulator(
         function,
         propagator,
-        comm=comm,
+        rng=rng,
+        island_comm=comm,
         generations=config.generations,
         checkpoint_path=config.checkpoint,
-        rng=rng,
     )
-    propulator.propulate(config.logging_int, config.verbosity)
+    propulator.propulate(
+        logging_interval=config.logging_interval, debug=config.verbosity
+    )
     propulator.summarize(top_n=config.top_n, debug=config.verbosity)
diff --git a/tutorials/torch_example.py b/tutorials/torch_example.py
index 9d671b7c..d6ceb16d 100755
--- a/tutorials/torch_example.py
+++ b/tutorials/torch_example.py
@@ -5,6 +5,7 @@
 This script was tested on a single compute node with 4 GPUs. Note that you need to adapt ``GPUS_PER_NODE`` (see ll. 25).
 """
 import logging
+import pathlib
 import random
 from typing import Union, Dict, Tuple
 
@@ -337,16 +338,16 @@ def ind_loss(params: Dict[str, Union[int, float, str]]) -> float:
     propagator = get_default_propagator(  # Get default evolutionary operator.
         pop_size=pop_size,  # Breeding population size
         limits=limits,  # Search space
-        mate_prob=0.7,  # Crossover probability
-        mut_prob=0.4,  # Mutation probability
-        random_prob=0.1,  # Random-initialization probability
-        rng=rng,  # Random number generator for evolutionary optimizer
+        crossover_prob=0.7,  # Crossover probability
+        mutation_prob=0.4,  # Mutation probability
+        random_init_prob=0.1,  # Random-initialization probability
+        rng=rng,  # Separate random number generator for Propulate optimization
     )
 
     # Set up separate logger for Propulate optimization.
     set_logger_config(
-        level=logging.INFO,  # logging level
-        log_file="./propulator.log",  # logging path
+        level=logging.INFO,  # Logging level
+        log_file=f"{log_path}/{pathlib.Path(__file__).stem}.log",  # Logging path
         log_to_stdout=True,  # Print log on stdout.
         log_rank=False,  # Do not prepend MPI rank to logging messages.
         colors=True,  # Use colors.
@@ -356,10 +357,10 @@ def ind_loss(params: Dict[str, Union[int, float, str]]) -> float:
     propulator = Propulator(
         loss_fn=ind_loss,  # Loss function to optimize
         propagator=propagator,  # Evolutionary operator
-        comm=comm,  # Communicator
+        rng=rng,  # Random number generator
+        island_comm=comm,  # Communicator
         generations=num_generations,  # Number of generations per worker
         checkpoint_path=log_path,  # Path to save checkpoints to
-        rng=rng,  # Random number generator
     )
 
     # Run optimization and print summary of results.


Test patch:
diff --git a/tests/conftest.py b/tests/conftest.py
deleted file mode 100644
index 573c292b..00000000
--- a/tests/conftest.py
+++ /dev/null
@@ -1,10 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-    Dummy conftest.py for propulate.
-
-    If you don't know what this is for, just leave it empty.
-    Read more about conftest.py under:
-    https://pytest.org/latest/plugins.html
-"""
-
-# import pytest
diff --git a/tests/test_cmaes.py b/tests/test_cmaes.py
index 2285ef0d..bb0ac39b 100644
--- a/tests/test_cmaes.py
+++ b/tests/test_cmaes.py
@@ -1,7 +1,6 @@
 import random
 import tempfile
 from typing import Dict
-from operator import attrgetter
 
 import numpy as np
 
@@ -18,20 +17,19 @@ def sphere(params: Dict[str, float]) -> float:
 
     Parameters
     ----------
-    params: dict[str, float]
-            function parameters
+    params : Dict[str, float]
+        The function parameters.
+
     Returns
     -------
     float
-        function value
+        The function value.
     """
-    return np.sum(np.array(list(params.values())) ** 2)
+    return np.sum(np.array(list(params.values())) ** 2).item()
 
 
-def test_PSO():
-    """
-    Test single worker using Propulator to optimize sphere using a PSO propagator.
-    """
+def test_cmaes():
+    """Test single worker using Propulator to optimize sphere using a CMA-ES propagator."""
     rng = random.Random(42)  # Separate random number generator for optimization.
     limits = {
         "a": (-5.12, 5.12),
@@ -39,22 +37,18 @@ def test_PSO():
     }
     with tempfile.TemporaryDirectory() as checkpoint_path:
         # Set up evolutionary operator.
-
         adapter = BasicCMA()
         propagator = CMAPropagator(adapter, limits, rng=rng)
 
-        # Set up propulator performing actual optimization.
+        # Set up Propulator performing actual optimization.
         propulator = Propulator(
             loss_fn=sphere,
             propagator=propagator,
-            generations=10,
-            checkpoint_path=checkpoint_path,
             rng=rng,
+            generations=100,
+            checkpoint_path=checkpoint_path,
         )
-
         # Run optimization and print summary of results.
         propulator.propulate()
-        propulator.summarize()
-        best = min(propulator.population, key=attrgetter("loss"))
-
-        assert best.loss < 10.0
+        best = propulator.summarize(top_n=1, debug=2)
+        assert best[0][0].loss < 10**-1
diff --git a/tests/test_propulator_sphere.py b/tests/test_propulator_sphere.py
index d6fae54e..d843d5f0 100644
--- a/tests/test_propulator_sphere.py
+++ b/tests/test_propulator_sphere.py
@@ -1,7 +1,6 @@
 import random
 import tempfile
 from typing import Dict
-from operator import attrgetter
 import logging
 
 import numpy as np
@@ -19,20 +18,19 @@ def sphere(params: Dict[str, float]) -> float:
 
     Parameters
     ----------
-    params: dict[str, float]
-            function parameters
+    params : Dict[str, float]
+        The function parameters.
+
     Returns
     -------
     float
-        function value
+        The function value.
     """
-    return np.sum(np.array(list(params.values())) ** 2)
+    return np.sum(np.array(list(params.values())) ** 2).item()
 
 
-def test_Propulator():
-    """
-    Test single worker using Propulator to optimize sphere.
-    """
+def test_propulator():
+    """Test single worker using Propulator to optimize sphere."""
     rng = random.Random(42)  # Separate random number generator for optimization.
     limits = {
         "a": (-5.12, 5.12),
@@ -47,27 +45,26 @@ def test_Propulator():
             colors=True,
         )
         # Set up evolutionary operator.
-        propagator = get_default_propagator(  # Get default evolutionary operator.
-            pop_size=4,  # Breeding pool size
-            limits=limits,  # Search-space limits
-            mate_prob=0.7,  # Crossover probability
-            mut_prob=9.0,  # Mutation probability
-            random_prob=0.1,  # Random-initialization probability
-            rng=rng,  # Random number generator
+        propagator = get_default_propagator(
+            pop_size=4,
+            limits=limits,
+            crossover_prob=0.7,
+            mutation_prob=9.0,
+            random_init_prob=0.1,
+            rng=rng,
         )
 
         # Set up propulator performing actual optimization.
         propulator = Propulator(
             loss_fn=sphere,
             propagator=propagator,
-            generations=10,
-            checkpoint_path=checkpoint_path,
             rng=rng,
+            generations=100,
+            checkpoint_path=checkpoint_path,
         )
 
         # Run optimization and print summary of results.
         propulator.propulate()
-        propulator.summarize()
-        best = min(propulator.population, key=attrgetter("loss"))
+        best = propulator.summarize()
 
-        assert best.loss < 0.8
+        assert best[0][0].loss < 0.8
diff --git a/tests/test_pso.py b/tests/test_pso.py
index 437384c3..e89289f9 100644
--- a/tests/test_pso.py
+++ b/tests/test_pso.py
@@ -1,8 +1,8 @@
 import random
 import tempfile
 from typing import Dict
-from operator import attrgetter
 
+from mpi4py import MPI
 import numpy as np
 
 from propulate import Propulator
@@ -19,20 +19,18 @@ def sphere(params: Dict[str, float]) -> float:
 
     Parameters
     ----------
-    params: dict[str, float]
-            function parameters
+    params : Dict[str, float]
+        The function parameters.
     Returns
     -------
     float
-        function value
+        The function value.
     """
-    return np.sum(np.array(list(params.values())) ** 2)
+    return np.sum(np.array(list(params.values())) ** 2).item()
 
 
-def test_PSO():
-    """
-    Test single worker using Propulator to optimize sphere using a PSO propagator.
-    """
+def test_pso():
+    """Test single worker using Propulator to optimize sphere using a PSO propagator."""
     rng = random.Random(42)  # Separate random number generator for optimization.
     limits = {
         "a": (-5.12, 5.12),
@@ -40,30 +38,27 @@ def test_PSO():
     }
     with tempfile.TemporaryDirectory() as checkpoint_path:
         # Set up evolutionary operator.
-
         pso_propagator = BasicPSO(
-            0.729,
-            1.49334,
-            1.49445,
-            0,  # MPI rank TODO fix when implemented proper MPI parallel tests
-            limits,
-            rng,
+            inertia=0.729,
+            c_cognitive=1.49334,
+            c_social=1.49445,
+            rank=MPI.COMM_WORLD.rank,  # MPI rank TODO fix when implemented proper MPI parallel tests
+            limits=limits,
+            rng=rng,
         )
-        init = InitUniformPSO(limits, rng=rng, rank=0)
+        init = InitUniformPSO(limits, rng=rng, rank=MPI.COMM_WORLD.rank)
         propagator = Conditional(1, pso_propagator, init)  # TODO MPIify
 
         # Set up propulator performing actual optimization.
         propulator = Propulator(
             loss_fn=sphere,
             propagator=propagator,
-            generations=10,
-            checkpoint_path=checkpoint_path,
             rng=rng,
+            generations=100,
+            checkpoint_path=checkpoint_path,
         )
 
         # Run optimization and print summary of results.
         propulator.propulate()
-        propulator.summarize()
-        best = min(propulator.population, key=attrgetter("loss"))
-
-        assert best.loss < 30.0
+        best = propulator.summarize(top_n=1, debug=2)
+        assert best[0][0].loss < 30.0